10/7/25About 17 min
第15章：数值线性代数基础
学习目标
- 理解数值计算中的误差问题
- 掌握矩阵分解的数值方法
- 学习迭代法求解线性方程组
- 理解条件数和稳定性
- 掌握特征值的数值计算方法
数值计算中的误差与稳定性
浮点数表示与舍入误差
计算机中的实数表示存在精度限制，导致舍入误差的累积。
条件数
矩阵的条件数衡量了线性系统对输入扰动的敏感性：
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import norm, solve, lu, qr, svd
from scipy.sparse import diags
import time
import warnings
warnings.filterwarnings('ignore')
def numerical_stability_demo():
"""数值稳定性和误差分析演示"""
print("数值线性代数：稳定性和误差分析")
print("=" * 60)
# 1. 机器精度演示
print("1. 机器精度分析")
print("-" * 30)
machine_epsilon = np.finfo(float).eps
print(f"机器精度 ε: {machine_epsilon:.2e}")
# 演示浮点数精度问题
a = 0.1 + 0.1 + 0.1
b = 0.3
print(f"0.1 + 0.1 + 0.1 = {a:.17f}")
print(f"0.3 = {b:.17f}")
print(f"相等? {a == b}")
print(f"差值: {abs(a - b):.2e}")
# 2. 病态矩阵示例
print(f"\n2. 病态矩阵分析")
print("-" * 30)
# Hilbert矩阵（著名的病态矩阵）
def hilbert_matrix(n):
"""生成n阶Hilbert矩阵"""
H = np.zeros((n, n))
for i in range(n):
for j in range(n):
H[i, j] = 1.0 / (i + j + 1)
return H
# 不同大小的Hilbert矩阵
sizes = [3, 5, 8, 10]
condition_numbers = []
print("Hilbert矩阵条件数分析:")
for n in sizes:
H = hilbert_matrix(n)
cond_num = np.linalg.cond(H)
condition_numbers.append(cond_num)
print(f"H_{n}: 条件数 = {cond_num:.2e}")
# 3. 条件数对求解精度的影响
print(f"\n3. 条件数对求解精度的影响")
print("-" * 40)
n = 6
H = hilbert_matrix(n)
x_exact = np.ones(n) # 精确解
b_exact = H @ x_exact # 精确右端项
# 添加小扰动
perturbation_levels = [1e-12, 1e-10, 1e-8, 1e-6]
results = []
for eps in perturbation_levels:
# 扰动右端项
b_perturbed = b_exact + eps * np.random.randn(n)
# 求解扰动系统
try:
x_computed = solve(H, b_perturbed)
relative_error = norm(x_computed - x_exact) / norm(x_exact)
results.append((eps, relative_error))
print(f"扰动水平 {eps:.0e}: 相对误差 = {relative_error:.2e}")
except:
print(f"扰动水平 {eps:.0e}: 求解失败")
# 可视化
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# 条件数增长
ax1 = axes[0, 0]
ax1.semilogy(sizes, condition_numbers, 'o-', linewidth=2, markersize=8)
ax1.set_xlabel('矩阵大小')
ax1.set_ylabel('条件数')
ax1.set_title('Hilbert矩阵条件数增长')
ax1.grid(True, alpha=0.3)
# Hilbert矩阵可视化
ax2 = axes[0, 1]
H_display = hilbert_matrix(8)
im = ax2.imshow(H_display, cmap='Blues')
ax2.set_title('8×8 Hilbert矩阵')
plt.colorbar(im, ax=ax2)
# 误差放大效应
ax3 = axes[0, 2]
if results:
perturbations, errors = zip(*results)
ax3.loglog(perturbations, errors, 'ro-', linewidth=2, markersize=8, label='实际误差')
# 理论误差界
cond_H = np.linalg.cond(H)
theoretical_errors = [cond_H * eps / norm(b_exact) for eps in perturbations]
ax3.loglog(perturbations, theoretical_errors, 'b--', linewidth=2, label='理论误差界')
ax3.set_xlabel('输入扰动水平')
ax3.set_ylabel('输出相对误差')
ax3.set_title('误差放大效应')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 不同矩阵类型的条件数比较
ax4 = axes[1, 0]
matrix_types = {
'单位矩阵': np.eye(5),
'对角优势': np.diag([10, 5, 3, 2, 1]) + 0.1*np.random.randn(5, 5),
'Hilbert-5': hilbert_matrix(5),
'随机矩阵': np.random.randn(5, 5)
}
names = []
conds = []
for name, A in matrix_types.items():
if name != 'Hilbert-5':
A = A @ A.T # 确保正定
try:
cond_A = np.linalg.cond(A)
names.append(name)
conds.append(cond_A)
except:
pass
ax4.bar(range(len(names)), conds, alpha=0.7)
ax4.set_yscale('log')
ax4.set_xlabel('矩阵类型')
ax4.set_ylabel('条件数（对数尺度）')
ax4.set_title('不同矩阵类型的条件数')
ax4.set_xticks(range(len(names)))
ax4.set_xticklabels(names, rotation=45)
ax4.grid(True, alpha=0.3)
# 舍入误差累积演示
ax5 = axes[1, 1]
def sum_series(n, method='forward'):
"""计算级数和的不同方法"""
if method == 'forward':
# 正向累加
result = 0.0
for i in range(1, n+1):
result += 1.0 / i**2
else: # backward
# 反向累加
result = 0.0
for i in range(n, 0, -1):
result += 1.0 / i**2
return result
n_values = np.logspace(2, 6, 20).astype(int)
forward_results = []
backward_results = []
for n in n_values:
forward_results.append(sum_series(n, 'forward'))
backward_results.append(sum_series(n, 'backward'))
ax5.semilogx(n_values, forward_results, 'r-', label='正向累加')
ax5.semilogx(n_values, backward_results, 'b-', label='反向累加')
ax5.axhline(y=np.pi**2/6, color='black', linestyle='--', label='理论值 π²/6')
ax5.set_xlabel('项数')
ax5.set_ylabel('级数和')
ax5.set_title('舍入误差累积：Σ(1/k²)')
ax5.legend()
ax5.grid(True, alpha=0.3)
# 数值稳定性比较
ax6 = axes[1, 2]
# 比较不同算法的数值稳定性
def unstable_algorithm(x):
"""数值不稳定算法示例"""
return (x**2 - 1) / (x - 1)
def stable_algorithm(x):
"""数值稳定算法"""
return x + 1
x_values = np.linspace(0.99, 1.01, 100)
x_values = x_values[x_values != 1.0] # 避免除零
unstable_results = []
stable_results = []
for x in x_values:
try:
unstable_results.append(unstable_algorithm(x))
stable_results.append(stable_algorithm(x))
except:
unstable_results.append(np.nan)
stable_results.append(np.nan)
ax6.plot(x_values, unstable_results, 'r-', label='不稳定: (x²-1)/(x-1)', linewidth=2)
ax6.plot(x_values, stable_results, 'b-', label='稳定: x+1', linewidth=2)
ax6.set_xlabel('x')
ax6.set_ylabel('函数值')
ax6.set_title('数值稳定性比较')
ax6.legend()
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return condition_numbers
def matrix_factorizations_numerical():
"""矩阵分解的数值实现"""
print(f"\n矩阵分解的数值方法")
print("=" * 40)
# 生成测试矩阵
np.random.seed(42)
n = 6
A_random = np.random.randn(n, n)
A_spd = A_random @ A_random.T # 对称正定矩阵
A_general = np.random.randn(n, n) + 0.1 * np.eye(n) # 一般矩阵
print("测试矩阵信息:")
print(f"随机对称正定矩阵条件数: {np.linalg.cond(A_spd):.2e}")
print(f"一般矩阵条件数: {np.linalg.cond(A_general):.2e}")
# 1. LU分解
print(f"\n1. LU分解")
print("-" * 20)
P, L, U = lu(A_general)
print("LU分解验证:")
print(f"||PA - LU||: {norm(P @ A_general - L @ U):.2e}")
# 检查L和U的性质
print(f"L是下三角: {np.allclose(L, np.tril(L))}")
print(f"U是上三角: {np.allclose(U, np.triu(U))}")
print(f"P是置换矩阵: {np.allclose(P @ P.T, np.eye(n))}")
# 2. QR分解
print(f"\n2. QR分解")
print("-" * 20)
Q, R = qr(A_general)
print("QR分解验证:")
print(f"||A - QR||: {norm(A_general - Q @ R):.2e}")
print(f"Q是正交矩阵: {np.allclose(Q @ Q.T, np.eye(n))}")
print(f"R是上三角: {np.allclose(R, np.triu(R))}")
# 3. SVD分解
print(f"\n3. 奇异值分解")
print("-" * 20)
U_svd, s, Vt = svd(A_general)
print("SVD分解验证:")
print(f"||A - UΣVᵀ||: {norm(A_general - U_svd @ np.diag(s) @ Vt):.2e}")
print(f"奇异值: {s}")
print(f"条件数 (σ_max/σ_min): {s[0]/s[-1]:.2e}")
# 4. Cholesky分解（仅对正定矩阵）
print(f"\n4. Cholesky分解")
print("-" * 20)
try:
from scipy.linalg import cholesky
L_chol = cholesky(A_spd, lower=True)
print("Cholesky分解验证:")
print(f"||A - LLᵀ||: {norm(A_spd - L_chol @ L_chol.T):.2e}")
print(f"L是下三角: {np.allclose(L_chol, np.tril(L_chol))}")
# Cholesky vs LU计算复杂度比较
print(f"Cholesky分解节省约50%计算量（对正定矩阵）")
except Exception as e:
print(f"Cholesky分解失败: {e}")
# 可视化矩阵分解
fig, axes = plt.subplots(3, 4, figsize=(16, 12))
# 原始矩阵
ax = axes[0, 0]
im = ax.imshow(A_general, cmap='RdBu')
ax.set_title('原矩阵 A')
plt.colorbar(im, ax=ax)
# LU分解
ax = axes[0, 1]
im = ax.imshow(L, cmap='Blues')
ax.set_title('L (下三角)')
plt.colorbar(im, ax=ax)
ax = axes[0, 2]
im = ax.imshow(U, cmap='Reds')
ax.set_title('U (上三角)')
plt.colorbar(im, ax=ax)
ax = axes[0, 3]
im = ax.imshow(P, cmap='Greys')
ax.set_title('P (置换)')
plt.colorbar(im, ax=ax)
# QR分解
ax = axes[1, 0]
im = ax.imshow(Q, cmap='RdBu')
ax.set_title('Q (正交)')
plt.colorbar(im, ax=ax)
ax = axes[1, 1]
im = ax.imshow(R, cmap='Reds')
ax.set_title('R (上三角)')
plt.colorbar(im, ax=ax)
# SVD分解
ax = axes[1, 2]
im = ax.imshow(U_svd, cmap='Blues')
ax.set_title('U (左奇异向量)')
plt.colorbar(im, ax=ax)
ax = axes[1, 3]
ax.bar(range(len(s)), s)
ax.set_title('奇异值')
ax.set_xlabel('索引')
ax.set_ylabel('值')
ax.grid(True, alpha=0.3)
# 分解精度比较
ax = axes[2, 0]
methods = ['LU', 'QR', 'SVD']
errors = [
norm(P @ A_general - L @ U),
norm(A_general - Q @ R),
norm(A_general - U_svd @ np.diag(s) @ Vt)
]
bars = ax.bar(methods, errors)
ax.set_yscale('log')
ax.set_ylabel('重构误差')
ax.set_title('分解精度比较')
ax.grid(True, alpha=0.3)
# 计算时间比较
ax = axes[2, 1]
def time_factorization(A, method, n_trials=10):
"""测试分解方法的计算时间"""
times = []
for _ in range(n_trials):
start_time = time.time()
if method == 'LU':
lu(A)
elif method == 'QR':
qr(A)
elif method == 'SVD':
svd(A)
end_time = time.time()
times.append(end_time - start_time)
return np.mean(times)
time_results = {}
for method in ['LU', 'QR', 'SVD']:
time_results[method] = time_factorization(A_general, method)
methods = list(time_results.keys())
times = list(time_results.values())
ax.bar(methods, times, alpha=0.7)
ax.set_ylabel('平均时间 (秒)')
ax.set_title('计算时间比较')
ax.grid(True, alpha=0.3)
# 条件数与分解稳定性
ax = axes[2, 2]
# 生成不同条件数的矩阵
condition_numbers = [1e2, 1e4, 1e6, 1e8]
lu_errors = []
qr_errors = []
for cond in condition_numbers:
# 构造指定条件数的矩阵
U_temp, _, Vt_temp = svd(np.random.randn(5, 5))
singular_vals = np.logspace(0, -np.log10(cond), 5)
A_temp = U_temp @ np.diag(singular_vals) @ Vt_temp
# 测试分解精度
P_temp, L_temp, U_temp = lu(A_temp)
Q_temp, R_temp = qr(A_temp)
lu_error = norm(P_temp @ A_temp - L_temp @ U_temp) / norm(A_temp)
qr_error = norm(A_temp - Q_temp @ R_temp) / norm(A_temp)
lu_errors.append(lu_error)
qr_errors.append(qr_error)
ax.loglog(condition_numbers, lu_errors, 'o-', label='LU分解', linewidth=2)
ax.loglog(condition_numbers, qr_errors, 's-', label='QR分解', linewidth=2)
ax.set_xlabel('条件数')
ax.set_ylabel('相对误差')
ax.set_title('分解稳定性 vs 条件数')
ax.legend()
ax.grid(True, alpha=0.3)
# 空白子图用于说明
ax = axes[2, 3]
ax.text(0.1, 0.8, '分解方法特点:', fontsize=12, weight='bold', transform=ax.transAxes)
ax.text(0.1, 0.7, '• LU: 一般矩阵，O(n³)', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.6, '• QR: 数值稳定，正交化', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.5, '• SVD: 最稳定，最昂贵', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.4, '• Cholesky: 正定矩阵专用', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.2, '选择原则:', fontsize=12, weight='bold', transform=ax.transAxes)
ax.text(0.1, 0.1, '• 速度: LU > QR > SVD', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.0, '• 稳定性: SVD > QR > LU', fontsize=10, transform=ax.transAxes)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.axis('off')
plt.tight_layout()
plt.show()
def iterative_methods():
"""迭代法求解线性方程组"""
print(f"\n迭代法求解线性方程组")
print("=" * 40)
# 生成对角优势矩阵（保证收敛）
n = 8
A = np.random.randn(n, n)
A = A + (n + 1) * np.eye(n) # 使其对角占优
x_true = np.random.randn(n)
b = A @ x_true
print(f"系统规模: {n} × {n}")
print(f"矩阵条件数: {np.linalg.cond(A):.2e}")
# 直接法求解（作为参考）
x_direct = solve(A, b)
print(f"直接法误差: {norm(x_direct - x_true):.2e}")
def jacobi_method(A, b, x0, max_iter=100, tol=1e-10):
"""雅可比迭代法"""
n = len(b)
x = x0.copy()
D = np.diag(np.diag(A))
R = A - D
residuals = []
errors = []
for k in range(max_iter):
x_new = solve(D, b - R @ x)
residual = norm(A @ x_new - b)
error = norm(x_new - x_true)
residuals.append(residual)
errors.append(error)
if residual < tol:
return x_new, k+1, residuals, errors
x = x_new
return x, max_iter, residuals, errors
def gauss_seidel_method(A, b, x0, max_iter=100, tol=1e-10):
"""高斯-赛德尔迭代法"""
n = len(b)
x = x0.copy()
residuals = []
errors = []
for k in range(max_iter):
x_new = x.copy()
for i in range(n):
sum1 = np.dot(A[i, :i], x_new[:i])
sum2 = np.dot(A[i, i+1:], x[i+1:])
x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
residual = norm(A @ x_new - b)
error = norm(x_new - x_true)
residuals.append(residual)
errors.append(error)
if residual < tol:
return x_new, k+1, residuals, errors
x = x_new
return x, max_iter, residuals, errors
def sor_method(A, b, x0, omega=1.2, max_iter=100, tol=1e-10):
"""超松弛迭代法"""
n = len(b)
x = x0.copy()
residuals = []
errors = []
for k in range(max_iter):
x_new = x.copy()
for i in range(n):
sum1 = np.dot(A[i, :i], x_new[:i])
sum2 = np.dot(A[i, i+1:], x[i+1:])
x_gs = (b[i] - sum1 - sum2) / A[i, i]
x_new[i] = (1 - omega) * x[i] + omega * x_gs
residual = norm(A @ x_new - b)
error = norm(x_new - x_true)
residuals.append(residual)
errors.append(error)
if residual < tol:
return x_new, k+1, residuals, errors
x = x_new
return x, max_iter, residuals, errors
# 运行各种迭代方法
x0 = np.zeros(n) # 初始猜测
print(f"\n迭代方法比较:")
print("-" * 30)
methods = {
'Jacobi': jacobi_method,
'Gauss-Seidel': gauss_seidel_method,
'SOR (ω=1.2)': lambda A, b, x0, **kwargs: sor_method(A, b, x0, omega=1.2, **kwargs)
}
results = {}
for name, method in methods.items():
x_sol, iterations, residuals, errors = method(A, b, x0)
final_error = norm(x_sol - x_true)
results[name] = {
'solution': x_sol,
'iterations': iterations,
'residuals': residuals,
'errors': errors,
'final_error': final_error
}
print(f"{name:12}: {iterations:3d} 次迭代，最终误差 = {final_error:.2e}")
# 可视化
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# 收敛历史
ax1 = axes[0, 0]
for name, result in results.items():
iterations_list = range(1, len(result['residuals']) + 1)
ax1.semilogy(iterations_list, result['residuals'], 'o-',
linewidth=2, markersize=4, label=name)
ax1.set_xlabel('迭代次数')
ax1.set_ylabel('残差 ||Ax - b||')
ax1.set_title('残差收敛历史')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 误差收敛历史
ax2 = axes[0, 1]
for name, result in results.items():
iterations_list = range(1, len(result['errors']) + 1)
ax2.semilogy(iterations_list, result['errors'], 'o-',
linewidth=2, markersize=4, label=name)
ax2.set_xlabel('迭代次数')
ax2.set_ylabel('误差 ||x - x_true||')
ax2.set_title('误差收敛历史')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 收敛速度比较
ax3 = axes[0, 2]
method_names = list(results.keys())
iterations_count = [results[name]['iterations'] for name in method_names]
final_errors = [results[name]['final_error'] for name in method_names]
bars = ax3.bar(range(len(method_names)), iterations_count, alpha=0.7)
ax3.set_xlabel('方法')
ax3.set_ylabel('收敛所需迭代次数')
ax3.set_title('收敛速度比较')
ax3.set_xticks(range(len(method_names)))
ax3.set_xticklabels(method_names, rotation=45)
ax3.grid(True, alpha=0.3)
# 添加数值标签
for bar, count in zip(bars, iterations_count):
ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
f'{count}', ha='center', va='bottom')
# SOR参数优化
ax4 = axes[1, 0]
omega_values = np.linspace(0.5, 1.95, 30)
sor_iterations = []
for omega in omega_values:
try:
_, iters, _, _ = sor_method(A, b, x0, omega=omega, max_iter=200)
sor_iterations.append(iters)
except:
sor_iterations.append(200) # 未收敛
ax4.plot(omega_values, sor_iterations, 'b-', linewidth=2)
optimal_idx = np.argmin(sor_iterations)
optimal_omega = omega_values[optimal_idx]
ax4.plot(optimal_omega, sor_iterations[optimal_idx], 'ro', markersize=10,
label=f'最优 ω = {optimal_omega:.2f}')
ax4.set_xlabel('松弛因子 ω')
ax4.set_ylabel('收敛迭代次数')
ax4.set_title('SOR方法参数优化')
ax4.legend()
ax4.grid(True, alpha=0.3)
# 谱半径分析
ax5 = axes[1, 1]
# 计算迭代矩阵的谱半径
D = np.diag(np.diag(A))
L = np.tril(A, k=-1)
U = np.triu(A, k=1)
# Jacobi迭代矩阵
T_jacobi = -solve(D, L + U)
rho_jacobi = max(abs(np.linalg.eigvals(T_jacobi)))
# Gauss-Seidel迭代矩阵
T_gs = -solve(D + L, U)
rho_gs = max(abs(np.linalg.eigvals(T_gs)))
methods_spectral = ['Jacobi', 'Gauss-Seidel']
spectral_radii = [rho_jacobi, rho_gs]
bars = ax5.bar(methods_spectral, spectral_radii, alpha=0.7)
ax5.axhline(y=1, color='red', linestyle='--', label='收敛界限')
ax5.set_ylabel('谱半径')
ax5.set_title('迭代矩阵谱半径')
ax5.legend()
ax5.grid(True, alpha=0.3)
# 添加数值标签
for bar, rho in zip(bars, spectral_radii):
ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01,
f'{rho:.3f}', ha='center', va='bottom')
# 解的收敛轨迹（选择前两个分量）
ax6 = axes[1, 2]
for name, result in results.items():
if name == 'Jacobi': # 只显示一种方法的轨迹，避免图形混乱
# 重新运行以获取所有中间解
x = x0.copy()
trajectory_x = [x[0]]
trajectory_y = [x[1]]
D = np.diag(np.diag(A))
R = A - D
for k in range(min(50, result['iterations'])):
x = solve(D, b - R @ x)
trajectory_x.append(x[0])
trajectory_y.append(x[1])
ax6.plot(trajectory_x, trajectory_y, 'b-o', markersize=4,
alpha=0.7, label='Jacobi轨迹')
ax6.plot(x_true[0], x_true[1], 'rs', markersize=10, label='真解')
ax6.plot(x0[0], x0[1], 'go', markersize=8, label='初始值')
ax6.set_xlabel('x₁')
ax6.set_ylabel('x₂')
ax6.set_title('解的收敛轨迹')
ax6.legend()
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return results
def eigenvalue_algorithms():
"""特征值计算的数值方法"""
print(f"\n特征值计算的数值方法")
print("=" * 40)
# 生成测试矩阵
np.random.seed(42)
n = 6
A_random = np.random.randn(n, n)
A_symmetric = A_random + A_random.T # 对称矩阵
A_general = A_random + 0.1 * np.eye(n) # 一般矩阵
print("测试矩阵:")
print(f"对称矩阵 A_sym 的条件数: {np.linalg.cond(A_symmetric):.2e}")
print(f"一般矩阵 A_gen 的条件数: {np.linalg.cond(A_general):.2e}")
# 真实特征值（用作参考）
true_eigenvals_sym = np.linalg.eigvals(A_symmetric)
true_eigenvals_gen = np.linalg.eigvals(A_general)
# 按模长排序
true_eigenvals_sym = true_eigenvals_sym[np.argsort(np.abs(true_eigenvals_sym))[::-1]]
true_eigenvals_gen = true_eigenvals_gen[np.argsort(np.abs(true_eigenvals_gen))[::-1]]
print(f"\n对称矩阵真实特征值: {true_eigenvals_sym}")
print(f"一般矩阵真实特征值: {true_eigenvals_gen}")
def power_method(A, max_iter=100, tol=1e-10):
"""幂法求主特征值"""
n = A.shape[0]
x = np.random.randn(n)
x = x / norm(x)
eigenval_history = []
eigenvec_history = []
for k in range(max_iter):
y = A @ x
eigenval = x.T @ y
x_new = y / norm(y)
eigenval_history.append(eigenval)
eigenvec_history.append(x_new.copy())
if k > 0 and abs(eigenval_history[-1] - eigenval_history[-2]) < tol:
return eigenval, x_new, k+1, eigenval_history
x = x_new
return eigenval_history[-1], x, max_iter, eigenval_history
def inverse_power_method(A, sigma, max_iter=100, tol=1e-10):
"""反幂法求接近σ的特征值"""
n = A.shape[0]
x = np.random.randn(n)
x = x / norm(x)
A_shifted = A - sigma * np.eye(n)
eigenval_history = []
for k in range(max_iter):
try:
y = solve(A_shifted, x)
mu = x.T @ y
x_new = y / norm(y)
eigenval = 1/mu + sigma # 逆变换回原特征值
eigenval_history.append(eigenval)
if k > 0 and abs(eigenval_history[-1] - eigenval_history[-2]) < tol:
return eigenval, x_new, k+1, eigenval_history
x = x_new
except:
break
return eigenval_history[-1] if eigenval_history else sigma, x, k+1, eigenval_history
def qr_algorithm_basic(A, max_iter=100, tol=1e-10):
"""基础QR算法"""
A_k = A.copy()
eigenval_history = []
for k in range(max_iter):
Q, R = qr(A_k)
A_k = R @ Q
# 记录对角元素（特征值的近似）
eigenvals = np.diag(A_k)
eigenval_history.append(eigenvals.copy())
# 检查收敛（下三角部分接近零）
if k > 5:
lower_triangular_norm = norm(np.tril(A_k, k=-1))
if lower_triangular_norm < tol:
break
return np.diag(A_k), A_k, k+1, eigenval_history
# 运行各种方法
print(f"\n幂法求主特征值:")
print("-" * 25)
# 对称矩阵
eigenval_pm_sym, eigenvec_pm_sym, iter_pm_sym, hist_pm_sym = power_method(A_symmetric)
error_pm_sym = abs(eigenval_pm_sym - true_eigenvals_sym[0])
print(f"对称矩阵:")
print(f" 幂法结果: {eigenval_pm_sym:.6f} ({iter_pm_sym} 次迭代)")
print(f" 真实值: {true_eigenvals_sym[0]:.6f}")
print(f" 误差: {error_pm_sym:.2e}")
# 一般矩阵
eigenval_pm_gen, eigenvec_pm_gen, iter_pm_gen, hist_pm_gen = power_method(A_general)
error_pm_gen = abs(eigenval_pm_gen - true_eigenvals_gen[0])
print(f"一般矩阵:")
print(f" 幂法结果: {eigenval_pm_gen:.6f} ({iter_pm_gen} 次迭代)")
print(f" 真实值: {true_eigenvals_gen[0]:.6f}")
print(f" 误差: {error_pm_gen:.2e}")
# 反幂法
print(f"\n反幂法求特定特征值:")
print("-" * 30)
# 选择一个目标值
target = 0.5
eigenval_ipm, eigenvec_ipm, iter_ipm, hist_ipm = inverse_power_method(A_symmetric, target)
# 找到最接近的真实特征值
closest_true = true_eigenvals_sym[np.argmin(np.abs(true_eigenvals_sym - target))]
error_ipm = abs(eigenval_ipm - closest_true)
print(f"目标值: {target}")
print(f"反幂法结果: {eigenval_ipm:.6f} ({iter_imp} 次迭代)" if 'iter_ipm' in locals() else f"反幂法结果: {eigenval_ipm:.6f}")
print(f"最接近的真实特征值: {closest_true:.6f}")
print(f"误差: {error_ipm:.2e}")
# QR算法
print(f"\nQR算法求所有特征值:")
print("-" * 30)
eigenvals_qr_sym, A_final_sym, iter_qr_sym, hist_qr_sym = qr_algorithm_basic(A_symmetric)
eigenvals_qr_gen, A_final_gen, iter_qr_gen, hist_qr_gen = qr_algorithm_basic(A_general)
print(f"对称矩阵 ({iter_qr_sym} 次迭代):")
print(f" QR结果: {np.sort(eigenvals_qr_sym)[::-1]}")
print(f" 真实值: {np.sort(true_eigenvals_sym)[::-1]}")
print(f" 误差: {norm(np.sort(eigenvals_qr_sym) - np.sort(true_eigenvals_sym)):.2e}")
print(f"一般矩阵 ({iter_qr_gen} 次迭代):")
print(f" QR结果: {np.sort(np.real(eigenvals_qr_gen))[::-1]}")
print(f" 真实值: {np.sort(np.real(true_eigenvals_gen))[::-1]}")
# 可视化
fig, axes = plt.subplots(3, 3, figsize=(18, 18))
# 幂法收敛历史
ax1 = axes[0, 0]
ax1.semilogy(range(1, len(hist_pm_sym)+1),
[abs(val - true_eigenvals_sym[0]) for val in hist_pm_sym],
'b-o', linewidth=2, markersize=4, label='对称矩阵')
ax1.semilogy(range(1, len(hist_pm_gen)+1),
[abs(val - true_eigenvals_gen[0]) for val in hist_pm_gen],
'r-s', linewidth=2, markersize=4, label='一般矩阵')
ax1.set_xlabel('迭代次数')
ax1.set_ylabel('误差')
ax1.set_title('幂法收敛历史')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 反幂法收敛历史
ax2 = axes[0, 1]
if hist_ipm:
ax2.semilogy(range(1, len(hist_ipm)+1),
[abs(val - closest_true) for val in hist_ipm],
'g-o', linewidth=2, markersize=4)
ax2.set_xlabel('迭代次数')
ax2.set_ylabel('误差')
ax2.set_title(f'反幂法收敛历史 (目标={target})')
ax2.grid(True, alpha=0.3)
# QR算法收敛历史（显示最大非对角元素）
ax3 = axes[0, 2]
def off_diagonal_norm(eigenval_matrix):
"""计算矩阵的非对角元素范数"""
return norm(eigenval_matrix - np.diag(np.diag(eigenval_matrix)))
if len(hist_qr_sym) > 1:
# 这里需要重新运行QR算法来获取矩阵历史
# 简化显示：只显示特征值的变化
max_changes_sym = []
for i in range(1, len(hist_qr_sym)):
change = max(abs(hist_qr_sym[i] - hist_qr_sym[i-1]))
max_changes_sym.append(change)
if max_changes_sym:
ax3.semilogy(range(1, len(max_changes_sym)+1), max_changes_sym,
'b-o', linewidth=2, markersize=4, label='对称矩阵')
ax3.set_xlabel('迭代次数')
ax3.set_ylabel('最大特征值变化')
ax3.set_title('QR算法收敛历史')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 特征值比较
ax4 = axes[1, 0]
x_pos = range(len(true_eigenvals_sym))
width = 0.35
ax4.bar([x - width/2 for x in x_pos], np.sort(true_eigenvals_sym)[::-1], width,
label='真实值', alpha=0.7)
ax4.bar([x + width/2 for x in x_pos], np.sort(eigenvals_qr_sym)[::-1], width,
label='QR算法', alpha=0.7)
ax4.set_xlabel('特征值索引')
ax4.set_ylabel('特征值')
ax4.set_title('对称矩阵特征值比较')
ax4.legend()
ax4.grid(True, alpha=0.3)
# 计算时间比较
ax5 = axes[1, 1]
# 简化的时间测试
methods = ['Power Method', 'QR Algorithm', 'NumPy (eig)']
times = []
# 幂法时间（估算）
times.append(iter_pm_sym * 1e-4) # 估算每次迭代时间
# QR算法时间（估算）
times.append(iter_qr_sym * 1e-3) # QR每次迭代更昂贵
# NumPy时间
start_time = time.time()
for _ in range(10):
np.linalg.eig(A_symmetric)
numpy_time = (time.time() - start_time) / 10
times.append(numpy_time)
bars = ax5.bar(methods, times, alpha=0.7)
ax5.set_ylabel('时间 (秒)')
ax5.set_title('计算时间比较 (估算)')
ax5.set_xticklabels(methods, rotation=45)
for bar, time_val in zip(bars, times):
ax5.text(bar.get_x() + bar.get_width()/2, bar.get_height(),
f'{time_val:.2e}', ha='center', va='bottom')
# 特征向量可视化（前两个分量）
ax6 = axes[1, 2]
# 比较主特征向量
true_eigenvec = np.linalg.eig(A_symmetric)[1][:, 0]
# 确保符号一致
if np.dot(eigenvec_pm_sym, true_eigenvec) < 0:
eigenvec_pm_sym = -eigenvec_pm_sym
ax6.arrow(0, 0, true_eigenvec[0], true_eigenvec[1],
head_width=0.05, head_length=0.05, fc='blue', ec='blue',
linewidth=3, label='真实特征向量')
ax6.arrow(0, 0, eigenvec_pm_sym[0], eigenvec_pm_sym[1],
head_width=0.05, head_length=0.05, fc='red', ec='red',
linewidth=3, alpha=0.7, label='幂法结果')
ax6.set_xlim(-1.2, 1.2)
ax6.set_ylim(-1.2, 1.2)
ax6.set_xlabel('v₁')
ax6.set_ylabel('v₂')
ax6.set_title('主特征向量比较')
ax6.legend()
ax6.grid(True, alpha=0.3)
ax6.set_aspect('equal')
# 算法复杂度比较
ax7 = axes[2, 0]
matrix_sizes = [10, 20, 50, 100]
complexities = {
'Power Method': [n for n in matrix_sizes], # O(n²) per iteration
'QR Algorithm': [n**3 for n in matrix_sizes], # O(n³)
'Direct (eig)': [n**3 for n in matrix_sizes] # O(n³)
}
for method, complexity in complexities.items():
ax7.loglog(matrix_sizes, complexity, 'o-', linewidth=2, label=method)
ax7.set_xlabel('矩阵大小 n')
ax7.set_ylabel('计算复杂度')
ax7.set_title('算法复杂度比较')
ax7.legend()
ax7.grid(True, alpha=0.3)
# 收敛速度比较
ax8 = axes[2, 1]
# 模拟不同特征值间隙对幂法收敛的影响
gaps = [0.1, 0.3, 0.5, 0.7, 0.9] # λ₂/λ₁ 比值
convergence_rates = []
for gap in gaps:
# 理论收敛率 ≈ |λ₂/λ₁|
rate = gap
convergence_rates.append(rate)
ax8.plot(gaps, convergence_rates, 'bo-', linewidth=2, markersize=8)
ax8.set_xlabel('特征值间隙 (λ₂/λ₁)')
ax8.set_ylabel('收敛率')
ax8.set_title('幂法收敛率 vs 特征值间隙')
ax8.grid(True, alpha=0.3)
# 方法选择指南
ax9 = axes[2, 2]
ax9.text(0.1, 0.9, '特征值算法选择指南', fontsize=14, weight='bold', transform=ax9.transAxes)
guidelines = [
'单个主特征值：',
'• 幂法 - 简单、内存效率高',
'• 适用于大型稀疏矩阵',
'',
'特定特征值：',
'• 反幂法 - 收敛快',
'• 需要近似值',
'',
'所有特征值：',
'• QR算法 - 数值稳定',
'• 对称矩阵优化版本',
'',
'大规模问题：',
'• Lanczos/Arnoldi方法',
'• Krylov子空间方法'
]
for i, text in enumerate(guidelines):
y_pos = 0.8 - i * 0.05
if text.startswith('•'):
ax9.text(0.15, y_pos, text, fontsize=10, transform=ax9.transAxes)
elif text == '':
continue
else:
ax9.text(0.1, y_pos, text, fontsize=11, weight='bold', transform=ax9.transAxes)
ax9.set_xlim(0, 1)
ax9.set_ylim(0, 1)
ax9.axis('off')
plt.tight_layout()
plt.show()
# 运行所有演示
condition_nums = numerical_stability_demo()
matrix_factorizations_numerical()
iter_results = iterative_methods()
eigenvalue_algorithms()
# 最终总结
def numerical_linear_algebra_summary():
"""数值线性代数总结"""
print(f"\n数值线性代数总结")
print("=" * 50)
summary = {
"误差分析": {
"舍入误差": "浮点数精度限制",
"条件数": "系统对输入扰动的敏感性",
"数值稳定性": "算法对舍入误差的鲁棒性",
"后向误差": "实际解对应的扰动问题"
},
"矩阵分解": {
"LU分解": "A = PLU，高效求解",
"QR分解": "A = QR，数值稳定",
"SVD分解": "A = UΣVᵀ，最稳定最通用",
"Cholesky": "A = LLᵀ，正定矩阵专用"
},
"线性方程组": {
"直接法": "分解求解，精确但昂贵",
"Jacobi": "简单并行，收敛慢",
"Gauss-Seidel": "收敛更快，串行",
"SOR": "加速收敛，需调参"
},
"特征值问题": {
"幂法": "主特征值，简单高效",
"反幂法": "特定特征值，快速收敛",
"QR算法": "所有特征值，工业标准",
"Krylov方法": "大规模稀疏矩阵"
}
}
for category, methods in summary.items():
print(f"\n{category}:")
print("-" * 30)
for method, description in methods.items():
print(f"• {method:15}: {description}")
print(f"\n实践建议:")
print("-" * 20)
recommendations = [
"• 优先考虑数值稳定性而非速度",
"• 利用矩阵的特殊结构（对称、稀疏等）",
"• 监控条件数，警惕病态问题",
"• 选择合适的分解方法",
"• 使用成熟的数值库（LAPACK, BLAS）",
"• 验证结果的数值精度"
]
for rec in recommendations:
print(rec)
numerical_linear_algebra_summary()
# 完成所有章节
print(f"\n" + "="*60)
print("线性代数完整课程学习笔记生成完成！")
print("="*60)
print(f"已生成章节：")
for i in range(8, 16):
print(f"• 第{i}章：la-{i:03d}.md")
# 最后更新todo列表
from datetime import datetime
print(f"\n生成完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")