Implementación SVD incremental en MATLAB

¿Hay alguna biblioteca / caja de herramientas que tenga implementación de SVD incremental en MATLAB? He implementado este documento , es rápido pero no funciona bien. Intenté esto, pero en este también el error se propaga rápidamente (dentro de la actualización, el error de 5-10 puntos es alto).

Si. Christopher Baker ha implementado su método SVD incremental en un paquete MATLAB llamado IncPACK ( archivado en GitHub, dentro del proyecto imtsl ). Implementa métodos que se describen en su tesis de maestría . Una breve discusión de por qué el algoritmo de Brand tiende a acumular errores se puede encontrar en un artículo de 2012 de Baker, et al . Un método relacionado por Chahlaoui, et al analiza los límites de error en el subespacio singular izquierdo y los valores singulares.

Ya he mencionado estos puntos en los comentarios sobre la respuesta de Stephen, pero la pena repetir que los métodos tanto por Baker y por escala Chahlaoui como para un truncado Rank- k SVD de un m por n matriz. Para aproximaciones de rango bajo, el término m n k domina y, dependiendo de la variante del algoritmo, tiene una constante principal que generalmente está entre 8 y 12.$O(mnk + nk^{3})$$k$$m$$n$$mnk$

Al igual que la respuesta de Stephen, el algoritmo de Chahlaoui comienza con una factorización QR. La respuesta de Stephen trabajará para el cálculo de vectores singulares izquierda, pero uno SVD densa del matriz tendría complejidad superlinear de m y n antes de truncamiento (sería O ( m n 2 ) ), lo que probablemente reducen la eficiencia, sino ser más preciso.$R$$m$$n$$O(mn^{2})$

Para lo que vale, he implementado el algoritmo de Brand yo mismo, y es algo sensible a la tolerancia interna del producto utilizada para el truncamiento de rango. No he usado el paquete de Baker, pero creo que sería mejor, porque existen estimaciones de error para el algoritmo de Baker (o uno estrechamente relacionado) y no para el algoritmo de Brand, y porque la tolerancia de truncamiento de rango para el algoritmo de Baker está en valores singulares, no internos productos

Un método para calcular el svd de una matriz Xes factorizar primero X=QRusando la descomposición QR (para estabilidad, use pivotar, entonces esto está [Q,R,E] = qr(X,0)en Matlab), y luego calcular el svd de R. Si la matriz es muy rectangular en cualquiera de los dos, entonces el cálculo más costoso es la factorización QR.

Por lo tanto, si incrementa su matriz Xcon otra fila o columna (esto es lo que quiso decir, ¿verdad?), Puede actualizar la factorización QR con la qrinsertfunción de Matlab y luego volver a hacer el cálculo SVD de R.

Si tiene una matriz cuadrada grande, este método no sería tan útil, ya que volver a hacer el SVD de Rllevará mucho tiempo.

Aquí hay un método que puede manejar adiciones de columna: http://pcc.byu.edu/resources.html . Lo actualicé para manejar las adiciones de fila:

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;

Pruébalo con

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

Verá que los resultados U, S, V en ambos lados son iguales.

También la versión de Python,

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

Una alternativa al SVD incremental es el HAPOD de descomposición ortogonal adecuada aproximada jerárquica , de la cual se puede encontrar una implementación en github: http://git.io/hapod . El HAPOD tiene límites de error rigurosos y un caso especial es una variante incremental.