文章

35 · 线性代数实战：特征值、SVD 与最小二乘法

#044 · 2026-04-17 · Python
PythonPython教程高级线性代数特征值SVD

🔗 知识图谱导航：阅读本文前，建议先掌握《31 · NumPy 实战》中的矩阵运算基础——本文是 NumPy 线性代数的深入篇，np.linalg 模块是机器学习算法的数学底层。

运行环境：pip install numpy（NumPy 内置 linalg 模块）。

极客解析：线性代数是机器学习的数学语言。特征值分解是 PCA 的核心，SVD 是推荐系统的基础，最小二乘法是线性回归的解析解。理解这三个操作，机器学习的大多数算法都有了数学直觉。

np.linalg 核心函数速查

np.linalg.eig(A)          特征值分解，返回 (eigenvalues, eigenvectors)
np.linalg.svd(A)          奇异值分解，返回 (U, s, Vt)
np.linalg.lstsq(A, b)     最小二乘法，求 Ax=b 的最优解
np.polyfit(x, y, deg)     多项式拟合，返回系数
np.polyval(coeffs, x)     用多项式系数预测
np.linalg.norm(v)         向量/矩阵范数
np.linalg.det(A)          行列式
np.linalg.inv(A)          矩阵求逆

步步为营：核心逻辑自适应拆解

这一篇不机械追求 8 步，而是按完整源码的真实结构拆成 6 个台阶：先准备终端输出工具，再依次进入 EVD、SVD、多项式拟合、最小二乘回归，最后用 CLI 统一调度。每个演示都带 Mock 数据和 print() 反馈。

Step 1：用 ASCII 工具把矩阵和条形图打印清楚

痛点与机制：

这一步先讲输出工具。线性代数结果经常是矩阵和一串数，如果直接 print(array)，新手很容易看懵。ascii_matrix 像把矩阵放进整齐表格，ascii_bar 像给奇异值画小柱状图，让抽象数字先变得可见。

核心源码（逐字来自文末完整源码）：

def section(title: str) -> None:
    print(f"\n{'='*60}\n  {title}\n{'='*60}")


def ascii_matrix(M: np.ndarray, fmt: str = ".3f", max_rows: int = 6) -> None:
    """打印矩阵的 ASCII 表示。"""
    rows = min(M.shape[0], max_rows)
    for i in range(rows):
        row_str = "  ".join(f"{v:{fmt}}" for v in M[i])
        print(f"  [{row_str}]")
    if M.shape[0] > max_rows:
        print(f"  ... ({M.shape[0]-max_rows} 行省略)")


def ascii_bar(values: np.ndarray, label: str = "", width: int = 30) -> None:
    """打印归一化条形图。"""
    max_v = values.max()
    for i, v in enumerate(values[:10]):
        bar = "█" * int(v / max_v * width)
        print(f"  {label}{i:>2}: {bar:<{width}} {v:.4f}")

可运行演示（补齐 Mock 数据与 print 反馈）：

import numpy as np


def section(title: str) -> None:
    print(f"\n{'='*60}\n  {title}\n{'='*60}")


def ascii_matrix(M: np.ndarray, fmt: str = ".3f", max_rows: int = 6) -> None:
    """打印矩阵的 ASCII 表示。"""
    rows = min(M.shape[0], max_rows)
    for i in range(rows):
        row_str = "  ".join(f"{v:{fmt}}" for v in M[i])
        print(f"  [{row_str}]")
    if M.shape[0] > max_rows:
        print(f"  ... ({M.shape[0]-max_rows} 行省略)")


def ascii_bar(values: np.ndarray, label: str = "", width: int = 30) -> None:
    """打印归一化条形图。"""
    max_v = values.max()
    for i, v in enumerate(values[:10]):
        bar = "█" * int(v / max_v * width)
        print(f"  {label}{i:>2}: {bar:<{width}} {v:.4f}")


section("工具函数演示")
ascii_matrix(np.array([[1.2345, 2.0], [3.0, 4.5678]]), fmt=".2f")
ascii_bar(np.array([3.0, 2.0, 1.0]), label="σ")

Step 2：用 mode_eig 找出协方差矩阵的主方向

痛点与机制：

特征值分解回答的是：“这个矩阵最重要的方向在哪里？” 对协方差矩阵来说，最大特征值对应方差最大的方向，也就是 PCA 的第一主成分。可以把它想成在一团数据云里找最长的影子。

