1 | function NULLSPACE = create_nullspace(M,INC,r,n) |
---|
2 | |
---|
3 | % Fit a rank r matrix to a matrix with missing data. This uses my |
---|
4 | % algorithm of taking the most orthogonal space to the null space generated |
---|
5 | % by r-tuples of the columns of M. In principle, it would be nice |
---|
6 | % to use all r-tuples, but instead we'll pick n random ones. |
---|
7 | |
---|
8 | %global NULLSPACEexternal |
---|
9 | |
---|
10 | % 10 is the default for the simulated motion experiments. |
---|
11 | % however, for the real intensity images, may want to change |
---|
12 | % since matrix dimensions are very different, eg, 9 rows |
---|
13 | % instead of 40 in M. |
---|
14 | |
---|
15 | %MAXNULLSIZE = 100; |
---|
16 | MAXNULLSIZE = 10; |
---|
17 | % if the nullspace gets MAXNULLSIZE times the number of columns in M, we say |
---|
18 | % that that's enough, and stop generating elements. |
---|
19 | |
---|
20 | Mdim = size(M); |
---|
21 | numrows = Mdim(1); |
---|
22 | numcols = Mdim(2); |
---|
23 | incsum = sum(INC); |
---|
24 | NULLSPACE = []; |
---|
25 | current_row = 1; |
---|
26 | for i = 1:n |
---|
27 | if size(NULLSPACE,2) >= numrows*MAXNULLSIZE, break, end |
---|
28 | |
---|
29 | col_nums = find(INC(current_row,:) ~= 0); |
---|
30 | |
---|
31 | % first check that there actually is enough data to recover that row. |
---|
32 | % note that although we simply choose columns at random if the check fails, |
---|
33 | % it will be impossible to recover the row at any point in the process, |
---|
34 | % since it is under-constrained. |
---|
35 | |
---|
36 | if (length(col_nums) < r) |
---|
37 | col_nums = diff_rand_ints([],r,1,num_cols); |
---|
38 | else |
---|
39 | col_nums = col_nums(diff_rand_ints([],r,1,length(col_nums))); |
---|
40 | end |
---|
41 | |
---|
42 | incsub = INC(:,col_nums); |
---|
43 | rowsum = sum(incsub')'; |
---|
44 | fullrows = rowsum==r; |
---|
45 | if sum(fullrows) > r |
---|
46 | |
---|
47 | submatrix = M(find(fullrows),col_nums); |
---|
48 | subnull = nulleps(submatrix,.001,-1); |
---|
49 | % .1 is a total hack. Also, I'm not sure if its good to check sum(fullrows) > r. |
---|
50 | if size(submatrix,1) == size(submatrix,2) + size(subnull,2) |
---|
51 | %% The null space of the submatrix, combined with the submatrix should be a full |
---|
52 | %% rank square matrix. However, if not, it means that the submatrix was rank |
---|
53 | %% deficient (to within tolerance) and the null matrix is too big. We can't know |
---|
54 | %% which of these components are the correct ones, and which just arise because the |
---|
55 | %% submatrix was rank deficient. This test will make the whole method not work if |
---|
56 | %% M really does have rank less than r (in which case the whole null space is right. |
---|
57 | nullTEMP = zeros(numrows,size(subnull,2)); |
---|
58 | nullTEMP(fullrows,:) = subnull; |
---|
59 | NULLSPACE = [NULLSPACE, nullTEMP]; |
---|
60 | end |
---|
61 | end |
---|
62 | |
---|
63 | %% Below is a hack added to try to focus attention on the rows for which |
---|
64 | %% we have no information. Half the time, it makes the current row to be |
---|
65 | %% one for which we don't have any information yet in the nullspace. The |
---|
66 | %% other half of the time it's just incremented. I'm not sure whether this |
---|
67 | %% helps significantly. |
---|
68 | if isempty(NULLSPACE) |
---|
69 | sumnonzeros = []; |
---|
70 | nonsamprows = []; |
---|
71 | else |
---|
72 | sumnonzeros = (sum(NULLSPACE' ~= 0)')'; |
---|
73 | nonsamprows = find(0 == sumnonzeros); |
---|
74 | end |
---|
75 | if (rand(1) > .5) & (length(nonsamprows) > 0) |
---|
76 | current_row = nonsamprows(random_int(1,length(nonsamprows))); |
---|
77 | else |
---|
78 | current_row = 1+rem(current_row,numrows); |
---|
79 | end |
---|
80 | end |
---|
81 | %NULLSPACEexternal = NULLSPACE; |
---|