机器学习数学推导（二）：线性代数与矩阵论

为什么写这一章，有什么不同#

如果你学过标准的线性代数课程，大部分内容你可能已经见过。但本章不是对传统课程的简单复述，而是面向机器学习实践者，聚焦于实际场景中高频使用的线性代数核心概念，如实现梯度下降、运行 PCA、训练神经网络或研读论文时所需的内容。

具体来说，我的目标有三个：帮你建立矩阵操作的几何直觉，如旋转、拉伸、投影和压缩；讲清楚四个无处不在的矩阵分解方法——谱分解、SVD、QR 和 Cholesky，并告诉你何时使用哪种方法；把矩阵求导讲明白，让你能在草稿纸上轻松推导出任何神经网络的梯度公式。

我们省略了行化简、基于余子式展开的行列式计算以及抽象向量空间的相关证明。如果需要这些内容，文末的参考文献中有经典教材推荐。每个概念都通过一张示意图或一行 NumPy 代码直观呈现。

先修知识： 矩阵乘法、转置、行列式的概念、偏导数。

向量空间、子空间、秩——几何核心#

暂不引入向量空间的八条公理。在机器学习中，真正有用的心智模型是：向量空间是一张过原点的无限延伸的「平面」，子空间则是嵌套在其中的更小的过原点的「平面」。

$\mathbb{R}^2$ 中过原点的一条直线是一个 1 维子空间； $\mathbb{R}^3$ 中过原点的一个平面是一个 2 维子空间。关键是「过原点」这三个字——一旦平面发生平移，向量加法与数乘运算将不再满足封闭性。

张成、线性无关、基——一口气讲完#

随便挑几个向量 $v_1, \ldots, v_k$ 。它们所有线性组合 $\sum \alpha_i v_i$ 的集合叫张成 (span)。这个张成一定是个子空间。

\alpha_1 v_1 + \cdots + \alpha_k v_k = 0 \;\implies\; \alpha_1 = \cdots = \alpha_k = 0. \tag{1}

基就是线性无关的张成集。基的大小就是子空间的维数。

图已经说明了一切：线性无关的向量构成非退化平行四边形（面积非零，可张成二维子空间），而线性相关的向量共线，对应的平行四边形退化为线段。

矩阵就是线性映射#

每个矩阵 $A \in \mathbb{R}^{m \times n}$ 都对应一个线性映射 $T : \mathbb{R}^n \to \mathbb{R}^m$ ，反过来也一样。两者之间的关系很简单：

$$A$$ 的第 $$j$$ 列就是 $$T(e_j)$$ ——也就是第 $$j$$ 个标准基向量被映射到的位置。

所以当我看到神经网络某一层的权重矩阵 $W \in \mathbb{R}^{h \times d}$ 时，我看到的就是一个函数：「输入 $$d$$ 维，输出 $$h$$ 维」。 $$W$$ 的每一列告诉我某个输入特征对输出的贡献是什么。

秩与四个基本子空间#

对于一个秩为 $$r$$ 的 $m \times n$ 矩阵 $$A$$ ，它的四个基本子空间如下：

子空间	定义	所在空间	维数	直观含义
列空间 $\text{Col}(A)$	$\{Ax : x \in \mathbb{R}^n\}$	$\mathbb{R}^m$	$$r$$	输出能到达的地方
零空间 $\text{Null}(A)$	$\{x : A x = 0\}$	$\mathbb{R}^n$	$$n-r$$	被抹掉的输入
行空间 $\text{Row}(A)$	$\text{Col}(A^\top)$	$\mathbb{R}^n$	$$r$$	$$A$$ 「真正看到」的输入
左零空间	$\{y : A^\top y = 0\}$	$\mathbb{R}^m$	$$m-r$$	永远到不了的输出

正交性（基本定理）： $\text{Null}(A) \perp \text{Row}(A)$ 。

证明： 如果 $$Ax = 0$$ ，行空间中任一向量可以写成 $w = A^\top z$ ，那么 $\langle x, w\rangle = z^\top (Ax) = 0$ 。 $\square$

所以 $\mathbb{R}^n$ 被分成两块正交的子空间：行空间（ $$A$$ 在这里干活）和零空间（ $$A$$ 在这里抹掉信息）。当模型参数数量超过方程数量（这在现代机器学习中极为常见）时，零空间非平凡，即存在多个不同参数向量产生完全相同的预测结果。