核心源码（逐字来自文末完整源码）：

def mode_eig() -> None:
    section("特征值分解（EVD）")

    # 构造一个对称正定矩阵（协方差矩阵形式）
    rng = np.random.RandomState(42)
    X = rng.randn(100, 4)
    cov = X.T @ X / 100   # 4×4 协方差矩阵

    print("\n  协方差矩阵 (4×4):")
    ascii_matrix(cov, fmt=".4f")

    eigvals, eigvecs = np.linalg.eig(cov)
    # 按特征值降序排列
    order = np.argsort(eigvals)[::-1]
    eigvals = eigvals[order].real
    eigvecs = eigvecs[:, order].real

    print(f"\n  特征值（方差解释量）:")
    total = eigvals.sum()
    cumsum = 0.0
    for i, v in enumerate(eigvals):
        cumsum += v
        bar = "█" * int(v / eigvals[0] * 25)
        print(f"  λ{i+1}: {v:.4f}  ({v/total*100:.1f}%)  累计: {cumsum/total*100:.1f}%  {bar}")

    print(f"\n  第1主成分方向（特征向量）: {eigvecs[:, 0].round(4)}")
    print(f"\n  验证 A·v = λ·v:")
    v0 = eigvecs[:, 0]
    lhs = cov @ v0
    rhs = eigvals[0] * v0
    print(f"  A·v₁ = {lhs.round(4)}")
    print(f"  λ₁·v₁= {rhs.round(4)}")
    print(f"  误差: {np.linalg.norm(lhs - rhs):.2e}  ✓")

    print(f"\n  💡 PCA 就是对协方差矩阵做特征值分解，取前 k 个特征向量作为主成分")

可运行演示（补齐 Mock 数据与 print 反馈）：

import numpy as np


def section(title: str) -> None:
    print(f"\n{'='*60}\n  {title}\n{'='*60}")


def ascii_matrix(M: np.ndarray, fmt: str = ".3f", max_rows: int = 6) -> None:
    rows = min(M.shape[0], max_rows)
    for i in range(rows):
        print("  [" + "  ".join(f"{v:{fmt}}" for v in M[i]) + "]")


def mode_eig() -> None:
    section("特征值分解（EVD）")
    rng = np.random.RandomState(42)
    X = rng.randn(20, 3)
    cov = X.T @ X / 20
    print("\n  协方差矩阵:")
    ascii_matrix(cov, fmt=".4f")
    eigvals, eigvecs = np.linalg.eig(cov)
    order = np.argsort(eigvals)[::-1]
    eigvals = eigvals[order].real
    eigvecs = eigvecs[:, order].real
    print("\n  特征值:", eigvals.round(4).tolist())
    v0 = eigvecs[:, 0]
    print("  验证 A·v = λ·v:", np.allclose(cov @ v0, eigvals[0] * v0))
    print("  第一主成分方向:", v0.round(4).tolist())


mode_eig()

Step 3：用 mode_svd 做矩阵低秩近似

痛点与机制：

SVD 把矩阵拆成若干“重要层”。奇异值越大，说明这一层携带的信息越多。低秩近似就像压缩照片：不保留每个细节，只保留最重要的几层结构，文件小了，但主体仍然看得出来。

核心源码（逐字来自文末完整源码）：

