source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/missing-data/convergence.m @ 37

% Here is code -- I have added a few comments to help you.
% N is dimension of square matrix, M is rank we want.  alpha is number
% between 0 and 1 that selects entries that are don't care.
% Program generates random [0,1] NxN matrix A1.  Z is generated from another
% similar random matrix whose elements are then set to one where entries 
% exceed alpha, and zero otherwise.  One entries determine the don't cares.
% Hence, about 1-alpha entries are don't care.  A is set to A1 in the care
% states, zero elsewhere.  Initial random entries are chosen for the
% don't care states in the outer for loop in range -50,+50.
% 
% Inner for loop is really to prevent an infinite loop (crummy programming).
% The real test is the break whent he error change is small.  That's also a
% hack for now.
% 
% XX accumulates sets of all final values of don't care values.

N = 5; M = 2;  alpha = 0.9; NUMSTEPS = 1000;
A1 = rand(N,N);
Z = max(0,sign(rand(N,N) - alpha));
A = A1 - A1.*Z;
Zs = reshape(Z,1,N^2);
index = find(Zs == 1)
count = size(index,2);
X = zeros(N,N);
XX = [];
STEPS = [];

for i = 1:20

 Xs = reshape(X,1,N^2);
 Xf = Xs(index);
 Xf = 100*(rand(1,count)-0.5);
 Xs(index) = Xf;
 X = reshape(Xs,N,N);
% X holds random initial starting conditions for unknown values.
 XfA = [];

 for j = 1:NUMSTEPS
  B = A + X;
% B is the matrix
  [U,D,V] = svd(B);
  BA = V(:,1:M)*D(1:M,1:M)*(V(:,1:M))';
% BA is the best fit rank M matrix to B.
  X = BA.*Z;
% X is updated with missing values from this fit.

  Xs = reshape(X,1,N^2);
  Xfold =Xf;
  Xf = Xs(index);
  XfA = [XfA; Xf];
  Error = sqrt(sum((Xf-Xfold).^2));
% Error is change in magnitude of X.

  A - (BA - X)
  if  Error < 1e-10 
    break
  end
 end
 STEPS = [STEPS; j];
 XX = [XX;Xf];
end
% XX
STEPS
XX(find(STEPS ~= NUMSTEPS), :)