数组和矩阵操作

在科学计算中，数组和矩阵操作是最基础也是最重要的技能。SciPy 在 NumPy 的基础上，通过 scipy.linalg 模块提供了更加丰富和高效的线性代数功能。本章将详细介绍这些功能的使用方法。

scipy.linalg 模块概述

scipy.linalg 模块包含了所有在 numpy.linalg 中的函数，并且提供了更多高级功能：

更完整的线性代数函数集合
更好的性能优化
更多的矩阵分解方法
特殊矩阵的处理
矩阵函数计算

python

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

# 查看 scipy.linalg 的主要功能
print("scipy.linalg 主要功能:")
functions = [attr for attr in dir(linalg) if not attr.startswith('_')]
print(f"总共 {len(functions)} 个函数和类")
print("主要函数:", functions[:10])

基本矩阵操作

1. 矩阵创建和基本属性

python

# 创建各种类型的矩阵
n = 4

# 随机矩阵
A = np.random.random((n, n))
print(f"随机矩阵 A:\n{A}")

# 对称矩阵
S = A + A.T
print(f"\n对称矩阵 S:\n{S}")

# 正定矩阵
P = A @ A.T
print(f"\n正定矩阵 P:\n{P}")

# 单位矩阵
I = np.eye(n)
print(f"\n单位矩阵 I:\n{I}")

# 零矩阵
Z = np.zeros((n, n))
print(f"\n零矩阵 Z:\n{Z}")

# 对角矩阵
D = np.diag([1, 2, 3, 4])
print(f"\n对角矩阵 D:\n{D}")

2. 矩阵基本运算

python

# 矩阵加法和减法
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print(f"矩阵 A:\n{A}")
print(f"矩阵 B:\n{B}")
print(f"A + B:\n{A + B}")
print(f"A - B:\n{A - B}")

# 矩阵乘法
print(f"\n矩阵乘法 A @ B:\n{A @ B}")
print(f"元素乘法 A * B:\n{A * B}")

# 矩阵幂
print(f"\n矩阵平方 A²:\n{linalg.matrix_power(A, 2)}")
print(f"矩阵立方 A³:\n{linalg.matrix_power(A, 3)}")

# 矩阵转置
print(f"\n矩阵转置 A.T:\n{A.T}")

# 矩阵迹（对角线元素之和）
print(f"\n矩阵的迹 tr(A): {np.trace(A)}")

3. 矩阵范数

python

# 不同类型的矩阵范数
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

print(f"矩阵 A:\n{A}")
print(f"\n各种范数:")
print(f"Frobenius 范数: {linalg.norm(A, 'fro'):.4f}")
print(f"1-范数 (列和的最大值): {linalg.norm(A, 1):.4f}")
print(f"∞-范数 (行和的最大值): {linalg.norm(A, np.inf):.4f}")
print(f"2-范数 (最大奇异值): {linalg.norm(A, 2):.4f}")
print(f"核范数 (奇异值之和): {linalg.norm(A, 'nuc'):.4f}")

# 向量范数
v = np.array([3, 4, 5])
print(f"\n向量 v: {v}")
print(f"L1 范数: {linalg.norm(v, 1):.4f}")
print(f"L2 范数: {linalg.norm(v, 2):.4f}")
print(f"L∞ 范数: {linalg.norm(v, np.inf):.4f}")

线性方程组求解

1. 基本线性方程组

python

# 求解线性方程组 Ax = b
A = np.array([[3, 2, -1], [2, -2, 4], [-1, 0.5, -1]])
b = np.array([1, -2, 0])

print(f"系数矩阵 A:\n{A}")
print(f"常数向量 b: {b}")

# 使用 scipy.linalg.solve
x = linalg.solve(A, b)
print(f"\n解向量 x: {x}")

# 验证解的正确性
residual = A @ x - b
print(f"残差 ||Ax - b||: {linalg.norm(residual):.2e}")

