本文介绍一种基于秩-1 更新的滑动窗口 svd 算法，利用前一窗口的 svd 结果，以 $o(n^2)$ 时间复杂度增量更新下一窗口的完整 svd，避免重复执行 $o(n^3)$ 的全量分解。
本文介绍一种基于秩-1 更新的滑动窗口 svd 算法，利用前一窗口的 svd 结果，以 $o(n^2)$ 时间复杂度增量更新下一窗口的完整 svd，避免重复执行 $o(n^3)$ 的全量分解。
在时间序列分析、流式矩阵处理、动态主成分分析（PCA）等场景中，常需对按固定长度滑动的子矩阵持续进行奇异值分解（SVD）。若对每个窗口独立调用 np.linalg.svd，当窗口数量大、维度高时，计算开销急剧上升——标准 SVD 复杂度为 $O(\min(p n^2, n p^2))$，对 $n \times n$ 方阵即为 $O(n^3)$。而“滑动 SVD”（Sliding SVD）的核心思想是：当前窗口 $A_{m+1} = A[:, m+1:m+1+n]$ 与前一窗口 $A_m = A[:, m:m+n]$ 仅相差一列（删左列、增右列），因此可建模为对已知 SVD 的矩阵施加低秩校正，进而复用原有分解结构。
关键数学洞察在于：设 $A_0 = A[:, m:m+n]$, $A_1 = A[:, m+1:m+1+n]$，则存在一个列循环置换矩阵 $P$（将所有列左移一位，末列补至首列），使得
$$
A_1 = A_0 P + u v^\top,
$$
其中
- $P = \operatorname{roll}(I_n, -1, \text{axis}=1)$ 是 $n \times n$ 置换矩阵；
- $u = \text{last_col}(A_1) - \text{first_col}(A_0) \in \mathbb{R}^n$；
- $v = [0,\dots,0,1]^\top \in \mathbb{R}^n$（单位向量，仅末位为 1）。
由于 $A_0 = U S V^\top$，则 $A_0 P = U S (V^\top P)$，且 $V^\top P$ 仍正交，故 $(U, S, V^\top P)$ 即为 $A_0 P$ 的有效 SVD。此时问题转化为：对已有 SVD 的矩阵 $A_0 P$ 执行秩-1 更新 $+u v^\top$，快速获得 $A_1$ 的新 SVD。这正是 Brand 提出的经典算法（见参考文献［1］），其核心是构造一个仅含 $O(n)$ 维的“小矩阵”，对其做 SVD 后映射回原空间，从而将总复杂度降至 $O(n^2)$。
以下为完整可运行的 NumPy 实现（支持方阵；非方阵可依类似思路扩展行/列空间）：
import numpy as np
def svd_rank1_update(U, S, Vt, u, v):
"""对 A = U @ diag(S) @ Vt 执行 A 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)
# 拼接小矩阵：尺寸最多为 (k+1) x (k+1)
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
# 对小矩阵做 SVD
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
# 示例：滑动窗口 SVD 计算
np.random.seed(42)
n = 5
m_start = 3
A = np.random.randn(n, m_start + n + 1) # 预留边界
A0 = A[:, m_start : m_start + n] # shape (n, n)
A1 = A[:, m_start + 1 : m_start + 1 + n]
# 构造循环左移置换矩阵 P
P = np.roll(np.eye(n), -1, axis=1)
# 计算 A0P 的 SVD（重用 A0 的 U,S,Vt）
U0, S0, Vt0 = np.linalg.svd(A0, full_matrices=False)
Vt0P = Vt0 @ P # V^T P 仍正交
# 定义 rank-1 更新项
v = np.zeros(n); v[-1] = 1
u = A1[:, -1] - A0[:, 0]
# 增量更新得到 A1 的 SVD
U1, S1, Vt1 = svd_rank1_update(U0, S0, Vt0P, u, v)
# 验证精度（与直接 SVD 对比）
U_ref, S_ref, Vt_ref = np.linalg.svd(A1, full_matrices=False)
def align_signs(X, Y):
"""修正奇异向量符号歧义（SVD 不唯一）"""
signs = np.sign(np.sum(X * Y, axis=0))
return X * signs
U1_aligned = align_signs(U1, U_ref)
Vt1_aligned = align_signs(Vt1.T, Vt_ref.T).T
print("误差验证（应 <p>✅ <strong>注意事项与实践建议</strong>： </p>
-
数值稳定性：代码中使用
1e-10作为秩判定阈值，实际应用中可根据数据尺度和精度需求调整（如1e-12或基于np.finfo(float).eps * max(S)）； -
非方阵扩展：若 $A$ 为 $(p,n)$ 矩阵（$p \neq n$），需分别处理行空间与列空间的增/删逻辑，并维护不同维度的 $U$ 和 $V^\top$；推荐先对 $A_m$ 做 thin SVD（
full_matrices=False），再按 Brand 方法推导广义更新公式； - 内存优化：该方法天然适合流式场景——无需缓存全部历史窗口，仅需保留上一窗口的 $(U,S,V^\top)$ 及最新两列（用于计算 $u$）；
-
替代方案权衡：若允许低秩近似（如 $r \ll \min(p,n)$），可考虑
scikit-learn的IncrementalPCA或fbpca库；但本文方案严格保持满秩，适用于需精确特征值/向量的场景（如稳定性判据、精确子空间跟踪）。
综上，滑动窗口 SVD 并非简单重复计算，而是通过代数建模将滑动操作解耦为“置换 + 秩-1 更新”，再借由 Brand 算法实现高效重构。该策略在保证数学严谨性的同时，显著提升大规模滑动分析的可行性，是动态矩阵分解领域的重要工程实践范式。
