奇异值分解、低秩近似与主成分分析

特征向量在同一空间比较方形映射输入输出。矩形设计、图像矩阵、嵌入表连接不同空间；即使方形非对称也可能没有标准正交特征基。奇异值分解为每个有限实矩阵提供标准正交输入、输出方向，用非负强度连接二者，推广对称谱几何而不声称任意特征向量正交。

本课推导分解、伪逆、最佳秩近似，再用于中心化主成分分析。你将区分完整、约化、紧致形状，计算舍弃能量，并只在训练记录拟合预处理。高解释方差衡量特定几何目标，不证明预测用途。

检索检查：写实对称谱定理，解释 AᵀA 为何半正定。复习模块 13。投影、最小范数、四空间见模块 12及11。本课 PCA 散布计算为代数，概率统计解释留到模块 23。实验见 NumPy 准备。

m×n 的 A 使 AᵀA 为 n×n 对称半正定。谱定理给标准正交基 v₁,…,v_n，非负值从大到小排序。令 σᵢ 为其非负平方根，正 σᵢ 时定义 uᵢ=Avᵢ/σᵢ。uᵢᵀuⱼ=vᵢᵀAᵀAvⱼ/(σᵢσⱼ)，i=j 为一，其他零，故非零像方向标准正交。零值使 ‖Avᵢ‖²=0，所以 Avᵢ=0。

用完整 v 基展开输入，Ax=Σ正值σᵢuᵢ(vᵢᵀx)，证明 A=Σ_{i=1}^rσᵢuᵢvᵢᵀ，r 为正值数。u 张成列空间，对应 v 张成行空间，其余 v 张成核。把 u 扩成 Rᵐ 正交基，得到奇异值分解 A=UΣVᵀ，即 SVD。存在来自对称定理，A 本身可以矩形或亏损。

完整 SVD：U 为 m×m，V 为 n×n，矩形对角 Σ 为 m×n。前 q=min(m,n) 对角放奇异值，秩小于 q 时含零；其他项零。完整 U,V 描述全部输出输入基，包括未使用方向。m≠n 时不能拿 n×n 对角中因子与 m×m U 相乘，中间矩形形状不可忽略。

约化 SVD 保留 q 列：U_q 为 m×q，V_q 为 n×q，中对角 q×q，r<q 时仍含零。秩 r 紧致 SVD 只保留正方向：U_r 为 m×r,V_r 为 n×r,Σ_r 为 r×r。秩亏时约化不等于紧致，库用语可不同，应写形状及是否含零。

秩零紧致正列表为空，和为零矩阵，完整或约化仍可有任意零方向正交基。重复正值允许对应 u,v 同步正交旋转；同步反转一对符号不改外积。不能直接要求有符号列相同。

通过 AᵀA 的构造证明存在，不要求数值算法形成格拉姆矩阵；平方弱强度可能损失分辨率。库 SVD 用专门算法，应查形状、正交、重构，决定有效零时记录 dtype 和容差。

正交输入输出换基使 A 变 PᵀAQ，格拉姆变 QᵀAᵀAQ，特征值相同，因此奇异值、秩、两矩阵范数不变，方向坐标变。全映射乘 α，使强度乘 |α|；只缩放一特征一般不均匀改变，因为不是正交换基。这区分保欧几里得几何的重定向与改变度量的单位换算。

完整公式三步：Vᵀ 将输入表达为正交坐标并保长度，Σ 对匹配坐标非负缩放并按矩形需要插入或移除坐标，U 将输出坐标表达为正交输出基并保长度。正交因子可为旋转或反射，不能只说旋转而遗漏行列式负一。中强度非负，不同于有符号或复特征值。

正对满足 Avᵢ=σᵢuᵢ、Aᵀuᵢ=σᵢvᵢ，第二来自 AᵀAvᵢ=σᵢ²vᵢ。所以 u 也是 AAᵀ 的值 σᵢ² 特征向量。两个格拉姆共享正值，零数可因维度不同而不同；是同一映射输入、输出方向强度，不说 A 普通特征值相同。

r<n 时输入方向消失，完整 V 仍保其坐标；r<m 时输出像为真子空间，部分方向不可达，完整 U 后面方向张成左零。前 r 右向量张成行、后面张成零；前 r 左张成列、后面张成左零。四空间同时可见。

单位 x=Vc，‖Ax‖²=Σ正值σᵢ²cᵢ²≤σ₁²，第一右方向达到，故算子范数 ‖A‖₂=σ₁，零映射范数零。方阵可逆时最小值为最小伸缩，逆范数 1/σ_n，所以 κ₂=σ₁/σ_n，不需对称，修正无条件特征值比。