范数与内积——几何层#

内积 $\langle \cdot, \cdot\rangle$ 提供角度信息，范数 $\|\cdot\|$ 提供长度信息。在 $\mathbb{R}^n$ 上，标准内积定义为 $\langle x, y\rangle = x^\top y$ ，对应的诱导范数是 $\|x\|_2 = \sqrt{x^\top x}$ 。

正则化问题中最常用的工具是 $\ell_p$ 范数族：

范数	公式	单位球形状	机器学习用途
$\ell_1$	$\sum_i \lvert x_i\rvert$	菱形	Lasso、稀疏性
$\ell_2$	$\sqrt{\sum_i x_i^2}$	圆形 / 球形	Ridge、权重衰减
$\ell_\infty$	$\max_i \lvert x_i\rvert$	方形 / 立方体	对抗扰动边界

Cauchy-Schwarz 与余弦#

定理（Cauchy-Schwarz）。 $|\langle u, v\rangle| \le \|u\| \cdot \|v\|$ ，等号成立当且仅当 $$u$$ 和 $$v$$ 平行。

证明： 考虑二次函数 $\|u + tv\|^2 = \|u\|^2 + 2t\langle u, v\rangle + t^2\|v\|^2$ ，对任意 $$t$$ 都非负，因此判别式必须 $\le 0$ 。 $\square$

将两边同时除以 $\|u\|\|v\|$ ，就得到夹角的余弦公式 $\cos\theta = \langle u, v\rangle / (\|u\|\|v\|)$ ，同时也证明了它的取值范围是 $$[-1, 1]$$ 。这个公式是检索系统中余弦相似度、注意力机制和对比学习的核心原理。

实际用到的矩阵范数#

范数	公式	衡量内容
Frobenius $\mid A\mid_F$	$\sqrt{\sum_{ij} a_{ij}^2}$	把矩阵 $$A$$ 当作向量时的欧几里得长度
谱范数 $\mid A\mid_2$	$\sigma_1(A)$	在所有单位输入上的最大拉伸
核范数 $\mid A\mid_*$	$\sum_i \sigma_i(A)$	矩阵秩 $\text{rank}(A)$ 的凸替代

谱范数是 Lipschitz 分析中的正确选择，比如 GAN 中的谱归一化。核范数则是低秩矩阵补全问题（如 Netflix 推荐系统）中真正需要最小化的目标。

特征值分解——「变换中方向不变的那些向量」#

定义与直观理解#

Av = \lambda v. \tag{2}

从几何上看，矩阵 $$A$$ 通常会旋转输入向量，但沿着特征向量的方向， $$A$$ 只会对向量进行缩放。换句话说，这些方向在变换后依然保持不变。

特征值是通过求解特征多项式 $\det(A - \lambda I) = 0$ 得到的。这里有两个重要的恒等式：

$\sum_i \lambda_i = \operatorname{tr}(A)$
$\prod_i \lambda_i = \det(A)$

谱定理（每天都会用到的那个）#

定理： 如果矩阵 $$A$$ 是实对称矩阵（即 $A = A^\top$ ），那么：

所有特征值都是实数；
不同特征值对应的特征向量彼此正交；
$$A$$ 可以正交对角化： $A = Q \Lambda Q^\top$ ，其中 $Q^\top Q = I$ ， $\Lambda$ 是对角矩阵。

证明思路（特征值为实数）。 假设 $Av = \lambda v$ 且 $v \neq 0$ 。分别用两种方法计算 $\bar v^\top A v$ ，并利用 $A = A^\top$ ，可以得到 $\lambda = \bar\lambda$ 。 $\square$

A = \sum_{i=1}^n \lambda_i\, q_i q_i^\top. \tag{3}

仔细看这个公式—— $$A$$ 实际上就是一组沿正交方向的投影，按特征值加权求和。这就是机器学习中看待对称矩阵的标准视角：协方差矩阵、 Gram 矩阵、 Hessian、图 Laplacian 都可以用这种方式理解。

正定与半正定矩阵#

