NumPy 线性代数

线性代数是数据科学和机器学习的数学基础。NumPy 的 numpy.linalg 模块提供了完整的线性代数运算功能，底层调用高效的 BLAS 和 LAPACK 库实现。本章将系统介绍 NumPy 中的线性代数操作。

矩阵运算基础

在 NumPy 中，矩阵就是二维数组。理解矩阵运算的基本概念是使用线性代数功能的前提。

矩阵乘法

矩阵乘法是最基本的线性代数运算。NumPy 提供了多种方式进行矩阵乘法：

import numpy as np

A = np.array([[1, 2],
              [3, 4]])
B = np.array([[5, 6],
              [7, 8]])

# 方法1: 使用 @ 运算符（推荐）
C1 = A @ B
print(f"@ 运算符:\n{C1}")

# 方法2: 使用 np.matmul
C2 = np.matmul(A, B)
print(f"\nnp.matmul:\n{C2}")

# 方法3: 使用 np.dot（较老的写法）
C3 = np.dot(A, B)
print(f"\nnp.dot:\n{C3}")

# 注意：* 是元素级乘法，不是矩阵乘法！
element_wise = A * B
print(f"\n元素级乘法 A * B:\n{element_wise}")

输出结果：

@ 运算符:
[[19 22]
 [43 50]]

元素级乘法 A * B:
[[ 5 12]
 [21 32]]

关键区别：@ 是矩阵乘法，* 是元素级乘法。这是初学者最容易混淆的地方。

点积与内积

import numpy as np

# 向量点积
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])

# 使用 np.dot
dot_product = np.dot(v1, v2)
print(f"点积: {dot_product}")  # 1*4 + 2*5 + 3*6 = 32

# 使用 @ 运算符
dot_product_v2 = v1 @ v2
print(f"@ 运算符: {dot_product_v2}")

# np.vdot 会先展平数组再计算点积
m1 = np.array([[1, 2], [3, 4]])
m2 = np.array([[5, 6], [7, 8]])
print(f"vdot: {np.vdot(m1, m2)}")  # 1*5 + 2*6 + 3*7 + 4*8

外积

外积将两个向量扩展为矩阵：

import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5])

# 外积
outer = np.outer(a, b)
print(f"外积 (a × b):\n{outer}")
# [[ 4,  5],
#  [ 8, 10],
#  [12, 15]]

# 等价于使用广播
outer_broadcast = a[:, np.newaxis] * b
print(f"\n广播实现:\n{outer_broadcast}")

张量积

对于高维数组的运算，NumPy 提供了 np.tensordot：

import numpy as np

# 两个三维数组
A = np.random.rand(2, 3, 4)
B = np.random.rand(4, 5, 6)

# 在 A 的最后一个轴和 B 的第一个轴上求和
C = np.tensordot(A, B, axes=([2], [0]))
print(f"tensordot 结果形状: {C.shape}")  # (2, 3, 5, 6)

# 矩阵乘法等价于 axes=1 的情况
M1 = np.random.rand(3, 4)
M2 = np.random.rand(4, 5)
C = np.tensordot(M1, M2, axes=1)
print(f"矩阵乘法结果形状: {C.shape}")  # (3, 5)

Kronecker 积

import numpy as np

A = np.array([[1, 2],
              [3, 4]])
B = np.array([[0, 5],
              [6, 7]])

# Kronecker 积
kron = np.kron(A, B)
print(f"Kronecker 积:\n{kron}")
# [[ 0,  5,  0, 10],
#  [ 6,  7, 12, 14],
#  [ 0, 15,  0, 20],
#  [18, 21, 24, 28]]

矩阵分解

矩阵分解是线性代数的核心内容，广泛应用于数据压缩、降维、推荐系统等领域。

奇异值分解（SVD）

奇异值分解是最重要和最通用的矩阵分解方法。任何实数矩阵都可以进行 SVD 分解。

对于矩阵 $A$，SVD 将其分解为：

$A = U \Sigma V^T$

其中：

• $U$ 是 $m \times m$ 的正交矩阵（左奇异向量）
• $\Sigma$ 是 $m \times n$ 的对角矩阵（奇异值）
• $V^T$ 是 $n \times n$ 的正交矩阵（右奇异向量的转置）

import numpy as np

# 创建一个矩阵
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9],
              [10, 11, 12]])

# SVD 分解
U, S, Vt = np.linalg.svd(A)

