Question: 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).

Answer: 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