如果对称矩阵 $$A$$ 满足对所有 $x \neq 0$ 都有 $x^\top A x > 0$ ，就称它是正定矩阵（记作 $A \succ 0$ ）。等价地说，所有特征值都为正。如果把条件放宽为 $x^\top A x \ge 0$ ，就称为半正定矩阵（记作 $A \succeq 0$ ）。

从几何上看，集合 $\{x : x^\top A x = 1\}$ 是一个椭球，主轴方向沿着特征向量，半轴长度为 $1/\sqrt{\lambda_i}$ 。如果某个特征值为负，结果就会变成双曲面（鞍面），并且不存在最小值。这就是为什么凸优化要求 Hessian 必须正定。

正定矩阵在机器学习中随处可见：

协方差矩阵 $\Sigma = \mathbb{E}[(x - \mu)(x - \mu)^\top]$ 总是半正定的；
Gram 矩阵 $X^\top X$ 半正定；当 $$X$$ 的列线性无关时，它是正定的；
核函数必须满足半正定性（Mercer 定理），否则无法在任何特征空间中复现；
凸损失函数的 Hessian 半正定；严格凸 ⇔ 正定。

奇异值分解 (SVD)——万能分解工具#

定理#

A = U \Sigma V^\top, \tag{4}

其中 $U \in \mathbb{R}^{m \times m}$ 和 $V \in \mathbb{R}^{n \times n}$ 是正交矩阵， $\Sigma$ 是对角矩阵，对角线上是奇异值 $\sigma_1 \ge \cdots \ge \sigma_r > 0$ ，其余位置为零。

为什么 SVD 总存在？ $A^\top A$ 是对称半正定矩阵，谱定理保证它有标准正交特征向量基 $$V$$ ，对应的特征值是 $\sigma_i^2$ 。对于 $\sigma_i > 0$ ，定义 $u_i = A v_i / \sigma_i$ ，然后补全成一组正交基 $$U$$ 。于是 $A v_i = \sigma_i u_i$ ，即 $AV = U\Sigma$ ，也就是 $A = U\Sigma V^\top$ 。 $\square$

几何直观：旋转、拉伸、再旋转#

任何线性变换都可以拆成三步：

$V^\top$ ：把输入旋转一下，让主轴对齐坐标轴；
$\Sigma$ ：沿每根坐标轴按 $\sigma_i$ 拉伸（这是唯一改变形状的一步）；
$$U$$ ：把结果旋转到输出空间。

这是线性代数里最简洁的模型：任何矩阵都可以看成「旋转 → 沿坐标轴拉伸 → 再旋转」。

秩不足会发生什么：可视化#

如果某个 $\sigma_i = 0$ ，对应方向上的输入就会被抹掉。输出椭球会塌缩成更低维的形状：

秩为 1 的矩阵 $A = u v^\top$ 把整个 $\mathbb{R}^n$ 都压到一条直线 $\text{span}(u)$ 上。所有与 $$v$$ 正交的向量都进入零空间，变成零。这就是低秩逼近能压缩数据的核心原因。

Eckart-Young 定理：最佳低秩逼近#

A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top, \tag{5}

误差分别是 $\|A - A_k\|_F = \sqrt{\sum_{i>k} \sigma_i^2}$ 和 $\|A - A_k\|_2 = \sigma_{k+1}$ 。

这是 PCA 的数学核心。把数据矩阵 $$X$$ 中心化后做 SVD，前 $$k$$ 个右奇异向量就是主成分方向，奇异值平方除以 $$n$$ 就是各方向的方差。这个方法也是图像压缩、潜在语义分析、推荐系统、 LoRA 等技术的基础。

伪逆与条件数#

Moore-Penrose 伪逆可以用 SVD 计算： $A^+ = V \Sigma^+ U^\top$ ，其中 $\Sigma^+$ 把非零奇异值取倒数。这样 $x^\star = A^+ b$ 就是 $$Ax = b$$ 的最小范数最小二乘解，也就是 numpy.linalg.lstsq 返回的结果。

条件数 $\kappa(A) = \sigma_1 / \sigma_r$ 表示 $$A$$ 对相对扰动的放大倍数。粗略来说，如果 $\kappa(A) \approx 10^k$ ，解 $$Ax = b$$ 时会丢失大约 $$k$$ 位精度。每次用最小二乘拟合之前，我都会先检查条件数。