def mode_svd() -> None:
    section("奇异值分解（SVD）— 图像压缩模拟")

    # 生成一个模拟"图像"矩阵（32×32，有结构的数据）
    rng = np.random.RandomState(42)
    x = np.linspace(0, 2 * np.pi, 32)
    # 构造低秩矩阵（模拟图像的低频结构）
    img = (
        np.outer(np.sin(x), np.cos(x)) * 128 +
        np.outer(np.cos(2*x), np.sin(2*x)) * 64 +
        rng.randn(32, 32) * 10   # 噪声
    )
    img = np.clip(img + 128, 0, 255)

    print(f"\n  原始矩阵: {img.shape}  秩: {np.linalg.matrix_rank(img)}")

    # SVD 分解
    U, s, Vt = np.linalg.svd(img, full_matrices=False)
    print(f"\n  SVD 分解: U{U.shape} · Σ({len(s)},) · Vᵀ{Vt.shape}")

    print(f"\n  奇异值分布（前10个）:")
    ascii_bar(s[:10], "σ")

    # 用不同 k 重建矩阵，计算压缩比和误差
    print(f"\n  {'k（保留奇异值数）':<20} {'压缩比':<10} {'重建误差(RMSE)':<18} {'信息保留率'}")
    print(f"  {'─'*65}")
    total_params = img.shape[0] * img.shape[1]
    for k in [1, 2, 4, 8, 16, 32]:
        # 重建：取前 k 个奇异值
        img_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
        compressed_params = k * (img.shape[0] + img.shape[1] + 1)
        ratio = total_params / compressed_params
        rmse = float(np.sqrt(np.mean((img - img_k) ** 2)))
        info = float((s[:k] ** 2).sum() / (s ** 2).sum() * 100)
        bar = "█" * int(info / 100 * 20)
        print(f"  k={k:<18} {ratio:.1f}x       {rmse:.2f}{'':12} {info:.1f}%  {bar}")

    print(f"\n  💡 k=4 时保留 ~95% 信息，压缩比 ~4x——这就是 JPEG 压缩的数学基础")

    # 验证 SVD 重建
    img_full = U @ np.diag(s) @ Vt
    print(f"\n  验证 U·Σ·Vᵀ = 原矩阵，误差: {np.linalg.norm(img - img_full):.2e}  ✓")

# ─── 模式3：多项式拟合 ─────────────────────────────────────────────────────────

可运行演示（补齐 Mock 数据与 print 反馈）：

import numpy as np


def ascii_bar(values: np.ndarray, label: str = "", width: int = 30) -> None:
    max_v = values.max()
    for i, v in enumerate(values[:10]):
        bar = "█" * int(v / max_v * width)
        print(f"  {label}{i:>2}: {bar:<{width}} {v:.4f}")


def mode_svd() -> None:
    print("\n=== 奇异值分解（SVD）— 低秩近似 ===")
    rng = np.random.RandomState(42)
    A = rng.randn(8, 6)
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    print(f"  A{A.shape} = U{U.shape} · Σ({len(s)},) · Vᵀ{Vt.shape}")
    print("\n  奇异值分布:")
    ascii_bar(s, "σ")
    for k in [1, 2, 4, 6]:
        A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
        info = (s[:k] ** 2).sum() / (s ** 2).sum() * 100
        rmse = np.sqrt(np.mean((A - A_k) ** 2))
        print(f"  k={k}: 信息保留 {info:5.1f}% | RMSE={rmse:.4f}")


mode_svd()

Step 4：用 mode_poly 拟合趋势线并预测未来

痛点与机制：

多项式拟合是在一堆带噪声的点里找趋势线。np.polyfit 负责找系数，np.polyval 负责拿系数去预测。它像用尺子沿着散点摆出一条“最顺眼”的趋势路径，但阶数太高也容易把噪声当规律。

核心源码（逐字来自文末完整源码）：

def mode_poly() -> None:
    section("多项式拟合 — 股价趋势线")

    # 生成模拟股价数据（带趋势+周期+噪声）
    rng = np.random.RandomState(42)
    n = 60
    t = np.arange(n, dtype=float)
    price = (
        100 +
        0.5 * t +                          # 上升趋势
        5 * np.sin(2 * np.pi * t / 20) +  # 周期波动
        rng.randn(n) * 3                   # 噪声
    )

    print(f"\n  模拟股价数据: {n} 个交易日，均价 {price.mean():.2f}")

    # 不同阶数的多项式拟合
    print(f"\n  {'阶数':<8} {'拟合误差(RMSE)':<18} {'系数'}")
    print(f"  {'─'*60}")
    for deg in [1, 2, 3, 5]:
        coeffs = np.polyfit(t, price, deg)
        fitted = np.polyval(coeffs, t)
        rmse = float(np.sqrt(np.mean((price - fitted) ** 2)))
        coeff_str = "  ".join(f"{c:.3f}" for c in coeffs[:3])
        print(f"  {deg}阶      {rmse:.4f}{'':12} [{coeff_str}...]")

    # 用1阶（线性趋势线）做预测
    coeffs_1 = np.polyfit(t, price, 1)
    slope, intercept = coeffs_1
    print(f"\n  线性趋势线: price = {slope:.4f}·t + {intercept:.4f}")
    print(f"  日均涨幅: {slope:.4f}  （原始设定: 0.5）")

    # 预测未来5天
    future_t = np.arange(n, n + 5, dtype=float)
    future_price = np.polyval(coeffs_1, future_t)
    print(f"\n  未来5天预测（线性外推）:")
    for i, (day, p) in enumerate(zip(future_t, future_price)):
        print(f"  第{int(day)+1}天: {p:.2f}")

    # ASCII 价格走势图
    print(f"\n  股价走势（ASCII，每5天采样）:")
    min_p, max_p = price.min(), price.max()
    H = 8
    for row in range(H, -1, -1):
        threshold = min_p + (max_p - min_p) * row / H
        line = ""
        for i in range(0, n, 5):
            if price[i] >= threshold:
                line += "█"
            else:
                line += "·"
        label = f"{threshold:6.1f} │"
        print(f"  {label}{line}")
    print(f"  {'':8}└{'─'*12}")