满列秩矩形 X 的 σ₁/σ_d 衡量输入最大、最小伸缩及伪逆范数；XᵀX 值 σᵢ²，所以 κ₂(XᵀX)=(σ₁/σ_d)²，证明先前平方机制。秩亏有零输入方向，在全输入空间没有有限逆敏感界；若只谈正方向比，必须声明限制。

实 NumPy s 一维，U 左向量在列，Vh 右向量在行。约化形状匹配时 U*s 广播缩放 U 列，缩放行是另一矩阵。复输出 Vh 为共轭转置，本课主要实例用普通转置。预期形状的重构能先发现方向错误。

紧致 A=U_rΣ_rV_rᵀ，定义 A⁺=V_rΣ_r⁻¹U_rᵀ，只逆正强度，形状 n×m。先取 u 坐标投影目标到列，再除强度，构造行空间输入，恰为模块 12 几何最小范数映射。普通逆是满秩方阵特例，不是公式前提。

乘得 AA⁺=U_rU_rᵀ、A⁺A=V_rV_rᵀ，为列和行投影，对称幂等。因此 AA⁺A=A、A⁺AA⁺=A⁺，得到全部 Moore–Penrose 条件。w_min=A⁺y 的预测是唯一最近列投影，全部最小化者 w_min+n，核与行垂直保证最小范数唯一。

保留方向 σᵢ 时，εuᵢ 目标扰动导致 (ε/σᵢ)vᵢ 参数变化，长度分别 |ε| 与 |ε|/σᵢ。这是精确方向敏感，说明给定数据解得准仍可对测量扰动不稳定。

diag(1,10⁻¹⁰) 对 (1,10⁻¹⁰) 给 (1,1)，加 (0,10⁻⁸) 后为 (1,101)。丢弃弱方向的数值伪逆给 (1,0)，留下第二残差。它主动改变有效逆以限制放大，不证明精确强度零。分辨率、拟合、参数响应之间有取舍，须声明容差目的。

固定伪逆接口按最大奇异值相对阈值。绝对秩阈值、相对伪逆阈值、测量噪声策略不自动可互换。全局缩放对绝对与相对不同，单列缩放又改变几何。报告保留方向、单位、dtype、容差；临近阈值小变化可改变选秩与代表。

约化列表含零，不能盲算 1/s；之后乘零也可能由无穷中间量产生 NaN。先建保留掩码或用明确库策略。全零伪逆为零，没有正方向可除。应验证选定投影与最小范数，而不只看调用无异常。

保留前 k 强项得 A_k=Σ_{i=1}^kσᵢuᵢvᵢᵀ，0≤k≤r，保留正值时秩 k，A₀=0。形状仍原矩阵。存因子可在 k 小时节省，标量数 k(m+n+1)，应与 mn 及实际编码开销比较；低秩不自动使文件小。

弗罗贝尼乌斯范数 ‖M‖F=√ΣᵢⱼMᵢⱼ²，等于各列长度平方和根。左正交乘保各列，右正交可对转置用同理保总和。因此在 SVD 坐标，差只含舍弃对角，‖A−A_k‖F²=Σ{i>k}σᵢ²，算子误差是最大舍弃 σ{k+1}，全部正项保留后零。

最佳秩 k 定理：秩至多 k 的任意 B 都不能使弗罗贝尼乌斯或算子误差小于 A_k。完整平方证明：P 投影到 C(B)，维 t≤k。逐列误差分内外成分，‖A−B‖_F²≥‖(I−P)A‖_F²=‖A‖_F²−‖PA‖_F²。左奇异基给后项 Σσᵢ²wᵢ，其中 wᵢ=‖Puᵢ‖²∈[0,1]。包含完整左基，对未用方向赋零强度，权重和 t，因为正交基投影平方和计 t 个单位方向。

加权强度和不超过前 k 平方和。k>0 时尾强度≤σ_k²，剩余预算≤k−Σ前k wᵢ，因此至多 Σ前k[σᵢ²wᵢ+σ_k²(1−wᵢ)]≤Σ前kσᵢ²。k=0 投影零。每候选误差至少尾平方，A_k 达到。论证覆盖任意候选像，而不是预设已经对齐。

算子证明：k<r，将 B 限制到前 k+1 右方向张成，其秩至多 k，由秩–零化度有单位 z 且 Bz=0。于是 ‖(A−B)z‖=‖Az‖≥σ_{k+1}，因为包含强度都不小于此值，截断达到。边界重复值时多个子空间可最优，最佳误差不保证唯一近似。

