HW5–6：从稳定 SVD 到并行优化

Ex.1 — Numerical stability

1.1 大数据中提高精度是否值得

题意：题目要求权衡精度、存储、网络和吞吐，不接受“精度越高越好”的单向结论。

知识点：float 通常约 7 位有效十进制数字；double 约 15–16 位。更高精度降低累计、相消、病态分解和长迭代中的舍入误差，但增加内存、磁盘与网络流量，并可能降低缓存、SIMD 或 GPU 吞吐。

解法结构：

1. 先说误差收益：敏感 reduction 与 factorization 更可靠。
2. 再说系统代价：double 的每元素容量通常约为 float 的两倍。
3. 指出边界：高精度不能修复噪声数据、错误模型或不稳定算法。
4. 给条件化结论：依据 condition number、误差预算与失败证据选择；可采用混合精度。

1.2 四种操作、100 个随机 1000×100 矩阵

题面事实：分别累计测量 svd(X)、svd(X')、eig(X*X') 与 eig(X'*X)。

知识点：XX^T 是 1000×1000，X^TX 是 100×100；两者的非零特征值相同。X 与 X^T 的非零奇异值也相同，但库实现、full/economy 模式与 workspace 会影响计时。

rng(0);
trials = 100;
elapsed = zeros(trials, 4);

% 可先用一份矩阵 warm-up，避免把首次库初始化混入主计时
X = randn(1000, 100);
svd(X, "econ"); svd(X', "econ"); eig(X*X'); eig(X'*X);

for t = 1:trials
    X = randn(1000, 100);
    tic; svd(X, "econ"); elapsed(t,1) = toc;
    tic; svd(X', "econ"); elapsed(t,2) = toc;
    tic; eig(X*X'); elapsed(t,3) = toc;
    tic; eig(X'*X); elapsed(t,4) = toc;
end

total = sum(elapsed, 1);
average = mean(elapsed, 1);
spread = std(elapsed, 0, 1);

结果解释框架：报告 MATLAB 版本、CPU/BLAS、随机种子、econ 约定、总时间、均值和离散度。小 Gram 通常更快，但显式形成 X^TX 会令 2-norm condition number 从 κ(X) 变为 κ(X)²；速度不自动代表稳定。

常见扣分点：只生成一个矩阵；只测一次；把本机秒数写成算法定律；不说明 full/econ；把形成 Gram 的时间排除却不声明。

1.3 1000 次随机扰动

题面事实：对给定 5×5 非对称矩阵 X，分别研究 X+δX 的 eigenvalues 与 singular values 在 1000 次小随机扰动下的变化。题目提示可参考 MATLAB eps。

rng(0);
trials = 1000;
scale = eps(norm(X, "fro")) * norm(X, "fro");
baseEig = eig(X);
baseS = sort(svd(X), "descend");

eigSamples = complex(zeros(5, trials));
sSamples = zeros(5, trials);
for t = 1:trials
    dX = scale * randn(size(X));
    e = eig(X + dX);
    [~, order] = sortrows([real(e), imag(e)], [1 2]);
    eigSamples(:,t) = e(order);
    sSamples(:,t) = sort(svd(X + dX), "descend");
end

% 报告 absolute error 的 median/std/quantiles，避免 signed error 抵消

知识点：非对称矩阵的 eigenvalues 可能为复数，跨 trial 比较前必须固定排序或匹配规则。singular values 非负且可降序匹配，并满足扰动界 |σ_i(X+E)−σ_i(X)|≤‖E‖₂。

常见扣分点：把确定性的逐元素 eps(X) 称为随机扰动；直接比较未匹配 eigenvalues；累加 signed difference 使正负误差抵消；不写扰动尺度。

1.4 用 Chapter 4 解释实验

一般特征分解 M=PDP⁻¹ 的局部扰动放大与 ‖P‖‖P⁻¹‖ 有关。SVD 的左右因子正交，因此对应 2-norm 条件因子为 1。答题时还必须说明：SVD 较稳不表示所有 SVD 路线都最快，且通过 Gram 求 SVD 会损失一部分稳定性优势。

Ex.2 — 不用计算机求完整 singular decomposition

题面矩阵：

题意：题面写的是 singular decomposition，不是只求 singular values。答案必须交出 U,Σ,V^T，说明 full/thin 约定，并详写步骤。

最省算力路线：因为 3<4，先算 3×3 的 XX^T：

调研材料记录其正特征值约为 235.6958、53.5435、0.7607，因此奇异值约为 15.3524、7.3173、0.8722。先求对应的正交归一左特征向量 u₁,u₂,u₃，再用 v_i=X^Tu_i/σ_i 恢复前三个右奇异向量。最后解 Xv₄=0 并归一，补出 full V∈ℝ^4×4。

完整手算合同与三个可算尽的小例见下一专题。对 HW5 大矩阵，最终必须核验 U^TU=I₃、V^TV=I₄ 与 UΣV^T≈X。

Ex.3 — PCA in Spark

题面：两个无 header 的 sensor CSV 中至少一个可能包含 1001 个电路传感器输出，最后一列 y 是每小时用电量。需要解释 PCA 的帮助、找出解释 90% 数据的最小列数 n、用前 n 个 principal components 建模 y，并判断 sensors2 是否同类。

3.1 PCA 如何帮助

PCA 把强相关传感器列旋转为正交主方向，按方差排序。它可以揭示低维结构、减少后续回归维数，并帮助观察异常或分布差异。PCA 不知道列的物理名称，也不能单独证明两个文件来自同一系统。

3.2 解释 90% 所需列数

先把最后一列 y 分离。只在 feature matrix 上拟合 scaler/PCA。若解释方差比为 e₁≥e₂≥…，取最小 n 使：

常见扣分点：把“每个 component 的 EVR 大于 0.01”当成累计 90% 规则；把 y 放进 PCA，造成 label leakage。

3.3 用 principal-component scores 回归

import numpy as np
from pyspark.ml.feature import PCA, StandardScaler, VectorAssembler
from pyspark.ml.regression import LinearRegression

# 无 header CSV：最后一列是 y，其余列是传感器特征
raw = spark.read.option("inferSchema", True).csv(path)
feature_cols = raw.columns[:-1]
label_col = raw.columns[-1]
data = raw.withColumnRenamed(label_col, "label")
train_raw, test_raw = data.randomSplit([0.8, 0.2], seed=42)

assembler = VectorAssembler(
    inputCols=feature_cols,
    outputCol="features"
)
train_assembled = assembler.transform(train_raw)
test_assembled = assembler.transform(test_raw)

scaler = StandardScaler(
    inputCol="features",
    outputCol="scaled",
    withMean=True,
    withStd=True
).fit(train_assembled)                 # 只在训练集拟合
train_scaled = scaler.transform(train_assembled)
test_scaled = scaler.transform(test_assembled)

# 只用训练集方差谱选择累计解释方差达到 90% 的最小 n
full_pca = PCA(
    k=len(feature_cols),
    inputCol="scaled",
    outputCol="all_pc_scores"
).fit(train_scaled)
cum_evr = np.cumsum(full_pca.explainedVariance.toArray())
n = int(np.searchsorted(cum_evr, 0.90) + 1)

pca_model = PCA(
    k=n,
    inputCol="scaled",
    outputCol="pc_scores"
).fit(train_scaled)                    # 仍只在训练集拟合
train_scores = pca_model.transform(train_scaled)
test_scores = pca_model.transform(test_scaled)

model = LinearRegression(
    featuresCol="pc_scores",
    labelCol="label"
).fit(train_scores)
test_predictions = model.transform(test_scores)

管线合同：先切分，再由训练集拟合 scaler、选择 n、拟合 PCA 和回归；测试集只调用 transform。test_predictions 才用于 held-out R²/RMSE，因而不会把测试分布泄漏进特征空间。

模型写为 y=β₀+Σβ_ip_i+ε。题面给 ε∼Normal(0,1) 是模型假设；应检查残差，而不能把它当作已证事实。报告 held-out R²/RMSE、split/seed 与 residual diagnostics。

3.4 判断 sensors2

对两份文件采用相同预处理、split 与评价口径，比较 cumulative spectrum、PC-score regression 的 held-out 表现、残差结构和基线。仅凭一个高 R² 或二维散点相似，不能证明来源相同。

1. 先写 shape：避免一开始就丢分

若 A∈ℝ^m×n、r=rank(A)：

题面要求 decomposition 时，若未声明 thin，应优先给 full 或明确说明所用约定。

2. 七步通用流程

1. 写形状与秩上界。r≤min(m,n)，Σ 的位置由形状决定。
2. 选较小 Gram。宽矩阵优先 AA^T；高矩阵优先 A^TA。
3. 求对称 Gram 的特征对。特征向量可取为正交归一；特征值理论上非负。
4. 排序并开方。σ_i=√λ_i，按降序排列；向量列必须同步重排。
5. 恢复另一侧。若先得 v_i，用 u_i=Av_i/σ_i；若先得 u_i，用 v_i=A^Tu_i/σ_i。只对 σ_i>0 使用。
6. 补 full bases。null(A) 补 V；null(A^T) 补 U。对补向量正交归一。
7. 三重验算。检查列单位正交、奇异方程与重构。

形状 → 小 Gram → eigenpairs → σ 降序
     → nonzero modes 恢复另一侧
     → nullspace completion
     → UᵀU / VᵀV / UΣVᵀ 三重验算

3. 左右奇异向量到底属于哪一边

• A^TA∈ℝ^n×n 的 eigenvectors 是 V 的列，位于输入空间。
• AA^T∈ℝ^m×m 的 eigenvectors 是 U 的列，位于输出空间。
• 正的非零 eigenvalues 相同，均为 σ_i²。

记忆法：Gram 的外侧字母决定向量：A^TA 最外是输入维度 n，所以给 V；AA^T 是输出维度 m，所以给 U。

Example 1：方阵、两个不同正奇异值

Step 1 — Gram：

Step 2 — 特征多项式：

故 λ₁=16、λ₂=4，已降序。

Step 3 — 右奇异向量：

Step 4 — 奇异值与左向量：

Step 5 — 组装：

Step 6 — 验算：

陷阱：本例 U=V 是因为 A 对称正定，不是所有方阵的规律。

Example 2：重复奇异值与旋转不唯一

Step 1：B^TB=4I₂，所以 λ₁=λ₂=4，σ₁=σ₂=2。

Step 2：最简单的选择是：

Step 3 — 展示全部不唯一性：对任意二维正交矩阵 Q，例如旋转矩阵：

结论：重复奇异值对应的左右 singular subspaces 可作同一个正交旋转。不同同学得到不同 U/V 不代表其中一个错误；应检查正交性和重构。

排序边界：相等奇异值之间没有唯一顺序。只要相同 block 的 U/V 变换同步即可。

Example 3：宽、秩亏、零奇异值与符号不唯一

Step 1 — 选小 Gram：算 2×2 的 CC^T：

Step 2 — 特征对：

Step 3 — 奇异值：σ₁=√3，σ₂=1。full Σ 为：

Step 4 — 恢复 V 的非零 modes：

Step 5 — 零空间补 V：解 Cv₃=0：

不能使用 C^Tu/σ 计算这个 mode，因为对应 σ=0。

Step 6 — 组装：

Step 7 — 用 rank-one sum 验算：

符号不唯一：把任一非零 mode 的 u_i 与 v_i 同时乘 −1，外积不变。只改一侧则会改变 A。

4. 零与重复奇异值：统一规则

5. 三重验算模板

1. 正交：U^TU=I，V^TV=I。
2. 奇异方程：对正 mode，Av_i=σ_iu_i 与 A^Tu_i=σ_iv_i。
3. 重构：UΣV^T=A，或逐个 rank-one term 相加。

数值答案还可报告 ||U'*U-I||、||V'*V-I|| 与 ||U*S*V'-A||。手算题应展示至少一个非平凡元素或 rank-one sum，不能只写“verified”。

6. 常见扣分点

• 把 A^TA 的 eigenvectors 放进 U，而不是 V。
• 把 λ 直接当奇异值，忘记 σ=√λ。
• 排序 σ 后没有同步重排 U/V 列。
• 把矩形 Σ 写成方阵而不声明 thin SVD。
• 只求 singular values，未给 singular vectors。
• 对零奇异值套恢复公式，发生除零。
• 忘记补 null(A) 或 null(A^T)。
• 把 UΣV^T 写成 UΣV。
• 因同步换号或重复 singular subspace 的旋转误判答案。
• 只检查 dimensions，不检查 orthogonality 与 reconstruction。

7. 计时训练：从 8 分钟到 25 分钟

建议节奏：前 10% 时间写 shape 与路线；中间 70% 算 eigenpairs 与恢复；最后 20% 做正交/重构检查。若特征多项式耗时过长，先保住结构分：明确哪一 Gram、哪个向量属于 U/V、如何补零空间。

Ex.1 — 题面称“steepest descent”的二阶方法

术语边界：题面明确导出的更新是 Newton step。对非线性 least squares 用 J^TJ 近似 Hessian 通常称 Gauss–Newton。标准 steepest descent 只沿 −∇f。作答可沿题面措辞，但 discussion 应区分术语。

1. Jacobian 与 Hessian

若 g:ℝⁿ→ℝ^m，Jacobian 的元素为 (J_g)_ij=∂g_i/∂x_j。若 f:ℝⁿ→ℝ，Hessian 为 (H_f)_ij=∂²f/(∂x_i∂x_j)，充分光滑时对称。

二阶 Taylor model 对步长 s=x−x₀ 求导：

H 非正定时，stationary point 不保证是最小点；实际方法需要 damping、line search 或 trust region。

2. 应用与速度权衡

更新为 w_k+1=w_k+s_k，其中通过解 H_ks_k=−g_k 得到步长。不要显式求逆。对 d 个参数，dense factor/solve 典型代价为 Θ(d³)，存储为 Θ(d²)；局部收敛可能快，但每步昂贵。

3. Least squares 的 gradient 与 Hessian

令 r(w)=Xw−y，f(w)=(1/2)r(w)^Tr(w)，J=J_r：

题面说明 residual 是线性函数，因此每个 H_ri=0，故 H_f=J^TJ 是精确 Hessian，不只是近似。

4–6. 更新、算法与性能

在 full-column-rank 线性 least squares 中：

initialize w
repeat:
    r := residual(w)
    J := residual Jacobian
    g := Jᵀr
    choose stable solve:
        QR / SVD / damped normal equations
    solve for Δw; do not form an explicit inverse
    w := w + Δw
until gradient norm, loss change, or max iterations stops

线性 full-rank least squares 在精确算术中可一步到最优解；一般非线性问题每轮重新线性化。小 singular values 会放大噪声，可用 pseudoinverse、truncated SVD 或 damping。

Ex.2 — PRAM XOR reduction

1. Sequential

s := A[1]
for i := 2..n:
    s := s XOR A[i]
return s

使用 n−1 次 XOR，所以 work 与 depth 均为 Θ(n)。

2. Parallel binary tree

for round r := 0 .. log2(n)-1 in sequence:
    parallel for each block start i:
        A[i] := A[i] XOR A[i + 2^r]
return A[1]

第 r 轮执行 n/2^r+1 次 XOR。总 work：

轮数为 log₂n，故 depth 为 Θ(log n)。

3. Brent 质量评价

算法 work-optimal。当 p≲n/log n 时可接近线性 speedup；processor 再增加后，span log n 成为下界。

Ex.3 — Spark 中六种梯度下降变体

未提交：当前没有上述六种实现的可审计代码、日志或结果。本页不提供虚构的 benchmark 结论。