# 检查矩阵条件数
cond_num = linalg.cond(A)
print(f"条件数: {cond_num:.2f}")
if cond_num > 1e12:
    print("警告: 矩阵接近奇异，解可能不稳定")

2. 特殊类型的线性方程组

python

# 对称正定矩阵的求解（使用 Cholesky 分解）
n = 5
A_spd = np.random.random((n, n))
A_spd = A_spd @ A_spd.T  # 确保正定
b = np.random.random(n)

# 通用求解
start_time = time.time()
x1 = linalg.solve(A_spd, b)
general_time = time.time() - start_time

# 利用正定性质求解
start_time = time.time()
x2 = linalg.solve(A_spd, b, assume_a='pos')
cholesky_time = time.time() - start_time

print(f"通用求解耗时: {general_time:.6f} 秒")
print(f"Cholesky 求解耗时: {cholesky_time:.6f} 秒")
print(f"性能提升: {general_time / cholesky_time:.2f} 倍")
print(f"解的差异: {linalg.norm(x1 - x2):.2e}")

# 三对角矩阵求解
from scipy.linalg import solve_banded

# 创建三对角矩阵
n = 100
diag_main = 2 * np.ones(n)
diag_upper = -1 * np.ones(n-1)
diag_lower = -1 * np.ones(n-1)

# 构造带状矩阵格式
ab = np.zeros((3, n))
ab[0, 1:] = diag_upper  # 上对角线
ab[1, :] = diag_main    # 主对角线
ab[2, :-1] = diag_lower # 下对角线

b_tri = np.random.random(n)

# 使用带状矩阵求解器
x_banded = solve_banded((1, 1), ab, b_tri)
print(f"\n三对角矩阵求解完成，解向量长度: {len(x_banded)}")

3. 最小二乘问题

python

# 超定线性方程组的最小二乘解
m, n = 100, 50  # m > n，超定系统
A_over = np.random.random((m, n))
b_over = np.random.random(m)

# 最小二乘解
x_lstsq, residuals, rank, s = linalg.lstsq(A_over, b_over)

print(f"超定系统: {m} 个方程，{n} 个未知数")
print(f"矩阵秩: {rank}")
print(f"残差平方和: {residuals[0]:.4f}")
print(f"解向量范数: {linalg.norm(x_lstsq):.4f}")

# 正则化最小二乘（岭回归）
alpha = 0.1  # 正则化参数
A_reg = np.vstack([A_over, alpha * np.eye(n)])
b_reg = np.hstack([b_over, np.zeros(n)])
x_ridge = linalg.lstsq(A_reg, b_reg)[0]

print(f"\n正则化解向量范数: {linalg.norm(x_ridge):.4f}")
print(f"正则化效果: 范数减小了 {(1 - linalg.norm(x_ridge)/linalg.norm(x_lstsq))*100:.1f}%")

矩阵分解

1. LU 分解

python

# LU 分解：A = PLU
A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])

P, L, U = linalg.lu(A)

print(f"原矩阵 A:\n{A}")
print(f"\n置换矩阵 P:\n{P}")
print(f"\n下三角矩阵 L:\n{L}")
print(f"\n上三角矩阵 U:\n{U}")

# 验证分解
reconstructed = P @ L @ U
print(f"\n重构矩阵 PLU:\n{reconstructed}")
print(f"重构误差: {linalg.norm(A - reconstructed):.2e}")

# 使用 LU 分解求解线性方程组
lu, piv = linalg.lu_factor(A)
b = np.array([4, 5, 6])
x = linalg.lu_solve((lu, piv), b)
print(f"\n使用 LU 分解求解 Ax = b:")
print(f"解: {x}")
print(f"验证: Ax = {A @ x}")

2. QR 分解

python

# QR 分解：A = QR
A = np.random.random((5, 3))

Q, R = linalg.qr(A)

print(f"原矩阵 A 形状: {A.shape}")
print(f"正交矩阵 Q 形状: {Q.shape}")
print(f"上三角矩阵 R 形状: {R.shape}")

