Implementação incremental de SVD no MATLAB

Existe alguma biblioteca / caixa de ferramentas que tenha implementação de SVD incremental no MATLAB. Eu implementei este documento , é rápido, mas não funciona bem. Eu tentei isso, mas neste erro também se propaga rapidamente (dentro de atualização de 5 a 10 pontos, o erro é alto).

Sim. Christopher Baker implementou seu método incremental de SVD em um pacote MATLAB chamado IncPACK ( arquivado no GitHub, dentro do projeto imtsl ). Implementa métodos descritos na tese de mestrado . Uma breve discussão sobre por que o algoritmo de Brand tende a acumular erros pode ser encontrada em um artigo de 2012 de Baker et al . Um método relacionado de Chahlaoui et al. Discute limites de erro no subespaço singular esquerdo e os valores singulares.

Eu já mencionei esses pontos nos comentários sobre a resposta de Stephen, mas vale a pena repetir que os métodos por tanto Baker e por escala Chahlaoui como para um truncada rank SVD de um por matriz . Para aproximações de baixa classificação, o termo domina e, dependendo da variante do algoritmo, possui uma constante líder que geralmente está entre 8 e 12.k m n m n $O(mnk + nk^{3})$$k$$m$$n$$mnk$

Como a resposta de Stephen, o algoritmo de Chahlaoui começa com uma fatoração QR. A resposta de Stephen irá funcionar para calcular os vectores singulares esquerda, mas uma densa SVD do matriz teria complexidade superlinear em e antes da truncagem (que seria ), o que provavelmente iria reduzir a eficiência, mas seja mais preciso.m n O ( m n 2$R$$m$$n$$O(mn^{2})$

Pelo que vale a pena, eu mesmo implementei o algoritmo de Brand e é um pouco sensível à tolerância interna do produto usada para o truncamento de classificação. Eu não usei o pacote de Baker, mas acredito que seria melhor, porque existem estimativas de erro para o algoritmo de Baker (ou um relacionado a ele) e não para o algoritmo de Brand, e porque a tolerância de truncamento de classificação para o algoritmo de Baker está em valores singulares, não internos produtos.

Um método para calcular o svd de uma matriz Xé o primeiro fator X=QRusando a decomposição QR (para estabilidade, use pivotamento, portanto, isso está [Q,R,E] = qr(X,0)no Matlab) e depois calcule o svd de R. Se a matriz for muito retangular em qualquer um deles, o cálculo mais caro é a fatoração QR.

Portanto, se você incrementar sua matriz Xcom outra linha ou coluna (é isso que você quis dizer, certo?), Basta atualizar a fatoração QR com a qrinsertfunção Matlab e refazer o cálculo SVD de R.

Se você tiver uma matriz quadrada grande, esse método não seria tão útil, pois refazer o SVD de Rconsumirá muito tempo.

Aqui está um método que pode lidar com adições de coluna: http://pcc.byu.edu/resources.html . Eu atualizei para lidar com adições de linha:

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;

Teste com

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

Você verá os resultados U, S, V de ambos os lados.

Também a versão 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

Uma alternativa ao SVD Incremental é o HAPOD de Decomposição Ortogonal Aproximada Aproximada Hierárquica , cuja implementação pode ser encontrada no github: http://git.io/hapod . O HAPOD possui limites de erro rigorosos e um caso especial é uma variante incremental.