X 有 m 行记录、d 列特征，训练行求均值 μ，Z=X−1μᵀ，每列和零。右奇异向量为主特征方向。保留 V_k，得分 T=ZV_k 为 m×k，重构 X_hat=TV_kᵀ+1μᵀ。均值减加使原记录重构为仿射，其中心位移投影线性。

散布 ZᵀZ 的方向 V，值 σᵢ²。m>1 时除 m−1 为惯常样本协方差，值 σᵢ²/(m−1)。此处仅代数归一约定，分布估计后讲。TᵀT=diag(σ₁²,…,σ_k²) 解释中心坐标不相关，不声称统计独立。

解释比例 σᵢ²/Σσⱼ²，保留和是中心训练平方散布的保留份额。全记录相同使 Z=0，分母零，比例未定义，不是自动完美。单记录样本协方差分母也零，应声明边界而非悄悄 NaN 或改定义。

PCA 选正交 k 维特征子空间最小化中心重构。任意正交 W 的投影 ZWWᵀ，保留能量为所选 ZᵀZ 瑞利值和。最佳秩的投影权重证明界为前 k 平方强度，V_k 达到。得分重构恰 Z_k=U_kΣ_kV_kᵀ，所以最大方差、最小平方误差是同一目标两种代数视角。

中心化 1ᵀZ=0，使秩≤min(d,m−1)，不只是 min(d,m)。小训练集尤重要，宽数据中心后不超过 m−1 正主成分。重复散布值使单方向不唯一，整体空间仍明确；符号反转同步反转得分，重构不变。

新数学列 x 用 t=V_kᵀ(x−μ)，长 k；存储行用 (x−μ)V_k。复用训练参数。新残差不必零、留出得分均值不必零，均不表示失败，而是在训练固定坐标中测量。每批测试重拟合会改变变换与评估含义。

对一个新位移 z，重构 V_kV_kᵀz，舍弃成分垂直，误差平方为舍弃坐标平方和。训练尾强度给训练总误差，不自动预测新点误差；新记录可主要在保留空间外，即使训练重构优秀。任务关心留出重构时，应单独评估。

PCA 依赖单位，因为欧几里得重构加坐标平方。某特征乘百，其散布贡献乘万，交叉项也变；同一物理测量可能改变最大方向。实验第二尺度改变使首方向近第二轴并改变解释比例，是度量改变，不自动是计算错。

按训练均值尺度标准化是另一建模选择，零方差不可除零，应移除或声明策略。尺度、特征选择也属拟合变换，须只从训练估计；中心、尺度、基要一起保存以计算新分数并逆到原单位。

泄漏指保留用于评估的信息影响拟合模型或预处理。PCA 不用标签，但训练和留出特征合拟均值、方向仍改变过程。通常问题要求未用这些记录的训练变换。实验合入两条远测试使均值 (13,.5) 变约 (58.666667,5.333333)，尚未重算方向就已是不同变换。

操作约定：先分记录，训练拟合中心、尺度、成分，保存后应用验证测试；验证选择 k 等，最终保留测试评价所选过程。交叉验证每折需自己的训练变换，不能把全部数据拟合组件复用于各折，否则访问留出特征。若有意采用允许未标记评估特征的传导式任务，应明确不同设置，不当普通归纳评估。

未中心化 SVD 也不同，围绕原点近似原向量，混合均值位置与变化。离零很远时首方向可能主要表达共同位置而不是均值周围偏差。有明确应用时可用，但不默默叫中心 PCA。实验比较线投影以区分真实方向变化与任意符号。

低秩支持图像、嵌入压缩和紧凑参数更新。定理保证指定范数最近矩阵，不保证下游最佳、最快实现、最小文件。两因子各 k 列的 UVᵀ 秩至多 k，但学习选择这种更新还引入优化统计问题，应分开代数秩与任务等价。

最后须声明分辨率。重复或近重复强度可使单方向不稳定，而整体子空间重构较稳定；弱逆放噪，丢弃改变有效算子。记录 k、保留强度、训练中心尺度、dtype、容差及误差，使近似可复现、边界可判断，不将成功分解当所有用途证据。

采用数组准备的 Python 3.11 和 NumPy。输出为运行可下载脚本捕获。先预测形状、方向与预算；符号可有无害差别，应比重构和子空间。

实验 A 生成矩阵压缩

由三组正交方向及六、三、一生成 8×6 矩阵。验证生成方向，比较约化尺寸与精确秩，每个 k 先预测舍弃误差及标量数。

下载 lab1_truncation.py

