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