# ─── 模式4：最小二乘 ───────────────────────────────────────────────────────────

可运行演示（补齐 Mock 数据与 print 反馈）：

import numpy as np


def mode_poly() -> None:
    print("\n=== 多项式拟合 — 趋势线 ===")
    rng = np.random.RandomState(42)
    t = np.arange(20, dtype=float)
    price = 100 + 0.5 * t + rng.randn(20) * 2
    for deg in [1, 2, 3]:
        coeffs = np.polyfit(t, price, deg)
        fitted = np.polyval(coeffs, t)
        rmse = float(np.sqrt(np.mean((price - fitted) ** 2)))
        print(f"  {deg}阶拟合 RMSE={rmse:.3f} | 前两个系数={coeffs[:2].round(3).tolist()}")
    coeffs_1 = np.polyfit(t, price, 1)
    future_t = np.arange(20, 23, dtype=float)
    future_price = np.polyval(coeffs_1, future_t)
    print("  未来3天预测:", future_price.round(2).tolist())


mode_poly()

Step 5：用 mode_lstsq 求多变量线性回归系数

痛点与机制：

最小二乘法是线性回归的数学核心：当方程太多、无法每条都完美满足时，它找一组系数，让整体误差平方和最小。可以把它想成调音台，三个特征像三个旋钮，lstsq 帮你调到最贴近真实声音的位置。

核心源码（逐字来自文末完整源码）：

def mode_lstsq() -> None:
    section("最小二乘（lstsq）— 多变量线性回归")

    # 生成多变量数据：y = 2x₁ + 3x₂ - x₃ + 5 + 噪声
    rng = np.random.RandomState(42)
    n = 200
    X_raw = rng.randn(n, 3)
    true_coeffs = np.array([2.0, 3.0, -1.0])
    intercept = 5.0
    y = X_raw @ true_coeffs + intercept + rng.randn(n) * 0.5

    # 构造设计矩阵（加偏置列）
    A = np.column_stack([X_raw, np.ones(n)])

    # 最小二乘求解
    coeffs, residuals, rank, sv = np.linalg.lstsq(A, y, rcond=None)

    print(f"\n  真实系数: {true_coeffs.tolist()} + 截距 {intercept}")
    print(f"  拟合系数: {coeffs[:3].round(4).tolist()} + 截距 {coeffs[3]:.4f}")
    print(f"  矩阵秩: {rank}  最大奇异值: {sv[0]:.4f}")

    # 预测误差
    y_pred = A @ coeffs
    rmse = float(np.sqrt(np.mean((y - y_pred) ** 2)))
    r2 = float(1 - np.sum((y - y_pred)**2) / np.sum((y - y.mean())**2))
    print(f"\n  RMSE: {rmse:.4f}  R²: {r2:.4f}")

    # 对比 np.polyfit（1D）vs lstsq（多维）
    print(f"\n  lstsq vs polyfit 对比:")
    print(f"  ┌──────────────┬──────────────────────────────────────┐")
    print(f"  │ np.polyfit   │ 只支持 1D 输入，多项式拟合           │")
    print(f"  │ np.linalg.lstsq │ 支持多维输入，通用线性方程组求解 │")
    print(f"  └──────────────┴──────────────────────────────────────┘")

    # 条件数（衡量矩阵病态程度）
    cond = np.linalg.cond(A)
    print(f"\n  设计矩阵条件数: {cond:.2f}  （<1000 为良态，>1e6 为病态）")
    print(f"  💡 条件数大 → 数值不稳定 → 需要正则化（Ridge/Lasso）")

