Joshua Carmichael

You are doing something like C * ((C'*C) \ C').
The latter part is (C'*C)\C' which is the normal eqn's form of just C\I. Isn't that prefered for accuracy? It is slower than the (C'*C \ C') form, a bit.

That's odd. Yes, I did download and review the version with "inv" in it. The newer version looks much better. For this updated version, I would suggest writing

projv = (I-P)*W(:,k) ;

and remove the Pc and v variables, in the interest of speed and simplicity, but that's a minor point. (For that matter, the variable P can be removed too).

Looks great.

John is right ... no one should ever multiply by the inverse. You should replace the line

    P=(A*U)*inv((A*U)'*(A*U))*(A*U)';

with

    P=(A*U) * (((A*U)'*(A*U)) \ (A*U)') ;

Multiplying by the inverse is shockingly common in files posted at the File Exchange. It should never be done. Even if A is diagonal, A\b is better than inv(A)*b.

Also, if you'd like to work with sparse W and A, you should replace the line:

   I=eye(size(A));

with

 if (issparse (A))
    I = speye (size (A)) ;
 else
    I = eye (size (A)) ;
 end

My guess is that if W and A are sparse then U will probably be all nonzero anyway, though, so the above probably isn't helpful. If you pass sparse W and A, then U will be returned full, with the current version of the code.

What is the statement v=W(:,k) doing inside the for loop? You don't use v later on at all, but you do use W(:,k) itself.

I would also recommend replacing the "disp" with "warning".