线性代数深度指南

文档概述

线性代数是现代数学与机器学习的基石。本指南从向量空间出发，系统涵盖矩阵运算与分解、特征值理论、最小二乘法以及矩阵微积分等核心内容，并深入探讨其在机器学习中的广泛应用。

关键词

序号	关键词	英文	核心概念
1	向量空间	Vector Space	$({\mathbb{R}}^n,+,\cdot )$
2	基	Basis	$\text{span}\{v_1,\dots ,v_n\}=V$
3	矩阵分解	Matrix Decomposition	$A=U\Sigma V^T$
4	特征值	Eigenvalue	$A\mathbf{v}=\lambda \mathbf{v}$
5	特征向量	Eigenvector	满足 $A\mathbf{v}=\lambda \mathbf{v}$ 的非零向量
6	正交性	Orthogonality	$⟨\mathbf{u},\mathbf{v}⟩=0$
7	最小二乘	Least Squares	$\min \Vert A\mathbf{x}-\mathbf{b}{\Vert }^2$
8	矩阵微积分	Matrix Calculus	$\frac{\partial f(X)}{\partial X}$
9	SVD	Singular Value Decomposition	$A=U\Sigma V^T$
10	正定矩阵	Positive Definite	${\mathbf{x}}^{T}A\mathbf{x}>0,\forall \mathbf{x}\ne 0$
11	迹	Trace	$\text{tr}(A)={\sum }_i{A}_{ii}$
12	行列式	Determinant	$\det (A)$

---

一、用几何思维理解向量：一切的开始

1.1 向量是什么？

很多人第一次学向量觉得它就是一个有方向的箭头，或者一组数字。但对于机器学习来说，向量其实是数据的表示方式。

拿一个具体的例子来说：假设你在收集一只猫的特征——体重4公斤、体长50厘米、毛发长度3厘米、年龄2岁。如果你把这些数字排成一列：

$\mathbf{x}=\left(\begin{array}{c}4.0\\ 50.0\\ 3.0\\ 2.0\end{array}\right)$

恭喜，你创造了一个4维向量。在现实世界你没法”看见”4维空间，但在数学里这完全没问题。

几何直觉：把向量想象成一个从原点指向某点的箭头。在2维空间，这个箭头可以在纸上画出来；在100维空间，这个箭头依然存在，只是我们没法画出来了。

1.2 向量的基本运算

加法：两个向量相加，就是把它们的”箭头”首尾相接。或者说，把对应坐标加起来：

$\left(\begin{array}{c}1\\ 2\end{array}\right)+\left(\begin{array}{c}3\\ 1\end{array}\right)=\left(\begin{array}{c}4\\ 3\end{array}\right)$

数乘：用一个标量（普通数字）乘以向量，就是把箭头拉长、缩短，或者反向：

$3\times \left(\begin{array}{c}1\\ 2\end{array}\right)=\left(\begin{array}{c}3\\ 6\end{array}\right)$

点积（内积）：这是最有趣的操作。两个向量的点积就是把对应坐标相乘后相加：

$\left(\begin{array}{c}1\\ 2\\ 3\end{array}\right)\cdot \left(\begin{array}{c}4\\ 5\\ 6\end{array}\right)=1\times 4+2\times 5+3\times 6=32$

点积的几何意义是：衡量两个向量有多”相似”。如果两个向量方向相同，点积为正；如果垂直，点积为0；如果方向相反，点积为负。

import numpy as np
 
# 向量运算示例
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
 
# 加法
print("a + b =", a + b)  # [5 7 9]
 
# 数乘
print("3 * a =", 3 * a)  # [3 6 9]
 
# 点积
print("a dot b =", np.dot(a, b))  # 32
# 或者 a @ b
print("a @ b =", a @ b)  # 32
 
# 向量长度（范数）
print("||a|| =", np.linalg.norm(a))  # sqrt(1+4+9) ≈ 3.741

1.3 线性组合：创造一切的基础

线性组合听起来高大上，其实就是一个简单的公式：

$c_1{\mathbf{v}}_1+c_2{\mathbf{v}}_2+\cdots +c_k{\mathbf{v}}_k$

就是用一些系数（标量）乘以向量，然后加起来。

为什么这很重要？因为神经网络里的每一层计算，本质上就是在做线性组合。输入特征通过权重做线性组合，加上偏置，再通过激活函数——整个过程就是这么简单操作的反复叠加。

1.4 向量空间与基

向量空间就是所有可能的向量组成的”世界”。${\mathbb{R}}^n$ 是 n 维实数向量空间。

基是一组特殊的向量，它们有两个性质：

1. 线性无关（你不能用其中一些拼出另一些）
2. 能张成整个空间（空间里任何向量都能用它们表示）

标准基就是 ${\mathbf{e}}_1=(1,0,\dots ,0),{\mathbf{e}}_2=(0,1,\dots ,0)$ 等等。

关键洞察：同一个向量，在不同基下有不同的坐标表示。就像北京同一个地方，可以用经纬度表示，也可以用距离天安门的距离和方向来表示。

---

二、矩阵：线性变换的数学语言

2.1 矩阵的几何意义

矩阵本质上是一个线性变换的表示。给定一个矩阵 $A$，当你用它乘以一个向量 $\mathbf{x}$ 时，你实际上是在对 $\mathbf{x}$ 做某种变换——旋转、拉伸、压缩，或者这三者的组合。

举一个具体的 $2\times 2$ 矩阵例子：

$A=\left(\begin{array}{cc}2& 0\\ 0& 0.5\end{array}\right)$

