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