在大规模数据流和在线学习场景中，我们经常需要维护矩阵的低秩近似，同时适应数据的动态变化。本章深入探讨流式环境下的矩阵分解更新、自适应秩选择、在线补全算法，以及这些技术与现代深度学习中模型压缩的深刻联系。我们将特别关注时间-空间-精度的三方权衡，以及在非平稳环境中的理论保证。
考虑数据流设定：时刻$t$观察到向量$\mathbf{x}_t \in \mathbb{R}^n$或矩阵更新$\Delta\mathbf{A}_t$。目标是维护当前数据矩阵$\mathbf{A}_t$的秩$k$近似$\tilde{\mathbf{A}}_t$：
\[\mathbf{A}_t = \sum_{i=1}^t \mathbf{x}_i\mathbf{x}_i^T \quad \text{或} \quad \mathbf{A}_t = \mathbf{A}_0 + \sum_{i=1}^t \Delta\mathbf{A}_i\]核心挑战：
数学形式化：定义流式低秩近似问题为 \(\min_{\tilde{\mathbf{A}}_t: \text{rank}(\tilde{\mathbf{A}}_t) \leq k} \|\mathbf{A}_t - \tilde{\mathbf{A}}_t\|_F^2\)
在约束条件下：
实际应用场景：
定义累积regret： \(R_T = \sum_{t=1}^T \|\mathbf{A}_t - \tilde{\mathbf{A}}_t\|_F^2 - \sum_{t=1}^T \|\mathbf{A}_t - \mathbf{A}_t^{(k)}\|_F^2\)
其中$\mathbf{A}_t^{(k)}$是$\mathbf{A}_t$的最优秩$k$近似。
定理13.1（流式SVD的regret界）：存在算法使得 \(R_T \leq O(k\log T \cdot \sum_{j=k+1}^n \lambda_j(\mathbf{A}_T))\)
这表明额外误差仅对数增长于时间$T$。
证明要点：
详细证明思路： 设$\mathbf{E}_t = \tilde{\mathbf{A}}_t - \mathbf{A}_t^{(k)}$为时刻$t$的额外误差。关键观察是： \(\|\mathbf{E}_t\|_F^2 \leq \|\mathbf{E}_{t-1}\|_F^2 + 2\langle\mathbf{E}_{t-1}, \Delta_t\rangle + \|\Delta_t\|_F^2\)
其中$\Delta_t$是时刻$t$的更新。通过精心选择更新策略，可以控制交叉项$\langle\mathbf{E}_{t-1}, \Delta_t\rangle$。
更细致的性能度量：
相对误差： \(\text{RelErr}_t = \frac{\|\mathbf{A}_t - \tilde{\mathbf{A}}_t\|_F}{\|\mathbf{A}_t\|_F}\)
子空间距离（更敏感的度量）： \(d_{\text{sub}}(\mathbf{U}_t, \tilde{\mathbf{U}}_t) = \|\mathbf{U}_t\mathbf{U}_t^T - \tilde{\mathbf{U}}_t\tilde{\mathbf{U}}_t^T\|_F\)
投影误差（应用导向）： \(\text{ProjErr}_t = \mathbb{E}_{\mathbf{x}}[\|\mathbf{x} - \tilde{\mathbf{U}}_t\tilde{\mathbf{U}}_t^T\mathbf{x}\|^2]\)
Grassmannian距离（几何视角）： \(d_{\text{Grass}}(\mathbf{U}_t, \tilde{\mathbf{U}}_t) = \left(\sum_{i=1}^k \theta_i^2\right)^{1/2}\) 其中$\theta_i$是两个子空间之间的主角。
定理13.2（子空间追踪）：在温和条件下， \(d_{\text{sub}}(\mathbf{U}_t, \tilde{\mathbf{U}}_t) \leq \frac{C}{\lambda_k(\mathbf{A}_t) - \lambda_{k+1}(\mathbf{A}_t)} \cdot \text{err}_t\)
其中分母是特征值间隙，决定了子空间分离的难度。
实际含义：
高阶误差分析：
对于有限精度计算，总误差由三部分组成：
定理13.3（总误差界）：在IEEE双精度下， \(\|\mathbf{A}_t - \tilde{\mathbf{A}}_t\|_F \leq \epsilon_{\text{approx}} + O(\sqrt{k\log t})\sigma_{k+1} + O(t \cdot \epsilon_{\text{machine}})\|\mathbf{A}_t\|_F\)
这给出了精度、秩选择和运行时间的三方权衡。
实用性能基准：
| 度量 | 优秀 | 良好 | 需改进 |
|---|---|---|---|
| 相对误差 | < 1% | 1-5% | > 5% |
| 子空间距离 | < 0.1 | 0.1-0.3 | > 0.3 |
| 每样本时间 | < 10μs | 10-100μs | > 100μs |
| 内存开销 | < 2nk | 2-5nk | > 5nk |
算法13.1：通用增量低秩近似框架
输入：初始分解 U₀Σ₀V₀ᵀ，流数据 {xₜ}
维护：当前近似 UₜΣₜVₜᵀ，秩为k
For t = 1, 2, ...
1. 接收新数据xₜ（或ΔAₜ）
2. 计算残差：r = xₜ - UₜUₜᵀxₜ
3. 更新子空间：
- 若‖r‖ > θ：扩展U空间
- 否则：仅更新系数
4. 可选：秩截断保持k
5. 可选：应用遗忘因子
关键设计选择：
详细的扩展策略分析：
策略1：固定阈值
if ‖r‖ > θ:
扩展子空间
else:
仅更新系数
优点：简单直观 缺点：阈值选择困难，对数据尺度敏感
策略2：相对阈值
if ‖r‖/‖xₜ‖ > θ_rel:
扩展子空间
优点：尺度不变性 适用：数据范数变化大的场景
策略3：累积误差
维护累积残差能量E_res
if E_res > θ_energy:
扩展并重置E_res
优点：考虑历史信息 适用：缓慢变化的数据流
截断策略比较：
将低秩近似视为约束优化： \(\min_{\mathbf{U}\in\mathbb{R}^{n\times k}} \sum_{t=1}^T \ell_t(\mathbf{U}) = \sum_{t=1}^T \|\mathbf{x}_t - \mathbf{U}\mathbf{U}^T\mathbf{x}_t\|^2\)
这建立了与在线梯度下降、镜像下降等经典算法的联系。特别地，Oja’s rule可视为此问题的随机梯度解法。
Oja’s算法详解： 更新规则： \(\mathbf{U}_{t+1} = \mathbf{U}_t + \eta_t(\mathbf{x}_t\mathbf{x}_t^T - \mathbf{U}_t\mathbf{U}_t^T\mathbf{x}_t\mathbf{x}_t^T)\mathbf{U}_t\)
收敛性质：
Oja算法的几何解释： Oja更新可分解为两步：
这种”增强-正交化”模式在神经科学中有对应：Hebbian学习与侧向抑制。
推广：在线矩阵流形优化
考虑Stiefel流形约束：$\mathbf{U}^T\mathbf{U} = \mathbf{I}_k$
Riemannian SGD：
1. 计算欧氏梯度：∇f(Uₜ)
2. 投影到切空间：grad = ∇f - Uₜ(Uₜᵀ∇f)
3. 沿测地线更新：Uₜ₊₁ = Retr(Uₜ, -ηₜ·grad)
其中Retraction可选：
在线镜像下降视角：
定义Bregman散度： \(D_\phi(\mathbf{U}, \mathbf{V}) = \phi(\mathbf{U}) - \phi(\mathbf{V}) - \langle\nabla\phi(\mathbf{V}), \mathbf{U} - \mathbf{V}\rangle\)
选择合适的$\phi$可得到不同算法：
加速方法：
动量加速Oja： \(\begin{aligned} \mathbf{V}_{t+1} &= \beta\mathbf{V}_t + (1-\beta)\nabla_{\mathbf{U}}\ell_t(\mathbf{U}_t) \\ \mathbf{U}_{t+1} &= \mathbf{U}_t - \eta_t\mathbf{V}_{t+1} \end{aligned}\)
Nesterov加速变体： \(\begin{aligned} \mathbf{Y}_t &= \mathbf{U}_t + \frac{t-1}{t+2}(\mathbf{U}_t - \mathbf{U}_{t-1}) \\ \mathbf{U}_{t+1} &= \mathbf{Y}_t - \eta_t\nabla\ell_t(\mathbf{Y}_t) \end{aligned}\)
与经典在线算法的对比：
| 算法 | 更新复杂度 | 收敛率 | 约束处理 | 内存需求 |
|---|---|---|---|---|
| Oja’s rule | $O(nk)$ | $O(1/t)$ | 渐近满足 | $O(nk)$ |
| 在线梯度下降 | $O(nk)$ | $O(\sqrt{T})$ | 投影步 | $O(nk)$ |
| Riemannian SGD | $O(nk^2)$ | $O(1/t)$ | 精确满足 | $O(nk)$ |
| Block power | $O(nkb)$ | $O(1/t^2)$ | 渐近满足 | $O(nkb)$ |
| SVRG-style | $O(nk)$ | $O(1/t^2)$ | 投影步 | $O(nk)$ |
在线核PCA扩展：
映射到RKHS空间$\phi: \mathbb{R}^n \rightarrow \mathcal{H}$： \(\min_{\mathbf{U}} \sum_{t=1}^T \|\phi(\mathbf{x}_t) - \mathbf{U}\mathbf{U}^T\phi(\mathbf{x}_t)\|_{\mathcal{H}}^2\)
使用核技巧避免显式映射：
研究方向：
给定SVD分解$\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T$，考虑秩一更新$\mathbf{A}’ = \mathbf{A} + \mathbf{c}\mathbf{d}^T$。
引理13.1（Bunch-Nielsen）：更新后的SVD可通过以下步骤计算：
计算复杂度：$O(nk + k^3)$，远优于重新计算的$O(n^2k)$。
详细推导：
将$\mathbf{A}’ = \mathbf{A} + \mathbf{c}\mathbf{d}^T$重写为： \(\mathbf{A}' = [\mathbf{U}, \mathbf{c}] \begin{bmatrix} \boldsymbol{\Sigma} & \mathbf{0} \\ \mathbf{0}^T & 1 \end{bmatrix} [\mathbf{V}, \mathbf{d}]^T\)
但$[\mathbf{U}, \mathbf{c}]$和$[\mathbf{V}, \mathbf{d}]$可能不正交。通过正交化：
\[[\mathbf{U}, \mathbf{c}] = [\mathbf{U}, \tilde{\mathbf{u}}]\begin{bmatrix} \mathbf{I} & \mathbf{p} \\ \mathbf{0}^T & \rho_c \end{bmatrix}\]其中$\tilde{\mathbf{u}} = \mathbf{r}_c/|\mathbf{r}_c|$，$\rho_c = |\mathbf{r}_c|$。
稳定性分析：
定理13.3：设$\kappa = |\mathbf{c}||\mathbf{d}|/\sigma_{\min}(\mathbf{A})$，则相对误差满足： \(\frac{\|\tilde{\mathbf{A}}' - \mathbf{A}'\|_F}{\|\mathbf{A}'\|_F} \leq \epsilon_{\text{machine}} \cdot O(\kappa^2)\)
这表明当更新幅度相对于最小奇异值过大时，算法可能不稳定。
特殊情况处理：
算法13.2：Brand’s增量SVD
输入：当前SVD UΣVᵀ (秩k)，新列c
1. 计算投影系数：m = Uᵀc
2. 计算正交分量：p = c - Um, ρ = ‖p‖
3. 若ρ > ε：
构造扩展矩阵并SVD分解
4. 否则：
仅更新Σ中的系数
5. 定期重正交化U, V
数值稳定性改进：
稳定性分析的深入考察：
浮点误差传播模型： 设$\epsilon_m$为机器精度，第$t$步后的误差界： \(E_t \leq E_0 \cdot \kappa^t + \epsilon_m \cdot \frac{\kappa^t - 1}{\kappa - 1}\)
其中$\kappa = \sigma_1/\sigma_k$是条件数。这表明：
自适应重正交化策略：
定义正交性损失：δ = ‖UᵀU - I‖_F
触发条件：
1. δ > τ_orth (典型值：10⁻⁸)
2. 更新次数达到N_max = ⌊1/(10·ε_machine)⌋
3. 条件数激增：κ(t)/κ(0) > 100
混合精度技术：
完整的Brand算法实现细节：
function updateSVD(U, Σ, V, c):
// 步骤1: 投影到当前子空间
m = Uᵀc
p = c - Um
ρ = ‖p‖
// 步骤2: 判断是否需要扩展秩
if ρ > ε_threshold:
// 步骤3a: 构造扩展矩阵
P = p/ρ // 归一化
Q = zeros(k+1, 1)
Q[k+1] = 1
// 构造(k+1)×(k+1)矩阵
K = [Σ m]
[0 ρ]
// 步骤4a: 对K进行SVD
[U', Σ', V'] = svd(K)
// 步骤5a: 更新大矩阵
U_new = [U P] @ U'
V_new = [V Q] @ V'
Σ_new = Σ'
else:
// 步骤3b: 仅更新系数
K = Σ + m*1ᵀ // 秩一更新
[U', Σ', V'] = svd(K)
U_new = U @ U'
V_new = V @ V'
Σ_new = Σ'
return U_new, Σ_new, V_new
高级实现技巧：
1. 延迟更新策略：
维护缓冲区B存储待处理更新
if size(B) < batch_size:
B.append(c)
else:
// 批量更新
[Q_B, R_B] = qr(B)
K = [Σ UᵀQ_B·R_B]
[0 Q_B'ᵀQ_B·R_B]
执行块更新
清空B
2. 自适应阈值选择：
ε_threshold = max(
ε_abs, // 绝对阈值
ε_rel · ‖c‖, // 相对阈值
√(ε_machine) · σ_k // 数值阈值
)
3. 增量正交化： 使用修正Gram-Schmidt的增量版本：
for j = 1:k:
α_j = u_j'p
p = p - α_j·u_j
// 二次修正提高精度
β_j = u_j'p
p = p - β_j·u_j
m[j] = m[j] + β_j
正交化策略比较：
自适应重正交化：
function checkOrthogonality(U):
orth_error = ‖UᵀU - I‖_F
if orth_error > ε_reorth:
U = reorthogonalize(U)
return U
重正交化时机：
对于批量更新$\mathbf{A}’ = \mathbf{A} + \mathbf{C}\mathbf{D}^T$（$\mathbf{C}, \mathbf{D} \in \mathbb{R}^{n \times b}$）：
定理13.2：块更新可分解为 \(\mathbf{A}' = \begin{bmatrix} \mathbf{U} & \tilde{\mathbf{U}} \end{bmatrix} \tilde{\boldsymbol{\Sigma}} \begin{bmatrix} \mathbf{V} & \tilde{\mathbf{V}} \end{bmatrix}^T\)
其中$\tilde{\mathbf{U}}, \tilde{\mathbf{V}}$来自$\mathbf{C}, \mathbf{D}$的QR分解。
这在处理mini-batch更新时特别高效。
块更新算法详解：
function blockUpdateSVD(U, Σ, V, C, D):
// 步骤1: QR分解输入矩阵
[Q_C, R_C] = qr(C - U(UᵀC)) // 正交化C
[Q_D, R_D] = qr(D - V(VᵀD)) // 正交化D
// 步骤2: 构造扩展矩阵
K = [Σ UᵀC]
[VᵀD R_C @ R_Dᵀ]
// 步骤3: 对K进行SVD (size: (k+b)×(k+b))
[U_K, Σ_K, V_K] = svd(K)
// 步骤4: 更新原始矩阵
U_new = [U, Q_C] @ U_K
V_new = [V, Q_D] @ V_K
Σ_new = Σ_K
// 步骤5: 可选秩截断
if rank(Σ_new) > k_max:
[U_new, Σ_new, V_new] = truncate(U_new, Σ_new, V_new, k_max)
return U_new, Σ_new, V_new
效率分析：
复杂度分解：
总复杂度：$O(nb^2 + n(k+b)^2 + (k+b)^3)$
与逐个更新的对比：
特殊块结构的优化：
在分布式设置下，不同节点观察数据流的不同部分：
算法13.3：异步分布式SVD更新
节点i维护局部近似Uᵢ
1. 局部更新：处理本地数据流
2. 周期通信：交换子空间信息
3. 共识步骤：对齐不同节点的子空间
4. 全局更新：合并得到全局U
关键挑战：
理论分析：
定理13.4（分布式收敛性）：设$m$个节点，通信周期$\tau$，则： \(\mathbb{E}[\|\mathbf{U}_T - \mathbf{U}_*\|_F] \leq O\left(\frac{\sqrt{m\tau}}{T} + \frac{1}{\sqrt{mT}}\right)\)
第一项来自通信延迟，第二项来自分布式平均。
通信复杂度下界： 定理13.5：任何$\epsilon$-近似分布式PCA算法需要至少$\Omega(mk\log(1/\epsilon))$比特通信。
证明基于信息论：每个节点需传递足够信息以重构全局主成分。
分布式架构详解：
// 节点i的本地算法
function nodeUpdate(local_data_stream):
// 本地状态
U_local = [] // n×k_local
Σ_local = [] // k_local×k_local
buffer = [] // 数据缓冲
while true:
// 阶段1: 本地处理
batch = collect_batch(local_data_stream)
U_local, Σ_local = incrementalSVD(U_local, Σ_local, batch)
// 阶段2: 周期通信
if time_to_communicate():
// 发送本地子空间的紧凑表示
sketch = compute_sketch(U_local, Σ_local)
broadcast(sketch)
// 接收其他节点的sketch
sketches = receive_all()
// 阶段3: 子空间对齐
U_global = align_subspaces(sketches)
// 阶段4: 本地投影更新
U_local = project_to_global(U_local, U_global)
子空间对齐算法：
方法1：Grassmannian平均
function align_subspaces(U_list):
// 计算Grassmannian重心
U_mean = grassmannian_mean(U_list)
return U_mean
function grassmannian_mean(U_list):
// 迭代算法
U_avg = U_list[0]
for iter in 1:max_iters:
// 计算到各子空间的测地线
tangents = []
for U in U_list:
v = log_map(U_avg, U)
tangents.append(v)
// 平均切向量
v_mean = mean(tangents)
// 沿测地线更新
U_avg = exp_map(U_avg, v_mean)
return U_avg
方法2：分布式特征分解
function distributed_eigen_alignment():
// 构造分布式协方差矩阵
C_global = sum_all_nodes(U_local @ Σ_local² @ U_localᵀ)
// 分布式特征分解
[V, Λ] = distributed_eig(C_global, k)
return V[:, 1:k]
通信优化技术：
function quantize_subspace(U, bits):
// 随机旋转
R = random_orthogonal(k)
U_rot = U @ R
// 量化
U_quant = quantize(U_rot, bits)
return U_quant, R
function sketch_subspace(U, Σ, sketch_size):
// Johnson-Lindenstrauss投影
S = randn(n, sketch_size) / sqrt(sketch_size)
Y = Sᵀ @ U @ Σ
return Y // sketch_size × k
拜占庭鲁棒性：
function byzantine_robust_aggregation(sketches, f):
// f = 拜占庭节点数量上界
// 方法1: 中位数聚合
median_sketch = geometric_median(sketches)
// 方法2: 修剪均值
distances = compute_distances(sketches)
good_nodes = trim_outliers(distances, f)
robust_mean = mean(sketches[good_nodes])
return robust_mean
高级通信协议：
1. 渐进式精度协议：
初始阶段：交换低精度sketch (1-2 bits)
中期阶段：提高到中等精度 (8-16 bits)
最终阶段：全精度交换关键成分
优势：早期快速收敛，后期精确对齐
2. 异步Gossip SVD：
每个节点维护邻居列表N(i)
随机选择邻居j ∈ N(i)
交换并平均子空间：
U_i^new = geodesic_mean(U_i, U_j)
U_j^new = U_i^new
收敛速度：O(log m / λ₂(G))
其中$\lambda_2(G)$是通信图的第二大特征值。
3. 分层聚合树：
叶节点：本地SVD更新
中间节点：合并子树结果
根节点：全局SVD
优势：通信复杂度O(log m)
缺点：单点故障风险
实际系统考虑：
带宽感知调度：
容错机制：
function robustAggregation(sketches):
// 检测异常值
median = geometric_median(sketches)
deviations = [distance(s, median) for s in sketches]
// 移除异常
threshold = median(deviations) * 3
valid = [s for s,d in zip(sketches, deviations) if d < threshold]
// 加权聚合
weights = [1/d for d in valid_deviations]
return weighted_mean(valid, weights)
研究方向：
流式环境下的秩选择面临独特挑战：需在观察完整数据前做出决策。核心权衡：
定义13.1（在线秩选择regret）： \(R_k(T) = \sum_{t=1}^T \|\mathbf{A}_t - \tilde{\mathbf{A}}_{t,k_t}\|_F^2 - \min_{k^*} \sum_{t=1}^T \|\mathbf{A}_t - \tilde{\mathbf{A}}_{t,k^*}\|_F^2\)
其中$k_t$是算法在时刻$t$选择的秩。
算法13.4：自适应能量阈值
参数：能量保留比例τ ∈ (0,1)
维护：奇异值{σᵢ}及对应向量
1. 更新奇异值分解
2. 计算累积能量：Eⱼ = Σᵢ₌₁ʲ σᵢ² / Σᵢ σᵢ²
3. 选择秩：k* = min{j : Eⱼ ≥ τ}
4. 自适应调整τ基于历史性能
定理13.3：在平稳假设下，自适应能量阈值达到 \(\mathbb{E}[R_k(T)] \leq O(\sqrt{T\log K})\) 其中$K$是最大允许秩。
证明思路： 使用在线学习的expert advice框架：
实用增强版本：
动态阈值调整：
function adaptThreshold(τ, history):
// 计算近期重构误差
recent_error = mean(history[-window:])
// 误差趋势分析
if increasing_trend(history):
τ = min(τ + Δτ, 0.99) // 需要更多能量
elif stable_low_error(history):
τ = max(τ - Δτ, 0.90) // 可以减少能量
return τ
多尺度能量分析：
function multiScaleEnergy(σ_values):
scales = [1, 5, 10, 20] // 不同观察尺度
energy_profiles = []
for s in scales:
// 计算前s个奇异值的能量比
E_s = sum(σ²[1:s]) / sum(σ²)
energy_profiles.append(E_s)
// 检测“肘部”位置
k_elbow = detectElbow(energy_profiles)
return k_elbow
噪声鲁棒性考虑：
当数据包含噪声时，需要避免过拟合：
Marcenko-Pastur律应用： 对于随机矩阵，奇异值分布遍循： \(\rho(\sigma) = \frac{1}{2\pi\sigma c}\sqrt{(\sigma_+ - \sigma)(\sigma - \sigma_-)}\)
其中$\sigma_\pm = (1 \pm \sqrt{c})^2$，$c = n/T$。
噪声阈值策略：
function noiseAwareRank(σ_values, n, T):
c = n/T
σ_threshold = (1 + √c)² * noise_level
// 只保留显著大于噪声的奇异值
k = sum(σ_values > σ_threshold)
return max(k, 1) // 至少保留一个
利用奇异值衰减模式预测未来：
模型假设：
算法13.5：预测性秩选择
1. 在线估计衰减参数α或β
2. 预测未来L步的奇异值分布
3. 选择秩k最小化预期误差：
k* = argmin_k E[Σₜ₊₁ᵗ⁺ᴸ ‖Aₜ - Ãₜ,ₖ‖²]
4. 使用滑动窗口更新参数估计
详细的参数估计方法：
1. 幂律模型的在线估计：
function estimatePowerLaw(σ_history):
// 取对数变换
log_i = log(1:length(σ))
log_σ = log(σ)
// 加权最小二乘
weights = 1./sqrt(1:length(σ)) // 前面的奇异值更重要
α = -weightedLeastSquares(log_i, log_σ, weights)
// 鲁棒性检查
residuals = log_σ - (-α * log_i)
if max(abs(residuals)) > threshold:
// 使用Huber回归
α = huberRegression(log_i, log_σ)
return α
2. 指数模型的自适应估计：
function estimateExponential(σ_history, window):
// 使用最近的window个点
recent_σ = σ_history[-window:]
// MLE估计
β_mle = 1 / mean(recent_σ)
// 贝叶斯更新（结合先验）
prior_mean = 1 / median(σ_history)
prior_weight = 0.1
β = (1-prior_weight) * β_mle + prior_weight * prior_mean
return β
3. 混合模型识别：
function identifyDecayPattern(σ):
// 计算不同模型的拟合度
power_fit = fitPowerLaw(σ)
exp_fit = fitExponential(σ)
mixed_fit = fitMixedModel(σ)
// AIC准则选择
aic_power = AIC(power_fit, σ)
aic_exp = AIC(exp_fit, σ)
aic_mixed = AIC(mixed_fit, σ)
return argmin([aic_power, aic_exp, aic_mixed])
预测误差的理论分析：
定理13.6（预测误差界）：设真实衰减率为$\alpha_*$，估计值为$\hat{\alpha}_t$，则： \(\mathbb{E}[\text{PredErr}_t] \leq C_1 \cdot \mathbb{E}[|\hat{\alpha}_t - \alpha_*|] + C_2 \cdot \frac{\log t}{t}\)
第一项来自参数估计误差，第二项来自有限样本效应。
实用策略：保守选择：
function conservativeRankSelection(predicted_k, confidence):
// 加入安全边界
safety_margin = ceil(predicted_k * (1 - confidence))
// 考虑预测不确定性
uncertainty = estimatePredictionUncertainty()
extra_ranks = ceil(2 * sqrt(uncertainty))
return predicted_k + safety_margin + extra_ranks
将秩视为隐变量，使用在线贝叶斯推断：
概率模型： \(p(\mathbf{A}|k) = \mathcal{MN}(\mathbf{0}, \mathbf{I}_n, \boldsymbol{\Sigma}_k)\) 其中$\boldsymbol{\Sigma}_k$是秩$k$协方差。
在线推断：
| 维护秩的后验分布$p(k | D_{1:t})$ |
详细的贝叶斯框架：
先验设计：
// 基于领域知识的先验
π(k) ∝ k^{-γ} * exp(-λ*k) // 偏好低秩
// 参数选择
γ: 控制幂律衰减 (typical: 1-2)
λ: 控制指数惩罚 (typical: 0.01-0.1)
似然模型：
// 概率性PCA模型
p(σ²_i | k) = {
InverseGamma(α, β), if i ≤ k // 信号奇异值
InverseGamma(α_0, β_0), if i > k // 噪声奇异值
}
// 参数学习
α, β: 从数据中学习
α_0, β_0: 噪声先验
在线更新算法：
function bayesianRankUpdate(posterior_t, new_data):
// 计算每个秩的似然
for k in 1:K_max:
likelihood[k] = computeLikelihood(new_data, k)
// 贝叶斯更新
posterior_t+1 = posterior_t .* likelihood
posterior_t+1 = posterior_t+1 / sum(posterior_t+1)
// 防止数值下溢
posterior_t+1 = renormalize(posterior_t+1)
return posterior_t+1
变分推断加速：
为避免计算所有$K$个可能秩的似然：
function variationalRankInference(data):
// 变分分布: q(k) = Categorical(π)
π = uniform(K_max)
for iter in 1:max_iters:
// E-step: 更新秩分布
for k in 1:K_max:
π[k] ∝ exp(ELBO(data, k))
π = normalize(π)
// M-step: 更新模型参数
updateModelParams(data, π)
// 收敛检查
if converged(π):
break
return π
实时决策策略：
function selectRank(posterior, risk_aversion):
// MAP估计（最激进）
k_map = argmax(posterior)
// 后验均值（平衡）
k_mean = sum(k * posterior[k] for k in 1:K_max)
// 风险敏感决策
// 最小化预期损失: L(k,k') = |k-k'|^α
k_bayes = argmin_k sum(L(k,k') * posterior[k'] for k' in 1:K_max)
// 根据风险偏好选择
if risk_aversion == "low":
return k_map
elif risk_aversion == "medium":
return k_mean
else:
return k_bayes
研究方向：
考虑部分观测的流式设定：
形式化：在线学习协议
For t = 1, 2, ..., T:
1. 算法预测M̂ₜ（秩≤r）
2. 对手选择(iₜ, jₜ)
3. 观测真实值Mᵢₜ,ⱼₜ
4. 遭受损失ℓₜ = (M̂ₜ[iₜ,jₜ] - Mᵢₜ,ⱼₜ)²
矩阵补全的因子分解形式：$\mathbf{M} \approx \mathbf{U}\mathbf{V}^T$
算法13.6：流式矩阵分解
初始化：U₀, V₀ ∈ ℝⁿˣʳ
For each观测(i,j,Mᵢⱼ):
1. 预测：M̂ᵢⱼ = ⟨uᵢ, vⱼ⟩
2. 梯度计算：
∇ᵤᵢ = -2(Mᵢⱼ - M̂ᵢⱼ)vⱼ
∇ᵥⱼ = -2(Mᵢⱼ - M̂ᵢⱼ)uᵢ
3. 更新：
uᵢ ← uᵢ - ηₜ∇ᵤᵢ
vⱼ ← vⱼ - ηₜ∇ᵥⱼ
4. 可选正则化/投影步骤
收敛性分析：
将低秩约束视为流形约束： \(\mathcal{M}_r = \{\mathbf{X} \in \mathbb{R}^{m \times n} : \text{rank}(\mathbf{X}) = r\}\)
算法13.7：在线Riemannian梯度下降
维护：当前点Xₜ ∈ 𝓜ᵣ
1. 计算欧氏梯度∇f(Xₜ)
2. 投影到切空间：gradₜ = Proj_{TXₜ}(∇f)
3. 沿测地线更新：
Xₜ₊₁ = Retr_Xₜ(-ηₜ gradₜ)
关键优势：
实际应用常有辅助信息（如用户特征、时间戳）：
模型扩展： \(M_{ij} = \langle \mathbf{u}_i, \mathbf{v}_j \rangle + f(\mathbf{x}_i, \mathbf{y}_j; \boldsymbol{\theta})\)
其中$f$可以是：
在线学习策略：
研究方向：
处理非平稳数据流的经典方法：
更新规则： \(\mathbf{A}_t = \gamma \mathbf{A}_{t-1} + \mathbf{x}_t\mathbf{x}_t^T\)
其中$\gamma \in (0,1)$是遗忘因子。
性质分析：
算法13.8：基于变化检测的自适应遗忘
维护：多个时间尺度的统计量
1. 短期统计：Sₛ = 最近k个样本
2. 长期统计：Sₗ = 指数平滑历史
3. 检测变化：D = d(Sₛ, Sₗ)
4. 调整遗忘：
if D > threshold:
γ ← max(γ - δ, γₘᵢₙ) # 加快遗忘
else:
γ ← min(γ + δ, γₘₐₓ) # 减慢遗忘
理论保证：在分段平稳假设下，自适应算法达到 \(\text{Regret} \leq O(\sqrt{T(1 + N_c)})\) 其中$N_c$是变化点数量。
保持固定大小$W$的数据窗口：
精确滑窗SVD：
数据结构：循环缓冲区存储最近W个样本
更新策略：
1. 加入新样本：秩1更新
2. 移除旧样本：秩1降级（downdating）
3. 周期性从头计算以控制误差累积
近似滑窗方法：
不同奇异向量可能有不同的时间尺度：
模型： \(\mathbf{u}_i(t) = \gamma_i \mathbf{u}_i(t-1) + \text{update}\)
其中$\gamma_i$依赖于第$i$个模式的稳定性。
自动尺度发现：
研究方向：
深度网络中的权重矩阵$\mathbf{W} \in \mathbb{R}^{m \times n}$常呈现低秩结构：
在线压缩框架：
训练过程中：
1. 监控权重矩阵的有效秩
2. 当检测到低秩结构时分解：W ≈ UV^T
3. 继续以因子形式训练
4. 周期性评估是否需要调整秩
与剪枝的对比：
算法13.9：训练时动态压缩
每E个epoch:
1. 对当前权重W计算SVD
2. 评估能量分布确定有效秩k
3. 压缩：W ← UₖΣₖVₖᵀ
4. 可选：切换到因子形式继续训练
5. 监控压缩对性能的影响
关键考虑：
彩票假说：随机初始化的网络包含能独立训练到同等精度的稀疏子网络。
低秩重述：
实验观察：
在序列任务学习中防止遗忘：
策略13.1：正交子空间分配
对每个新任务t:
1. 识别之前任务使用的子空间Uₚᵣₑᵥ
2. 在正交补空间中学习：U⊥
3. 合并：Uₜₒₜₐₗ = [Uₚᵣₑᵥ | Uₙₑw]
4. 当空间耗尽时压缩或合并
策略13.2：弹性权重巩固(EWC)的低秩版本
研究方向：
本章系统介绍了动态低秩近似的理论与算法：
核心概念：
关键公式：
实践要点：
习题13.1：推导秩一更新的计算复杂度 给定$n \times n$矩阵的SVD和秩一更新$\mathbf{u}\mathbf{v}^T$，证明更新后的SVD可在$O(n^2)$时间内计算。
习题13.2：设计滑动窗口SVD算法 实现保持最近$W$个样本的精确SVD，分析添加新样本和删除旧样本的复杂度。
习题13.3：比较不同遗忘因子 对于$\gamma \in {0.9, 0.95, 0.99}$，计算有效窗口长度，并分析其对突变检测的影响。
习题13.4：在线矩阵补全的收敛性 证明在强凸假设下，在线梯度下降达到$O(\log T)$ regret。
习题13.5：非均匀采样的偏差校正 在矩阵补全中，若位置$(i,j)$以概率$p_{ij}$被观测，设计无偏估计算法。
习题13.6：多尺度遗忘的最优性 证明对于分层时间序列数据，多尺度遗忘优于单一遗忘因子。构造具体反例。
习题13.7：分布式流式PCA的通信下界 $m$个节点观测数据流，证明达到$\epsilon$-近似需要$\Omega(mk/\epsilon)$通信量。
习题13.8：神经网络压缩的理论保证 给定预训练网络，证明存在秩$k$分解使得精度损失不超过$\epsilon$，其中$k$如何依赖于$\epsilon$？