function [Up,Sp,Vp] = svd_update( U, S, V, A, B, force_orth )
% function [Up,Sp,Vp] = svd_update( U, S, V, A, B, force_orth )
%
% Given the SVD of
%
% X = U*S*V'
%
% update it to be the SVD of
%
% X + AB' = Up*Sp*Vp'
%
% Depending on A,B there may be considerable structure that could
% be exploited, but which this code does not.
%
% The subspace rotations involved may not preserve orthogonality due
% to numerical round-off errors. To compensate, you can set the
% "force_orth" flag, which will force orthogonality via a QR plus
% another SVD. In a long loop, you may want to force orthogonality
% every so often.
%
% See Matthew Brand, "Fast low-rank modifications of the thin
% singular value decomposition".
%
% D. Wingate 8/17/2007
%
current_rank = size( U, 2 );
% P is an orthogonal basis of the column-space
% of (I-UU')A, which is the component of "A" that is
% orthogonal to U.
m = U' * A;
p = A - U*m;
P = orth( p );
% p may not have full rank. If not, P will be too small. Pad
% with zeros.
P = [ P zeros(size(P,1), size(p,2)-size(P,2)) ];
Ra = P' * p;
% Q is an orthogonal basis of the column-space
% of (I-VV')b.
n = V' * B;
q = B - V*n;
Q = orth( q );
% q may not have full rank. If not, Q will be too small. Pad
% with zeros.
Q = [ Q zeros(size(Q,1), size(q,2)-size(Q,2)) ];
Rb = Q' * q;
%
% Diagonalize K, maintaining rank
%
z = zeros( size(m) );
z2 = zeros( size(m,2), size(m,2) );
K = [ S z ; z' z2 ] + [ m; Ra ] * [ n; Rb ]';
[tUp,tSp,tVp] = svds( K, current_rank );
%
% Now update our matrices!
%
Sp = tSp;
Up = [ U P ] * tUp;
Vp = [ V Q ] * tVp;
% The above rotations may not preserve orthogonality, so we explicitly
% deal with that via a QR plus another SVD. In a long loop, you may
% want to force orthogonality every so often.
if ( force_orth )
[UQ,UR] = qr( Up, 0 );
[VQ,VR] = qr( Vp, 0 );
[tUp,tSp,tVp] = svds( UR * Sp * VR', current_rank );
Up = UQ * tUp;
Vp = VQ * tVp;
Sp = tSp;
end;
return;