print(f"U 的形状: {U.shape}")   # (4, 4)
print(f"S 的形状: {S.shape}")   # (3,) - 只存储对角线上的非零元素
print(f"Vt 的形状: {Vt.shape}") # (3, 3)

print(f"\n奇异值: {S}")

# 验证分解是否正确
# 需要将 S 转换为完整的对角矩阵
Sigma = np.zeros((A.shape[0], A.shape[1]))
Sigma[:len(S), :len(S)] = np.diag(S)

# 重构矩阵
A_reconstructed = U @ Sigma @ Vt
print(f"\n重构矩阵:\n{A_reconstructed}")
print(f"\n重构误差: {np.allclose(A, A_reconstructed)}")

SVD 应用：低秩近似

SVD 的重要应用之一是矩阵的低秩近似，可以用于图像压缩：

import numpy as np

# 创建一个有结构的矩阵
A = np.random.rand(100, 100)
A = A @ A.T  # 使其具有较低的秩

# 完整 SVD
U, S, Vt = np.linalg.svd(A)

# 使用前 k 个奇异值进行近似
def low_rank_approximation(U, S, Vt, k):
    """使用前 k 个奇异值进行低秩近似"""
    Sigma_k = np.diag(S[:k])
    return U[:, :k] @ Sigma_k @ Vt[:k, :]

# 比较不同 k 值的近似效果
for k in [5, 10, 20, 50]:
    A_approx = low_rank_approximation(U, S, Vt, k)
    error = np.linalg.norm(A - A_approx) / np.linalg.norm(A)
    print(f"k={k:2d}: 相对误差 = {error:.4f}")

# 计算原矩阵的有效秩（使用奇异值衰减）
print(f"\n奇异值衰减:")
print(f"前 5 个奇异值占比: {sum(S[:5]) / sum(S):.4f}")
print(f"前 10 个奇异值占比: {sum(S[:10]) / sum(S):.4f}")

特征值分解

特征值分解适用于方阵。对于矩阵 $A$，如果存在标量 $\lambda$ 和非零向量 $v$ 满足：

$A v = \lambda v$

则 $\lambda$ 是特征值，$v$ 是对应的特征向量。

import numpy as np

# 创建一个对称矩阵（保证特征值为实数）
A = np.array([[4, -2],
              [-2, 4]])

# 特征值分解
eigenvalues, eigenvectors = np.linalg.eig(A)

print(f"特征值: {eigenvalues}")
print(f"特征向量:\n{eigenvectors}")

# 验证: A @ v = λ @ v
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lambda_v = eigenvalues[i]
    print(f"\n验证第 {i+1} 个特征值:")
    print(f"A @ v = {A @ v}")
    print(f"λ @ v = {lambda_v * v}")
    print(f"相等: {np.allclose(A @ v, lambda_v * v)}")

对称矩阵的特征值分解

对于对称矩阵，使用 eigh 更高效，且保证返回实数：

import numpy as np

# 创建对称矩阵
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]])

# 对于对称矩阵，使用 eigh（更高效且保证实数特征值）
eigenvalues, eigenvectors = np.linalg.eigh(A)

print(f"特征值（升序）: {eigenvalues}")
print(f"特征向量（按特征值排序）:\n{eigenvectors}")

# 验证特征向量的正交性
print(f"\n特征向量正交性验证:")
print(f"V @ V.T = \n{eigenvectors @ eigenvectors.T}")
# 应该接近单位矩阵

QR 分解

QR 分解将矩阵分解为正交矩阵 $Q$ 和上三角矩阵 $R$ 的乘积：

$A = QR$

import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 10]])

# QR 分解
Q, R = np.linalg.qr(A)

print(f"Q (正交矩阵):\n{Q}")
print(f"\nR (上三角矩阵):\n{R}")

# 验证
print(f"\n验证 Q 是正交矩阵:")
print(f"Q @ Q.T:\n{Q @ Q.T}")  # 应该是单位矩阵

print(f"\n验证 A = Q @ R:")
print(f"重构 A:\n{Q @ R}")

Cholesky 分解

对于正定对称矩阵，Cholesky 分解将其分解为下三角矩阵和其转置的乘积：

$A = LL^T$

import numpy as np

# 创建一个正定对称矩阵
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]])

# Cholesky 分解
L = np.linalg.cholesky(A)

print(f"L (下三角矩阵):\n{L}")
print(f"\n验证 A = L @ L.T:")
print(f"重构 A:\n{L @ L.T}")

求解线性方程组

求解线性方程组 $Ax = b$ 是线性代数的基本问题。

直接求解

import numpy as np

