[37] | 1 | %approximate Compute r-rank approximation of MM using null-space. |
---|
| 2 | % |
---|
| 3 | % The approximated matrix is P*X. |
---|
| 4 | |
---|
| 5 | function [P,X, u1,u2] = approximate(M, r, P, opt) |
---|
| 6 | |
---|
| 7 | [m n] = size(M); |
---|
| 8 | rows = find(mean(abs(P')) > opt.tol); %find(sum(P'))'; |
---|
| 9 | nonrows = setdiff(1:m, rows); |
---|
| 10 | |
---|
| 11 | if opt.verbose, fprintf(1,'Immersing columns of MM into the basis...'); tic;end |
---|
| 12 | [Mapp, noncols, X] = approx_matrix(M(rows,:), P(rows,:),r, opt); |
---|
| 13 | if opt.verbose, disp(['(' num2str(toc) ' sec)']); end |
---|
| 14 | |
---|
| 15 | % The rows of NULLSPACE now contain all the nullspaces of the crossproduct |
---|
| 16 | % spaces. Use the nullspace to approximate M, taking r least principal |
---|
| 17 | % components. |
---|
| 18 | cols = setdiffJ(1:n, noncols); |
---|
| 19 | if length(cols) < r; % it will not be able to extend the matrix correctly |
---|
| 20 | u1 = 1:m; u2 = 1:n; |
---|
| 21 | else |
---|
| 22 | if isempty(nonrows), u1 = []; r1 = 1:m; else |
---|
| 23 | if opt.verbose, |
---|
| 24 | disp('Filling rows of MM which have not been computed yet...'); end |
---|
| 25 | [Mapp1, u1, Pnonrows] = extend_matrix(M(:,cols), Mapp(:,cols), ... |
---|
| 26 | X(:,cols), rows, nonrows, r,opt.tol); |
---|
| 27 | Mapp = []; Mapp(:,cols) = Mapp1; |
---|
| 28 | if length(u1) < length(nonrows), P(setdiff(nonrows,u1),:) = Pnonrows; end |
---|
| 29 | r1 = union(rows, setdiff(nonrows,u1)); P = P(r1,:); |
---|
| 30 | end |
---|
| 31 | if isempty(noncols), u2 = []; else |
---|
| 32 | if opt.verbose |
---|
| 33 | disp('Filling columns of MM which have not been computed yet...'); end |
---|
| 34 | [Mapp_tr, u2, Xnoncols_tr] = extend_matrix(M(r1,:)', Mapp(r1,cols)', P',... |
---|
| 35 | cols, noncols, r, opt.tol); |
---|
| 36 | Mapp2 = Mapp_tr'; Xnoncols = Xnoncols_tr'; |
---|
| 37 | if length(u2) < length(noncols), X(:,setdiff(noncols,u2)) = Xnoncols; end |
---|
| 38 | X = X(:, union(cols, setdiff(noncols, u2))); |
---|
| 39 | if sum(~sum(X)), keyboard; end |
---|
| 40 | end |
---|
| 41 | % These extensions are necessary because the nullspace might not allow us to |
---|
| 42 | % compute every row of the approximating rank r linear space, and then this |
---|
| 43 | % linear space might not allow us to fill in some columns, if they are |
---|
| 44 | % missing too much data. |
---|
| 45 | end |
---|
| 46 | if sum(~sum(X)), keyboard; end |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | %approx_matrix Immerse columns of MM into the basis. |
---|
| 50 | |
---|
| 51 | function [Mapp, misscols, X] = approx_matrix(M, P, r, opt) |
---|
| 52 | |
---|
| 53 | [m n] = size(M); misscols = []; |
---|
| 54 | Mapp(m, n) = NaN; X(r,n) = NaN; % memory allocation |
---|
| 55 | |
---|
| 56 | if ~isempty(P) |
---|
| 57 | % P has r columns, which span the r-D space that gives the best linear |
---|
| 58 | % surface to approximate M. |
---|
| 59 | for j = 1:n |
---|
| 60 | rows = find(~isnan(M(:,j))); |
---|
| 61 | if rank(P(rows,:), opt.tol) == r |
---|
| 62 | X(:,j) = P(rows,:)\M(rows,j); Mapp(:,j) = P*X(:,j); |
---|
| 63 | else |
---|
| 64 | misscols = [misscols, j]; |
---|
| 65 | end |
---|
| 66 | end |
---|
| 67 | end |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | %extend_matrix Fill rows of MM which have not been computed yet. |
---|
| 71 | % |
---|
| 72 | % [E, unrecovered] = extend_matrix(M,INC,subM, rows, nonrows, r) |
---|
| 73 | % |
---|
| 74 | % Parameters: |
---|
| 75 | % rows ... rows of M which have been used to find the solution |
---|
| 76 | % subM ... is a fit to just these rows |
---|
| 77 | % nonrows ... indicate rows of that still need to be fit |
---|
| 78 | % |
---|
| 79 | % Output parameters: |
---|
| 80 | % Pnonrows ... in terms of unrecovered rows |
---|
| 81 | |
---|
| 82 | function [E, unrecovered, Pnonrows] = extend_matrix(M, subM, X, rows, ... |
---|
| 83 | nonrows, r, tol) |
---|
| 84 | |
---|
| 85 | E(size(M,1),size(M,2)) = NaN; |
---|
| 86 | E(rows,:) = subM; unrecovered = []; row = 0; Pnonrows = []; |
---|
| 87 | for i = nonrows |
---|
| 88 | cols = find(~isnan(M(i,:))); |
---|
| 89 | if rank(X(:,cols)) == r |
---|
| 90 | if max(abs(M(i,cols)/X(:,cols) * X)) > 1/tol |
---|
| 91 | unrecovered = [unrecovered i]; |
---|
| 92 | else, row = row + 1; |
---|
| 93 | Pnonrows(row,:) = M(i,cols)/X(:,cols); E(i,:) = Pnonrows(row,:)*X; |
---|
| 94 | end |
---|
| 95 | else |
---|
| 96 | unrecovered = [unrecovered i]; |
---|
| 97 | end |
---|
| 98 | end |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | function c = setdiffJ(a,b) %J as Jacobs so that matlab's setdiff isn't damaged |
---|
| 102 | c = []; |
---|
| 103 | for i = a |
---|
| 104 | if ~member(i,b) |
---|
| 105 | c = [c,i]; |
---|
| 106 | end |
---|
| 107 | end |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | %member(e,s) Return 1 if element e is a member of set s. Otherwise return 0. |
---|
| 111 | |
---|
| 112 | function c = member(e,s) |
---|
| 113 | c = 0; |
---|
| 114 | for i = s |
---|
| 115 | if i == e |
---|
| 116 | c = 1; |
---|
| 117 | break |
---|
| 118 | end |
---|
| 119 | end |
---|