# 验证正交性
print(f"\nQ 的正交性验证:")
print(f"Q.T @ Q:\n{Q.T @ Q}")
print(f"是否接近单位矩阵: {np.allclose(Q.T @ Q, np.eye(Q.shape[1]))}")

# 验证分解
reconstructed = Q @ R
print(f"\n重构误差: {linalg.norm(A - reconstructed):.2e}")

# QR 分解的应用：最小二乘问题
m, n = 10, 5
A_qr = np.random.random((m, n))
b_qr = np.random.random(m)

Q, R = linalg.qr(A_qr)
# 解 Rx = Q.T @ b
x_qr = linalg.solve_triangular(R, Q.T @ b_qr)

print(f"\n使用 QR 分解求解最小二乘问题:")
print(f"解向量: {x_qr}")
print(f"残差: {linalg.norm(A_qr @ x_qr - b_qr):.4f}")

3. 奇异值分解 (SVD)

python

# SVD 分解：A = UΣV^T
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])

U, s, Vt = linalg.svd(A)

print(f"原矩阵 A 形状: {A.shape}")
print(f"U 矩阵形状: {U.shape}")
print(f"奇异值向量 s: {s}")
print(f"V^T 矩阵形状: {Vt.shape}")

# 重构矩阵
S = np.zeros(A.shape)
S[:len(s), :len(s)] = np.diag(s)
reconstructed = U @ S @ Vt

print(f"\n重构误差: {linalg.norm(A - reconstructed):.2e}")

# SVD 的应用：矩阵的秩和条件数
rank = np.sum(s > 1e-10)
cond_num = s[0] / s[-1] if s[-1] > 1e-15 else np.inf

print(f"\n矩阵的数值秩: {rank}")
print(f"条件数: {cond_num:.2f}")

# 低秩近似
for k in [1, 2, 3]:
    A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    error = linalg.norm(A - A_k, 'fro')
    print(f"秩-{k} 近似误差: {error:.4f}")

4. Cholesky 分解

python

# Cholesky 分解：A = L @ L.T（仅适用于正定矩阵）
n = 4
A_random = np.random.random((n, n))
A_pd = A_random @ A_random.T  # 构造正定矩阵

try:
    L = linalg.cholesky(A_pd, lower=True)
    print(f"正定矩阵 A:\n{A_pd}")
    print(f"\nCholesky 因子 L:\n{L}")
    
    # 验证分解
    reconstructed = L @ L.T
    print(f"\n重构误差: {linalg.norm(A_pd - reconstructed):.2e}")
    
    # 使用 Cholesky 分解求解线性方程组
    b = np.random.random(n)
    
    # 两步求解：Ly = b, L^T x = y
    y = linalg.solve_triangular(L, b, lower=True)
    x = linalg.solve_triangular(L.T, y, lower=False)
    
    print(f"\n使用 Cholesky 分解求解:")
    print(f"解: {x}")
    print(f"验证: {linalg.norm(A_pd @ x - b):.2e}")
    
except linalg.LinAlgError:
    print("矩阵不是正定的，无法进行 Cholesky 分解")

特征值和特征向量

1. 一般矩阵的特征值分解

python

# 计算特征值和特征向量
A = np.array([[1, 2], [2, 1]])

eigenvalues, eigenvectors = linalg.eig(A)

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

# 验证特征值方程 Av = λv
for i in range(len(eigenvalues)):
    lambda_i = eigenvalues[i]
    v_i = eigenvectors[:, i]
    
    Av = A @ v_i
    lambda_v = lambda_i * v_i
    
    print(f"\n特征值 {i+1}: λ = {lambda_i:.4f}")
    print(f"Av = {Av}")
    print(f"λv = {lambda_v}")
    print(f"误差: {linalg.norm(Av - lambda_v):.2e}")

# 矩阵对角化
D = np.diag(eigenvalues)
P = eigenvectors
A_reconstructed = P @ D @ linalg.inv(P)