# ─── 入口 ─────────────────────────────────────────────────────────────────────

可运行演示（补齐 Mock 数据与 print 反馈）：

import numpy as np


def mode_lstsq() -> None:
    print("\n=== 最小二乘（lstsq）— 多变量线性回归 ===")
    rng = np.random.RandomState(42)
    n = 80
    X_raw = rng.randn(n, 3)
    true_coeffs = np.array([2.0, 3.0, -1.0])
    intercept = 5.0
    y = X_raw @ true_coeffs + intercept + rng.randn(n) * 0.5
    A = np.column_stack([X_raw, np.ones(n)])
    coeffs, residuals, rank, sv = np.linalg.lstsq(A, y, rcond=None)
    y_pred = A @ coeffs
    rmse = float(np.sqrt(np.mean((y - y_pred) ** 2)))
    r2 = float(1 - np.sum((y - y_pred)**2) / np.sum((y - y.mean())**2))
    print("  真实系数:", true_coeffs.tolist(), "+ 截距", intercept)
    print("  拟合系数:", coeffs[:3].round(3).tolist(), "+ 截距", round(float(coeffs[3]), 3))
    print(f"  矩阵秩: {rank} | RMSE={rmse:.4f} | R²={r2:.4f}")


mode_lstsq()

Step 6：用 main 做 eig/svd/poly/lstsq/all 脚本遥控器

痛点与机制：

main 是脚本遥控器。用户不用改代码，只要传 --mode eig/svd/poly/lstsq/all 就能切换功能。这个结构让教程代码从“演示片段”升级成真正可运行的小工具。

核心源码（逐字来自文末完整源码）：

def main() -> None:
    parser = argparse.ArgumentParser(description="NumPy 线性代数深度实战")
    parser.add_argument(
        "--mode",
        choices=["eig", "svd", "poly", "lstsq", "all"],
        default="all",
    )
    args = parser.parse_args()
    dispatch = {
        "eig":   mode_eig,
        "svd":   mode_svd,
        "poly":  mode_poly,
        "lstsq": mode_lstsq,
        "all":   lambda: [mode_eig(), mode_svd(), mode_poly(), mode_lstsq()],
    }
    dispatch[args.mode]()


if __name__ == "__main__":
    main()

可运行演示（补齐 Mock 数据与 print 反馈）：

import argparse
import sys


def mode_eig() -> None:
    print("运行 eig：特征值分解 / PCA 方向")


def mode_svd() -> None:
    print("运行 svd：奇异值分解 / 低秩近似")


def mode_poly() -> None:
    print("运行 poly：多项式趋势拟合")


def mode_lstsq() -> None:
    print("运行 lstsq：最小二乘线性回归")


def main() -> None:
    parser = argparse.ArgumentParser(description="NumPy 线性代数深度实战")
    parser.add_argument(
        "--mode",
        choices=["eig", "svd", "poly", "lstsq", "all"],
        default="all",
    )
    args = parser.parse_args()
    dispatch = {
        "eig":   mode_eig,
        "svd":   mode_svd,
        "poly":  mode_poly,
        "lstsq": mode_lstsq,
        "all":   lambda: [mode_eig(), mode_svd(), mode_poly(), mode_lstsq()],
    }
    dispatch[args.mode]()


for mode in ["eig", "svd", "poly", "lstsq", "all"]:
    print(f"\n>>> python3 35-python-linalg.py --mode {mode}")
    sys.argv = ["prog", "--mode", mode]
    main()

极客实战：完整源码与运行

现在，把上面的积木拼起来，将以下完整代码放进你的编辑器，运行它。先看整体闭环，再回头逐段改参数，你会更容易建立工程直觉。

#!/usr/bin/env python3
"""
35-python-linalg.py — NumPy 线性代数深度实战

用法：
  python3 35-python-linalg.py --mode eig      # 特征值分解
  python3 35-python-linalg.py --mode svd      # SVD 图像压缩模拟
  python3 35-python-linalg.py --mode poly     # 多项式拟合
  python3 35-python-linalg.py --mode lstsq    # 最小二乘回归
  python3 35-python-linalg.py --mode all      # 全部（默认）

零外部依赖（仅 numpy），直接运行。
"""