# 方程组:
# 3x + y = 9
# x + 2y = 8
A = np.array([[3, 1],
              [1, 2]])
b = np.array([9, 8])

# 求解
x = np.linalg.solve(A, b)
print(f"解: x = {x[0]:.2f}, y = {x[1]:.2f}")

# 验证
print(f"验证 A @ x = {A @ x}")
print(f"原向量 b = {b}")

逆矩阵求解

也可以使用逆矩阵求解，但效率较低：

import numpy as np

A = np.array([[3, 1],
              [1, 2]])
b = np.array([9, 8])

# 使用逆矩阵求解（不推荐，效率低）
A_inv = np.linalg.inv(A)
x = A_inv @ b
print(f"解: {x}")

# 对于数值稳定性，优先使用 solve
x_solve = np.linalg.solve(A, b)
print(f"solve 解: {x_solve}")

最小二乘解

当方程组无精确解（超定方程组）时，使用最小二乘法：

import numpy as np

# 超定方程组（方程数 > 未知数）
# 拟合直线 y = mx + b
x_data = np.array([0, 1, 2, 3, 4])
y_data = np.array([1, 2.2, 2.8, 4.1, 5.0])

# 构造矩阵形式: A @ [m, b] = y
A = np.column_stack([x_data, np.ones_like(x_data)])
b = y_data

# 最小二乘解
result = np.linalg.lstsq(A, b, rcond=None)
m, b_intercept = result[0]
residual = result[1]

print(f"拟合直线: y = {m:.3f}x + {b_intercept:.3f}")
print(f"残差平方和: {residual[0]:.6f}")

# 使用拟合参数计算预测值
y_pred = m * x_data + b_intercept
print(f"预测值: {y_pred}")

齐次线性方程组

求解 $Ax = 0$ 的非零解（零空间）：

import numpy as np

# 创建一个秩不满的矩阵
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 计算矩阵的秩
rank = np.linalg.matrix_rank(A)
print(f"矩阵的秩: {rank}")  # 2 < 3，存在非零解

# 使用 SVD 找零空间
U, S, Vt = np.linalg.svd(A)

# 零空间是 Vt 中对应奇异值为 0 的行向量
# 对于这个例子，最小奇异值对应的向量就是零空间的一个基
null_space = Vt[-1]
print(f"零空间向量: {null_space}")
print(f"验证 A @ null_space ≈ 0: {A @ null_space}")

矩阵属性

矩阵的秩

矩阵的秩是其线性独立的行（或列）的最大数目：

import numpy as np

# 满秩矩阵
A_full = np.array([[1, 2],
                   [3, 4]])
print(f"满秩矩阵的秩: {np.linalg.matrix_rank(A_full)}")  # 2

# 秩亏矩阵
A_rank_deficient = np.array([[1, 2],
                              [2, 4]])  # 第二行是第一行的两倍
print(f"秩亏矩阵的秩: {np.linalg.matrix_rank(A_rank_deficient)}")  # 1

# 实际应用：检查矩阵是否可逆
def is_invertible(A):
    return A.shape[0] == A.shape[1] and np.linalg.matrix_rank(A) == A.shape[0]

print(f"\n满秩矩阵可逆: {is_invertible(A_full)}")
print(f"秩亏矩阵可逆: {is_invertible(A_rank_deficient)}")

行列式

行列式是方阵的一个标量值，对于判断矩阵是否可逆很有用：

import numpy as np

A = np.array([[1, 2],
              [3, 4]])

det = np.linalg.det(A)
print(f"行列式: {det}")

# 行列式为 0 的矩阵不可逆
B = np.array([[1, 2],
              [2, 4]])  # 奇异矩阵
det_B = np.linalg.det(B)
print(f"\n奇异矩阵的行列式: {det_B} (接近 0)")

# 对于高维数组，可以批量计算行列式
C = np.random.rand(3, 4, 4)  # 3 个 4x4 矩阵
dets = np.linalg.det(C)
print(f"\n批量行列式: {dets}")

条件数

条件数衡量矩阵对输入扰动的敏感程度，条件数越大，数值计算越不稳定：

import numpy as np

# 良好条件的矩阵
A_good = np.array([[1, 0],
                   [0, 1]])
cond_good = np.linalg.cond(A_good)
print(f"单位矩阵的条件数: {cond_good}")  # 接近 1

# 病态矩阵
A_ill = np.array([[1, 1],
                  [1, 1.0001]])
cond_ill = np.linalg.cond(A_ill)
print(f"病态矩阵的条件数: {cond_ill}")  # 很大