特征分解 vs SVD#

	特征分解	SVD
适用矩阵	方阵	任意 $m \times n$
形式	$A = P \Lambda P^{-1}$	$A = U \Sigma V^\top$
是否总存在	仅当可对角化	总存在
「值」	可能复数	非负实数
左右因子互逆？	不一定	一定（ $$U, V$$ 正交）

对于对称半正定矩阵，特征分解和 SVD 是一样的（ $$U = V = Q$$ ， $\Sigma = \Lambda$ ）。其他情况下，我更喜欢用 SVD——它稳定、通用，还能处理长方矩阵。

矩阵求导——优化的核心#

大多数机器学习读者都希望看到一份清晰的参考。我用的是分子布局：导数的形状和输出一致。比如，标量对向量的梯度是列向量，和输入保持相同形状。

每天都会用到的四个公式#

#	$$f$$	$\nabla f$	备注
1	$a^\top x$	$$a$$	线性函数
2	$x^\top A x$	$(A + A^\top) x$	如果 $$A$$ 对称，则结果为 $$2Ax$$
3	$\mid A x - b\mid_2^2$	$2 A^\top (A x - b)$	最小二乘问题
4	$\ln \det X$	$X^{-\top}$	log-det，常见于高斯最大似然估计

公式 (2) 的证明
设 $f = \sum_{i,j} A_{ij} x_i x_j$ ，那么 $\partial f / \partial x_k = \sum_j A_{kj} x_j + \sum_i A_{ik} x_i = [(A + A^\top) x]_k$ 。 $\square$

A = \begin{pmatrix} 1 & 4 \\ 0 & 1 \end{pmatrix},\qquad x = \begin{pmatrix} 1 \\ 1 \end{pmatrix}.

计算 $f(x) = x^\top A x = 1\cdot 1\cdot 1 + 1\cdot 4\cdot 1 + 1\cdot 0\cdot 1 + 1\cdot 1\cdot 1 = 6$ 。如果误用 $\nabla f = 2Ax$ ，得 $2Ax = 2(5,\ 1)^\top = (10,\ 2)^\top$ 。但正确公式 $(A+A^\top)x = \begin{pmatrix}2 & 4\\ 4 & 2\end{pmatrix}\begin{pmatrix}1\\1\end{pmatrix} = (6,\ 6)^\top$ 。差距大到方向都不同。

数值验证： $f(x + (\varepsilon, 0)^\top) - f(x) \approx 6\varepsilon$ ，所以 $\partial f/\partial x_1 = 6$ ——正确公式给的也是 6，错误公式给的是 10，错。机器学习里 Hessian 几乎总是对称的（因二阶混合偏导对称），所以 ML 课本可以直接写 $$2Ax$$ ；一旦你在 RL、控制、二阶随机优化里碰到非对称 $$A$$ （如 $x^\top A x$ 表示有向博弈的成本），就必须老老实实用 $A + A^\top$ 。

公式 (3) 的证明
应用链式法则： $\nabla_x \tfrac12 \|Ax - b\|^2 = A^\top (Ax - b)$ 。乘以 2 就能抵消掉机器学习论文中通常省略的 $\tfrac12$ 。最优解满足正规方程 $A^\top A x = A^\top b$ 。

链式法则与反向传播#

\nabla_x L = J_g^\top \cdot \nabla_g f, \tag{6}

这里 $$J_g$$ 是 $$g$$ 的 $m \times n$ 雅可比矩阵。

\delta := \frac{\partial L}{\partial a} \odot \sigma'(z), \quad \frac{\partial L}{\partial W} = \delta\, x^\top, \quad \frac{\partial L}{\partial b} = \delta. \tag{7}

这就是整个算法的核心。 PyTorch 的 autograd 引擎所做的，就是在计算图上运行公式 (6)。

有限差分验证#

\hat g_i = \frac{f(x + \varepsilon e_i) - f(x - \varepsilon e_i)}{2 \varepsilon}.

取 $\varepsilon = 10^{-5}$ ，解析梯度和数值梯度的相对误差应接近 $10^{-6}$ 。第 7 节的代码正是这样做的。

数值分解——怎么选，什么时候用#

QR 分解解决最小二乘问题#