当你用这个矩阵乘以任何2维向量时，$x$ 坐标会被乘以2（拉长2倍），$y$ 坐标会被乘以0.5（压缩到一半）。

import numpy as np
import matplotlib.pyplot as plt
 
# 定义变换矩阵
A = np.array([[2, 0], 
              [0, 0.5]])
 
# 原始向量
v = np.array([1, 1])
 
# 应用变换
v_transformed = A @ v
 
print(f"原始向量: {v}")
print(f"变换后: {v_transformed}")
# 原始向量: [1 1]
# 变换后: [2.  0.5]

2.2 矩阵乘法的直观理解

很多人背公式 $(AB{)}_{ij}={\sum }_k{A}_{ik}{B}_{kj}$，但其实有更好的理解方式。

矩阵乘法 $AB$ 的本质是”先做 B 变换，再做 A 变换”。

这就像函数的复合：$(A\circ B)(\mathbf{x})=A(B(\mathbf{x}))$。矩阵乘法从右往左读。

# 矩阵乘法示例
A = np.array([[1, 2], 
              [3, 4]])
B = np.array([[5, 6], 
              [7, 8]])
 
# 标准乘法
C = A @ B
print("A @ B =")
print(C)
 
# 验证：(AB)x = A(Bx)
x = np.array([1, 1])
print("(AB)x =", C @ x)
print("A(Bx) =", A @ (B @ x))
# 两者相等

2.3 为什么神经网络里到处都是矩阵乘法？

这是个好问题。神经网络的核心运算是：

$\mathbf{y}=f(W\mathbf{x}+\mathbf{b})$

其中 $W$ 是权重矩阵，$\mathbf{x}$ 是输入，$\mathbf{b}$ 是偏置，$f$ 是激活函数。

原因一：高效计算。用矩阵乘法一次性处理整个层的所有神经元，比用循环逐个计算快几十上百倍。GPU的并行计算能力在这里发挥得淋漓尽致。

原因二：表达能力。一个 $W$ 矩阵可以表示旋转、拉伸、投影等任意线性变换的组合。足够多的参数就能逼近任意复杂的函数。

原因三：梯度计算简洁。反向传播本质上就是链式法则，而矩阵形式让梯度计算变得整齐划一。

# 神经网络前向传播的矩阵形式
class SimpleLayer:
    def __init__(self, input_dim, output_dim):
        # Xavier初始化
        self.W = np.random.randn(output_dim, input_dim) * np.sqrt(2.0 / input_dim)
        self.b = np.zeros(output_dim)
    
    def forward(self, x):
        # x: (batch_size, input_dim)
        # W: (output_dim, input_dim)
        # b: (output_dim,)
        # 输出: (batch_size, output_dim)
        return self.activation(self.W @ x.T).T + self.b
    
    def activation(self, x):
        return np.maximum(0, x)  # ReLU
 
# 批量处理示例
batch_size = 32
input_dim = 784  # MNIST图像维度
output_dim = 128
 
layer = SimpleLayer(input_dim, output_dim)
x_batch = np.random.randn(batch_size, input_dim)
y_batch = layer.forward(x_batch)
print(f"输入形状: {x_batch.shape}")
print(f"输出形状: {y_batch.shape}")

2.4 矩阵的转置、逆、行列式

转置 $A^T$：把矩阵的行和列互换。几何上，这相当于做一个变换然后”翻转”回来。

逆矩阵 $A^{-1}$：满足 $A{A}^{-1}=I$。几何上，这相当于”撤销”之前的变换。不是所有矩阵都有逆（奇异矩阵就没有）。

行列式 $\det (A)$：衡量矩阵”缩放”空间的程度。如果 $\det (A)=0$，说明变换把空间”压扁”了（降维了），这样的矩阵不可逆。

A = np.array([[1, 2], 
              [3, 4]])
 
# 转置
print("A.T =")
print(A.T)
 
# 行列式
print("det(A) =", np.linalg.det(A))  # -2.0
 
# 逆矩阵
A_inv = np.linalg.inv(A)
print("A^(-1) =")
print(A_inv)
print("A @ A^(-1) =")
print(A @ A_inv)  # 应该接近单位矩阵

---

三、特征值与特征向量：变换的本质

3.1 它们到底在说什么？

特征值和特征向量可能是线性代数里最”魔法”的概念了。让我用一个生活化的例子来解释。

想象你有一块橡皮泥，你对它做一个变换（比如拉伸、旋转、剪切）。特征向量就是那些在变换后仍然指向原来方向的向量。特征值则是这个向量被拉伸或压缩了多少倍。

数学语言：$A\mathbf{v}=\lambda \mathbf{v}$

• $A$ 是变换矩阵
• $\mathbf{v}$ 是特征向量
• $\lambda$ 是特征值

几何直觉：

• 如果 $\lambda =2$，特征向量被拉长2倍
• 如果 $\lambda =0.5$，特征向量被压缩到一半
• 如果 $\lambda =-1$，特征向量反向（转180度）
• 如果 $\lambda =0$，这个方向的信息完全丢失了

3.2 为什么特征值这么重要？

原因一：理解高维变换。一个 $n\times n$ 的矩阵变换可能很复杂，但沿着特征向量的方向，这个变换就是简单的拉伸/压缩。知道了特征值，你就知道在哪些方向上变化大，哪些方向上变化小。

原因二：降维。在PCA里，我们保留特征值最大的几个方向，因为它们包含最多的”信息”。

原因三：稳定性分析。如果所有特征值的绝对值都小于1，迭代 $A^k$ 会收敛到零。这在动力系统和差分方程里非常重要。