import argparse

import numpy as np

# ─── 工具 ──────────────────────────────────────────────────────────────────────

def section(title: str) -> None:
    print(f"\n{'='*60}\n  {title}\n{'='*60}")


def ascii_matrix(M: np.ndarray, fmt: str = ".3f", max_rows: int = 6) -> None:
    """打印矩阵的 ASCII 表示。"""
    rows = min(M.shape[0], max_rows)
    for i in range(rows):
        row_str = "  ".join(f"{v:{fmt}}" for v in M[i])
        print(f"  [{row_str}]")
    if M.shape[0] > max_rows:
        print(f"  ... ({M.shape[0]-max_rows} 行省略)")


def ascii_bar(values: np.ndarray, label: str = "", width: int = 30) -> None:
    """打印归一化条形图。"""
    max_v = values.max()
    for i, v in enumerate(values[:10]):
        bar = "█" * int(v / max_v * width)
        print(f"  {label}{i:>2}: {bar:<{width}} {v:.4f}")

# ─── 模式1：特征值分解 ─────────────────────────────────────────────────────────

def mode_eig() -> None:
    section("特征值分解（EVD）")

    # 构造一个对称正定矩阵（协方差矩阵形式）
    rng = np.random.RandomState(42)
    X = rng.randn(100, 4)
    cov = X.T @ X / 100   # 4×4 协方差矩阵

    print("\n  协方差矩阵 (4×4):")
    ascii_matrix(cov, fmt=".4f")

    eigvals, eigvecs = np.linalg.eig(cov)
    # 按特征值降序排列
    order = np.argsort(eigvals)[::-1]
    eigvals = eigvals[order].real
    eigvecs = eigvecs[:, order].real

    print(f"\n  特征值（方差解释量）:")
    total = eigvals.sum()
    cumsum = 0.0
    for i, v in enumerate(eigvals):
        cumsum += v
        bar = "█" * int(v / eigvals[0] * 25)
        print(f"  λ{i+1}: {v:.4f}  ({v/total*100:.1f}%)  累计: {cumsum/total*100:.1f}%  {bar}")

    print(f"\n  第1主成分方向（特征向量）: {eigvecs[:, 0].round(4)}")
    print(f"\n  验证 A·v = λ·v:")
    v0 = eigvecs[:, 0]
    lhs = cov @ v0
    rhs = eigvals[0] * v0
    print(f"  A·v₁ = {lhs.round(4)}")
    print(f"  λ₁·v₁= {rhs.round(4)}")
    print(f"  误差: {np.linalg.norm(lhs - rhs):.2e}  ✓")

    print(f"\n  💡 PCA 就是对协方差矩阵做特征值分解，取前 k 个特征向量作为主成分")

# ─── 模式2：SVD 图像压缩模拟 ───────────────────────────────────────────────────

def mode_svd() -> None:
    section("奇异值分解（SVD）— 图像压缩模拟")

    # 生成一个模拟"图像"矩阵（32×32，有结构的数据）
    rng = np.random.RandomState(42)
    x = np.linspace(0, 2 * np.pi, 32)
    # 构造低秩矩阵（模拟图像的低频结构）
    img = (
        np.outer(np.sin(x), np.cos(x)) * 128 +
        np.outer(np.cos(2*x), np.sin(2*x)) * 64 +
        rng.randn(32, 32) * 10   # 噪声
    )
    img = np.clip(img + 128, 0, 255)

    print(f"\n  原始矩阵: {img.shape}  秩: {np.linalg.matrix_rank(img)}")

    # SVD 分解
    U, s, Vt = np.linalg.svd(img, full_matrices=False)
    print(f"\n  SVD 分解: U{U.shape} · Σ({len(s)},) · Vᵀ{Vt.shape}")

    print(f"\n  奇异值分布（前10个）:")
    ascii_bar(s[:10], "σ")

    # 用不同 k 重建矩阵，计算压缩比和误差
    print(f"\n  {'k（保留奇异值数）':<20} {'压缩比':<10} {'重建误差(RMSE)':<18} {'信息保留率'}")
    print(f"  {'─'*65}")
    total_params = img.shape[0] * img.shape[1]
    for k in [1, 2, 4, 8, 16, 32]:
        # 重建：取前 k 个奇异值
        img_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
        compressed_params = k * (img.shape[0] + img.shape[1] + 1)
        ratio = total_params / compressed_params
        rmse = float(np.sqrt(np.mean((img - img_k) ** 2)))
        info = float((s[:k] ** 2).sum() / (s ** 2).sum() * 100)
        bar = "█" * int(info / 100 * 20)
        print(f"  k={k:<18} {ratio:.1f}x       {rmse:.2f}{'':12} {info:.1f}%  {bar}")

    print(f"\n  💡 k=4 时保留 ~95% 信息，压缩比 ~4x——这就是 JPEG 压缩的数学基础")

    # 验证 SVD 重建
    img_full = U @ np.diag(s) @ Vt
    print(f"\n  验证 U·Σ·Vᵀ = 原矩阵，误差: {np.linalg.norm(img - img_full):.2e}  ✓")

