% updatesvd.m
%
% Jason Blevins
%
% $Id: updatesvd.m,v 1.4 2004/12/08 21:23:47 jrblevin Exp $
%
% Returns the revised SVD corresponding to the rank-one
% update U*S*V' + A*B'. Based on the method outlined in
% M. E. Brand, Incremental singular value decomposition of
% uncertain data with missing values, European Conference on
% Computer Vision (ECCV), 2350:707--720, 2002.
%
%------------------------------------------------------
% Usage example:
%
% % Let X be some random data matrix and find its "economy
% % size" SVD
% X=randn(15,10);
% [U,S,V]=svd(X,0);
%
% % Suppose we want to append a random row r' to X
% r = randn(10,1);
%
% % First add a row of zeros to X (equivalently add a row of
% % zeros to U).
% U = [U; zeros(1,10)];
% X0 = [X; zeros(1,10)];
% % X0 will be equal to U*S*V'
%
% % In this case, we want to find the SVD of
% %[X; r'] = X0 + A*B' where
% % A = (0,...,0,1)' and B = r
% A = zeros(16,1); A(16,1) = 1;
% B = r;
%
% % Call the program like this, we want to find the SVD of
% [U2,S2,V2,time] = updatesvd(U,S,V,A,B);
%
% % Check the result, norm should be near machine epsilon
% norm((X0 + A*B') - U2*S2*V2')
%------------------------------------------------------
function [U,S,V,time]=updatesvd(U,S,V,A,B);
[m_u, n_u] = size(U);
[m_s, n_s] = size(S);
[m_v, n_v] = size(V);
[m_a, n_a] = size(A);
[m_b, n_b] = size(B);
m = m_u; % Number of rows in X=U*S*V'
n = m_v; % Number of columns in X
r = n_u; % Rank of the SVD we are given ( r = n_u = n_v )
c = n_a; % Number of columns in A and B
%%% Step 1: Find components of A in terms of left singular vectors, U.
% We decompose A into two components: M = U'*A, the component of A in U and
% Ra, the component of A orthogonal to U. We also find P, an orthogonal basis
% of Ra. We can find these via a QR decomposition.
%
% [U,A] = Q*R = [U,P] * [I, M; 0, Ra]
%
% The dimensions are as follows:
% U = [m_u, n_u]
% A = [m_a, n_a]
%
% These steps come from Appendix A and Equation (1.1) in Brand
time = clock;
[Q,R] = qr([U,A]);
%Q,R
[m_q, n_q] = size(Q);
[m_r, n_r] = size(R);
% If Q and U have the same number of columns, A lies completely
% within U and there is no orthogonal component and thus we don't need
% a basis P.
%
% Otherwise, the first r = n_u columns of Q will be the same as U, and the
% remaining m_q - r columns will be the orthogonal basis P.
if (n_q == r)
P = [];
Ra = [];
dimRa = 0;
else
P = Q(:, r+1:n_q);
Ra = R(r+1:m_r, r+1:n_r);
dimRa = m_r - r;
end
% I encounter round-off error when I use the result from the QR, why?
M = R(1:r, r+1:r+c);
%M = U'*A;
%%% Step 2: Find components of B in terms of right singular vectors, V.
% This step is very similar to step 1. We find two components of V:
% N, the component of B in V and Rb, the component of B orthogonal to V.
% Then we find a basis Q for Ra.
[Q,R] = qr([V,B]);
[m_q, n_q] = size(Q);
[m_r, n_r] = size(R);
if (n_q == r)
Q = [];
Rb = [];
dimRb = 0;
else
Q = Q(:, r+1:n_q); % Sorry for confusion, but yes, I have used Q for two
% different things. From now on, Q is the basis for Rb.
Rb = R(r+1:m_r, r+1:n_r);
dimRb = m_r - r;
end
% I encounter round-off error when I use the result from the QR, why?
%r,c
%N = R(1:r, r+1:r+c);
N = V'*B;
%%% Step 3: Construct the temporary new S and rediagonalize it
% Equation (1.4) in Brand
Saug = [M; Ra] * [N; Rb]';
S = [S, zeros(r, dimRb); zeros(dimRa, r), zeros(dimRa, dimRb)] + Saug;
% Now we rediagonalize via an SVD
[U2,S2,V2] = svd(S,0);
% And compute the new revised U, S, and V as in equation (1.5)
U = [U,P]*U2;
S = S2;
V = [V,Q]*V2;
time = etime(clock,time);