import numpy as np
 
# 定义一个变换矩阵
A = np.array([[3, 1], 
              [0, 2]])
 
# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)
 
print("特征值:")
print(eigenvalues)  # [3. 2.]
 
print("\n特征向量（列）:")
print(eigenvectors)
 
# 验证：A @ v = lambda * v
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lam = eigenvalues[i]
    print(f"\n验证特征对 {i+1}:")
    print(f"A @ v = {A @ v}")
    print(f"lambda * v = {lam * v}")

3.3 特征值的性质

几个关键性质：

• 迹（Trace） $=\text{tr}(A)=\sum {\lambda }_i$：特征值的和等于矩阵对角线元素之和
• 行列式 $=\det (A)=\prod {\lambda }_i$：特征值的乘积等于矩阵行列式
• 秩 $=\text{rank}(A)=\text{非零特征值的个数}$

这些性质在分析矩阵时非常有用。比如你知道了一个矩阵的迹和行列式，你就知道了特征值的和与乘积。

3.4 对称矩阵的特征值

对称矩阵（$A=A^T$）的特征值有特别好的性质：都是实数。

非对称矩阵的特征值可能是复数，比如旋转矩阵：

$R=\left(\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)$

它的特征值是 $e^{i\theta }$ 和 $e^{-i\theta }$，复数形式对应着纯粹的旋转（没有拉伸）。

# 对称矩阵 vs 非对称矩阵
A_sym = np.array([[2, 1], 
                  [1, 2]])
                  
A_rot = np.array([[0, -1], 
                  [1, 0]])  # 旋转90度
 
print("对称矩阵特征值（实数）:")
print(np.linalg.eigvals(A_sym))  # [3. 1.]
 
print("\n旋转矩阵特征值（复数）:")
print(np.linalg.eigvals(A_rot))  # [0.+1.j 0.-1.j]

---

四、SVD：最通用的矩阵分解

4.1 SVD的几何直觉

奇异值分解（SVD）是线性代数里最”全能”的分解。对于任意 $m\times n$ 矩阵，都存在分解：

$A=U\Sigma V^T$

其中：

• $U$ 是 $m\times m$ 正交矩阵（左奇异向量）
• $\Sigma$ 是 $m\times n$ 对角矩阵，对角线是奇异值 ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$
• $V$ 是 $n\times n$ 正交矩阵（右奇异向量）

几何意义：任何线性变换 $A$ 都可以分解为三步：

1. 旋转/反射 $V^T$：先旋转到”容易处理”的方向
2. 缩放 $\Sigma$：在各个方向做不同程度的拉伸
3. 旋转/反射 $U$：再旋转回来

这就像你有一块形状奇怪的布料（$A$），SVD告诉你：先旋转、再拉伸、再旋转，就能得到这块布料。

4.2 为什么SVD比特征分解更通用？

特征分解 $A=Q\Lambda Q^{-1}$ 有个前提：$A$ 必须有 $n$ 个线性无关的特征向量。这要求 $A$ 是可对角化的。

但SVD不一样。任何矩阵都有SVD，不管它是方的还是长方的、奇异的还是满秩的。

SVD与特征分解的关系：

• $A{A}^T=U{\Sigma }^2{U}^T$（$A{A}^T$ 的特征值是 ${\sigma }_i^2$）
• $A^{T}A=V{\Sigma }^2{V}^T$（$A^{T}A$ 的特征值也是 ${\sigma }_i^2$）

import numpy as np
 
# 任意矩阵的SVD
A = np.array([[1, 2, 3], 
              [4, 5, 6]])
 
U, S, Vt = np.linalg.svd(A)
 
print("U (左奇异向量):")
print(U)
print(f"\n奇异值: {S}")
print("\nVt (右奇异向量的转置):")
print(Vt)
 
# 验证分解
Sigma = np.zeros_like(A, dtype=float)
Sigma[:min(A.shape), :min(A.shape)] = np.diag(S)
reconstructed = U @ Sigma @ Vt
print("\n重构矩阵 (应该接近A):")
print(reconstructed)

4.3 SVD与降维

截断SVD是SVD最强大的应用之一。保留最大的 $k$ 个奇异值：

$A_k=U_k{\Sigma }_k{V}_k^T$

这就是矩阵的最优 $k$ 秩逼近。

# 图像压缩示例
from skimage import data
import matplotlib.pyplot as plt
 
# 加载一张图片作为矩阵
image = data.camera()  # 512x512 灰度图
print(f"原始图像大小: {image.shape}")
print(f"原始存储: {512 * 512} 个像素值")
 
# SVD分解
U, S, Vt = np.linalg.svd(image.astype(float))
 
# 用前k个奇异值重建
def reconstruct_svd(U, S, Vt, k):
    Sigma = np.zeros_like(image, dtype=float)
    Sigma[:k, :k] = np.diag(S[:k])
    return U @ Sigma @ Vt
 
# 比较不同k的压缩效果
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
k_values = [1, 5, 10, 30, 100, 512]
 
for ax, k in zip(axes.flat, k_values):
    compressed = reconstruct_svd(U, S, Vt, k)
    compression_ratio = (k * (512 + 512 + 1)) / (512 * 512)
    ax.imshow(compressed, cmap='gray')
    ax.set_title(f'k={k}, 压缩比={compression_ratio:.1%}')
    ax.axis('off')
 
plt.tight_layout()
plt.show()
 
# 计算重构误差
for k in [10, 30, 50, 100]:
    compressed = reconstruct_svd(U, S, Vt, k)
    error = np.linalg.norm(image - compressed) / np.linalg.norm(image)
    print(f"k={k}: 重构误差 = {error:.4f}")