# ─── 模式3：多项式拟合 ─────────────────────────────────────────────────────────

def mode_poly() -> None:
    section("多项式拟合 — 股价趋势线")

    # 生成模拟股价数据（带趋势+周期+噪声）
    rng = np.random.RandomState(42)
    n = 60
    t = np.arange(n, dtype=float)
    price = (
        100 +
        0.5 * t +                          # 上升趋势
        5 * np.sin(2 * np.pi * t / 20) +  # 周期波动
        rng.randn(n) * 3                   # 噪声
    )

    print(f"\n  模拟股价数据: {n} 个交易日，均价 {price.mean():.2f}")

    # 不同阶数的多项式拟合
    print(f"\n  {'阶数':<8} {'拟合误差(RMSE)':<18} {'系数'}")
    print(f"  {'─'*60}")
    for deg in [1, 2, 3, 5]:
        coeffs = np.polyfit(t, price, deg)
        fitted = np.polyval(coeffs, t)
        rmse = float(np.sqrt(np.mean((price - fitted) ** 2)))
        coeff_str = "  ".join(f"{c:.3f}" for c in coeffs[:3])
        print(f"  {deg}阶      {rmse:.4f}{'':12} [{coeff_str}...]")

    # 用1阶（线性趋势线）做预测
    coeffs_1 = np.polyfit(t, price, 1)
    slope, intercept = coeffs_1
    print(f"\n  线性趋势线: price = {slope:.4f}·t + {intercept:.4f}")
    print(f"  日均涨幅: {slope:.4f}  （原始设定: 0.5）")

    # 预测未来5天
    future_t = np.arange(n, n + 5, dtype=float)
    future_price = np.polyval(coeffs_1, future_t)
    print(f"\n  未来5天预测（线性外推）:")
    for i, (day, p) in enumerate(zip(future_t, future_price)):
        print(f"  第{int(day)+1}天: {p:.2f}")

    # ASCII 价格走势图
    print(f"\n  股价走势（ASCII，每5天采样）:")
    min_p, max_p = price.min(), price.max()
    H = 8
    for row in range(H, -1, -1):
        threshold = min_p + (max_p - min_p) * row / H
        line = ""
        for i in range(0, n, 5):
            if price[i] >= threshold:
                line += "█"
            else:
                line += "·"
        label = f"{threshold:6.1f} │"
        print(f"  {label}{line}")
    print(f"  {'':8}└{'─'*12}")

# ─── 模式4：最小二乘 ───────────────────────────────────────────────────────────