print(f"\n对角化验证:")
print(f"重构矩阵:\n{A_reconstructed}")
print(f"重构误差: {linalg.norm(A - A_reconstructed):.2e}")

2. 对称矩阵的特征值分解

python

# 对称矩阵有实特征值和正交特征向量
A_sym = np.array([[4, 2, 1], [2, 3, 0], [1, 0, 2]])

# 使用专门的对称矩阵算法
eigenvalues_sym, eigenvectors_sym = linalg.eigh(A_sym)

print(f"对称矩阵 A:\n{A_sym}")
print(f"\n特征值（按升序排列）: {eigenvalues_sym}")
print(f"特征向量矩阵:\n{eigenvectors_sym}")

# 验证正交性
print(f"\n特征向量的正交性:")
print(f"V^T @ V:\n{eigenvectors_sym.T @ eigenvectors_sym}")
print(f"是否正交: {np.allclose(eigenvectors_sym.T @ eigenvectors_sym, np.eye(3))}")

# 谱分解
spectral_decomp = eigenvectors_sym @ np.diag(eigenvalues_sym) @ eigenvectors_sym.T
print(f"\n谱分解重构误差: {linalg.norm(A_sym - spectral_decomp):.2e}")

3. 广义特征值问题

python

# 广义特征值问题：Av = λBv
A = np.array([[1, 2], [3, 4]])
B = np.array([[2, 1], [1, 2]])

eigenvalues_gen, eigenvectors_gen = linalg.eig(A, B)

print(f"矩阵 A:\n{A}")
print(f"矩阵 B:\n{B}")
print(f"\n广义特征值: {eigenvalues_gen}")
print(f"广义特征向量:\n{eigenvectors_gen}")

# 验证广义特征值方程
for i in range(len(eigenvalues_gen)):
    lambda_i = eigenvalues_gen[i]
    v_i = eigenvectors_gen[:, i]
    
    Av = A @ v_i
    lambda_Bv = lambda_i * (B @ v_i)
    
    print(f"\n广义特征值 {i+1}: λ = {lambda_i:.4f}")
    print(f"误差: {linalg.norm(Av - lambda_Bv):.2e}")

矩阵函数

1. 矩阵指数

python

# 矩阵指数 exp(A)
A = np.array([[1, 2], [0, 1]])

# 计算矩阵指数
exp_A = linalg.expm(A)

print(f"矩阵 A:\n{A}")
print(f"\n矩阵指数 exp(A):\n{exp_A}")

# 验证性质：exp(A + B) = exp(A) * exp(B)（当 AB = BA 时）
B = np.array([[0, 1], [0, 0]])
if np.allclose(A @ B, B @ A):
    exp_A_plus_B = linalg.expm(A + B)
    exp_A_exp_B = linalg.expm(A) @ linalg.expm(B)
    
    print(f"\nexp(A + B):\n{exp_A_plus_B}")
    print(f"exp(A) * exp(B):\n{exp_A_exp_B}")
    print(f"差异: {linalg.norm(exp_A_plus_B - exp_A_exp_B):.2e}")

# 矩阵指数的应用：求解线性微分方程组
# dx/dt = Ax, x(0) = x0
# 解为 x(t) = exp(At) * x0
t = 1.0
x0 = np.array([1, 0])
exp_At = linalg.expm(A * t)
x_t = exp_At @ x0

print(f"\n线性微分方程组的解:")
print(f"t = {t} 时的解: {x_t}")

2. 矩阵对数和幂函数

python

# 矩阵对数
A_pos = np.array([[2, 1], [0, 3]])  # 确保特征值为正

log_A = linalg.logm(A_pos)
print(f"矩阵 A:\n{A_pos}")
print(f"\n矩阵对数 log(A):\n{log_A}")

# 验证：exp(log(A)) = A
exp_log_A = linalg.expm(log_A)
print(f"\nexp(log(A)):\n{exp_log_A}")
print(f"重构误差: {linalg.norm(A_pos - exp_log_A):.2e}")