4.4 SVD与伪逆

当 $A$ 不是方阵时，普通的逆不存在，但伪逆存在：

$A^{+}=V{\Sigma }^{+}{U}^T$

其中 ${\Sigma }^{+}$ 是把非零奇异值取倒数后转置。

伪逆的性质：

• 如果 $A\mathbf{x}=\mathbf{b}$ 无解，$A^{+}\mathbf{b}$ 给出最小二乘解
• 如果有多解，$A^{+}\mathbf{b}$ 给出范数最小的解

# 伪逆计算
A = np.array([[1, 1], 
              [1, 2], 
              [2, 2]])
 
A_pinv = np.linalg.pinv(A)
print("伪逆:")
print(A_pinv)
 
# 验证性质
b = np.array([1, 2, 3])
x_ls = A_pinv @ b  # 最小二乘解
print(f"\n最小二乘解: {x_ls}")
print(f"残差: {np.linalg.norm(A @ x_ls - b)}")

---

五、最小二乘法：拟合数据的神器

5.1 几何意义

最小二乘法的目标是找到 $\mathbf{x}$ 使得 $\Vert A\mathbf{x}-\mathbf{b}{\Vert }^2$ 最小。

几何直觉：把 $\mathbf{b}$ 投影到 $A$ 的列空间，投影点就是 $A\hat{\mathbf{x}}$。我们要找的就是那个最优投影对应的 $\hat{\mathbf{x}}$。

如果 $\mathbf{b}$ 本身就在列空间里，完美！残差为0。 如果 $\mathbf{b}$ 不在列空间里，我们就找最”接近”的点。

5.2 正规方程

最小二乘解满足正规方程：

$A^{T}A\mathbf{x}=A^T\mathbf{b}$

解为 $\hat{\mathbf{x}}=(A^{T}A{)}^{-1}{A}^T\mathbf{b}$

import numpy as np
 
# 生成一些带噪声的数据
np.random.seed(42)
x = np.linspace(0, 10, 50)
y_true = 3 * x + 2  # 真实关系
noise = np.random.randn(50) * 2  # 噪声
y = y_true + noise
 
# 构建设计矩阵
X = np.column_stack([np.ones(50), x])  # 添加截距项
 
# 正规方程求解
beta = np.linalg.inv(X.T @ X) @ X.T @ y
print(f"最小二乘估计: 截距={beta[0]:.3f}, 斜率={beta[1]:.3f}")
print(f"真实值: 截距=2.000, 斜率=3.000")
 
# 使用NumPy内置函数
beta_np, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None)
print(f"NumPy lstsq: 截距={beta_np[0]:.3f}, 斜率={beta_np[1]:.3f}")

5.3 数值稳定性：为什么不用正规方程直接解？

直接算 $(A^{T}A{)}^{-1}$ 听起来简单，但有两个问题：

问题一：数值不稳定。$A^{T}A$ 的条件数是 $A$ 的平方。当 $A$ 病态时，$A^{T}A$ 更病态。

问题二：计算量大。求逆是 $O(n^3)$，而用 QR 分解或 SVD 求解更高效。

推荐方法：

1. QR分解：$A=QR$，解 $\hat{\mathbf{x}}=R^{-1}{Q}^T\mathbf{b}$
2. SVD：最稳定，尤其当 $A$ 接近奇异时

import numpy as np
 
# 病态矩阵示例
np.random.seed(123)
n = 10
A = np.random.randn(n, n)
# 让A接近奇异
A[n-1] = A[0] + 1e-10 * np.random.randn(n)
 
b = np.random.randn(n)
 
# 方法1：正规方程（不稳定）
try:
    x1 = np.linalg.inv(A.T @ A) @ A.T @ b
    print("正规方程解成功")
except:
    print("正规方程求解失败（病态）")
 
# 方法2：QR分解（推荐）
Q, R = np.linalg.qr(A)
x2 = np.linalg.solve(R, Q.T @ b)
print(f"QR分解误差: {np.linalg.norm(A @ x2 - b):.6f}")
 
# 方法3：SVD（最稳定）
U, S, Vt = np.linalg.svd(A, full_matrices=False)
x3 = Vt.T @ (U.T @ b / S)
print(f"SVD误差: {np.linalg.norm(A @ x3 - b):.6f}")

---

六、PCA：从数学到代码的完整演示

6.1 PCA在做什么？

PCA（主成分分析）本质上是找到一个坐标变换，让数据在新坐标下的不相关性最大化、方差最大化。

几何理解：想象你在看一群点云。PCA做的事情就是找到一条线，让所有点到这个线的垂距平方和最小——这条线就是第一主成分。然后找第二条线垂直于第一条，让剩余方差最大——这是第二主成分。

数学表述：

1. 数据中心化：$\tilde{X}=X-\bar{X}$
2. 计算协方差：$C=\frac{1}{n-1}{\tilde{X}}^T\tilde{X}$
3. SVD分解：$\tilde{X}=U\Sigma V^T$
4. 主成分方向就是 $V$ 的列（前 $k$ 列就是前 $k$ 个主成分）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
 
# 生成3D数据（不同方向方差不同）
np.random.seed(42)
n = 500
 
# 数据点分布在两个主方向附近
t = np.random.randn(n)
x = 3 * t + np.random.randn(n) * 0.5  # 主方向1，方差大
y = t + np.random.randn(n) * 0.3  # 主方向2
z = 0.5 * t + np.random.randn(n) * 0.2  # 主方向3，方差小
 
