Is there any library/toolbox which has implementation of incremental SVD in MATLAB. I have implemented this paper, it is fast but does not work well. I tried this but in this also error propagates fast (within updating 5-10 points error is high).
Add a comment
|
4 Answers
$\begingroup$
$\endgroup$
6
Yes. Christopher Baker has implemented his incremental SVD method in a MATLAB package called IncPACK (archived on GitHub, within the imtsl project). It implements methods that are described in his master's thesis. A brief discussion of why Brand's algorithm tends to accumulate error can be found in a 2012 paper by Baker, et al. A related method by Chahlaoui, et al discusses error bounds on the left singular subspace and the singular values.
I've already mentioned these points in the comments on Stephen's answer, but it bears repeating that the methods by both Baker and by Chahlaoui scale as $O(mnk + nk^{3})$ for a truncated rank-$k$ SVD of an $m$ by $n$ matrix. For low-rank approximations, the $mnk$ term dominates and, depending on the algorithm variant, has a leading constant that is usually between 8 and 12.
Like Stephen's answer, Chahlaoui's algorithm starts with a QR factorization. Stephen's answer will work for calculating left singular vectors, but a dense SVD of the $R$ matrix would have superlinear complexity in $m$ and $n$ prior to truncation (it would be $O(mn^{2})$), which would probably reduce efficiency, but be more accurate.
For what it's worth, I've implemented Brand's algorithm myself, and it's somewhat sensitive to the inner product tolerance used for rank truncation. I haven't used Baker's package, but I believe it would be better, because error estimates exist for Baker's algorithm (or one closely related) and not Brand's algorithm, and because the rank truncation tolerance for Baker's algorithm is on singular values, not inner products.
answered Mar 30, 2015 at 7:27
-
$\begingroup$ I have checked IncPACK package, although it has seqkl_update function it doesn't look like accepting any parameters for new rows and columns. Also from paper abstract(may not be correct I have to read it all) it looks like it is a multipass approach which they call it incremental. $\endgroup$– ParagCommented Mar 30, 2015 at 10:31
-
$\begingroup$ Baker discusses both single-pass and multipass approaches. The
seqklfunction looks to be the main function, and has options for single and multiple passes. A single pass is given byseqkl_stdpass, which callsseqkl_update, so you probably would want to useseqklfor an initial factorization, followed by calls toseqkl_updatefor column updates. $\endgroup$ Commented Mar 30, 2015 at 18:36 -
$\begingroup$ Yeah, till now what I found is it has column update only, new data Ai is stored in U(1:m,o:op) that is comment from seqkl_update file. But what about row update? $\endgroup$– ParagCommented Mar 30, 2015 at 18:44
-
$\begingroup$ @p.j From what I've read, most of the literature focuses on column updates. If it's possible, you could instead calculate the SVD of the transpose of your matrix; however, I recognize that having to transpose your data may not be an option. I presume that the algorithm could be re-expressed in terms of row updates, but that might require a bit of work. $\endgroup$ Commented Mar 30, 2015 at 23:32
-
$\begingroup$ The first two links are dead. $\endgroup$ Commented Apr 29, 2016 at 15:16
$\begingroup$
$\endgroup$
5
One method to compute the svd of a matrix X is to first factor X=QR using the QR decomposition (for stability, use pivoting, so this is [Q,R,E] = qr(X,0) in Matlab), and then compute the svd of R. If the matrix is very rectangular in either, then the most expensive computation is the QR factorization.
Thus if you increment your matrix X with another row or column (this is what you meant, right?), you can just update the QR factorization with Matlab's qrinsert function, and then re-do the SVD calculation of R.
If you have a large square matrix, this method would not be as useful, since re-doing the SVD of R will be time-consuming.
-
$\begingroup$ This is only fast for a low-rank svd, right? $\endgroup$– dranxoCommented Mar 30, 2015 at 3:53
-
$\begingroup$ @Stephen Yes that is what I want, add a column and a row(that is adding one data point for me). I think methods which I have tried also do QR initially and then update. My major concern is error. What I have noticed is that svd for old points which should not change much, after 5-6 incremental updates gets corrupted with huge error, same is the case with new points. $\endgroup$– ParagCommented Mar 30, 2015 at 6:59
-
$\begingroup$ The method you're proposing isn't quite an incremental SVD. This method vaguely resembles the incremental SVD algorithm of Chahlaoui, but recomputing the SVD of $R$ by the usual dense algorithm is going to kill efficiency. The work by Chahlaoui and the work by Baker, et al. are more efficient methods and scale linearly in the size of the matrix
Xand the rank $k$ of the truncated SVD desired. $\endgroup$ Commented Mar 30, 2015 at 7:13 -
$\begingroup$ @GeoffOxberry, that is what I said: if you have a square matrix, this is not efficient. I mention it because often you have a skinny/fat matrix, and this algorithm can be implemented trivially using existing functions in Matlab $\endgroup$– StephenCommented Mar 30, 2015 at 20:38
-
$\begingroup$ @dranxo, well, if you are very rectangular, it is relatively fast. It's not as fast as a true incremental SVD, but if you are very rectangular, the cost of the inner SVD is trivial, and you can balance this with the ease of implementation $\endgroup$– StephenCommented Mar 30, 2015 at 20:39
$\begingroup$
$\endgroup$
Here is a method that can handle column additions: http://pcc.byu.edu/resources.html. I updated it to handle row additions:
function [Up1,Sp,Vp1] = addblock_svd_update2( Uarg, Sarg, Varg, Aarg, force_orth )
U = Varg;
V = Uarg;
S = Sarg;
A = Aarg';
current_rank = size( U, 2 );
m = U' * A;
p = A - U*m;
P = orth( p );
P = [ P zeros(size(P,1), size(p,2)-size(P,2)) ];
Ra = P' * p;
z = zeros( size(m) );
K = [ S m ; z' Ra ];
[tUp,tSp,tVp] = svds( K, current_rank );
Sp = tSp;
Up = [ U P ] * tUp;
Vp = V * tVp( 1:current_rank, : );
Vp = [ Vp ; tVp( current_rank+1:size(tVp,1), : ) ];
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;
Up1 = Vp;
Vp1 = Up;
return;
Test it with
X = [[ 2.180116 2.493767 -0.047867;
-1.562426 2.292670 0.139761;
0.919099 -0.887082 -1.197149;
0.333190 -0.632542 -0.013330]];
A = [1 1 1];
X2 = [X; A];
[U,S,V] = svds(X);
[Up,Sp,Vp] = addblock_svd_update2(U, S, V, A, true);
Up
Sp
Vp
[U2,S2,V2] = svds(X2);
U2
S2
V2
You will see U,S,V results on both sides are the same.
Also the Python version,
import numpy as np
import scipy.linalg as lin
def addblock_svd_update( Uarg, Sarg, Varg, Aarg, force_orth = False):
U = Varg
V = Uarg
S = np.eye(len(Sarg),len(Sarg))*Sarg
A = Aarg.T
current_rank = U.shape[1]
m = np.dot(U.T,A)
p = A - np.dot(U,m)
P = lin.orth(p)
Ra = np.dot(P.T,p)
z = np.zeros(m.shape)
K = np.vstack(( np.hstack((S,m)), np.hstack((z.T,Ra)) ))
tUp,tSp,tVp = lin.svd(K);
tUp = tUp[:,:current_rank]
tSp = np.diag(tSp[:current_rank])
tVp = tVp[:,:current_rank]
Sp = tSp
Up = np.dot(np.hstack((U,P)),tUp)
Vp = np.dot(V,tVp[:current_rank,:])
Vp = np.vstack((Vp, tVp[current_rank:tVp.shape[0], :]))
if force_orth:
UQ,UR = lin.qr(Up,mode='economic')
VQ,VR = lin.qr(Vp,mode='economic')
tUp,tSp,tVp = lin.svd( np.dot(np.dot(UR,Sp),VR.T));
tSp = np.diag(tSp)
Up = np.dot(UQ,tUp)
Vp = np.dot(VQ,tVp)
Sp = tSp;
Up1 = Vp;
Vp1 = Up;
return Up1,Sp,Vp1
$\begingroup$
$\endgroup$
An alternative to the Incremental SVD is the Hierarchical Approximate Proper Orthogonal Decomposition HAPOD, of which an implementation can be found on github: http://git.io/hapod . The HAPOD has rigorous error bounds and a special case is an incremental variant.