# 矩阵幂函数
for p in [0.5, 2, -1]:
    A_power = linalg.fractional_matrix_power(A_pos, p)
    print(f"\nA^{p}:\n{A_power}")
    
    if p == -1:
        # 验证逆矩阵
        A_inv_check = linalg.inv(A_pos)
        print(f"使用 inv() 计算的逆矩阵:\n{A_inv_check}")
        print(f"差异: {linalg.norm(A_power - A_inv_check):.2e}")

3. 矩阵平方根

python

# 矩阵平方根
A_sqrt = linalg.sqrtm(A_pos)
print(f"矩阵平方根 sqrt(A):\n{A_sqrt}")

# 验证：sqrt(A) * sqrt(A) = A
A_reconstructed = A_sqrt @ A_sqrt
print(f"\nsqrt(A) * sqrt(A):\n{A_reconstructed}")
print(f"重构误差: {linalg.norm(A_pos - A_reconstructed):.2e}")

# 对于对称正定矩阵，可以使用特征值分解计算平方根
A_spd = np.array([[4, 2], [2, 3]])
evals, evecs = linalg.eigh(A_spd)
A_sqrt_eig = evecs @ np.diag(np.sqrt(evals)) @ evecs.T

print(f"\n使用特征值分解计算的平方根:\n{A_sqrt_eig}")
print(f"验证: {linalg.norm(A_sqrt_eig @ A_sqrt_eig - A_spd):.2e}")

稀疏矩阵操作

1. 稀疏矩阵的创建和基本操作

python

from scipy import sparse

# 创建稀疏矩阵
n = 1000
# 创建三对角矩阵
diagonals = [np.ones(n-1), 2*np.ones(n), np.ones(n-1)]
offsets = [-1, 0, 1]
A_sparse = sparse.diags(diagonals, offsets, format='csr')

print(f"稀疏矩阵形状: {A_sparse.shape}")
print(f"非零元素数量: {A_sparse.nnz}")
print(f"稀疏度: {1 - A_sparse.nnz / (n * n):.4f}")

# 稀疏矩阵的线性方程组求解
from scipy.sparse.linalg import spsolve

b_sparse = np.random.random(n)
x_sparse = spsolve(A_sparse, b_sparse)

print(f"\n稀疏线性方程组求解完成")
print(f"解向量范数: {linalg.norm(x_sparse):.4f}")
print(f"残差: {linalg.norm(A_sparse @ x_sparse - b_sparse):.2e}")

# 稀疏矩阵的特征值计算
from scipy.sparse.linalg import eigsh

# 计算最大的几个特征值
k = 5  # 计算前 5 个特征值
eigenvalues_sparse, eigenvectors_sparse = eigsh(A_sparse, k=k, which='LA')

print(f"\n最大的 {k} 个特征值: {eigenvalues_sparse}")
print(f"对应特征向量形状: {eigenvectors_sparse.shape}")

实际应用示例

示例 1：图像压缩（SVD 应用）

python

# 使用 SVD 进行图像压缩
import matplotlib.pyplot as plt

# 创建一个简单的"图像"（二维数组）
np.random.seed(42)
original_image = np.random.random((50, 50))
# 添加一些结构
for i in range(10, 40):
    for j in range(10, 40):
        original_image[i, j] = 0.8

# SVD 分解
U, s, Vt = linalg.svd(original_image)

# 不同压缩比的重构
compression_ratios = [5, 10, 20, 30]

plt.figure(figsize=(15, 3))

plt.subplot(1, len(compression_ratios) + 1, 1)
plt.imshow(original_image, cmap='gray')
plt.title('原始图像')
plt.axis('off')

for i, k in enumerate(compression_ratios):
    # 使用前 k 个奇异值重构
    compressed_image = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    
    # 计算压缩比和误差
    compression_ratio = k * (U.shape[0] + Vt.shape[1]) / (U.shape[0] * Vt.shape[1])
    mse = np.mean((original_image - compressed_image)**2)
    
    plt.subplot(1, len(compression_ratios) + 1, i + 2)
    plt.imshow(compressed_image, cmap='gray')
    plt.title(f'k={k}\n压缩比: {compression_ratio:.2f}\nMSE: {mse:.4f}')
    plt.axis('off')