如果矩阵 $$A$$ 列满秩，可以分解为 $$A = QR$$ ，其中 $$Q$$ 的列正交， $$R$$ 是上三角矩阵。要解 $\min_x \|Ax - b\|$ ，先计算 $Q^\top b$ ，然后回代求解 $Rx = Q^\top b$ 。

为什么用 QR 而不是正规方程？ 正规方程需要计算 $A^\top A$ ，条件数会变成 $\kappa(A)^2$ 。比如 $\kappa(A) = 10^4$ ，那么条件数就变成了 $$10^8$$ ，误差预算直接翻倍平方了。 QR 分解直接处理 $$A$$ ，条件数还是 $\kappa(A)$ ，更稳定。

Cholesky 分解解决正定系统#

对于正定矩阵 $$A$$ ，可以用 Cholesky 分解写成 $A = LL^\top$ ，其中 $$L$$ 是下三角矩阵且对角元为正。计算复杂度是 $$O(n^3 / 3)$$ ，比 LU 分解快一倍。解 $$Ax = b$$ 就是两次三角回代。这是高斯过程推断、核方法以及任何带正定 Hessian 的二次规划的核心工具。

决策表#

场景	方法
对称矩阵（协方差、 Hessian）	特征分解 / 谱分解
任意矩阵；需要准确的秩、 PCA、伪逆	SVD
长方形最小二乘问题 $\min \mid Ax - b\mid$	QR
正定线性方程组 $$Ax = b$$	Cholesky
检测秩亏 / 监控数值稳定性	通过 SVD 计算 $\kappa(A)$

代码：验证所有结论#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
from scipy.linalg import svd, qr, cholesky

rng = np.random.default_rng(42)

m, n, true_rank = 100, 80, 10
U_t = rng.standard_normal((m, true_rank))
S_t = np.diag(np.linspace(10, 1, true_rank))
V_t = rng.standard_normal((n, true_rank))
A = U_t @ S_t @ V_t.T + 0.5 * rng.standard_normal((m, n))

U, s, Vt = svd(A, full_matrices=False)
print("rank-k 逼近：实际误差 vs 理论误差")
for k in [1, 5, 10, 20]:
    A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    actual = np.linalg.norm(A - A_k, "fro")
    theory = np.sqrt(np.sum(s[k:] ** 2))
    print(f"  k={k:2d}: 实际误差={actual:6.2f}  理论误差={theory:6.2f}")

m, n = 100, 10
A = rng.standard_normal((m, n))
x_true = rng.standard_normal(n)
b = A @ x_true + 0.1 * rng.standard_normal(m)

x_normal = np.linalg.solve(A.T @ A, A.T @ b)
Q, R = qr(A, mode="economic")
x_qr = np.linalg.solve(R, Q.T @ b)
print(f"\nkappa(A) = {np.linalg.cond(A):.1f},  "
      f"kappa(A^T A) = {np.linalg.cond(A.T @ A):.1f}")
print(f"正规方程误差: {np.linalg.norm(x_normal - x_true):.2e}")
print(f"QR 误差     : {np.linalg.norm(x_qr - x_true):.2e}")

n, d = 100, 5
X = rng.standard_normal((n, d))
w_true = rng.standard_normal(d)
y = X @ w_true + 0.1 * rng.standard_normal(n)

loss = lambda w: 0.5 / n * np.sum((X @ w - y) ** 2)
grad = lambda w: 1 / n * X.T @ (X @ w - y)

def grad_numeric(w, eps=1e-7):
    g = np.zeros_like(w)
    for i in range(len(w)):
        wp, wm = w.copy(), w.copy()
        wp[i] += eps; wm[i] -= eps
        g[i] = (loss(wp) - loss(wm)) / (2 * eps)
    return g

w_test = rng.standard_normal(d)
print(f"\n梯度检查: ||解析梯度 - 数值梯度|| = "
      f"{np.linalg.norm(grad(w_test) - grad_numeric(w_test)):.2e}")

运行结果里有三点值得注意： Eckart-Young 的实际误差和理论值在机器精度范围内完全一致；即使矩阵条件数不算差， QR 方法的精度仍然比正规方程高出一到两个数量级；解析梯度和数值梯度的差距在 $10^{-9}$ 量级。

