NumPy: 线性代数

Q: @ 和 dot 有什么区别？

对于 2D 矩阵，三者完全等价。区别在高维： np.dot 不广播（将最后两维做矩阵乘）， np.matmul 和 @ 支持广播（前几维自动扩展）。推荐用 @ ——简洁、广播友好、可读性最佳。

Q: 什么矩阵没有逆？

行列式为 0 的方阵（奇异矩阵）没有逆，例如行/列线性相关的矩阵。非方阵也没有逆（但可以用 pinv 求伪逆）。 inv 会抛 LinAlgError ，此时用 pinv 或 solve 。

Q: SVD 和 PCA 有什么关系？

PCA 本质上就是对中心化数据的协方差矩阵做特征分解。而 SVD 直接对数据矩阵分解， X = UΣVᵀ ，其中 V 的列就是主成分方向，奇异值的平方正比于特征值。SVD 比 eig 更稳定，是 PCA 的推荐实现方式。

Q: solve 和 inv(A) @ b 有什么区别？

solve(A, b) 直接求解方程组，不计算逆矩阵； inv(A) @ b 先求逆再乘。solve 更快（避免冗余计算）、更稳定（避免逆矩阵的数值误差放大）。原则上：需要解方程时永远用 solve ，只在确实需要逆矩阵本身时才用 inv 。

线性代数

1. 你将学到

❶ 矩阵乘法 ❷ 转置与逆 ❸ 行列式与迹 ❹ 特征值 ❺ SVD

2. 故事

Charlie 给 1000 用户推荐电影：用户偏好(1000×50) × 电影特征(50×500) = 评分预测(1000×500)。"线性代数不是数学课，是推荐系统/图像压缩/PCA 的底层引擎。"

3. 矩阵乘法

(1) 三种写法：dot / matmul / @

NumPy 提供三种矩阵乘法写法，结果完全相同：

PYTHON

import numpy as np

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

B = np.array([[5, 6],
              [7, 8]])

# Three equivalent ways
r1 = np.dot(A, B)
r2 = np.matmul(A, B)
r3 = A @ B

print(r1)
# [[19 22]
#  [43 50]]

print(np.array_equal(r1, r2), np.array_equal(r2, r3))
# True True

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(2) dot vs matmul vs @ 对比

特性	`np.dot`	`np.matmul`	`@`
语法风格	函数调用	函数调用	中缀运算符
2D 矩阵乘	✅	✅	✅
1D 内积	✅ 返回标量	✅ 返回标量	✅
高维广播	❌ 不广播	✅ 广播	✅ 广播
1D×2D	自动提升	自动提升	自动提升
可读性	一般	一般	最佳

⚠️ 注意： 2D 矩阵乘法三者等价，推荐使用 @ 运算符——最简洁直观。

(3) 矩阵乘 vs 逐元素乘

这是初学者最常见的混淆：

运算	符号	含义	形状要求
逐元素乘	`*`	对应位置相乘	形状相同（或可广播）
矩阵乘	`@`	行列内积累加	`(m,n) @ (n,p) = (m,p)`

PYTHON

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

B = np.array([[5, 6],
              [7, 8]])

print("Element-wise (A * B):")
print(A * B)
# [[ 5 12]
#  [21 32]]

print("Matrix multiply (A @ B):")
print(A @ B)
# [[19 22]
#  [43 50]]

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：矩阵乘法与@运算符（难度⭐）

PYTHON

import numpy as np

# Movie recommendation: user preferences x movie features
users = np.random.rand(1000, 50)    # 1000 users, 50 features
movies = np.random.rand(50, 500)    # 50 features, 500 movies

# Predict ratings: (1000, 50) @ (50, 500) = (1000, 500)
ratings = users @ movies
print(f"Ratings shape: {ratings.shape}")  # (1000, 500)

# Each element: dot product of one user's preferences and one movie's features
print(f"User 0, Movie 0 rating: {ratings[0, 0]:.4f}")

