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