# SVD of Matrix w/ updated columns

**URL:** <https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048>\
**Category:** Numerics\
**Created:** [February 21, 2019, 1:36pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048 "2019-02-21T13:36:58Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [February 21, 2019, 1:36pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/1 "2019-02-21T13:36:58Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [February 21, 2019, 5:17pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/2 "2019-02-21T17:17:09Z")

</div>

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

```julia
# 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)

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 21, 2019, 6:08pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/3 "2019-02-21T18:08:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [February 21, 2019, 6:33pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/4 "2019-02-21T18:33:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [February 21, 2019, 9:07pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/5 "2019-02-21T21:07:31Z")

</div>

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://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/2678/updatable-svd-implementation-in-python-c-or-fortran)

[https://scicomp.stackexchange.com/questions/19253/incremental-svd-implementation-in-matlab](https://scicomp.stackexchange.com/questions/19253/incremental-svd-implementation-in-matlab)

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [February 21, 2019, 11:14pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/6 "2019-02-21T23:14:28Z")

</div>

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

> [@DrPapa](#):
>
> Would this be considered a rank-1 update?

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

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [February 22, 2019, 12:39pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/7 "2019-02-22T12:39:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [February 22, 2019, 5:26pm UTC](https://discourse.julialang.org/t/svd-of-matrix-w-updated-columns/21048/8 "2019-02-22T17:26:20Z")

</div>

Here is a working example of the SVD Rank-1 Update process outline in [[Brand 2006]](https://pdfs.semanticscholar.org/4e8f/82b0741c2151d36f2201fc11b0b148beab60.pdf). I also referenced [this](https://github.com/sioan/auto_xtc_hdf5_converter/blob/v0.2/user_created_analysis_functions/SVD_lib/manual_svd_update.py) 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**

```julia
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:

```julia
## 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

```julia
@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.