Inkrementelle SVD-Implementierung in MATLAB

Gibt es eine Bibliothek / Toolbox, in der inkrementelle SVD in MATLAB implementiert ist? Ich habe dieses Papier implementiert , es ist schnell, funktioniert aber nicht gut. Ich habe dies versucht, aber auch hier breitet sich der Fehler schnell aus (beim Aktualisieren von 5-10 Punkten ist der Fehler hoch).

Ja. Christopher Baker hat seine inkrementelle SVD-Methode in einem MATLAB-Paket namens IncPACK implementiert ( archiviert auf GitHub im imtsl-Projekt ). Es implementiert Methoden, die in seiner Masterarbeit beschrieben sind . Eine kurze Diskussion darüber, warum der Algorithmus von Brand dazu neigt, Fehler zu akkumulieren, findet sich in einem Artikel von Baker et al . Eine verwandte Methode von Chahlaoui et al. Diskutiert Fehlergrenzen im linken singulären Unterraum und die singulären Werte.

Ich habe bereits diese Punkte in den Kommentaren auf Stephen Antwort erwähnt, aber es trägt zu wiederholen , dass die Verfahren sowohl von Baker und durch Chahlaoui Skala als für ein verkürztes Rang SVD einen durch Matrix. Bei niedrigrangigen Näherungen dominiert der Term und hat je nach Algorithmusvariante eine führende Konstante, die normalerweise zwischen 8 und 12 liegt.$O(mnk + nk^{3})$$k$$m$$n$$mnk$

Wie Stephens Antwort beginnt Chahlaouis Algorithmus mit einer QR-Faktorisierung. Stephens Antwort wird für die Berechnung linker Singularvektoren funktionieren, aber eine dichte SVD der Matrix hätte vor dem Abschneiden eine superlineare Komplexität in und (es wäre ), was wahrscheinlich die Effizienz verringern würde, aber genauer sein.$R$$m$$n$$O(mn^{2})$

Für das, was es wert ist, habe ich den Brand-Algorithmus selbst implementiert und er reagiert etwas empfindlich auf die innere Produkttoleranz, die für die Rangabschneidung verwendet wird. Ich habe das Baker-Paket nicht verwendet, aber ich glaube, es wäre besser, da Fehlerschätzungen für den Baker-Algorithmus (oder einen eng verwandten) und nicht für den Brand-Algorithmus existieren und weil die Rangabschneidungstoleranz für den Baker-Algorithmus auf singulären Werten beruht, nicht auf inneren Produkte.

Eine Methode zum Berechnen der svd einer Matrix Xbesteht darin, zuerst X=QRdie QR-Zerlegung zu verwenden (aus Stabilitätsgründen [Q,R,E] = qr(X,0)Pivoting verwenden, also in Matlab) und dann die svd von zu berechnen R. Wenn die Matrix in beiden Fällen sehr rechteckig ist, ist die QR-Faktorisierung die teuerste Berechnung.

Wenn Sie also Ihre Matrix Xum eine andere Zeile oder Spalte erhöhen (das haben Sie gemeint, oder?), Können Sie einfach die QR-Faktorisierung mit der Matlab- qrinsertFunktion aktualisieren und dann die SVD-Berechnung von erneut durchführen R.

Wenn Sie eine große quadratische Matrix haben, ist diese Methode nicht so nützlich, da das erneute Ausführen der SVD von Rzeitaufwändig ist.

Hier ist eine Methode, die das Hinzufügen von Spalten verarbeiten kann: http://pcc.byu.edu/resources.html . Ich habe es aktualisiert, um Zeilenergänzungen zu verarbeiten:

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;

Testen Sie es mit

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

Sie werden sehen, dass die U-, S- und V-Ergebnisse auf beiden Seiten gleich sind.

Auch die 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

Eine Alternative zur inkrementellen SVD ist das HAPOD " Hierarchical Approximate Proper Orthogonal Decomposition " , dessen Implementierung auf github zu finden ist: http://git.io/hapod . Der HAPOD hat strenge Fehlergrenzen und ein Sonderfall ist eine inkrementelle Variante.