"""Generated low-rank matrix with known spectrum and measured truncation error."""
import numpy as np
left = np.column_stack((np.ones(8), [1,-1,1,-1,1,-1,1,-1], [1,1,-1,-1,1,1,-1,-1]))/np.sqrt(8)
right = np.column_stack((np.ones(6)/np.sqrt(6), np.array([1,-1,1,-1,1,-1])/np.sqrt(6), np.array([1,1,-2,1,1,-2])/np.sqrt(12)))
assert np.allclose(left.T@left, np.eye(3)) and np.allclose(right.T@right, np.eye(3))
A = (left*np.array([6., 3., 1.]))@right.T
U, s, Vh = np.linalg.svd(A, full_matrices=False)
print("A / reduced U / s / Vh shapes:", A.shape, U.shape, s.shape, Vh.shape)
print("Singular values:", [round(float(v), 9) for v in s])
print("Original scalar entries:", A.size)
assert np.allclose(s[:3], [6,3,1]) and np.allclose(s[3:], 0)
for k in range(4):
    approximation = (U[:, :k]*s[:k])@Vh[:k, :]
    error_squared = np.linalg.norm(A-approximation, "fro")**2
    tail_squared = s[k:]@s[k:]
    scalar_storage = k*(A.shape[0]+A.shape[1]+1)
    print("Rank / squared Frobenius error / discarded sum / factor scalars:", k, round(float(error_squared), 9), round(float(tail_squared), 9), scalar_storage)
    assert np.allclose(error_squared, tail_squared)
assert np.allclose((U*s)@Vh, A)
print("Scalar-count savings are distinct from file-format overhead or predictive usefulness.")

A / reduced U / s / Vh shapes: (8, 6) (8, 6) (6,) (6, 6)
Singular values: [6.0, 3.0, 1.0, 0.0, 0.0, 0.0]
Original scalar entries: 48
Rank / squared Frobenius error / discarded sum / factor scalars: 0 46.0 46.0 0
Rank / squared Frobenius error / discarded sum / factor scalars: 1 10.0 10.0 15
Rank / squared Frobenius error / discarded sum / factor scalars: 2 1.0 1.0 30
Rank / squared Frobenius error / discarded sum / factor scalars: 3 0.0 0.0 45
Scalar-count savings are distinct from file-format overhead or predictive usefulness.

交付：秩、形状、误差、存储表及尾公式证明。舍入打印零不证明精确秩，本例精确主张来自已知生成构造。

实验 B 训练中心化 PCA

手算均值、散布，比谱与平方奇异值，重构一成分训练数据。对两留出批次复用参数不重拟，核对形状及总尾误差。

下载 lab2_training_pca.py

"""Fit mean and principal directions on training records; reuse on held-out rows."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
training = np.array([[10.,0.], [12.,1.], [14.,0.], [16.,1.]])
test = np.array([[20.,2.], [8.,-1.]])
mean = training.mean(axis=0)
centered = training-mean
U, s, Vh = np.linalg.svd(centered, full_matrices=False)
directions = Vh.T
scores = centered@directions[:, :1]
reconstruction = scores@directions[:, :1].T+mean
test_scores = (test-mean)@directions[:, :1]
scatter = centered.T@centered
variances = s**2/(len(training)-1)
ratio = s**2/np.sum(s**2)
assert np.allclose(centered.mean(axis=0), 0)
assert np.allclose(scatter@directions, directions*s**2)
print("Training mean:", mean)
print("Principal directions in columns:", directions.tolist())
print("Sample covariance eigenvalues:", variances)
print("Explained variance ratios:", ratio)
print("Training rank-one reconstruction:", reconstruction.tolist())
print("Held-out scores using training transform:", test_scores.ravel())
error_squared = np.linalg.norm(training-reconstruction, "fro")**2
print("Rank-one reconstruction squared error:", error_squared)
assert np.allclose(error_squared, s[1]**2)
changed_test = np.array([[2000.,200.], [-800.,-100.]])
changed_scores = (changed_test-mean)@directions[:, :1]
print("Changed held-out scores using the same training transform:", changed_scores.ravel())
print("Mean reused for both held-out batches:", mean)
print("All-zero centred data variance ratio: undefined; guard division by zero.")

Training mean: [13.   0.5]
Principal directions in columns: [[0.9946229901291114, 0.10356209493161814], [0.10356209493161814, -0.9946229901291114]]
Sample covariance eigenvalues: [6.736081 0.263919]
Explained variance ratios: [0.962297 0.037703]
Training rank-one reconstruction: [[9.980672702257417, 0.18562172467192445], [12.062227727769086, 0.4023573132283918], [13.937772272230914, 0.5976426867716081], [16.019327297742585, 0.8143782753280755]]
Held-out scores using training transform: [ 7.117704 -5.128458]
Rank-one reconstruction squared error: 0.7917560805262006
Changed held-out scores using the same training transform: [1996.976519 -819.036482]
Mean reused for both held-out batches: [13.   0.5]
All-zero centred data variance ratio: undefined; guard division by zero.

交付：保存均值和方向、训练得分重构、留出分数、解释比例。说明零散布未定义策略及训练、新记录误差区别。

实验 C 中心、泄漏、尺度与弱逆

比较中心与未中心线投影、算泄漏均值，改变单位解释主线旋转，比较小扰动下两个相对伪逆阈值。

下载 lab3_centering_and_pinv_faults.py

"""Uncentred directions, evaluation leakage, feature units and weak inverse modes."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
training = np.array([[10.,0.], [12.,1.], [14.,0.], [16.,1.]])
test = np.array([[100.,10.], [200.,20.]])
mean = training.mean(axis=0)
centered = training-mean
_, _, Vh = np.linalg.svd(centered, full_matrices=False)
_, _, uncentered_Vh = np.linalg.svd(training, full_matrices=False)
leaked_mean = np.vstack((training, test)).mean(axis=0)
print("Correct training mean / leaked combined mean:", mean, leaked_mean)
print("Centred / uncentred leading line projectors:", np.outer(Vh[0], Vh[0]).tolist(), np.outer(uncentered_Vh[0], uncentered_Vh[0]).tolist())
assert not np.allclose(mean, leaked_mean)
assert not np.allclose(np.outer(Vh[0], Vh[0]), np.outer(uncentered_Vh[0], uncentered_Vh[0]))
rescaled = centered*np.array([1.,100.])
_, scaled_s, scaled_Vh = np.linalg.svd(rescaled, full_matrices=False)
print("After feature two multiplied by 100, leading projector:", np.outer(scaled_Vh[0], scaled_Vh[0]).tolist())
print("Scaled leading explained ratio:", scaled_s[0]**2/(scaled_s@scaled_s))

