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 |
---|