# 实际影响：求解线性方程组的精度
b = np.array([1, 1])

x_good = np.linalg.solve(A_good, b)
x_ill = np.linalg.solve(A_ill, b)

# 添加小扰动
b_perturbed = b + np.array([0.0001, 0])

x_good_p = np.linalg.solve(A_good, b_perturbed)
x_ill_p = np.linalg.solve(A_ill, b_perturbed)

print(f"\n良好矩阵解的变化: {np.linalg.norm(x_good_p - x_good)}")
print(f"病态矩阵解的变化: {np.linalg.norm(x_ill_p - x_ill)}")

范数

范数用于衡量向量或矩阵的"大小"：

import numpy as np

v = np.array([3, 4])

# 向量范数
print(f"L2 范数（默认）: {np.linalg.norm(v)}")  # sqrt(3^2 + 4^2) = 5
print(f"L1 范数: {np.linalg.norm(v, ord=1)}")   # |3| + |4| = 7
print(f"无穷范数: {np.linalg.norm(v, ord=np.inf)}")  # max(|3|, |4|) = 4

# 矩阵范数
A = np.array([[1, 2],
              [3, 4]])

print(f"\nFrobenius 范数: {np.linalg.norm(A)}")      # sqrt(sum(a_ij^2))
print(f"谱范数（2-范数）: {np.linalg.norm(A, ord=2)}")  # 最大奇异值
print(f"1-范数: {np.linalg.norm(A, ord=1)}")  # 列和的最大值
print(f"无穷范数: {np.linalg.norm(A, ord=np.inf)}")  # 行和的最大值

矩阵的迹

矩阵的迹是对角线元素之和：

import numpy as np

A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 使用 np.trace
trace = np.trace(A)
print(f"矩阵的迹: {trace}")  # 1 + 5 + 9 = 15

# 使用 np.linalg.trace（Array API 兼容）
trace_v2 = np.linalg.trace(A)
print(f"np.linalg.trace: {trace_v2}")

# 迹的性质验证
B = np.random.rand(3, 3)
print(f"\n验证 tr(AB) = tr(BA):")
print(f"tr(AB) = {np.trace(A @ B)}")
print(f"tr(BA) = {np.trace(B @ A)}")

矩阵的逆与伪逆

矩阵求逆

import numpy as np

A = np.array([[1, 2],
              [3, 4]])

# 计算逆矩阵
A_inv = np.linalg.inv(A)
print(f"逆矩阵:\n{A_inv}")

# 验证
print(f"\n验证 A @ A_inv = I:")
print(f"A @ A_inv:\n{A @ A_inv}")

# 使用逆矩阵求解线性方程组
b = np.array([5, 6])
x = A_inv @ b
print(f"\n解: {x}")
print(f"验证 A @ x = b: {A @ x}")

Moore-Penrose 伪逆

对于非方阵或奇异矩阵，使用伪逆：

import numpy as np

# 非方阵
A = np.array([[1, 2, 3],
              [4, 5, 6]])  # 2x3 矩阵

# 伪逆
A_pinv = np.linalg.pinv(A)
print(f"伪逆形状: {A_pinv.shape}")  # 3x2
print(f"伪逆:\n{A_pinv}")

# 验证伪逆的性质
print(f"\n验证 A @ A_pinv @ A ≈ A:")
print(f"结果:\n{A @ A_pinv @ A}")

# 使用伪逆求解最小二乘问题
b = np.array([7, 8])
x = A_pinv @ b
print(f"\n最小二乘解: {x}")

批量矩阵运算

NumPy 的线性代数函数支持批量处理多个矩阵：

import numpy as np

# 批量矩阵（3 个 2x2 矩阵）
A_batch = np.random.rand(3, 2, 2)
B_batch = np.random.rand(3, 2, 2)

# 批量矩阵乘法
C_batch = A_batch @ B_batch
print(f"批量矩阵乘法结果形状: {C_batch.shape}")

# 批量求逆
A_inv_batch = np.linalg.inv(A_batch)
print(f"批量求逆结果形状: {A_inv_batch.shape}")

# 批量求行列式
det_batch = np.linalg.det(A_batch)
print(f"批量行列式: {det_batch}")

# 批量求解线性方程组
b_batch = np.random.rand(3, 2)
x_batch = np.linalg.solve(A_batch, b_batch)
print(f"批量求解结果形状: {x_batch.shape}")

实战示例

示例1：主成分分析（PCA）

使用特征值分解实现 PCA：