# 组成数据矩阵
X = np.column_stack([x, y, z])
X_centered = X - X.mean(axis=0)
 
# 使用sklearn的PCA
pca = PCA(n_components=3)
pca.fit(X_centered)
 
print("主成分方向（V的行）:")
print(pca.components_)
print(f"\n各主成分解释的方差比例: {pca.explained_variance_ratio_}")
print(f"累计解释方差: {np.cumsum(pca.explained_variance_ratio_)}")
 
# 手动计算验证
U, S, Vt = np.linalg.svd(X_centered)
print("\n手动SVD验证奇异值:", S)
print("sklearn解释方差:", pca.explained_variance_ * (n - 1))

6.2 PCA的降维应用

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
 
# 人脸数据降维示例
np.random.seed(42)
 
# 模拟一些人脸特征（实际中用人脸数据集）
n_faces = 100
n_pixels = 64 * 64  # 64x64灰度图
 
# 生成一些"假"人脸数据
faces = np.random.randn(n_faces, n_pixels)
# 添加一些结构（让数据有低维结构）
latent = np.random.randn(n_faces, 50)  # 50维隐变量
W = np.random.randn(50, n_pixels)
faces = faces + latent @ W
 
print(f"原始数据维度: {faces.shape}")
print(f"数据内在维度估计: ~50")
 
# PCA降维到不同维度
fig, axes = plt.subplots(2, 4, figsize=(12, 6))
 
for i, n_components in enumerate([5, 10, 20, 50, 100, 200, 500, 4096]):
    pca = PCA(n_components=min(n_components, n_faces - 1))
    faces_reduced = pca.fit_transform(faces)
    faces_reconstructed = pca.inverse_transform(faces_reduced)
    
    # 计算重构误差
    error = np.mean((faces - faces_reconstructed) ** 2)
    
    ax = axes.flat[i]
    # 显示重构后的人脸（取第一张）
    ax.imshow(faces_reconstructed[0].reshape(64, 64), cmap='gray')
    ax.set_title(f'k={n_components}\nMSE={error:.2f}')
    ax.axis('off')
 
plt.tight_layout()
plt.show()

6.3 PCA与SVD的关系

PCA和SVD本质上是同一件事：

$X=U\Sigma V^T$

• $V$ 的列（行 $V^T$）就是主成分方向
• $\Sigma$ 的对角元素与方差有关
• $U\Sigma$ 就是降维后的数据

为什么SVD更常用？

1. SVD对奇异矩阵也有效
2. 计算更稳定
3. 可以增量计算（适合大数据）

---

七、NumPy代码实战

7.1 矩阵分解大全

import numpy as np
 
A = np.array([[1, 2, 3], 
              [4, 5, 6], 
              [7, 8, 0]])
 
# 1. LU分解
P, L, U = np.linalg.lu(A)
print("LU分解:")
print("P (置换矩阵):\n", P)
print("L (下三角):\n", L)
print("U (上三角):\n", U)
print("P @ L @ U =\n", P @ L @ U)
 
# 2. QR分解
Q, R = np.linalg.qr(A)
print("\nQR分解:")
print("Q (正交):\n", Q)
print("R (上三角):\n", R)
print("Q @ R =\n", Q @ R)
 
# 3. 特征分解
eigenvalues, eigenvectors = np.linalg.eig(A)
print("\n特征分解:")
print("特征值:", eigenvalues)
print("特征向量 (列):\n", eigenvectors)
 
# 4. SVD
U, S, Vt = np.linalg.svd(A)
print("\nSVD:")
print("U:\n", U)
print("奇异值:", S)
print("Vt:\n", Vt)
 
# 5. Cholesky分解（需要正定矩阵）
# 先构造一个正定矩阵
B = A @ A.T  # 这保证是正定的
L_chol = np.linalg.cholesky(B)
print("\nCholesky分解:")
print("L:\n", L_chol)
print("L @ L.T =\n", L_chol @ L_chol.T)

7.2 解线性方程组

import numpy as np
 
A = np.array([[3, 1, 2], 
              [1, 2, 1], 
              [2, 1, 3]])
b = np.array([10, 8, 13])
 
# 方法1：直接求解
x1 = np.linalg.solve(A, b)
print("直接求解:", x1)
 
# 方法2：逆矩阵
x2 = np.linalg.inv(A) @ b
print("逆矩阵求解:", x2)
 
# 方法3：LU分解
P, L, U = np.linalg.lu(A)
y = np.linalg.solve(L, P @ b)
x3 = np.linalg.solve(U, y)
print("LU分解求解:", x3)
 
# 验证
print("残差:", np.linalg.norm(A @ x1 - b))

7.3 稀疏矩阵操作

import numpy as np
from scipy import sparse
 
# 创建稀疏矩阵（CSR格式）
data = np.array([1, 2, 3, 4, 5])
row = np.array([0, 0, 1, 2, 2])
col = np.array([0, 2, 2, 0, 1])
 
A_sparse = sparse.csr_matrix((data, (row, col)), shape=(3, 3))
print("稀疏矩阵:\n", A_sparse.toarray())
 
# 稀疏矩阵乘法
x = np.array([1, 2, 3])
y = A_sparse @ x
print("\n稀疏矩阵乘法结果:", y)
 
# 稠密 vs 稀疏 存储对比
n = 10000
# 稠密矩阵
A_dense = np.random.randn(n, n)
# 稀疏矩阵（只有1%的非零元素）
A_sparse = sparse.random(n, n, density=0.01, format='csr')
 
