第5–6次作业：奇异值分解与并行优化

第1题：数值稳定性

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

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

知识点：IEEE 754单精度浮点数 float 约有7位有效十进制数字；双精度浮点数 double 约有15–16位。更高精度可以降低累加、相消、病态分解和长迭代中的舍入误差。更高精度也会增加内存、磁盘和网络用量，并可能降低缓存、单指令多数据（SIMD）或图形处理器（GPU）的吞吐量。

解题步骤：

1. 说明误差收益。敏感的归约和分解计算可以获得更可靠的结果。
2. 说明系统代价。在IEEE 754格式下，double 占64位，float 占32位，因此前者的单元素存储量是后者的两倍。
3. 说明适用边界。高精度不能修复噪声数据、错误模型或不稳定算法。
4. 根据条件数、误差预算和失败证据选择精度；若低精度计算达到吞吐目标但残差不达标，则用高精度执行残差校正。

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 的非零奇异值也相同。库实现、完整或紧致计算模式以及工作区大小会影响计时。

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-范数条件数从 κ(X) 变为 κ(X)²。计算更快不表示数值更稳定。

常见错误：只生成一个矩阵；只测量一次；把本机运行时间写成算法定律；不说明完整或紧致计算模式；排除形成Gram矩阵的时间但不声明。

1.3 1000次随机扰动

题面事实：对给定的5×5非对称矩阵 X，分别研究 X+δX 的特征值和奇异值在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 抵消

知识点：非对称矩阵的特征值可能是复数。比较不同试验的结果之前，必须固定排序规则或匹配规则。奇异值非负，因此可以按降序匹配。奇异值满足扰动界 |σ_i(X+E)−σ_i(X)|≤‖E‖₂。

常见错误：把确定性的逐元素 eps(X) 称为随机扰动；直接比较未匹配的特征值；累加带符号差值而使正负误差抵消；不说明扰动尺度。

1.4使用第4章解释实验

采用非正交基的特征分解 M=PDP⁻¹，其局部扰动放大与 ‖P‖‖P⁻¹‖ 有关。SVD的左右因子正交，因此对应的2-范数条件因子为1。SVD数值较稳定，不表示所有SVD计算方法都最快。通过Gram矩阵求SVD会损失一部分稳定性优势。

第2题：不用计算机求完整奇异值分解

题面矩阵：

题意：题目要求奇异值分解，不是只求奇异值。答案必须给出 U,Σ,V^T。答案必须说明使用完整SVD还是紧致SVD。答案还必须写出计算步骤。

计算步骤：因为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 并归一化，以补齐完整的 V∈ℝ^4×4。

下一节给出完整手算步骤和三个小矩阵例题。对于第5次作业中的矩阵，最后必须检查 U^TU=I₃、V^TV=I₄ 和 UΣV^T≈X。

第3题：Spark中的PCA

题面：两个没有表头的传感器CSV文件中，至少一个文件可能包含1001个电路传感器输出。最后一列 y 表示每小时用电量。需要说明PCA的作用。需要找出累计解释90% 方差的最小主成分数 n。需要使用前 n 个主成分建立 y 的模型。还需要判断 sensors2 是否属于同一类数据。

3.1 PCA如何帮助

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

3.2解释90% 所需列数

先分离最后一列 y。只在特征矩阵上拟合缩放器和PCA。若解释方差比为 e₁≥e₂≥…，取满足下式的最小 n：

常见错误：把“每个主成分的解释方差比（EVR）大于0.01”当成累计90% 规则；把 y 输入PCA，造成标签泄漏。

3.3使用主成分得分进行回归

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)

处理流程：先切分数据。只使用训练集拟合缩放器、选择 n、拟合PCA和拟合回归模型。测试集只调用 transform。使用 test_predictions 计算留出测试集的决定系数（R²）和均方根误差（RMSE）。该流程不会把测试集分布泄漏到特征空间。

模型写为 y=β₀+Σβ_ip_i+ε。题面给出的 ε∼Normal(0,1) 是模型假设。必须检查残差，不能把该假设当作已经证明的事实。报告留出测试集的R² 和RMSE。还要报告数据切分方式、随机种子和残差诊断结果。

3.4判断sensors2

对两份文件使用相同的预处理、数据切分和评价方法。比较累计方差谱、主成分得分回归模型在留出测试集上的表现、残差结构和基线结果。单个较高的R² 或相似的二维散点图不能证明两个文件来源相同。

1.先写矩阵形状

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

如果题目要求分解，但没有说明形式，则给出完整SVD。也可以明确声明使用紧致SVD。

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. 补齐完整正交基。使用 null(A) 补齐V。使用 null(A^T) 补齐U。对补充向量进行正交归一化。
7. 检查结果。检查列的单位正交性、奇异向量方程和矩阵重构结果。

形状 → 较小的 Gram 矩阵 → 特征对 → σ 降序
     → 根据正奇异值计算另一侧向量
     → 使用零空间补齐正交基
     → UᵀU / VᵀV / UΣVᵀ 三重验算

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

• A^TA∈ℝ^n×n 的特征向量是V的列，位于输入空间。
• AA^T∈ℝ^m×m 的特征向量是U的列，位于输出空间。
• 两个矩阵的正非零特征值相同，均为 σ_i²。

