[37] | 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;
|
---|