print(f"\n稠密矩阵存储: {A_dense.nbytes / 1e6:.2f} MB")
print(f"稀疏矩阵存储: {A_sparse.data.nbytes / 1e6:.2f} MB")

---

八、线性代数在神经网络中的应用

8.1 前向传播的矩阵形式

全连接层的前向传播：

$Z^{(l)}=W^{(l)}{A}^{(l-1)}+{\mathbf{b}}^{(l)}$

其中：

• $W^{(l)}$：权重矩阵，shape $(n^{(l)},n^{(l-1)})$
• $A^{(l-1)}$：上一层的激活值，shape $(n^{(l-1)},batch\_size)$
• $Z^{(l)}$：线性组合结果
• ${\mathbf{b}}^{(l)}$：偏置向量

import numpy as np
 
class DenseLayer:
    def __init__(self, input_dim, output_dim):
        # Xavier初始化
        self.W = np.random.randn(output_dim, input_dim) * np.sqrt(2.0 / input_dim)
        self.b = np.zeros(output_dim)
        
    def forward(self, x):
        """x: (input_dim, batch_size)"""
        self.x = x
        self.z = self.W @ x + self.b.reshape(-1, 1)
        return self.activation(self.z)
    
    def activation(self, z):
        return np.maximum(0, z)  # ReLU
    
    def activation_grad(self, z):
        return (z > 0).astype(float)  # ReLU导数
 
# 多层网络
np.random.seed(42)
layers = [
    DenseLayer(784, 256),   # 输入层
    DenseLayer(256, 128),   # 隐藏层1
    DenseLayer(128, 64),    # 隐藏层2
    DenseLayer(64, 10)      # 输出层
]
 
# 模拟输入（batch_size=32）
x = np.random.randn(784, 32)
 
# 前向传播
a = x
for layer in layers:
    a = layer.forward(a)
 
print(f"输入形状: {x.shape}")
print(f"最终输出形状: {a.shape}")

8.2 反向传播的矩阵形式

反向传播的核心是链式法则。梯度计算：

$\frac{\partial \mathcal{L}}{\partial W^{(l)}}=\frac{\partial \mathcal{L}}{\partial Z^{(l)}}\cdot (A^{(l-1)}{)}^T$

$\frac{\partial \mathcal{L}}{\partial A^{(l-1)}}=(W^{(l)}{)}^T\cdot \frac{\partial \mathcal{L}}{\partial Z^{(l)}}$

class DenseLayer:
    def __init__(self, input_dim, output_dim):
        self.W = np.random.randn(output_dim, input_dim) * np.sqrt(2.0 / input_dim)
        self.b = np.zeros(output_dim)
    
    def forward(self, x):
        self.x = x
        self.z = self.W @ x + self.b.reshape(-1, 1)
        self.a = self.activation(self.z)
        return self.a
    
    def backward(self, da):
        """da: 来自上一层的梯度 (output_dim, batch_size)"""
        # 激活函数梯度
        dz = da * self.activation_grad(self.z)
        
        # 权重梯度: dL/dW = dz @ x.T / batch_size
        batch_size = self.x.shape[1]
        self.dW = dz @ self.x.T / batch_size
        self.db = np.mean(dz, axis=1)
        
        # 传递到前一层的梯度
        dA = self.W.T @ dz
        return dA
 
# 梯度检验
np.random.seed(42)
layer = DenseLayer(10, 5)
x = np.random.randn(10, 3)
 
# 前向传播
a = layer.forward(x)
 
# 假设损失是MSE
y_true = np.random.randn(5, 3)
loss = np.mean((a - y_true) ** 2)
 
# 反向传播
da = 2 * (a - y_true) / 3  # MSE的梯度
dA = layer.backward(da)
 
# 数值梯度检验
epsilon = 1e-5
numerical_grad = np.zeros_like(layer.W)
for i in range(layer.W.shape[0]):
    for j in range(layer.W.shape[1]):
        W_eps = layer.W.copy()
        W_eps[i, j] += epsilon
        layer_eps = DenseLayer(10, 5)
        layer_eps.W = W_eps
        a_eps = layer_eps.forward(x)
        loss_eps = np.mean((a_eps - y_true) ** 2)
        numerical_grad[i, j] = (loss_eps - loss) / epsilon
 
print("解析梯度 vs 数值梯度:")
print(f"最大相对误差: {np.max(np.abs(layer.dW - numerical_grad) / (np.abs(numerical_grad) + 1e-8)):.2e}")

8.3 注意力机制中的矩阵运算

Transformer的核心是自注意力：

$\text{Attention}(Q,K,V)=\text{softmax}\left(\frac{Q{K}^T}{\sqrt{d_k}}\right)V$

线性代数视角：

• $Q=X{W}_Q,K=X{W}_K,V=X{W}_V$：三个线性变换
• $Q{K}^T$：计算所有位置之间的”相关性”
• softmax归一化：把相关性变成概率分布
• 与 $V$ 相乘：用概率分布对值加权求和

import numpy as np
 
def self_attention(Q, K, V, d_k):
    """
    自注意力机制
    
    Q, K, V: (seq_len, d_model)
    d_k: 缩放因子
    """
    # 计算注意力分数
    scores = Q @ K.T / np.sqrt(d_k)
    
    # softmax归一化
    exp_scores = np.exp(scores - np.max(scores, axis=-1, keepdims=True))
    attention_weights = exp_scores / np.sum(exp_scores, axis=-1, keepdims=True)
    
    # 加权求和
    output = attention_weights @ V
    
    return output, attention_weights
 
# 示例
np.random.seed(42)
seq_len = 10
d_model = 64
d_k = d_v = 64
 