import numpy as np

# 生成示例数据
np.random.seed(42)
X = np.random.randn(100, 5)  # 100 个样本，5 个特征

# 步骤1: 标准化
X_centered = X - X.mean(axis=0)

# 步骤2: 计算协方差矩阵
cov_matrix = np.cov(X_centered, rowvar=False)

# 步骤3: 特征值分解
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

# 步骤4: 按特征值排序（降序）
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# 步骤5: 选择前 k 个主成分
k = 2
principal_components = eigenvectors[:, :k]

# 步骤6: 投影数据
X_pca = X_centered @ principal_components

print(f"原始数据形状: {X.shape}")
print(f"PCA 后形状: {X_pca.shape}")
print(f"\n各主成分解释的方差比例:")
total_variance = sum(eigenvalues)
for i in range(k):
    print(f"PC{i+1}: {eigenvalues[i] / total_variance:.4f}")

示例2：线性回归

使用最小二乘法实现多元线性回归：

import numpy as np

# 生成数据
np.random.seed(42)
n_samples = 100
n_features = 3

X = np.random.randn(n_samples, n_features)
true_coef = np.array([1.5, -2.0, 0.5])
y = X @ true_coef + np.random.randn(n_samples) * 0.1

# 添加截距项
X_with_intercept = np.column_stack([np.ones(n_samples), X])

# 最小二乘求解
result = np.linalg.lstsq(X_with_intercept, y, rcond=None)
coef = result[0]

print(f"真实系数: {true_coef}")
print(f"估计系数: {coef[1:]}")
print(f"截距: {coef[0]:.4f}")

# 计算预测和 R²
y_pred = X_with_intercept @ coef
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - y.mean()) ** 2)
r_squared = 1 - ss_res / ss_tot

print(f"\nR²: {r_squared:.4f}")

示例3：协方差矩阵分析

import numpy as np

# 生成相关的多维数据
np.random.seed(42)
n = 1000

# 创建有相关性的数据
X = np.random.randn(n, 3)
X[:, 1] = 0.5 * X[:, 0] + np.sqrt(0.75) * X[:, 1]  # X1 与 X0 相关
X[:, 2] = 0.3 * X[:, 0] + 0.3 * X[:, 1] + np.sqrt(0.82) * X[:, 2]  # X2 与 X0, X1 相关

# 计算协方差矩阵
cov = np.cov(X, rowvar=False)
print(f"协方差矩阵:\n{cov}")

# 计算相关系数矩阵
corr = np.corrcoef(X, rowvar=False)
print(f"\n相关系数矩阵:\n{corr}")

# 特征值分解分析
eigenvalues, eigenvectors = np.linalg.eigh(corr)
print(f"\n特征值（升序）: {eigenvalues}")
print(f"总方差: {sum(eigenvalues):.4f}")

性能优化建议

使用专门的函数

import numpy as np

A = np.random.rand(100, 100)

# 不要这样做
det_via_eig = np.prod(np.linalg.eigvals(A))

# 应该这样做
det_direct = np.linalg.det(A)

利用矩阵的特殊结构

import numpy as np
import time

# 对称正定矩阵
A = np.random.rand(500, 500)
A = A @ A.T + np.eye(500) * 0.1

# 通用求逆
start = time.time()
inv_general = np.linalg.inv(A)
time_general = time.time() - start

# 使用 Cholesky 分解（更快）
start = time.time()
L = np.linalg.cholesky(A)
inv_cholesky = np.linalg.inv(L.T) @ np.linalg.inv(L)
time_cholesky = time.time() - start

print(f"通用求逆: {time_general:.4f}s")
print(f"Cholesky 方法: {time_cholesky:.4f}s")
print(f"结果相近: {np.allclose(inv_general, inv_cholesky)}")

小结

本章介绍了 NumPy 线性代数模块的主要功能：

类别	函数	用途
矩阵乘法	@, matmul, dot	矩阵乘法、点积
分解	svd, eig, eigh, qr, cholesky	矩阵分解
求解	solve, lstsq	解线性方程组
属性	det, matrix_rank, cond, norm, trace	矩阵特性
逆	inv, pinv	矩阵求逆和伪逆

理解这些线性代数操作，是进行数据分析、机器学习和科学计算的基础。

练习

1. 实现一个函数判断矩阵是否正定
2. 使用 SVD 实现图像压缩
3. 手动实现矩阵的 QR 分解
4. 使用最小二乘法拟合多项式曲线
5. 实现一个简单的 PCA 类