判断方法：A^TA 的形状对应输入维度 n，所以它给出V。AA^T 的形状对应输出维度 m，所以它给出U。

例1：有两个不同正奇异值的方阵

步骤1：计算Gram矩阵。

步骤2：计算特征多项式。

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

步骤3：计算右奇异向量。

步骤4：计算奇异值和左奇异向量。

步骤5：组装三个矩阵。

步骤6：检查重构结果。

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

例2：重复奇异值和非唯一旋转

步骤1：计算特征值和奇异值。B^TB=4I₂，所以 λ₁=λ₂=4，σ₁=σ₂=2。

步骤2：选择一组奇异向量。可直接采用以下构造：

步骤3：说明结果不唯一。对任意二维正交矩阵Q，例如旋转矩阵：

结论：重复奇异值对应的左右奇异子空间可以进行相同的正交旋转。不同的U和V可能都正确。使用正交性和重构结果判断答案。

排序规则：相等奇异值之间没有唯一顺序。同一个重复值分块内的U和V必须同步变换。

例3：宽矩阵、秩亏、零奇异值和符号不唯一

步骤1：选择较小的Gram矩阵。计算2×2的 CC^T：

步骤2：计算特征对。

步骤3：计算奇异值。σ₁=√3，σ₂=1。完整 Σ 为：

步骤4：计算V中对应正奇异值的向量。

步骤5：使用零空间补齐V。求解 Cv₃=0：

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

步骤6：组装三个矩阵。

步骤7：使用秩一矩阵之和检查结果。

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

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

5.三重验算模板

1. 正交：U^TU=I，V^TV=I。
2. 奇异向量方程：对每个正奇异值，检查 Av_i=σ_iu_i 和 A^Tu_i=σ_iv_i。
3. 重构：检查 UΣV^T=A，或逐个相加秩一矩阵。

数值答案还可以报告 ||U'*U-I||、||V'*V-I|| 和 ||U*S*V'-A||。手算题应展示至少一个非平凡元素或一组秩一矩阵之和。不能只写“已验证”。

6.常见错误

• 把 A^TA 的特征向量放进U，而不是V。
• 把 λ 直接当奇异值，忘记 σ=√λ。
• 排序 σ 后没有同步重排U/V列。
• 把矩形 Σ 写成方阵，但不声明使用紧致SVD。
• 只求奇异值，未给出奇异向量。
• 对零奇异值套恢复公式，发生除零。
• 忘记补 null(A) 或 null(A^T)。
• 把 UΣV^T 写成 UΣV。
• 因为同步换号或重复奇异子空间的旋转而误判答案。
• 只检查矩阵维度，不检查正交性和重构结果。

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

时间分配：使用前10% 的时间写矩阵形状和计算步骤。使用中间70% 的时间计算特征对和另一侧向量。使用最后20% 的时间检查正交性和重构结果。如果特征多项式计算耗时过长，先说明选择哪个Gram矩阵、每个向量属于U还是V，以及如何补齐零空间。

第1题：题面称为“最速下降”的二阶方法

术语说明：题面推导的更新是Newton步长。非线性最小二乘中，以 J^TJ 近似Hessian矩阵的方法称为Gauss–Newton方法。标准最速下降只沿 −∇f 方向更新。作答时可以沿用题面措辞，但讨论部分必须区分这些术语。

1. Jacobian矩阵与Hessian矩阵

Jacobian矩阵表示向量值函数的一阶偏导数。若 g:ℝⁿ→ℝ^m，其元素为 (J_g)_ij=∂g_i/∂x_j。Hessian矩阵表示标量值函数的二阶偏导数。若 f:ℝⁿ→ℝ，其元素为 (H_f)_ij=∂²f/(∂x_i∂x_j)。函数充分光滑时，Hessian矩阵对称。

对二阶Taylor模型关于步长 s=x−x₀ 求导：

如果H不是正定矩阵，驻点不一定是最小点。实际算法可以使用阻尼、线搜索或信赖域处理该情况。

2.应用与速度权衡

更新式为 w_k+1=w_k+s_k。通过求解 H_ks_k=−g_k 得到步长。不得显式计算逆矩阵。对于 d 个参数，稠密分解和求解的典型时间复杂度为 Θ(d³)，存储复杂度为 Θ(d²)。该方法可能具有较快的局部收敛速度，但每一步的成本较高。

3.最小二乘的梯度与Hessian矩阵

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

题面说明残差是线性函数，因此每个 H_ri=0。所以 H_f=J^TJ 是精确的Hessian矩阵，不是近似值。

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

对于满列秩的线性最小二乘问题：

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

满秩线性最小二乘在精确算术中可以一步到达最优解。非线性问题需要在每一轮按当前迭代点重新线性化。较小的奇异值会放大噪声。可以使用伪逆、截断SVD或阻尼处理该问题。

第2题：PRAM异或归约

1.顺序算法

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

算法执行 n−1 次异或（XOR）操作，所以工作量和计算深度均为 Θ(n)。

2.并行二叉树

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。总工作量为：

轮数为 log₂n，所以计算深度为 Θ(log n)。

3.使用Brent定理评价算法

该算法达到渐近最优工作量。当 p≲n/log n 时，加速比可以接近处理器数量。继续增加处理器后，log n 的计算跨度成为时间下界。

第3题：Spark中的六种梯度下降变体

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