Throughout the question, I will use Python notation.
Suppose I have a matrix A of shape (p, nb) and I create a sliding window, taking the submatrix of p rows and n columns Am = A[:, m : m + n].
Now I want to compute it singular value decomposition (SVD):
U_m, S_m, Vh_m = svd(Am) = svd(A[:, m:m+n]
and then go to the next window m+1 and compute its SVD:
U_m1, S_m1, Vh_m1 = svd(Am1) = svd(A[:, m+1:m+1+n]
Computing full SVD from scratch has the complexity of o(min(m*n**2, n*m**2)).
I want to compute the SVD of the m+1 window using the SVD of the m window without computing the full SVD from scratch, so it will be more efficient (with less complexity) than doing a full SVD from scratch. Also I prefer it won't resort to low-rank approximation but will assume that rank(Am) = min(p, n).
U_m1, S_m1, Vh_m1 = sliding_svd(U_m, S_m, Vh_m, A[:,m+1:m+1+n], A[:,m],A[:,m+n])
A similar problem is called incremental SVD. To my problem, I call sliding SVD or moving SVD or sliding window SVD.
I asked a similar question in Math Stack Exchange: https://math.stackexchange.com/questions/5046319/sliding-singular-value-decomposition-sliding-svd-moving-svd
I am looking for code or a paper that can be implemented in Python using NumPy and SciPy to solve this problem, or at least some guiding (e.g. related papers that deal with a similar problem).
I'm going to assume square matrices, I leave the extension to non-square matrices to you (it may require further experimentation followed by a proof, but the general approach remains the same.)
The comments by a user in the post you made in Mathematics stack exchange pretty much point you at the right answer.
Lets call A0 = A[:,m:m+n] and A1 = A[:, m+1:m+1+n] following your nomenclature.
Furthermore, to simplify things, let L(A) be the leftmost column of matrix A and R(A) be the rightmost column of matrix A.
Then A1 = A0 P + uv^T for some suitable permutation matrix P and vectors u and v. In fact, P should be obvious (rotates all columns to the left by 1, and R(AP) = L(A), i.e. a cyclic rotation of columns.) The values of u and v are easy to work out, they are u = R(A1) - L(A) and v=[0,...,0,1]^T, that is v is a vector of all zeros save for 1 in the last entry.
Now, as pointed out in the comments in Math stack exchange, if A0=USV^T then A0P=US(V^TP) where V^TP remains orthogonal, therefore is the SVD of A0P.
Putting this altogether, as a coherent python example:
import numpy as np
def svd_rank1_update(U, S, Vt, u, v):
Su = U.T @ u
Sv = Vt @ v.T
pu = u - U @ Su
pv = v - Vt.T @ Sv
norm_pu = np.linalg.norm(pu)
norm_pv = np.linalg.norm(pv)
if norm_pu > 1e-10:
u_hat = pu / norm_pu
U_aug = np.column_stack([U, u_hat])
else:
U_aug = U
if norm_pv > 1e-10:
v_hat = pv / norm_pv
V_aug = np.column_stack([Vt.T, v_hat])
else:
V_aug = Vt.T
k = len(S)
S_mat = np.diag(S)
top_left = S_mat + np.outer(Su, Sv)
if norm_pu > 1e-10:
top = np.column_stack([top_left, norm_pv * Su])
else:
top = top_left
if norm_pv > 1e-10:
bottom = np.append(norm_pu * Sv, norm_pu * norm_pv)
small_matrix = np.vstack([top, bottom])
else:
small_matrix = top
U_small, S_new, Vt_small = np.linalg.svd(small_matrix, full_matrices=False)
U_new = U_aug @ U_small
V_new = V_aug @ Vt_small.T
Vt_new = V_new.T
return U_new, S_new, Vt_new
# Step 1: Create A and extract A0, A1
n = 5
m = 3
A = np.random.randn(n, m + n + 1) # Extra +1 for safe slicing
A0 = A[:, m:m+n] # A0 shape (n x n)
A1 = A[:, m+1:m+1+n] # A1 is next "sliding window" of size (n x n)
# Step 2: Create cyclic left permutation matrix P
P = np.roll(np.eye(n), -1, axis=1) # left shift columns
# Step 3: Compute A0 P
A0P = A0 @ P
# Step 4: Compute u and v
v = np.zeros(n)
v[-1] = 1 # [0, 0, ..., 0, 1]^T
u = A1[:, -1] - A0[:, 0] # u = R(A1) - L(A0)
# Step 5: Compute SVD of A0 and use it for A0 P
U, S, Vt = np.linalg.svd(A0, full_matrices=False)
Vt = Vt @ P # Apply the column permutation to get V^T P (still orthogonal)
# Step 6: Rank-1 update to get A1 SVD
U1, S1, Vt1 = svd_rank1_update(U, S, Vt, u, v)
# Step 7: Compare to actual SVD of A1
U_ref, S_ref, Vt_ref = np.linalg.svd(A1, full_matrices=False)
def align_signs(U1, U2):
signs = np.sign(np.sum(U1 * U2, axis=0))
return U1 * signs
U1_aligned = align_signs(U1, U_ref)
Vt1_aligned = align_signs(Vt1.T, Vt_ref.T).T
print("||U1_aligned - U_ref|| =", np.linalg.norm(U1_aligned - U_ref))
print("||Vt1_aligned - Vt_ref|| =", np.linalg.norm(Vt1_aligned - Vt_ref))
print("||S1 - S_ref|| =", np.linalg.norm(S1 - S_ref))
A1_rebuilt = U1 @ np.diag(S1) @ Vt1
err = np.linalg.norm(A1 - A1_rebuilt)
print("||A1 - A1_rebuilt|| =", err)
Which outputs something like:
||U1_aligned - U_ref|| = 1.4339826964875264e-15
||Vt1_aligned - Vt_ref|| = 1.2477860054076786e-15
||S1 - S_ref|| = 4.447825579847493e-15
||A1 - A1_rebuilt|| = 4.406110208978076e-15
SVD computation is typically O(n^3) for square matrices of size n whereas this uses the Brand method described in [1] which is O(n^2).
[1] Matthew Brand, Fast low-rank modifications of the thin singular value decomposition, Linear Algebra and its Applications, Volume 415, Issue 1, 1 June 2006, Pages 20–30. https://doi.org/10.1016/j.laa.2005.07.021