SVD of Matrix w/ updated columns

Specific Domains Numerics
1

I have a M\times N matrix A where N>>M. I need to compute the SVD of sub-matrices using a sliding window over the columns, i.e. svd(A[:,1:k]), svd(A[:,2:k+1]), svd(A[:,3:k+2])…

Does anybody know of an efficient way to update an SVD when the 1st column is removed and a new last column is added? Would this be considered a rank-1 update?

2

It’s a rank-1 update and a permutation. Here is some explanatory code:

# Whole matrix
A = rand(3, 3)

# First 2 cols
B = A[:, 1:2]
s = svd(B)

# The matrix that replaces the first col by the third
C = zeros(size(B))
C[:, 1] .= -A[:, 1] + A[:, 3]

# Standard basis vector with 1 at the index of the vector to be replaced in B
# Length of e is the number of columns of B
e = zeros(2); e[1] = 1;

# Low rank representation of C
C == (-A[:, 1] + A[:, 3]) * e'

# Low rank update
B += C

# New SVD
s = svd(B)

# But we want the third col to be after the second col, not before
# Easy, just permute the columns of Vt
new_Vt = s.Vt[:, [2, 1]]
s = SVD(s.U, s.S, new_Vt)

# Sanity check
s2 = svd(A[:,2:3])
isapprox(s.U, s2.U)
isapprox(s.Vt, s2.Vt)
isapprox(s.S, s2.S)
3

I don’t think it’s possible in general. You can easily update lu and qr, but not svd. In fact, Galois theory tells you that you can’t find the eigenvalues of a 5x5 matrix without some kind of (in principle infinite) iteration, so there should be no straightforward algorithm for svd updating. Of course there can be non-straightforward ones, I just don’t know any.

1 Like
4

Thanks for the example, but unless I am overlooking something I don’t think this is quite what I am looking for. My goal is to only compute the SVD once upfront and then compute the new SVD decomposition based on the already computed one directly w/ knowledge of the rank-1 update.

5

After some more digging I did find some references to an algorithm including a few Python & Matlab examples. I will take a crack and implementing one of these approaches.

https://stats.stackexchange.com/questions/177007/updating-svd-decomposition-after-adding-one-new-row-to-the-matrix

https://scicomp.stackexchange.com/questions/2678/updatable-svd-implementation-in-python-c-or-fortran

https://scicomp.stackexchange.com/questions/19253/incremental-svd-implementation-in-matlab

1 Like
6

I know I am not fully answering your question. I was mostly answering this part:

Once you have the low rank update mechanism, I was assuming that you can find resources online to take you the extra step :wink: Should have probably clarified.

7

:+1: Confirming that it was a rank-1 update did help me find some resources. Thank you.

1 Like
8

Here is a working example of the SVD Rank-1 Update process outline in [Brand 2006]. I also referenced this Python code w/ a few modifications. Note, I added an isapprox helper to compare results.

EDIT: This still needs some work. Although the original matrix can be recovered, the singular values are not correct. I think there is an issue with the normality constraint for the left/right singular vectors

using LinearAlgebra
import Base: ≈

function ≈(A::SVD{T},B::SVD{T}) where T<:Number
	A.U*Diagonal(A.S)*A.Vt≈B.U*Diagonal(B.S)*B.Vt
end

function svd_rank_one_update(Q::SVD{T}, a::AbstractArray{T,1}) where T<:Number
	U = Q.U
	S = Q.S
	Vt = Q.Vt
	r = length(S)

    b = zeros(r+1)
    b[end] = 1

    m = U'*a
    p = a-U*m
    Ra = norm(p)
    P = p/Ra

    Vt = vcat(Vt, zeros(1,size(Vt,2)))

    n = Vt'*b
    q = b - Vt*n
    Rb = norm(q)
    Q = q/Rb

	K = zeros(r+1, r+1)  # could use spzeros

 	K[1:end-1,end]  = m
 	K[1:end-1,1:end-1] = Diagonal(S)
 	K[end,end] = Ra

    U_P = hcat(U,P)
    V_Q = hcat(Vt,Q)

    eig_vals , eig_vec = eigen(K)

	U_P_dot_Up = U_P*eig_vec
	V_Q_dot_Vq = inv(eig_vec)*V_Q

 	return SVD(U_P_dot_Up, eig_vals, V_Q_dot_Vq)
end

Test:

## Generate data and compute SVD
r = 10
data = rand(r,1000000)
Q = svd(data)

## Augment data w/ a vector & compute SVD
a = rand(r)
data_aug = hcat(data,a)
Q_aug = svd(data_aug)

## Rank-1 update of original SVD
Q_rank1 = svd_rank_one_update(Q,a)
Q_aug ≈ Q_rank1  #true

Benchmark rank-1 update vs. concatenate then svd

@btime svd_rank_one_update($Q,$a) #100.994 ms (43 allocations: 267.05 MiB)
@btime svd(hcat($data,$a))        #504.592 ms (14 allocations: 228.89 MiB)

I am a relative noob to Julia, so any suggested improvements are welcome.

2 Likes