# Manual verification
manual = np.dot(users[0], movies[:, 0])
print(f"Manual calculation: {manual:.4f}")
print(f"Match: {np.isclose(ratings[0, 0], manual)}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

4. 转置与逆

(1) 转置

转置交换行与列，.T 是最简洁的写法：

PYTHON

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

print(A.T)
# [[1 4]
#  [2 5]
#  [3 6]]

print(f"Original shape: {A.shape}")    # (2, 3)
print(f"Transpose shape: {A.T.shape}")  # (3, 2)

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： .T 返回视图（不复制数据），修改转置数组会影响原数组。

(2) 逆矩阵

np.linalg.inv 计算方阵的逆矩阵。若 A 可逆，则 A @ inv(A) = I：

PYTHON

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

A_inv = np.linalg.inv(A)
print(A_inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

# Verify: A @ A_inv should be identity matrix
I = A @ A_inv
print(np.round(I, 10))
# [[1. 0.]
#  [0. 1.]]

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(3) 伪逆 pinv

奇异矩阵（不可逆）没有逆矩阵，但有伪逆（Moore-Penrose 伪逆）：

PYTHON

# Singular matrix (not invertible)
A = np.array([[1, 2],
              [2, 4]])  # row 2 = 2 * row 1

try:
    np.linalg.inv(A)  # LinAlgError!
except np.linalg.LinAlgError as e:
    print(f"inv failed: {e}")

# Pseudo-inverse always works
A_pinv = np.linalg.pinv(A)
print(f"pinv shape: {A_pinv.shape}")
# (2, 2)

# A @ pinv(A) @ A ≈ A (not identity, but projection)
result = A @ A_pinv @ A
print(f"A @ pinv(A) @ A ≈ A? {np.allclose(result, A)}")
# True

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：逆矩阵与方程组求解（难度⭐⭐）

PYTHON

import numpy as np

# Solve linear system: Ax = b
# 2x + y = 5
#  x + 3y = 10

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

b = np.array([5, 10])

# Method 1: Using inverse (conceptually clear, numerically less stable)
x_inv = np.linalg.inv(A) @ b
print(f"Solution via inv: x = {x_inv}")

# Method 2: Using solve (preferred — faster and more stable)
x_solve = np.linalg.solve(A, b)
print(f"Solution via solve: x = {x_solve}")

# Verify
print(f"Ax = {A @ x_solve}")
print(f"Matches b? {np.allclose(A @ x_solve, b)}")

# Method 3: solve is preferred because:
# - inv computes full inverse (O(n³)), then multiply (O(n²))
# - solve directly finds solution (O(n³) but with smaller constant)
# - solve is numerically more stable

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

5. 行列式与迹

(1) 行列式 det

行列式衡量矩阵"体积缩放因子"——det 为 0 表示矩阵不可逆：

PYTHON

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

d = np.linalg.det(A)
print(f"det(A) = {d:.4f}")  # -2.0

B = np.array([[1, 2],
              [2, 4]])  # singular

d_b = np.linalg.det(B)
print(f"det(B) = {d_b:.10f}")  # ≈ 0 (singular)

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(2) 迹 trace

迹是主对角线元素之和：

PYTHON

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

t = np.trace(A)
print(f"trace(A) = {t}")  # 1 + 5 + 9 = 15

# Equivalent to:
t_manual = sum(A[i, i] for i in range(A.shape[0]))
print(f"Manual trace: {t_manual}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：行列式与迹（难度⭐）

PYTHON

import numpy as np

# 3x3 matrix
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])

det = np.linalg.det(A)
trace = np.trace(A)

print(f"A =\n{A}")
print(f"det(A) = {det:.4f}")    # -3.0
print(f"trace(A) = {trace}")    # 2 + (-1) + 2 = 3

# det != 0 means A is invertible
A_inv = np.linalg.inv(A)
print(f"A_inv @ A ≈ I? {np.allclose(A_inv @ A, np.eye(3))}")

# Geometric meaning: |det| = volume scaling factor
# det = -3 means: orientation reversed, volume scaled by 3

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

6. 特征值与特征向量

(1) eig 基本用法

若 Av = λv，则 λ 是特征值，v 是特征向量：

PYTHON

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

eigenvalues, eigenvectors = np.linalg.eig(A)

print(f"Eigenvalues: {eigenvalues}")
# [5. 2.]

print(f"Eigenvectors:\n{eigenvectors}")
# Each column is an eigenvector

# Verify: A @ v = λ * v
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lam = eigenvalues[i]
    Av = A @ v
    lv = lam * v
    print(f"λ={lam:.1f}: A@v ≈ λ*v? {np.allclose(Av, lv)}")
# Both True

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(2) 特征值的实际意义

主成分分析（PCA）：特征值越大，对应方向上的方差越大（信息量越多）
振动分析：特征值 = 自然频率的平方
网页排名（PageRank）：最大特征值对应的特征向量就是网页排名

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：特征值与PCA直觉（难度⭐⭐）

PYTHON

import numpy as np

# Simulate 2D data with different spread along axes
np.random.seed(42)
data = np.random.randn(500, 2)
# Stretch along one direction
transform = np.array([[3, 1],
                      [1, 1]])
data = data @ transform.T  # shape (500, 2)

# Covariance matrix
cov = np.cov(data, rowvar=False)
print(f"Covariance matrix:\n{cov}")

# Eigen decomposition of covariance matrix
eigenvalues, eigenvectors = np.linalg.eig(cov)
print(f"\nEigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# The largest eigenvalue's eigenvector = direction of max variance
idx = np.argsort(eigenvalues)[::-1]  # sort descending
print(f"\nPrincipal component direction: {eigenvectors[:, idx[0]]}")
print(f"Variance explained ratio: {eigenvalues[idx[0]] / eigenvalues.sum():.2%}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

7. SVD 奇异值分解

(1) SVD 原理

SVD 将任意矩阵 A 分解为 U·Σ·Vᵀ，无论方阵还是非方阵、可逆还是不可逆：

PYTHON

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

U, s, Vt = np.linalg.svd(A, full_matrices=False)

print(f"U shape: {U.shape}")    # (2, 2)
print(f"s shape: {s.shape}")    # (2,)  singular values
print(f"Vt shape: {Vt.shape}")  # (2, 3)

# Reconstruct: A = U @ diag(s) @ Vt
S = np.diag(s)
A_reconstructed = U @ S @ Vt
print(f"Reconstruction match: {np.allclose(A, A_reconstructed)}")
# True

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(2) SVD 流程

graph LR
    A["矩阵 A<br/>(m×n)"] --> U["U<br/>(m×k)"]
    A --> S["Σ<br/>(k×k)"]
    A --> V["Vᵀ<br/>(k×n)"]
    U --> R["U @ Σ @ Vᵀ<br/>= A"]
    S --> R
    V --> R

    S --> Low["取前 r 个奇异值"]
    Low --> Comp["压缩矩阵<br/>A ≈ U_r @ Σ_r @ Vᵀ_r"]
    Comp --> App["应用：图像压缩<br/>推荐系统<br/>PCA"]

    style A fill:#4CAF50,color:#fff
    style U fill:#2196F3,color:#fff
    style S fill:#FF9800,color:#fff
    style V fill:#9C27B0,color:#fff
    style Comp fill:#F44336,color:#fff
    style App fill:#00BCD4,color:#fff

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(3) 分解方法对比

方法	函数	适用矩阵	分解形式	典型用途
特征分解	`eig`	方阵	Av = λv	PCA、振动分析
SVD	`svd`	任意矩阵	A = UΣVᵀ	图像压缩、推荐、LSA
Cholesky	`cholesky`	正定方阵	A = LLᵀ	高效求逆、采样
QR	`qr`	任意矩阵	A = QR	最小二乘、特征值算法
LU	无直接函数	方阵	A = LU	求解线性方程组

(4) 矩阵运算速查

运算	函数	输入	输出
矩阵乘	`A @ B`	(m,n), (n,p)	(m,p)
转置	`A.T`	(m,n)	(n,m)
逆	`inv(A)`	(n,n)	(n,n)
伪逆	`pinv(A)`	(m,n)	(n,m)
行列式	`det(A)`	(n,n)	标量
迹	`trace(A)`	(n,n)	标量
特征值	`eig(A)`	(n,n)	λ, V
SVD	`svd(A)`	(m,n)	U, s, Vt
求解 Ax=b	`solve(A, b)`	(n,n), (n,)	(n,)
范数	`norm(A)`	任意	标量
条件数	`cond(A)`	(n,n)	标量

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：SVD 图像压缩原理（难度⭐⭐⭐）

PYTHON

import numpy as np

# Create a 100x100 grayscale image (gradient + pattern)
x = np.linspace(0, 2 * np.pi, 100)
y = np.linspace(0, 2 * np.pi, 100)
X, Y = np.meshgrid(x, y)
image = (np.sin(X) + np.cos(Y) + 2) / 4 * 255  # range [0, 255]
image = image.astype(np.float64)

print(f"Image shape: {image.shape}")

# SVD decomposition
U, s, Vt = np.linalg.svd(image, full_matrices=False)
print(f"U: {U.shape}, s: {s.shape}, Vt: {Vt.shape}")
print(f"Top 5 singular values: {s[:5]}")

# Total information = sum of all singular values squared
total_energy = np.sum(s ** 2)

# Compress with different ranks
for rank in [5, 10, 20, 50]:
    # Keep only top-rank singular values
    U_r = U[:, :rank]
    s_r = s[:rank]
    Vt_r = Vt[:rank, :]

    # Reconstruct
    compressed = U_r @ np.diag(s_r) @ Vt_r

    # Compression ratio
    original_size = 100 * 100
    compressed_size = 100 * rank + rank + rank * 100
    ratio = compressed_size / original_size

    # Energy retained
    energy = np.sum(s_r ** 2) / total_energy

    # Reconstruction error
    error = np.linalg.norm(image - compressed) / np.linalg.norm(image)

    print(f"rank={rank:3d}: ratio={ratio:.2%}, "
          f"energy={energy:.2%}, error={error:.4f}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

8. solve / norm / cond

(1) solve 求解线性方程组

solve 比 inv(A) @ b 更快更稳定：

PYTHON

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

x = np.linalg.solve(A, b)
print(f"Solution: {x}")         # [2. 3.]
print(f"Verify: A@x = {A @ x}") # [9. 8.]

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(2) norm 范数

PYTHON

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

print(np.linalg.norm(v))           # L2 norm: 5.0
print(np.linalg.norm(v, ord=1))    # L1 norm: 7.0
print(np.linalg.norm(v, ord=np.inf))  # L-inf norm: 4.0

A = np.array([[1, 2],
              [3, 4]])
print(np.linalg.norm(A))           # Frobenius norm: 5.477...
print(np.linalg.norm(A, 'fro'))    # Same as above

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

(3) cond 条件数

条件数衡量矩阵"对扰动的敏感度"——越大越不稳定：

PYTHON

# Well-conditioned matrix
A_good = np.array([[1, 0],
                   [0, 1]])
print(f"cond(I) = {np.linalg.cond(A_good):.1f}")  # 1.0

# Ill-conditioned matrix (Hilbert-like)
A_bad = np.array([[1, 1],
                  [1, 1.0001]])
print(f"cond(A_bad) = {np.linalg.cond(A_bad):.1f}")  # ~40000

# Small perturbation → large error in solution
b = np.array([1, 1])
x1 = np.linalg.solve(A_bad, b)

b_perturbed = np.array([1, 1.0001])
x2 = np.linalg.solve(A_bad, b_perturbed)

print(f"x1 = {x1}")
print(f"x2 = {x2}")
print(f"Relative change in b: {np.linalg.norm(b_perturbed - b) / np.linalg.norm(b):.6f}")
print(f"Relative change in x: {np.linalg.norm(x2 - x1) / np.linalg.norm(x1):.2f}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：条件数与数值稳定性（难度⭐⭐）

PYTHON

import numpy as np

# Compare solve vs inv for ill-conditioned system
n = 10
np.random.seed(0)

# Create a moderately ill-conditioned matrix
A = np.random.randn(n, n)
U, s, Vt = np.linalg.svd(A)
s[-1] = s[0] * 1e-8  # Make smallest singular value tiny → large cond
A = U @ np.diag(s) @ Vt

b = np.random.randn(n)

print(f"cond(A) = {np.linalg.cond(A):.2e}")

# Method 1: inv
x_inv = np.linalg.inv(A) @ b
res_inv = np.linalg.norm(A @ x_inv - b)

# Method 2: solve
x_solve = np.linalg.solve(A, b)
res_solve = np.linalg.norm(A @ x_solve - b)

# Method 3: pinv (more robust for ill-conditioned)
x_pinv = np.linalg.pinv(A) @ b
res_pinv = np.linalg.norm(A @ x_pinv - b)

print(f"Residual (inv):   {res_inv:.2e}")
print(f"Residual (solve): {res_solve:.2e}")
print(f"Residual (pinv):  {res_pinv:.2e}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

9. 综合示例

▶ 示例

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

：SVD 压缩 100×100 灰度图像（难度⭐⭐⭐）

Charlie 要把一张 100×100 灰度图像压缩传输。用 SVD 保留不同数量的奇异值，对比压缩率和重建质量：

PYTHON

import numpy as np

np.random.seed(42)

# Generate a 100x100 synthetic grayscale image
x = np.linspace(0, 2 * np.pi, 100)
y = np.linspace(0, 2 * np.pi, 100)
X, Y = np.meshgrid(x, y)

# Composite pattern: multiple frequency components
image = (np.sin(X) * np.cos(Y) +
         0.5 * np.sin(3 * X) * np.cos(2 * Y) +
         0.3 * np.cos(5 * X + Y) + 2) / 3 * 255
image = image.astype(np.float64)

print(f"Original image shape: {image.shape}")
print(f"Original data size: {image.nbytes / 1024:.1f} KB")

# Step 1: SVD decomposition
U, s, Vt = np.linalg.svd(image, full_matrices=False)
print(f"\nSingular values (top 10): {np.round(s[:10], 2)}")
print(f"Total singular values: {len(s)}")

# Step 2: Analyze energy distribution
total_energy = np.sum(s ** 2)
cumulative = np.cumsum(s ** 2) / total_energy
print(f"\nEnergy captured by top 5: {cumulative[4]:.2%}")
print(f"Energy captured by top 10: {cumulative[9]:.2%}")
print(f"Energy captured by top 20: {cumulative[19]:.2%}")
print(f"Energy captured by top 50: {cumulative[49]:.2%}")

# Step 3: Compress with different ranks
def svd_compress(U, s, Vt, rank):
    U_r = U[:, :rank]
    s_r = s[:rank]
    Vt_r = Vt[:rank, :]
    compressed = U_r @ np.diag(s_r) @ Vt_r
    original_size = U.shape[0] * Vt.shape[1]
    compressed_size = U.shape[0] * rank + rank + rank * Vt.shape[1]
    ratio = compressed_size / original_size
    return compressed, ratio

print(f"\n{'Rank':>5} {'Compression':>12} {'Energy':>8} {'Rel Error':>10}")
print("-" * 40)

for rank in [1, 5, 10, 20, 30, 50]:
    compressed, ratio = svd_compress(U, s, Vt, rank)
    energy = cumulative[rank - 1]
    rel_error = np.linalg.norm(image - compressed) / np.linalg.norm(image)
    print(f"{rank:5d} {ratio:11.2%} {energy:7.2%} {rel_error:10.4f}")

# Step 4: Demonstrate best rank-10 approximation
rank = 10
compressed_10, ratio_10 = svd_compress(U, s, Vt, rank)
print(f"\nRank-{rank} approximation:")
print(f"  Compression ratio: {ratio_10:.2%}")
print(f"  Data saved: {(1 - ratio_10):.2%}")
print(f"  Max pixel error: {np.max(np.abs(image - compressed_10)):.2f}")
print(f"  Mean pixel error: {np.mean(np.abs(image - compressed_10)):.2f}")

# Step 5: Verify SVD properties
# Singular values are non-negative and in descending order
print(f"\nSVD property checks:")
print(f"  All s >= 0: {np.all(s >= 0)}")
print(f"  s is descending: {np.all(np.diff(s) <= 1e-10)}")
print(f"  U columns orthonormal: {np.allclose(U.T @ U, np.eye(U.shape[1]), atol=1e-10)}")
print(f"  Vt rows orthonormal: {np.allclose(Vt @ Vt.T, np.eye(Vt.shape[0]), atol=1e-10)}")

TEXT 📖 仅展示

> **输出：** 在本地 Python 环境运行 NumPy 2.x，输出 ndarray 数组内容。Piston 服务器未预装 NumPy，请在本机安装（`pip install numpy`）后实操对照。实际数值可能因 NumPy 版本、随机种子略有差异。

⚠️ 注意： 以下代码需在本地 Python 环境中运行。

❓ 常见问题

Q @ 和 dot 有什么区别？

A 对于 2D 矩阵，三者完全等价。区别在高维：np.dot 不广播（将最后两维做矩阵乘），np.matmul 和 @ 支持广播（前几维自动扩展）。推荐用 @——简洁、广播友好、可读性最佳。

Q 什么矩阵没有逆？

A 行列式为 0 的方阵（奇异矩阵）没有逆，例如行/列线性相关的矩阵。非方阵也没有逆（但可以用 pinv 求伪逆）。inv 会抛 LinAlgError，此时用 pinv 或 solve。

Q 特征值有什么实际意义？

A 特征值衡量矩阵在对应特征向量方向上的"拉伸力度"。在 PCA 中，协方差矩阵的最大特征值对应最大方差方向；在 PageRank 中，最大特征值的特征向量给出网页排名；在振动分析中，特征值是自然频率的平方。

Q SVD 和 PCA 有什么关系？

A PCA 本质上就是对中心化数据的协方差矩阵做特征分解。而 SVD 直接对数据矩阵分解，X = UΣVᵀ，其中 V 的列就是主成分方向，奇异值的平方正比于特征值。SVD 比 eig 更稳定，是 PCA 的推荐实现方式。

Q 条件数大意味着什么？

A 条件数 = 最大奇异值 / 最小奇异值，衡量矩阵对扰动的敏感度。条件数越大，求解 Ax=b 时 b 的微小变化会导致 x 的巨大变化。实用建议：cond > 10¹⁵ 基本无法求解；cond > 10⁸ 需要谨慎；cond ≈ 1 最稳定。

Q solve 和 inv(A) @ b 有什么区别？

A solve(A, b) 直接求解方程组，不计算逆矩阵；inv(A) @ b 先求逆再乘。solve 更快（避免冗余计算）、更稳定（避免逆矩阵的数值误差放大）。原则上：需要解方程时永远用 solve，只在确实需要逆矩阵本身时才用 inv。

📖 小节

矩阵乘法推荐用 @ 运算符，与 np.dot/np.matmul 在 2D 下等价
* 是逐元素乘，@ 是矩阵乘——最常见混淆点
.T 返回视图（零拷贝），inv 返回新数组
不可逆矩阵用 pinv（伪逆），总是存在且稳定
solve(A, b) 比 inv(A) @ b 更快更稳定，优先使用
行列式为 0 → 矩阵不可逆；迹 = 对角线元素之和
特征分解 eig 只适用于方阵；SVD 适用于任意矩阵
SVD 压缩：保留前 r 个奇异值，压缩率 = (m×r + r + r×n) / (m×n)
条件数衡量数值稳定性，越大越不稳定
范数 norm 支持多种阶数：L1、L2（默认）、L-inf、Frobenius

📝 作业

基础题（难度⭐）：创建一个 3×3 矩阵表示 2D 坐标变换（旋转 30°），对 5 个点 [1,0], [0,1], [1,1], [2,0], [0,3] 应用变换。提示：旋转矩阵为 [[cosθ, -sinθ], [sinθ, cosθ]]，用 @ 运算。打印变换前后坐标对比。
进阶题（难度⭐⭐）：用 np.linalg.solve 求解线性方程组：3x + 2y - z = 1，2x - 2y + 4z = -2，-x + 0.5y - z = 0。验证解的正确性，并计算矩阵的条件数判断解的可靠性。
进阶题（难度⭐⭐）：生成一张 100×100 合成灰度图像（自选公式），用 SVD 分别取 rank=5、10、20、50 压缩，打印每个 rank 的压缩率、能量保留比例和相对误差。找出"能量保留 95% 以上"的最小 rank。

Web-Tutorial 技术团队

由多位开发者共同维护的编程教程平台。每篇教程由对应领域的开发者编写和审核，确保内容准确可靠。如发现任何问题，欢迎向我们反馈。

关于我们

NumPy: 线性代数

线性代数

1. 你将学到

2. 故事

3. 矩阵乘法

(1) 三种写法：dot / matmul / @

(2) dot vs matmul vs @ 对比

(3) 矩阵乘 vs 逐元素乘

▶ 示例

4. 转置与逆

(1) 转置

(2) 逆矩阵

(3) 伪逆 pinv

▶ 示例

5. 行列式与迹

(1) 行列式 det

(2) 迹 trace

▶ 示例

6. 特征值与特征向量

(1) eig 基本用法

(2) 特征值的实际意义

▶ 示例

7. SVD 奇异值分解

(1) SVD 原理

(2) SVD 流程

(3) 分解方法对比

(4) 矩阵运算速查

▶ 示例

8. solve / norm / cond

(1) solve 求解线性方程组

(2) norm 范数

(3) cond 条件数

▶ 示例

9. 综合示例

▶ 示例

❓ 常见问题

📖 小节

📝 作业

🙏 帮我们做得更好