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