A = np.diag([1., 1e-10])
target = np.array([1., 1e-10])
perturbed = target+np.array([0.,1e-8])
for cutoff in [1e-12,1e-8]:
    inverse = np.linalg.pinv(A, rcond=cutoff)
    base = inverse@target
    changed = inverse@perturbed
    print("Relative cutoff / base / changed weights / new residual:", cutoff, base, changed, format(np.linalg.norm(perturbed-A@changed), ".6e"))
assert np.allclose(np.linalg.pinv(A, rcond=1e-12)@perturbed, [1,101])
assert np.allclose(np.linalg.pinv(A, rcond=1e-8)@perturbed, [1,0])
print("Dropping a weak singular direction trades resolution and residual for reduced amplification; it changes the effective inverse.")

Correct training mean / leaked combined mean: [13.   0.5] [58.666667  5.333333]
Centred / uncentred leading line projectors: [[0.9892748924933744, 0.10300524052492092], [0.10300524052492092, 0.010725107506625487]] [[0.998380120869607, 0.04021511061772047], [0.04021511061772047, 0.0016198791303925755]]
After feature two multiplied by 100, leading projector: [[0.000401121600379904, 0.02002400314226935], [0.02002400314226935, 0.9995988783996206]]
Scaled leading explained ratio: 0.998403833102307
Relative cutoff / base / changed weights / new residual: 1e-12 [1. 1.] [  1. 101.] 1.654361e-24
Relative cutoff / base / changed weights / new residual: 1e-08 [1. 0.] [1. 0.] 1.010000e-08
Dropping a weak singular direction trades resolution and residual for reduced amplification; it changes the effective inverse.

交付：识别每种改变的目标或评估过程，报告逆分辨率取舍及实际残差，不把主动截断说成精确方向消失。

1–12 必做，13–14 额外 35 分钟。实因子应说明秩、保留维数及中心假设。

作者 Mathematics for Machine Learning 配套网站的 SVD、伪逆、PCA 选段，及官方 NumPy 1.26 SVD、伪逆形状与阈值约定。原创推导证明几何与近似。

不看笔记从 AᵀA 推方向，写三形状与四空间，证明一个最佳秩界，给伪逆与训练 PCA 流程及标签反例。

结业任务：diag(4,2,0) 秩二，秩一平方误差四、算子二，精确伪逆 diag(1/4,1/2,0)，目标第三项留残差。解释约化仍三项、零不可逆，以及阈值去掉二会改变什么。

进入下一课：区分精确几何、数值分辨率、重构与任务评估。下一块用线性代数发展极限、导数与优化。开放状态见课程概览。