SciPy - 线性代数

SciPy - 线性代数 (scipy.linalg)

scipy.linalg 模块提供了一整套线性代数运算。它建立在经过优化的库之上，如 ATLAS、LAPACK（线性代数包）和 BLAS（基础线性代数程序集），确保了矩阵计算的高性能。SciPy 的线性代数例程期望以 NumPy 数组（或可转换为 NumPy 数组的对象）作为输入，并且通常将结果作为 NumPy 数组返回。

scipy.linalg 与 numpy.linalg

scipy.linalg 包含了 numpy.linalg 中所有可用的函数。然而，scipy.linalg 提供了更多更高级的函数。一个主要优势是 scipy.linalg 总是使用 BLAS/LAPACK 支持进行编译，保证了优化的性能。对于 NumPy，BLAS/LAPACK 支持在安装时是可选的，因此 scipy.linalg 函数有时会更快或更健壮。

在基于 SciPy 的项目中，通常建议使用 scipy.linalg 进行线性代数任务。

求解线性方程组

一个常见任务是求解线性方程组，表示为 $A \mathbf{x} = \mathbf{b}$，其中 $A$ 是系数矩阵（matrix of coefficients），$\mathbf{x}$ 是未知数向量（vector of unknowns），$\mathbf{b}$ 是结果向量（result vector）。

scipy.linalg.solve(a, b) 函数计算 $\mathbf{x}$。

考虑以下系统：

• $3x_1 + 2x_2 = 2$
• $1x_1 - 1x_2 = 4$
• $5x_2 + 1x_3 = -1$

这可以写成 $A \mathbf{x} = \mathbf{b}$，其中：

$$ A = \begin{bmatrix} 3 & 2 & 0 \ 1 & -1 & 0 \ 0 & 5 & 1 \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} 2 \ 4 \ -1 \end{bmatrix} $$

import numpy as np
from scipy import linalg

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

# Result vector b
b_vec = np.array([2, 4, -1])

# Solve for x
x_solution = linalg.solve(A, b_vec)

print("Solution vector x:")
print(x_solution)

# Verify the solution: A @ x should be close to b
print("\nVerification (A @ x):")
print(A @ x_solution)
# Check if A @ x_solution is close to b_vec using np.allclose
print(f"Is solution correct? {np.allclose(A @ x_solution, b_vec)}")

输出：

Solution vector x:
[ 2. -2.  9.]

Verification (A @ x):
[ 2.  4. -1.]
Is solution correct? True

计算行列式

行列式是与方阵关联的一个标量值，表示为 $\det(A)$ 或 $|A|$。它使用 scipy.linalg.det(a) 计算。

import numpy as np
from scipy import linalg

# Define a square matrix
M = np.array([
    [1, 2],
    [3, 4]
])

# Calculate the determinant
det_M = linalg.det(M)

print(f"Matrix M:\n{M}")
print(f"Determinant of M: {det_M}") # Expected: 1*4 - 2*3 = -2

输出：

Matrix M:
[[1 2]
 [3 4]]
Determinant of M: -2.0

特征值和特征向量

特征值-特征向量问题涉及为方阵 $A$ 找到标量 $\lambda$（特征值，eigenvalues）和对应的非零向量 $\mathbf{v}$（特征向量，eigenvectors），使得 $A\mathbf{v} = \lambda\mathbf{v}$。

scipy.linalg.eig(a) 计算特征值和特征向量。它返回一个元组 (eigenvalues, eigenvectors)，其中 eigenvectors[:, i] 对应于 eigenvalues[i]。

import numpy as np
from scipy import linalg

# Define a square matrix
A_eig = np.array([
    [1, 5],
    [3, 4]
])

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = linalg.eig(A_eig)

print(f"Matrix A:\n{A_eig}")
print(f"\nEigenvalues (lambda):\n{eigenvalues}")
print(f"\nEigenvectors (v columns):\n{eigenvectors}")

# Verification for the first eigenvalue/eigenvector pair
lambda1 = eigenvalues[0]
v1 = eigenvectors[:, 0]
print(f"\nVerification for first pair (A @ v1):\n{A_eig @ v1}")
print(f"Verification for first pair (lambda1 * v1):\n{lambda1 * v1}")
print(f"Are they close? {np.allclose(A_eig @ v1, lambda1 * v1)}")

输出（对于非对称矩阵，特征值/向量可能是复数）：

Matrix A:
[[1 5]
 [3 4]]

Eigenvalues (lambda):
[-1.40512484+0.j  6.40512484+0.j]

Eigenvectors (v columns):
[[-0.84797312 -0.65708331]
 [ 0.53004982 -0.75381815]]

Verification for first pair (A @ v1):
[ 1.19165751+0.j -0.74484869+0.j]
Verification for first pair (lambda1 * v1):
[ 1.19165751+0.j -0.74484869+0.j]
Are they close? True

奇异值分解 (SVD)

奇异值分解（Singular Value Decomposition，SVD）是将实数或复数矩阵 $M$ 分解为 $M = U \Sigma V^H$，其中：

• $U$ 是一个酉矩阵（unitary matrix），其列向量是左奇异向量（left singular vectors）。
• $\Sigma$ 是一个矩形对角矩阵（rectangular diagonal matrix），其对角线上是非负实数（奇异值，singular values）。
• $V^H$（$V$ 的共轭转置，conjugate transpose）是一个酉矩阵，其行向量（或 $V$ 的列向量）是右奇异向量（right singular vectors）。

SVD 是一个强大的工具，用于降维（dimensionality reduction）、求解线性系统等。scipy.linalg.svd(a) 计算 SVD。

import numpy as np
from scipy import linalg

# Define a matrix (can be non-square)
# Using a random matrix for demonstration
np.random.seed(0) # for reproducibility
A_svd = np.random.randn(3, 2) + 1j * np.random.randn(3, 2) # A 3x2 complex matrix

# Perform SVD
U, s, Vh = linalg.svd(A_svd)

# s contains the singular values (1D array)
# To reconstruct Sigma, create a diagonal matrix of appropriate shape
Sigma = np.zeros(A_svd.shape, dtype=complex) # A_svd.shape is (3,2)
Sigma[:A_svd.shape[1], :A_svd.shape[1]] = np.diag(s) # s has min(3,2)=2 elements

# Reconstruct the original matrix
A_reconstructed = U @ Sigma @ Vh

print(f"Original Matrix A (shape {A_svd.shape}):\n{A_svd}")
print(f"\nU matrix (shape {U.shape}):\n{U}")
print(f"\nSingular values (s - 1D array, shape {s.shape}):\n{s}")
print(f"\nSigma matrix (shape {Sigma.shape}):\n{Sigma}")
print(f"\nVh matrix (conjugate transpose of V, shape {Vh.shape}):\n{Vh}")
print(f"\nReconstructed Matrix (U @ Sigma @ Vh):")
print(A_reconstructed)
print(f"\nIs reconstruction close to original? {np.allclose(A_svd, A_reconstructed)}")

输出将显示矩阵 U、Vh、奇异值 s、构造的 Sigma 以及重构的 A，重构的 A 应该非常接近原始矩阵 A_svd。

scipy.linalg 还提供了许多其他函数，包括矩阵分解（LU, QR, Cholesky）、矩阵指数（matrix exponentials）和特殊矩阵构造。有关完整列表和高级用法，请务必查阅 SciPy linalg 官方文档。