plt.tight_layout()
plt.show()

# 奇异值的衰减
plt.figure(figsize=(10, 6))
plt.semilogy(s, 'bo-')
plt.xlabel('奇异值索引')
plt.ylabel('奇异值大小')
plt.title('奇异值衰减')
plt.grid(True, alpha=0.3)
plt.show()

示例 2：主成分分析 (PCA)

python

# 使用 SVD 实现 PCA
from sklearn.datasets import make_blobs

# 生成二维数据
n_samples = 300
X, _ = make_blobs(n_samples=n_samples, centers=2, n_features=2, 
                 random_state=42, cluster_std=1.5)

# 数据中心化
X_centered = X - np.mean(X, axis=0)

# 使用 SVD 进行 PCA
U, s, Vt = linalg.svd(X_centered, full_matrices=False)

# 主成分（特征向量）
components = Vt.T
print(f"主成分方向:")
print(f"第一主成分: {components[:, 0]}")
print(f"第二主成分: {components[:, 1]}")

# 解释方差比
explained_variance_ratio = (s**2) / np.sum(s**2)
print(f"\n解释方差比:")
print(f"第一主成分: {explained_variance_ratio[0]:.3f}")
print(f"第二主成分: {explained_variance_ratio[1]:.3f}")

# 数据投影到主成分空间
X_pca = X_centered @ components

# 可视化
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(X[:, 0], X[:, 1], alpha=0.6)
plt.title('原始数据')
plt.xlabel('特征 1')
plt.ylabel('特征 2')
plt.axis('equal')
plt.grid(True, alpha=0.3)

# 绘制主成分方向
mean_point = np.mean(X, axis=0)
for i in range(2):
    direction = components[:, i] * s[i] * 0.3
    plt.arrow(mean_point[0], mean_point[1], direction[0], direction[1],
              head_width=0.3, head_length=0.2, fc=f'C{i}', ec=f'C{i}',
              label=f'PC{i+1}')
plt.legend()

plt.subplot(1, 2, 2)
plt.scatter(X_pca[:, 0], X_pca[:, 1], alpha=0.6)
plt.title('PCA 变换后的数据')
plt.xlabel('第一主成分')
plt.ylabel('第二主成分')
plt.axis('equal')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

示例 3：线性回归的矩阵形式

python

# 使用矩阵运算实现线性回归
from sklearn.datasets import make_regression

# 生成回归数据
X, y = make_regression(n_samples=100, n_features=3, noise=10, random_state=42)

# 添加偏置项
X_with_bias = np.column_stack([np.ones(X.shape[0]), X])

print(f"设计矩阵形状: {X_with_bias.shape}")
print(f"目标向量形状: {y.shape}")

# 方法 1：正规方程 θ = (X^T X)^(-1) X^T y
XTX = X_with_bias.T @ X_with_bias
XTy = X_with_bias.T @ y
theta_normal = linalg.solve(XTX, XTy)

print(f"\n正规方程解: {theta_normal}")

# 方法 2：使用 QR 分解
Q, R = linalg.qr(X_with_bias)
theta_qr = linalg.solve_triangular(R, Q.T @ y)

print(f"QR 分解解: {theta_qr}")

# 方法 3：使用 SVD（最稳定）
U, s, Vt = linalg.svd(X_with_bias, full_matrices=False)
theta_svd = Vt.T @ (np.diag(1/s) @ (U.T @ y))

print(f"SVD 解: {theta_svd}")

# 比较解的差异
print(f"\n解的差异:")
print(f"正规方程 vs QR: {linalg.norm(theta_normal - theta_qr):.2e}")
print(f"正规方程 vs SVD: {linalg.norm(theta_normal - theta_svd):.2e}")
print(f"QR vs SVD: {linalg.norm(theta_qr - theta_svd):.2e}")