常见问题#

Q1. 什么时候应该用正规方程？
当 $\kappa(A) < 10^3$ 时可以用，比如小规模问题且特征尺度均匀的情况。如果问题病态严重，就改用 QR 或 SVD。条件数平方后确实会大幅损失精度。

Q2. 为什么 $\ell_1$ 能产生稀疏解？
$\ell_1$ 的单位球在坐标轴上有尖角。优化过程中，等值线通常先碰到这些尖角，而尖角处某些坐标正好为零。相比之下， $\ell_2$ 的单位球是光滑的，最优解一般不会落在坐标轴上。

Q3. SVD 和 PCA 是一回事吗？
基本是一回事，只是计算方式不同。假设数据 $X \in \mathbb{R}^{n \times d}$ 已经中心化：

协方差法： $\Sigma = \tfrac{1}{n} X^\top X$ ，取前几个特征向量就是主方向；
SVD 法： $X = U \Sigma V^\top$ ， $$V$$ 的列就是主方向， $\sigma_i^2 / n$ 是各方向的解释方差。
对于高维数据（ $d \gg n$ ）， SVD 更高效，因为它不需要显式构造 $d \times d$ 的协方差矩阵。

Q4. 对称矩阵的特征值为什么特殊？
它们是实数，可以排序，而且特征向量构成正交基。这种正交性让我能将任意向量写成 $\sum \langle v, q_i\rangle q_i$ 的形式——这是描述 $$A$$ 行为的最简洁坐标系。机器学习里的 Hessian、协方差、 Gram、 Laplacian 都是对称矩阵，所以这一点非常重要。

Q5. 「秩」的实际意义是什么？
秩表示矩阵 $$A$$ 保留了多少输入信息。一个秩为 $$r$$ 的 $\mathbb{R}^{m \times n}$ 矩阵会丢掉 $$(n - r)$$ 维子空间（零空间），输出则限制在 $$r$$ 维子空间（列空间）内。在机器学习中，这就是自动编码器的瓶颈维度、矩阵补全的秩或嵌入向量的潜在维度。

总结#

工具	关键公式	什么时候用
谱定理	$A = Q \Lambda Q^\top$	对称矩阵：比如协方差、 Hessian
SVD	$A = U \Sigma V^\top$	任意矩阵： PCA、伪逆、低秩分解
QR	$$A = QR$$	病态条件下的最小二乘问题
Cholesky	$A = LL^\top$	正定线性系统（高斯 MLE、 GP）
二次型梯度	$\nabla (x^\top A x) = 2Ax$	所有二次目标函数
链式法则	$\nabla_x L = J_g^\top \nabla_g f$	一句话概括反向传播

下一步#

线性代数把"参数"和"梯度"具象成了向量与矩阵，但它本身只是确定性的代数。真实的数据是带噪声的，模型会犯错，估计需要置信区间——这些都是概率的语言。下一章我换一套词汇表：随机变量、期望、方差、似然、后验。

我会刻意把贝叶斯公式、大数定律、中心极限定理、MLE/MAP 这条主线从头到尾走一遍，因为后续每一章都默认你能在这几个概念之间自由切换。线性回归会用 MLE 解释为什么是最小二乘；逻辑回归会用 MLE 解释交叉熵；EM 会用 MLE 解释隐变量；正则化会用 MAP 解释惩罚项。所以这一章是"概率视角下的统计推断"——不是为了教概率，而是为了把每个 ML 算法都还原成"对某个似然函数做估计"的具体动作。

参考文献#

Strang, G. (2023). Introduction to Linear Algebra (6th ed.). Wellesley-Cambridge Press.
Trefethen, L. N. & Bau III, D. (1997). Numerical Linear Algebra. SIAM.——稳定性问题写得最干净。
Golub, G. H. & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.——算法圣经。
Petersen, K. B. & Pedersen, M. S. (2012). The Matrix Cookbook. Technical University of Denmark.——大家都在查的梯度速查表。
Eckart, C. & Young, G. (1936). The approximation of one matrix by another of lower rank. Psychometrika, 1(3), 211-218.
Boyd, S. & Vandenberghe, L. (2018). Introduction to Applied Linear Algebra. Cambridge University Press.——以数据为中心的优秀配套读物。