# 模拟输入
X = np.random.randn(seq_len, d_model)
 
# 生成Q, K, V
W_Q = np.random.randn(d_k, d_model) * np.sqrt(2.0 / d_model)
W_K = np.random.randn(d_k, d_model) * np.sqrt(2.0 / d_model)
W_V = np.random.randn(d_v, d_model) * np.sqrt(2.0 / d_model)
 
Q = X @ W_Q.T
K = X @ W_K.T
V = X @ W_V.T
 
# 计算注意力
output, attention = self_attention(Q, K, V, d_k)
 
print(f"输入形状: {X.shape}")
print(f"注意力权重形状: {attention.shape}")
print(f"输出形状: {output.shape}")
print(f"\n注意力权重（第一行，代表第一个位置对所有位置的注意力）:")
print(attention[0])

---

九、常见数值问题与解决方案

9.1 病态矩阵

问题：当矩阵接近奇异（行列式接近0）时，解对输入扰动非常敏感。

诊断：看条件数 $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$

• $\kappa \approx 1$：良态
• $\kappa >10^6$：严重病态

import numpy as np
 
# 构造病态矩阵
n = 10
A = np.random.randn(n, n)
# 让它病态：添加一个接近0的特征值
eigenvalues = np.random.randn(n)
eigenvalues[-1] = 1e-10  # 一个特征值接近0
Q = np.linalg.qr(np.random.randn(n, n))[0]
A = Q @ np.diag(eigenvalues) @ Q.T
 
# 计算条件数
cond = np.linalg.cond(A)
print(f"条件数: {cond:.2e}")
 
# 病态矩阵的影响
b = np.random.randn(n)
b_noise = b + 1e-8 * np.random.randn(n)  # 微小扰动
 
x = np.linalg.solve(A, b)
x_noise = np.linalg.solve(A, b_noise)
 
print(f"解的相对变化: {np.linalg.norm(x - x_noise) / np.linalg.norm(x):.2e}")
print(f"输入的相对变化: {np.linalg.norm(b - b_noise) / np.linalg.norm(b):.2e}")

解决方案：

1. 正则化：$A+\lambda I$
2. 截断SVD：扔掉小奇异值
3. 迭代优化：用共轭梯度法

# 解决方案1：正则化（岭回归）
lamda = 1e-6
A_reg = A + lamda * np.eye(n)
x_reg = np.linalg.solve(A_reg, b)
print(f"\n正则化后解的相对变化: {np.linalg.norm(x_reg - x_noise) / np.linalg.norm(x_reg):.2e}")
 
# 解决方案2：截断SVD
U, S, Vt = np.linalg.svd(A)
k = 9  # 保留前k个奇异值
S_truncated = np.copy(S)
S_truncated[k:] = 0  # 把小奇异值设为0
A_truncated = U @ np.diag(S_truncated) @ Vt
x_svd = Vt.T @ (U.T @ b / S_truncated)
print(f"SVD后解的相对变化: {np.linalg.norm(x_svd - x_noise) / np.linalg.norm(x_svd):.2e}")

9.2 秩亏矩阵

问题：矩阵的秩小于维度，存在无穷多解或无解。

诊断：SVD奇异值里有明显的”零值”。

# 秩亏矩阵示例
A = np.array([[1, 2, 3], 
              [2, 4, 6],  # 第三行是第一行的两倍
              [1, 1, 1]])
 
# 计算秩
rank = np.linalg.matrix_rank(A)
print(f"矩阵秩: {rank} (满秩应该是3)")
 
# SVD诊断
U, S, Vt = np.linalg.svd(A)
print(f"奇异值: {S}")
print(f"非零奇异值个数: {np.sum(S > 1e-10)}")

9.3 浮点精度问题

import numpy as np
 
# 大数吃小数
a = 1e9
b = 1e-9
c = 1e9 + b
print(f"1e9 + 1e-9 = {c}")  # 丢失了b
 
# 解决方案：使用Kahan求和
def kahan_sum(x):
    """高精度求和"""
    total = 0.0
    compensation = 0.0
    for xi in x:
        y = xi - compensation
        t = total + y
        compensation = (t - total) - y
        total = t
    return total
 
# 大矩阵运算的精度问题
n = 1000
A = np.random.randn(n, n)
A = A @ A.T / n  # 让特征值范围合理
 
# 检查正交性
Q, R = np.linalg.qr(A)
print(f"Q.T @ Q 是否接近单位阵: {np.max(np.abs(Q.T @ Q - np.eye(n))):.2e}")

---

十、NumPy/SciPy工具函数速查

10.1 常用NumPy函数

import numpy as np
import numpy.linalg as la
import numpy.random as npr
 
# 向量操作
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
 
np.dot(a, b)          # 点积
np.cross(a, b)        # 叉积（3D向量）
np.linalg.norm(a)     # 范数（默认L2）
np.linalg.norm(a, ord=1)  # L1范数
np.linalg.norm(a, ord=np.inf)  # 无穷范数
 
# 矩阵创建
np.eye(3)             # 单位矩阵
np.zeros((3, 4))      # 零矩阵
np.ones((3, 4))       # 全1矩阵
np.random.randn(3, 4)  # 随机矩阵（正态分布）
np.diag([1, 2, 3])    # 对角矩阵
 
# 矩阵属性
A.shape               # 形状
A.T                   # 转置
np.trace(A)           # 迹
np.linalg.det(A)      # 行列式
np.linalg.matrix_rank(A)  # 秩
 