# 预测和评估
y_pred = X_with_bias @ theta_svd
mse = np.mean((y - y_pred)**2)
r2 = 1 - np.sum((y - y_pred)**2) / np.sum((y - np.mean(y))**2)

print(f"\n模型性能:")
print(f"均方误差: {mse:.2f}")
print(f"R² 分数: {r2:.3f}")

# 可视化（选择一个特征进行可视化）
plt.figure(figsize=(10, 4))

for i in range(3):
    plt.subplot(1, 3, i+1)
    plt.scatter(X[:, i], y, alpha=0.6)
    
    # 绘制回归线
    x_range = np.linspace(X[:, i].min(), X[:, i].max(), 100)
    X_range = np.column_stack([np.ones(100), 
                              np.zeros((100, 3))])
    X_range[:, i+1] = x_range
    y_range = X_range @ theta_svd
    
    plt.plot(x_range, y_range, 'r-', linewidth=2)
    plt.xlabel(f'特征 {i+1}')
    plt.ylabel('目标值')
    plt.title(f'特征 {i+1} vs 目标值')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

性能优化技巧

1. 选择合适的算法

python

import time

# 比较不同算法的性能
n = 1000

# 创建不同类型的矩阵
A_general = np.random.random((n, n))
A_symmetric = A_general + A_general.T
A_positive_definite = A_general @ A_general.T
b = np.random.random(n)

# 通用求解
start = time.time()
x1 = linalg.solve(A_general, b)
general_time = time.time() - start

# 对称矩阵求解
start = time.time()
x2 = linalg.solve(A_symmetric, b, assume_a='sym')
symmetric_time = time.time() - start

# 正定矩阵求解
start = time.time()
x3 = linalg.solve(A_positive_definite, b, assume_a='pos')
positive_definite_time = time.time() - start

print(f"性能比较 (n={n}):")
print(f"通用算法: {general_time:.4f} 秒")
print(f"对称矩阵算法: {symmetric_time:.4f} 秒")
print(f"正定矩阵算法: {positive_definite_time:.4f} 秒")
print(f"\n性能提升:")
print(f"对称 vs 通用: {general_time / symmetric_time:.2f}x")
print(f"正定 vs 通用: {general_time / positive_definite_time:.2f}x")

2. 内存优化

python

# 就地操作和内存管理
n = 2000

# 避免不必要的拷贝
A = np.random.random((n, n))
B = np.random.random((n, n))

# 低效：创建新数组
start = time.time()
C1 = A + B
C1 = C1 @ A
copy_time = time.time() - start

# 高效：就地操作
start = time.time()
C2 = A.copy()  # 只拷贝一次
C2 += B        # 就地加法
C2 @= A        # 就地矩阵乘法（Python 3.5+）
inplace_time = time.time() - start

print(f"拷贝操作耗时: {copy_time:.4f} 秒")
print(f"就地操作耗时: {inplace_time:.4f} 秒")
print(f"性能提升: {copy_time / inplace_time:.2f}x")

# 使用合适的数据类型
A_float64 = np.random.random((n, n)).astype(np.float64)
A_float32 = A_float64.astype(np.float32)

print(f"\n内存使用:")
print(f"float64: {A_float64.nbytes / 1024**2:.1f} MB")
print(f"float32: {A_float32.nbytes / 1024**2:.1f} MB")
print(f"内存节省: {(1 - A_float32.nbytes / A_float64.nbytes) * 100:.1f}%")

小结

在本章中，我们深入学习了 SciPy 的数组和矩阵操作：

基本操作：矩阵创建、运算、范数计算
线性方程组：各种求解方法和特殊情况处理
矩阵分解：LU、QR、SVD、Cholesky 分解及其应用
特征值问题：一般和对称矩阵的特征值分解
矩阵函数：指数、对数、幂函数等
稀疏矩阵：大规模稀疏问题的高效处理
实际应用：图像压缩、PCA、线性回归等
性能优化：算法选择和内存管理