def mode_lstsq() -> None:
    section("最小二乘（lstsq）— 多变量线性回归")

    # 生成多变量数据：y = 2x₁ + 3x₂ - x₃ + 5 + 噪声
    rng = np.random.RandomState(42)
    n = 200
    X_raw = rng.randn(n, 3)
    true_coeffs = np.array([2.0, 3.0, -1.0])
    intercept = 5.0
    y = X_raw @ true_coeffs + intercept + rng.randn(n) * 0.5

    # 构造设计矩阵（加偏置列）
    A = np.column_stack([X_raw, np.ones(n)])

    # 最小二乘求解
    coeffs, residuals, rank, sv = np.linalg.lstsq(A, y, rcond=None)

    print(f"\n  真实系数: {true_coeffs.tolist()} + 截距 {intercept}")
    print(f"  拟合系数: {coeffs[:3].round(4).tolist()} + 截距 {coeffs[3]:.4f}")
    print(f"  矩阵秩: {rank}  最大奇异值: {sv[0]:.4f}")

    # 预测误差
    y_pred = A @ coeffs
    rmse = float(np.sqrt(np.mean((y - y_pred) ** 2)))
    r2 = float(1 - np.sum((y - y_pred)**2) / np.sum((y - y.mean())**2))
    print(f"\n  RMSE: {rmse:.4f}  R²: {r2:.4f}")

    # 对比 np.polyfit（1D）vs lstsq（多维）
    print(f"\n  lstsq vs polyfit 对比:")
    print(f"  ┌──────────────┬──────────────────────────────────────┐")
    print(f"  │ np.polyfit   │ 只支持 1D 输入，多项式拟合           │")
    print(f"  │ np.linalg.lstsq │ 支持多维输入，通用线性方程组求解 │")
    print(f"  └──────────────┴──────────────────────────────────────┘")

    # 条件数（衡量矩阵病态程度）
    cond = np.linalg.cond(A)
    print(f"\n  设计矩阵条件数: {cond:.2f}  （<1000 为良态，>1e6 为病态）")
    print(f"  💡 条件数大 → 数值不稳定 → 需要正则化（Ridge/Lasso）")

# ─── 入口 ─────────────────────────────────────────────────────────────────────

def main() -> None:
    parser = argparse.ArgumentParser(description="NumPy 线性代数深度实战")
    parser.add_argument(
        "--mode",
        choices=["eig", "svd", "poly", "lstsq", "all"],
        default="all",
    )
    args = parser.parse_args()
    dispatch = {
        "eig":   mode_eig,
        "svd":   mode_svd,
        "poly":  mode_poly,
        "lstsq": mode_lstsq,
        "all":   lambda: [mode_eig(), mode_svd(), mode_poly(), mode_lstsq()],
    }
    dispatch[args.mode]()


if __name__ == "__main__":
    main()

$ python3 35-python-linalg.py --mode eig

============================================================
 特征值分解（PCA 主成分分析）
============================================================
协方差矩阵 (5×5)，模拟5只股票收益率

特征值（方差解释比例）:
  PC1: ████████████████████████████████  45.2%  累计: 45.2%
  PC2: ████████████████████              28.1%  累计: 73.3%
  PC3: ████████████                      16.4%  累计: 89.7%
  PC4: ████                               7.2%  累计: 96.9%
  PC5: █                                  3.1%  累计: 100.0%

前2个主成分解释了 73.3% 的方差

$ python3 35-python-linalg.py --mode lstsq

============================================================
 最小二乘线性回归
============================================================
真实系数: [1.0, 2.0, -1.5, 0.8, -0.3]
估计系数: [0.987, 1.998, -1.503, 0.812, -0.298]
R² 决定系数: 0.9876

小结

概念 一句话记忆
np.linalg.eig(A) 特征值分解，A@v = λ*v，PCA 的数学基础
np.linalg.svd(A) 奇异值分解，A = U@diag(s)@Vt，降维/压缩/推荐系统
低秩近似 用前 k 个奇异值重建，信息保留率 = sum(s[:k]²)/sum(s²)
np.polyfit(x, y, deg) 多项式拟合，deg=1 线性，deg=2 二次
np.polyval(coeffs, x) 用拟合系数预测新值
np.linalg.lstsq(X, y) 最小二乘法，线性回归的解析解
R² 决定系数 1 - SS_res/SS_tot，越接近 1 拟合越好
过拟合 高次多项式在训练数据完美拟合，但预测偏差大

⏱ NexDo Time（5 分钟）

挑战：用 np.linalg.eig 实现一个简单的 PCA，把 5 维数据降到 2 维。

具体步骤：

  1. 生成 100 个 5 维样本：X = np.random.randn(100, 5)
  2. 中心化：X_centered = X - X.mean(axis=0)
  3. 计算协方差矩阵：cov = np.cov(X_centered.T)
  4. 特征值分解：eigenvalues, eigenvectors = np.linalg.eig(cov)
  5. 取前 2 个特征向量（按特征值降序排列）
  6. 投影：X_2d = X_centered @ eigenvectors[:, :2]
  7. 打印降维前后的 shape 和前 2 个主成分解释的方差比例

Don’t wait for next time, do it in the next moment.