# 矩阵分解
np.linalg.svd(A)           # SVD
np.linalg.eig(A)           # 特征分解
np.linalg.eigvals(A)       # 只求特征值
np.linalg.qr(A)            # QR分解
np.linalg.lu(A)            # LU分解
np.linalg.cholesky(A)      # Cholesky分解（正定）
np.linalg.solve(A, b)      # 解线性方程组 Ax = b
np.linalg.lstsq(A, b)      # 最小二乘解
 
# 范数与条件数
np.linalg.norm(A)           # 默认为Frobenius范数
np.linalg.norm(A, 'fro')   # Frobenius范数
np.linalg.norm(A, 2)        # 谱范数（最大奇异值）
np.linalg.cond(A)          # 条件数

10.2 SciPy额外功能

from scipy import linalg, sparse, stats
 
# 更多矩阵分解
U, s, Vh = linalg.svd(A)          # 更多选项的SVD
lu = linalg.lu_factor(A)          # LU分解（用于多次求解）
x = linalg.lu_solve(lu, b)        # 用LU分解求解
 
# 稀疏矩阵
B = sparse.csr_matrix(A)          # 从稠密转稀疏
B.todense()                       # 转回稠密
sparse.random(1000, 1000, density=0.01)  # 创建稀疏随机矩阵
 
# 迭代求解器
x, info = linalg.cg(A, b)         # 共轭梯度法
x, info = linalg.bicgstab(A, b)  # 双共轭梯度稳定版
 
# 特殊矩阵
linalg.hilbert(5)                 # Hilbert矩阵（病态）
linalg.toeplitz([1,2,3], [1,4,5])  # Toeplitz矩阵

10.3 实用代码片段

import numpy as np
 
# 批量计算余弦相似度
def cosine_similarity(X):
    """X: (n_samples, n_features)"""
    X_norm = X / np.linalg.norm(X, axis=1, keepdims=True)
    return X_norm @ X_norm.T
 
# 批量矩阵的迹
def batch_trace(X):
    """X: (batch, n, n)"""
    return np.trace(X, axis1=1, axis2=2)
 
# 矩阵的平方根（正定矩阵）
def matrix_sqrt(A):
    """A: 对称正定矩阵"""
    eigenvalues, eigenvectors = np.linalg.eigh(A)
    return eigenvectors @ np.diag(np.sqrt(eigenvalues)) @ eigenvectors.T
 
# 计算协方差矩阵
def covariance(X):
    """X: (n_samples, n_features)，假设已中心化"""
    return X.T @ X / (len(X) - 1)
 
# 批量矩阵乘法
def batch_mm(A, B):
    """A: (batch, m, n), B: (batch, n, k) -> (batch, m, k)"""
    return np.einsum('bij,bjk->bik', A, B)

---

十一、其他重要章节（来自原文）

11.1 矩阵微积分基础

标量对向量求导： $f:{\mathbb{R}}^n\to \mathbb{R},\frac{\partial f}{\partial \mathbf{x}}=\left(\begin{array}{c}\frac{\partial f}{\partial x_1}\\ ⋮\\ \frac{\partial f}{\partial x_n}\end{array}\right)$

常用公式速记：

函数	导数
${\mathbf{a}}^T\mathbf{x}$	$\mathbf{a}$
${\mathbf{x}}^{T}A\mathbf{x}$（$A$ 对称）	$2A\mathbf{x}$
${\mathbf{x}}^T\mathbf{x}$	$2\mathbf{x}$
$\Vert A\mathbf{x}-\mathbf{b}{\Vert }^2$	$2A^T(A\mathbf{x}-\mathbf{b})$

11.2 约束优化与KKT条件

约束优化问题： $\min f(\mathbf{x})\text{s.t.}g(\mathbf{x})=0,h(\mathbf{x})\le 0$

KKT条件给出最优解的必要条件，涉及拉格朗日乘子。神经网络的对抗训练、GAN的优化都涉及这类问题。

11.3 随机矩阵理论

当矩阵维度很大时，特征值的分布呈现出规律性：

• Wigner半圆定律：随机对称矩阵的特征值分布趋向半圆
• Marchenko-Pastur定律：样本协方差矩阵的特征值分布

这在金融风险管理和高维统计分析中有重要应用。

11.4 流形上的优化

很多机器学习问题涉及约束：

• 正交矩阵优化（Stiefel流形）
• 正定矩阵优化（对称正定流形）
• 概率分布优化（概率单纯形）

需要使用黎曼梯度下降等专门的优化方法。

---

参考文献

1. Strang, G. (2009). Introduction to Linear Algebra (4th ed.). Wellesley-Cambridge Press.
2. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
3. Horn, R. A., & Johnson, C. R. (2012). Matrix Analysis (2nd ed.). Cambridge University Press.
4. Petersen, K. B., & Pedersen, M. S. (2012). The Matrix Cookbook. Technical University of Denmark.
5. Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning. MIT Press.
6. Demmel, J. W. (1997). Applied Numerical Linear Algebra. SIAM.
7. Trefethen, L. N., & Bau III, D. (1997). Numerical Linear Algebra. SIAM.
8. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms (2nd ed.). SIAM.
9. Boyd, S., & Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press.
10. Absil, P.-A., Mahony, R., & Sepulchre, R. (2009). Optimization Algorithms on Matrix Manifolds. Princeton University Press.

---

相关文档

• 凸优化深度指南 - 矩阵优化问题的求解方法
• 概率论深度指南 - 多元正态分布与协方差矩阵
• 统计学深度指南 - PCA与因子分析的统计基础
• 信息论深度指南 - 熵与矩阵分析的交叉
• 图神经网络 - 图谱理论与拉普拉斯矩阵