[37] | 1 | function [error, Wapp,stable] = shum(W,INC,r,INIT,num_it) |
---|
| 2 | % shum produces an approximation to the W matrix, using method of Shum, Ikeuchi and Reddy |
---|
| 3 | % W is a matrix of data. INC is an incidence matrix of the same dimension. If |
---|
| 4 | % INC is 1, the corresponding value in W is a real value. Otherwise, it is really |
---|
| 5 | % missing data, and that entry in W should be ignored. |
---|
| 6 | % r is the rank of the matrix we want to approximate W with. |
---|
| 7 | % INIT is an initial guess about the matrix, of the same size as W. If INIT is 0, |
---|
| 8 | % we initialize with a random matrix, with elements from 0 to 1. |
---|
| 9 | % We repeat for either num_it iterations, or till method converges to nearly |
---|
| 10 | % 0 error. |
---|
| 11 | % We return the matrix, and the residual error. By convention, if the method |
---|
| 12 | % doesn't converge to 10 decimal places we return the negation of the residual error. |
---|
| 13 | |
---|
| 14 | % W is FxP. V is rxP, U is Fxr. UV approximates W. |
---|
| 15 | % ENTERING_SHUM = 1 |
---|
| 16 | MAXSTEPS = num_it; |
---|
| 17 | stable = 1; |
---|
| 18 | Wdim = size(W); |
---|
| 19 | F = Wdim(1); |
---|
| 20 | P = Wdim(2); |
---|
| 21 | |
---|
| 22 | if (sum(sum(INC)) < r*(F + P - r) + P) |
---|
| 23 | % Return -2 is there's not enough information for this iterative method to work. |
---|
| 24 | % See Shum et al, for an explanation of this formula. |
---|
| 25 | error=-2; |
---|
| 26 | Wapp = []; |
---|
| 27 | else |
---|
| 28 | if INIT == 0 |
---|
| 29 | V = rand(P,r); |
---|
| 30 | else |
---|
| 31 | [Ui,Si,Vi] = svd(INIT); |
---|
| 32 | V = (Si(1:r,:)*Vi')'; |
---|
| 33 | end |
---|
| 34 | |
---|
| 35 | Wm = W.*INC; |
---|
| 36 | |
---|
| 37 | old_error = 999999999.0; |
---|
| 38 | error = 77777777.0; |
---|
| 39 | j = 1; |
---|
| 40 | |
---|
| 41 | while (j <= MAXSTEPS) & (abs(error - old_error) > 1e-10) |
---|
| 42 | j = j+1; |
---|
| 43 | |
---|
| 44 | [U,s] = whole_invert(V, W, INC, r); |
---|
| 45 | stable = s & stable; |
---|
| 46 | [V,s] = whole_invert(U, W', INC', r); |
---|
| 47 | stable = s & stable; |
---|
| 48 | |
---|
| 49 | Wapp = U*V'; |
---|
| 50 | old_error = error; |
---|
| 51 | error = sum(sum((Wm - Wapp.*INC).^2)); |
---|
| 52 | |
---|
| 53 | end |
---|
| 54 | |
---|
| 55 | if j == MAXSTEPS+1 |
---|
| 56 | [j, error, old_error]; |
---|
| 57 | error = - error; |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | end % if on conditioning of matrix. |
---|