这些技能是科学计算的基础，为后续学习更高级的 SciPy 功能奠定了坚实基础。接下来，我们将在统计分析中学习 SciPy 强大的统计功能。

练习题

矩阵分解应用：实现一个函数，使用不同的矩阵分解方法求解线性方程组，并比较它们的性能
特征值应用：编写代码分析一个网络的连通性（使用图的邻接矩阵的特征值）
SVD 应用：实现一个简单的推荐系统，使用 SVD 进行矩阵分解
性能优化：比较稠密矩阵和稀疏矩阵操作的性能差异

python

# 练习题参考答案

# 1. 矩阵分解应用
def solve_comparison(A, b):
    """
    使用不同方法求解线性方程组并比较性能
    """
    methods = {}
    
    # 直接求解
    start = time.time()
    x_direct = linalg.solve(A, b)
    methods['直接求解'] = (time.time() - start, x_direct)
    
    # LU 分解
    start = time.time()
    lu, piv = linalg.lu_factor(A)
    x_lu = linalg.lu_solve((lu, piv), b)
    methods['LU分解'] = (time.time() - start, x_lu)
    
    # QR 分解
    start = time.time()
    Q, R = linalg.qr(A)
    x_qr = linalg.solve_triangular(R, Q.T @ b)
    methods['QR分解'] = (time.time() - start, x_qr)
    
    return methods

# 测试
n = 500
A_test = np.random.random((n, n))
b_test = np.random.random(n)
results = solve_comparison(A_test, b_test)

for method, (time_taken, solution) in results.items():
    error = linalg.norm(A_test @ solution - b_test)
    print(f"{method}: {time_taken:.4f}s, 误差: {error:.2e}")

数组和矩阵操作 ​

scipy.linalg 模块概述 ​

基本矩阵操作 ​

1. 矩阵创建和基本属性 ​

2. 矩阵基本运算 ​

3. 矩阵范数 ​

线性方程组求解 ​

1. 基本线性方程组 ​

2. 特殊类型的线性方程组 ​

3. 最小二乘问题 ​

矩阵分解 ​

1. LU 分解 ​

2. QR 分解 ​

3. 奇异值分解 (SVD) ​

4. Cholesky 分解 ​

特征值和特征向量 ​

1. 一般矩阵的特征值分解 ​

2. 对称矩阵的特征值分解 ​

3. 广义特征值问题 ​

矩阵函数 ​

1. 矩阵指数 ​

2. 矩阵对数和幂函数 ​

3. 矩阵平方根 ​

稀疏矩阵操作 ​

1. 稀疏矩阵的创建和基本操作 ​

实际应用示例 ​

示例 1：图像压缩（SVD 应用） ​

示例 2：主成分分析 (PCA) ​

示例 3：线性回归的矩阵形式 ​

性能优化技巧 ​

1. 选择合适的算法 ​

2. 内存优化 ​

小结 ​

练习题 ​

数组和矩阵操作

scipy.linalg 模块概述

基本矩阵操作

1. 矩阵创建和基本属性

2. 矩阵基本运算

3. 矩阵范数

线性方程组求解

1. 基本线性方程组

2. 特殊类型的线性方程组

3. 最小二乘问题

矩阵分解

1. LU 分解

2. QR 分解

3. 奇异值分解 (SVD)

4. Cholesky 分解

特征值和特征向量

1. 一般矩阵的特征值分解

2. 对称矩阵的特征值分解

3. 广义特征值问题

矩阵函数

1. 矩阵指数

2. 矩阵对数和幂函数

3. 矩阵平方根

稀疏矩阵操作

1. 稀疏矩阵的创建和基本操作

实际应用示例

示例 1：图像压缩（SVD 应用）

示例 2：主成分分析 (PCA)

示例 3：线性回归的矩阵形式

性能优化技巧

1. 选择合适的算法

2. 内存优化

小结

练习题