[37] | 1 | %fill_prmm Compute the null-space and fill PRMM. |
---|
| 2 | % |
---|
| 3 | % Parameters: |
---|
| 4 | |
---|
| 5 | function [P,X, u1,u2, lambda, info] = fill_prmm(M, Idepths, central, ... |
---|
| 6 | opt, info) |
---|
| 7 | |
---|
| 8 | [NULLSPACE, result] = create_nullspace(M, Idepths, central, ... |
---|
| 9 | opt.create_nullspace); |
---|
| 10 | |
---|
| 11 | info.create_nullspace = opt.create_nullspace; |
---|
| 12 | info.sequence{end}.tried = result.tried; |
---|
| 13 | info.sequence{end}.tried_perc = result.tried/comb(size(M,2), 4) *100; |
---|
| 14 | info.sequence{end}.used = result.used; |
---|
| 15 | info.sequence{end}.used_perc = result.used/result.tried *100; |
---|
| 16 | info.sequence{end}.failed = result.failed; |
---|
| 17 | info.sequence{end}.size_nullspace = size(NULLSPACE); |
---|
| 18 | |
---|
| 19 | if opt.verbose |
---|
| 20 | disp(sprintf('Tried/used: %d/%d (%.1e %%/ %2.1f %%)', result.tried, ... |
---|
| 21 | result.used, info.sequence{end}.tried_perc, ... |
---|
| 22 | info.sequence{end}.used_perc)); |
---|
| 23 | disp(sprintf('%d x %d'' is size of the nullspace', size(NULLSPACE,2), ... |
---|
| 24 | size(NULLSPACE,1))); end |
---|
| 25 | |
---|
| 26 | [m,n] = size(M); m = m/3; |
---|
| 27 | |
---|
| 28 | if size(NULLSPACE) == [0 0] |
---|
| 29 | P = []; X = []; u1 = 1:m; u2 = 1:n; lambda=[]; |
---|
| 30 | else |
---|
| 31 | r = 4; |
---|
| 32 | [L, S] = nullspace2L(NULLSPACE, r, opt); |
---|
| 33 | clear NULLSPACE; |
---|
| 34 | |
---|
| 35 | if isempty(opt.create_nullspace), threshold = .01; |
---|
| 36 | else, threshold = opt.create_nullspace.threshold; end |
---|
| 37 | if svd_suff_data(S, r, threshold) |
---|
| 38 | if opt.verbose |
---|
| 39 | dS = diag(S); %disp(diag(S)); |
---|
| 40 | fprintf(1,'Smallest 2r singular values:%s.\n', sprintf(' %f', ... |
---|
| 41 | dS(end-2*r:end))); end |
---|
| 42 | [Mdepths, lambda] = L2depths(L, M, Idepths, opt); info.Mdepths = Mdepths; |
---|
| 43 | [P,X, u1b,u2] = approximate(Mdepths, r, L, opt); |
---|
| 44 | u1 = union(ceil(u1b/3),[]); killb = setdiff(k2i(u1),u1b); |
---|
| 45 | if ~isempty(killb), r1b = setdiff(1:3*m,u1b); kill = killb; |
---|
| 46 | for ib = killb(1:end-1), lower = find(killb > ib); |
---|
| 47 | if kill(lower(1)-1) < kill(lower(1)) - 1, |
---|
| 48 | kill(lower) = kill(lower) -1; end; end |
---|
| 49 | P = P(setdiff(1:length(r1b),kill),:);end |
---|
| 50 | lambda = lambda(setdiff(1:m,u1),setdiff(1:n,u2)); % to fit P*X |
---|
| 51 | else P=[]; X=[]; u1=1:m; u2=1:n; lambda=[]; end |
---|
| 52 | end |
---|
| 53 | |
---|
| 54 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 55 | |
---|
| 56 | function [L,S] = nullspace2L(NULLSPACE, r, opt) |
---|
| 57 | % Compute the basis of MM (L) from the null-space. |
---|
| 58 | |
---|
| 59 | if opt.verbose, fprintf(1,'Computing the basis...'); tic; end |
---|
| 60 | if size(NULLSPACE,2) < 10*size(NULLSPACE,1) % orig:[A,S,U] = svd(NULLSPACE',0); |
---|
| 61 | [U,S,V] = svd(NULLSPACE); |
---|
| 62 | else |
---|
| 63 | [U,SS] = eig(NULLSPACE*NULLSPACE'); |
---|
| 64 | dSS = diag(SS); l = length(dSS); |
---|
| 65 | [sortdSS(l:-1:1),I] = sort(dSS); I(l:-1:1) = I; U = U(:,I); |
---|
| 66 | sortdSS = max(sortdSS,zeros(1,l)); sortdSS=sqrt(sortdSS); |
---|
| 67 | for i=1:l, S(i,i) = sortdSS(i); end |
---|
| 68 | end |
---|
| 69 | L = U(:,end+1-r:end); |
---|
| 70 | if opt.verbose, disp(['(' num2str(toc) ' sec)']); end |
---|
| 71 | |
---|
| 72 | function r = comb(n,k) |
---|
| 73 | |
---|
| 74 | % returns combination number n over k |
---|
| 75 | |
---|
| 76 | r=1; |
---|
| 77 | for i=1:k |
---|
| 78 | r=r*(n-i+1)/i; |
---|
| 79 | end |
---|
| 80 | |
---|
| 81 | function y = svd_suff_data(S,r, threshold) |
---|
| 82 | |
---|
| 83 | % S is the singular value part of the svd of the nullspaces of the column |
---|
| 84 | % r-tuples. We'll want to be able to take the r least significant columns |
---|
| 85 | % of U. This is right because the columns of U should span the whole space |
---|
| 86 | % that M's columns might span. That is, M is FxP. The columns of U should |
---|
| 87 | % span the F-dimensional Euclidean space, since U is FxF. However, we want |
---|
| 88 | % to make sure that the F-r-1'th singular value of S isn't tiny. If it is, |
---|
| 89 | % our answer is totally unreliable, because the nullspaces of the column |
---|
| 90 | % r-tuples don't have sufficient rank. If this happens, it means that the |
---|
| 91 | % intersection of the column cross-product spaces is bigger than r-dimensional, |
---|
| 92 | % and randomly choosing an r-dimensional subspace of that isn't likely to |
---|
| 93 | % give the right answer. |
---|
| 94 | |
---|
| 95 | Snumcols = size(S,2); |
---|
| 96 | Snumrows = size(S,1); |
---|
| 97 | if (Snumrows == 0 | Snumcols + r < Snumrows | Snumrows <= r) |
---|
| 98 | y = 0; |
---|
| 99 | else |
---|
| 100 | y = S(Snumrows-r,Snumrows-r) > threshold; |
---|
| 101 | end |
---|