Chapter 15: Introduction to Numerical Linear Algebra
Haiyue
39min
Chapter 15: Introduction to Numerical Linear Algebra
Learning Objectives
- Understand error issues in numerical computation
- Master numerical methods for matrix decomposition
- Learn iterative methods for solving linear equations
- Understand condition numbers and stability
- Master numerical methods for eigenvalue computation
Errors and Stability in Numerical Computation
Floating-Point Representation and Round-off Errors
The representation of real numbers in computers has precision limitations, leading to the accumulation of round-off errors.
Condition Number
The condition number of a matrix measures the sensitivity of a linear system to perturbations in the input:
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():
"""Numerical Stability and Error Analysis Demonstration"""
print("Numerical Linear Algebra: Stability and Error Analysis")
print("=" * 60)
# 1. Machine Epsilon Demonstration
print("1. Machine Epsilon Analysis")
print("-" * 30)
machine_epsilon = np.finfo(float).eps
print(f"Machine epsilon ε: {machine_epsilon:.2e}")
# Demonstrate floating-point precision issues
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"Equal? {a == b}")
print(f"Difference: {abs(a - b):.2e}")
# 2. Ill-conditioned Matrix Example
print(f"\n2. Ill-conditioned Matrix Analysis")
print("-" * 30)
# Hilbert matrix (a famous ill-conditioned matrix)
def hilbert_matrix(n):
"""Generates an n-order Hilbert matrix"""
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 matrices of different sizes
sizes = [3, 5, 8, 10]
condition_numbers = []
print("Hilbert Matrix Condition Number Analysis:")
for n in sizes:
H = hilbert_matrix(n)
cond_num = np.linalg.cond(H)
condition_numbers.append(cond_num)
print(f"H_{n}: Condition number = {cond_num:.2e}")
# 3. Effect of Condition Number on Solution Accuracy
print(f"\n3. Effect of Condition Number on Solution Accuracy")
print("-" * 40)
n = 6
H = hilbert_matrix(n)
x_exact = np.ones(n) # Exact solution
b_exact = H @ x_exact # Exact right-hand side
# Add small perturbations
perturbation_levels = [1e-12, 1e-10, 1e-8, 1e-6]
results = []
for eps in perturbation_levels:
# Perturb the right-hand side
b_perturbed = b_exact + eps * np.random.randn(n)
# Solve the perturbed system
try:
x_computed = solve(H, b_perturbed)
relative_error = norm(x_computed - x_exact) / norm(x_exact)
results.append((eps, relative_error))
print(f"Perturbation level {eps:.0e}: Relative error = {relative_error:.2e}")
except:
print(f"Perturbation level {eps:.0e}: Solution failed")
# Visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# Condition number growth
ax1 = axes[0, 0]
ax1.semilogy(sizes, condition_numbers, 'o-', linewidth=2, markersize=8)
ax1.set_xlabel('Matrix Size')
ax1.set_ylabel('Condition Number')
ax1.set_title('Growth of Hilbert Matrix Condition Number')
ax1.grid(True, alpha=0.3)
# Visualize Hilbert matrix
ax2 = axes[0, 1]
H_display = hilbert_matrix(8)
im = ax2.imshow(H_display, cmap='Blues')
ax2.set_title('8×8 Hilbert Matrix')
plt.colorbar(im, ax=ax2)
# Error amplification effect
ax3 = axes[0, 2]
if results:
perturbations, errors = zip(*results)
ax3.loglog(perturbations, errors, 'ro-', linewidth=2, markersize=8, label='Actual Error')
# Theoretical error bound
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='Theoretical Error Bound')
ax3.set_xlabel('Input Perturbation Level')
ax3.set_ylabel('Output Relative Error')
ax3.set_title('Error Amplification Effect')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Comparison of condition numbers for different matrix types
ax4 = axes[1, 0]
matrix_types = {
'Identity Matrix': np.eye(5),
'Diagonally Dominant': np.diag([10, 5, 3, 2, 1]) + 0.1*np.random.randn(5, 5),
'Hilbert-5': hilbert_matrix(5),
'Random Matrix': np.random.randn(5, 5)
}
names = []
conds = []
for name, A in matrix_types.items():
if name != 'Hilbert-5':
A = A @ A.T # Ensure positive definite
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('Matrix Type')
ax4.set_ylabel('Condition Number (Log Scale)')
ax4.set_title('Condition Numbers of Different Matrix Types')
ax4.set_xticks(range(len(names)))
ax4.set_xticklabels(names, rotation=45)
ax4.grid(True, alpha=0.3)
# Round-off error accumulation demonstration
ax5 = axes[1, 1]
def sum_series(n, method='forward'):
"""Different methods for summing a series"""
if method == 'forward':
# Forward summation
result = 0.0
for i in range(1, n+1):
result += 1.0 / i**2
else: # backward
# Backward summation
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='Forward Summation')
ax5.semilogx(n_values, backward_results, 'b-', label='Backward Summation')
ax5.axhline(y=np.pi**2/6, color='black', linestyle='--', label='Theoretical Value π²/6')
ax5.set_xlabel('Number of Terms')
ax5.set_ylabel('Series Sum')
ax5.set_title('Accumulation of Round-off Error: Σ(1/k²)')
ax5.legend()
ax5.grid(True, alpha=0.3)
# Numerical stability comparison
ax6 = axes[1, 2]
# Compare numerical stability of different algorithms
def unstable_algorithm(x):
"""Example of a numerically unstable algorithm"""
return (x**2 - 1) / (x - 1)
def stable_algorithm(x):
"""Numerically stable algorithm"""
return x + 1
x_values = np.linspace(0.99, 1.01, 100)
x_values = x_values[x_values != 1.0] # Avoid division by zero
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='Unstable: (x²-1)/(x-1)', linewidth=2)
ax6.plot(x_values, stable_results, 'b-', label='Stable: x+1', linewidth=2)
ax6.set_xlabel('x')
ax6.set_ylabel('Function Value')
ax6.set_title('Numerical Stability Comparison')
ax6.legend()
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return condition_numbers
def matrix_factorizations_numerical():
"""Numerical Implementation of Matrix Factorizations"""
print(f"\nNumerical Methods for Matrix Factorization")
print("=" * 40)
# Generate test matrix
np.random.seed(42)
n = 6
A_random = np.random.randn(n, n)
A_spd = A_random @ A_random.T # Symmetric positive definite matrix
A_general = np.random.randn(n, n) + 0.1 * np.eye(n) # General matrix
print("Test Matrix Information:")
print(f"Condition number of random symmetric positive definite matrix: {np.linalg.cond(A_spd):.2e}")
print(f"Condition number of general matrix: {np.linalg.cond(A_general):.2e}")
# 1. LU Decomposition
print(f"\n1. LU Decomposition")
print("-" * 20)
P, L, U = lu(A_general)
print("LU Decomposition Verification:")
print(f"||PA - LU||: {norm(P @ A_general - L @ U):.2e}")
# Check properties of L and U
print(f"L is lower triangular: {np.allclose(L, np.tril(L))}")
print(f"U is upper triangular: {np.allclose(U, np.triu(U))}")
print(f"P is a permutation matrix: {np.allclose(P @ P.T, np.eye(n))}")
# 2. QR Decomposition
print(f"\n2. QR Decomposition")
print("-" * 20)
Q, R = qr(A_general)
print("QR Decomposition Verification:")
print(f"||A - QR||: {norm(A_general - Q @ R):.2e}")
print(f"Q is orthogonal: {np.allclose(Q @ Q.T, np.eye(n))}")
print(f"R is upper triangular: {np.allclose(R, np.triu(R))}")
# 3. SVD Decomposition
print(f"\n3. Singular Value Decomposition")
print("-" * 20)
U_svd, s, Vt = svd(A_general)
print("SVD Decomposition Verification:")
print(f"||A - UΣVᵀ||: {norm(A_general - U_svd @ np.diag(s) @ Vt):.2e}")
print(f"Singular values: {s}")
print(f"Condition number (σ_max/σ_min): {s[0]/s[-1]:.2e}")
# 4. Cholesky Decomposition (only for positive definite matrices)
print(f"\n4. Cholesky Decomposition")
print("-" * 20)
try:
from scipy.linalg import cholesky
L_chol = cholesky(A_spd, lower=True)
print("Cholesky Decomposition Verification:")
print(f"||A - LLᵀ||: {norm(A_spd - L_chol @ L_chol.T):.2e}")
print(f"L is lower triangular: {np.allclose(L_chol, np.tril(L_chol))}")
# Cholesky vs LU computational complexity comparison
print(f"Cholesky decomposition saves approximately 50% computation (for positive definite matrices)")
except Exception as e:
print(f"Cholesky decomposition failed: {e}")
# Visualize matrix factorization
fig, axes = plt.subplots(3, 4, figsize=(16, 12))
# Original matrix
ax = axes[0, 0]
im = ax.imshow(A_general, cmap='RdBu')
ax.set_title('Original Matrix A')
plt.colorbar(im, ax=ax)
# LU Decomposition
ax = axes[0, 1]
im = ax.imshow(L, cmap='Blues')
ax.set_title('L (Lower Triangular)')
plt.colorbar(im, ax=ax)
ax = axes[0, 2]
im = ax.imshow(U, cmap='Reds')
ax.set_title('U (Upper Triangular)')
plt.colorbar(im, ax=ax)
ax = axes[0, 3]
im = ax.imshow(P, cmap='Greys')
ax.set_title('P (Permutation)')
plt.colorbar(im, ax=ax)
# QR Decomposition
ax = axes[1, 0]
im = ax.imshow(Q, cmap='RdBu')
ax.set_title('Q (Orthogonal)')
plt.colorbar(im, ax=ax)
ax = axes[1, 1]
im = ax.imshow(R, cmap='Reds')
ax.set_title('R (Upper Triangular)')
plt.colorbar(im, ax=ax)
# SVD Decomposition
ax = axes[1, 2]
im = ax.imshow(U_svd, cmap='Blues')
ax.set_title('U (Left Singular Vectors)')
plt.colorbar(im, ax=ax)
ax = axes[1, 3]
ax.bar(range(len(s)), s)
ax.set_title('Singular Values')
ax.set_xlabel('Index')
ax.set_ylabel('Value')
ax.grid(True, alpha=0.3)
# Decomposition accuracy comparison
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('Reconstruction Error')
ax.set_title('Comparison of Decomposition Accuracy')
ax.grid(True, alpha=0.3)
# Computation time comparison
ax = axes[2, 1]
def time_factorization(A, method, n_trials=10):
"""Tests the computation time of factorization methods"""
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('Average Time (seconds)')
ax.set_title('Computation Time Comparison')
ax.grid(True, alpha=0.3)
# Condition number and decomposition stability
ax = axes[2, 2]
# Generate matrices with different condition numbers
condition_numbers = [1e2, 1e4, 1e6, 1e8]
lu_errors = []
qr_errors = []
for cond in condition_numbers:
# Construct a matrix with a specified condition number
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
# Test decomposition accuracy
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 Decomposition', linewidth=2)
ax.loglog(condition_numbers, qr_errors, 's-', label='QR Decomposition', linewidth=2)
ax.set_xlabel('Condition Number')
ax.set_ylabel('Relative Error')
ax.set_title('Decomposition Stability vs. Condition Number')
ax.legend()
ax.grid(True, alpha=0.3)
# Blank subplot for illustration
ax = axes[2, 3]
ax.text(0.1, 0.8, 'Decomposition Method Features:', fontsize=12, weight='bold', transform=ax.transAxes)
ax.text(0.1, 0.7, '• LU: General matrix, O(n³)', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.6, '• QR: Numerically stable, orthogonalization', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.5, '• SVD: Most stable, most expensive', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.4, '• Cholesky: Specific for positive definite matrices', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.2, 'Selection Principle:', fontsize=12, weight='bold', transform=ax.transAxes)
ax.text(0.1, 0.1, '• Speed: LU > QR > SVD', fontsize=10, transform=ax.transAxes)
ax.text(0.1, 0.0, '• Stability: 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():
"""Iterative Methods for Solving Linear Equations"""
print(f"\nIterative Methods for Solving Linear Equations")
print("=" * 40)
# Generate diagonally dominant matrix (guarantees convergence)
n = 8
A = np.random.randn(n, n)
A = A + (n + 1) * np.eye(n) # Make it diagonally dominant
x_true = np.random.randn(n)
b = A @ x_true
print(f"System size: {n} × {n}")
print(f"Matrix condition number: {np.linalg.cond(A):.2e}")
# Direct method solution (for reference)
x_direct = solve(A, b)
print(f"Direct method error: {norm(x_direct - x_true):.2e}")
def jacobi_method(A, b, x0, max_iter=100, tol=1e-10):
"""Jacobi Iteration Method"""
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):
"""Gauss-Seidel Iteration Method"""
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):
"""Successive Over-Relaxation (SOR) Method"""
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
# Run various iterative methods
x0 = np.zeros(n) # Initial guess
print(f"\nComparison of Iterative Methods:")
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} iterations, final error = {final_error:.2e}")
# Visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
# Convergence history
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('Iterations')
ax1.set_ylabel('Residual ||Ax - b||')
ax1.set_title('Residual Convergence History')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Error convergence history
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('Iterations')
ax2.set_ylabel('Error ||x - x_true||')
ax2.set_title('Error Convergence History')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Convergence speed comparison
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('Method')
ax3.set_ylabel('Iterations Required for Convergence')
ax3.set_title('Convergence Speed Comparison')
ax3.set_xticks(range(len(method_names)))
ax3.set_xticklabels(method_names, rotation=45)
ax3.grid(True, alpha=0.3)
# Add numerical labels
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 parameter optimization
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) # Not converged
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 ω = {optimal_omega:.2f}')
ax4.set_xlabel('Relaxation Factor ω')
ax4.set_ylabel('Iterations for Convergence')
ax4.set_title('SOR Method Parameter Optimization')
ax4.legend()
ax4.grid(True, alpha=0.3)
# Spectral radius analysis
ax5 = axes[1, 1]
# Calculate the spectral radius of the iteration matrices
D = np.diag(np.diag(A))
L = np.tril(A, k=-1)
U = np.triu(A, k=1)
# Jacobi iteration matrix
T_jacobi = -solve(D, L + U)
rho_jacobi = max(abs(np.linalg.eigvals(T_jacobi)))
# Gauss-Seidel iteration matrix
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='Convergence Limit')
ax5.set_ylabel('Spectral Radius')
ax5.set_title('Spectral Radius of Iteration Matrices')
ax5.legend()
ax5.grid(True, alpha=0.3)
# Add numerical labels
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')
# Convergence trajectory of the solution (selecting the first two components)
ax6 = axes[1, 2]
for name, result in results.items():
if name == 'Jacobi': # Only show trajectory for one method to avoid clutter
# Rerun to get all intermediate solutions
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 Trajectory')
ax6.plot(x_true[0], x_true[1], 'rs', markersize=10, label='True Solution')
ax6.plot(x0[0], x0[1], 'go', markersize=8, label='Initial Value')
ax6.set_xlabel('x₁')
ax6.set_ylabel('x₂')
ax6.set_title('Convergence Trajectory of the Solution')
ax6.legend()
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return results
def eigenvalue_algorithms():
"""Numerical Methods for Eigenvalue Computation"""
print(f"\nNumerical Methods for Eigenvalue Computation")
print("=" * 40)
# Generate test matrix
np.random.seed(42)
n = 6
A_random = np.random.randn(n, n)
A_symmetric = A_random + A_random.T # Symmetric matrix
A_general = A_random + 0.1 * np.eye(n) # General matrix
print("Test Matrices:")
print(f"Condition number of symmetric matrix A_sym: {np.linalg.cond(A_symmetric):.2e}")
print(f"Condition number of general matrix A_gen: {np.linalg.cond(A_general):.2e}")
# True eigenvalues (for reference)
true_eigenvals_sym = np.linalg.eigvals(A_symmetric)
true_eigenvals_gen = np.linalg.eigvals(A_general)
# Sort by absolute value
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"\nTrue eigenvalues for symmetric matrix: {true_eigenvals_sym}")
print(f"True eigenvalues for general matrix: {true_eigenvals_gen}")
def power_method(A, max_iter=100, tol=1e-10):
"""Power method for dominant eigenvalue"""
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, x, max_iter, eigenval_history
def inverse_power_method(A, sigma, max_iter=100, tol=1e-10):
"""Inverse power method for eigenvalue closest to σ"""
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 # Transform back to original eigenvalue
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):
"""Basic QR Algorithm"""
A_k = A.copy()
eigenval_history = []
for k in range(max_iter):
Q, R = qr(A_k)
A_k = R @ Q
# Record diagonal elements (approximation of eigenvalues)
eigenvals = np.diag(A_k)
eigenval_history.append(eigenvals.copy())
# Check for convergence (lower triangular part approaches zero)
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
# Run various methods
print(f"\nPower Method for Dominant Eigenvalue:")
print("-" * 25)
# Symmetric matrix
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"Symmetric matrix:")
print(f" Power method result: {eigenval_pm_sym:.6f} ({iter_pm_sym} iterations)")
print(f" True value: {true_eigenvals_sym[0]:.6f}")
print(f" Error: {error_pm_sym:.2e}")
# General matrix
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"General matrix:")
print(f" Power method result: {eigenval_pm_gen:.6f} ({iter_pm_gen} iterations)")
print(f" True value: {true_eigenvals_gen[0]:.6f}")
print(f" Error: {error_pm_gen:.2e}")
# Inverse Power Method
print(f"\nInverse Power Method for Specific Eigenvalue:")
print("-" * 30)
# Choose a target value
target = 0.5
eigenval_ipm, eigenvec_ipm, iter_ipm, hist_ipm = inverse_power_method(A_symmetric, target)
# Find the closest true eigenvalue
closest_true = true_eigenvals_sym[np.argmin(np.abs(true_eigenvals_sym - target))]
error_ipm = abs(eigenval_ipm - closest_true)
print(f"Target value: {target}")
print(f"Inverse power method result: {eigenval_ipm:.6f} ({iter_ipm} iterations)" if 'iter_ipm' in locals() else f"Inverse power method result: {eigenval_ipm:.6f}")
print(f"Closest true eigenvalue: {closest_true:.6f}")
print(f"Error: {error_ipm:.2e}")
# QR Algorithm
print(f"\nQR Algorithm for All Eigenvalues:")
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"Symmetric matrix ({iter_qr_sym} iterations):")
print(f" QR result: {np.sort(eigenvals_qr_sym)[::-1]}")
print(f" True value: {np.sort(true_eigenvals_sym)[::-1]}")
print(f" Error: {norm(np.sort(eigenvals_qr_sym) - np.sort(true_eigenvals_sym)):.2e}")
print(f"General matrix ({iter_qr_gen} iterations):")
print(f" QR result: {np.sort(np.real(eigenvals_qr_gen))[::-1]}")
print(f" True value: {np.sort(np.real(true_eigenvals_gen))[::-1]}")
# Visualization
fig, axes = plt.subplots(3, 3, figsize=(18, 18))
# Power method convergence history
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='Symmetric Matrix')
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='General Matrix')
ax1.set_xlabel('Iterations')
ax1.set_ylabel('Error')
ax1.set_title('Power Method Convergence History')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Inverse power method convergence history
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('Iterations')
ax2.set_ylabel('Error')
ax2.set_title(f'Inverse Power Method Convergence History (Target={target})')
ax2.grid(True, alpha=0.3)
# QR algorithm convergence history (showing max off-diagonal element)
ax3 = axes[0, 2]
def off_diagonal_norm(eigenval_matrix):
"""Calculates the norm of off-diagonal elements of a matrix"""
return norm(eigenval_matrix - np.diag(np.diag(eigenval_matrix)))
if len(hist_qr_sym) > 1:
# Re-run QR algorithm to get matrix history for off-diagonal norms
# Simplified display: only show changes in eigenvalues
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='Symmetric Matrix')
ax3.set_xlabel('Iterations')
ax3.set_ylabel('Max Eigenvalue Change')
ax3.set_title('QR Algorithm Convergence History')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Eigenvalue comparison
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='True Values', alpha=0.7)
ax4.bar([x + width/2 for x in x_pos], np.sort(eigenvals_qr_sym)[::-1], width,
label='QR Algorithm', alpha=0.7)
ax4.set_xlabel('Eigenvalue Index')
ax4.set_ylabel('Eigenvalue')
ax4.set_title('Comparison of Symmetric Matrix Eigenvalues')
ax4.legend()
ax4.grid(True, alpha=0.3)
# Computation time comparison
ax5 = axes[1, 1]
# Simplified time testing
methods = ['Power Method', 'QR Algorithm', 'NumPy (eig)']
times = []
# Power method time (estimate)
times.append(iter_pm_sym * 1e-4) # Estimate time per iteration
# QR algorithm time (estimate)
times.append(iter_qr_sym * 1e-3) # QR is more expensive per iteration
# NumPy time
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('Time (seconds)')
ax5.set_title('Computation Time Comparison (Estimated)')
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')
# Eigenvector visualization (first two components)
ax6 = axes[1, 2]
# Compare dominant eigenvectors
true_eigenvec = np.linalg.eig(A_symmetric)[1][:, 0]
# Ensure consistent sign
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='True Eigenvector')
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='Power Method Result')
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('Dominant Eigenvector Comparison')
ax6.legend()
ax6.grid(True, alpha=0.3)
ax6.set_aspect('equal')
# Algorithm complexity comparison
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('Matrix Size n')
ax7.set_ylabel('Computational Complexity')
ax7.set_title('Algorithm Complexity Comparison')
ax7.legend()
ax7.grid(True, alpha=0.3)
# Convergence speed comparison
ax8 = axes[2, 1]
# Simulate effect of eigenvalue gap on power method convergence
gaps = [0.1, 0.3, 0.5, 0.7, 0.9] # λ₂/λ₁ ratio
convergence_rates = []
for gap in gaps:
# Theoretical convergence rate ≈ |λ₂/λ₁|
rate = gap
convergence_rates.append(rate)
ax8.plot(gaps, convergence_rates, 'bo-', linewidth=2, markersize=8)
ax8.set_xlabel('Eigenvalue Gap (λ₂/λ₁)')
ax8.set_ylabel('Convergence Rate')
ax8.set_title('Power Method Convergence Rate vs. Eigenvalue Gap')
ax8.grid(True, alpha=0.3)
# Method selection guidelines
ax9 = axes[2, 2]
ax9.text(0.1, 0.9, 'Eigenvalue Algorithm Selection Guide', fontsize=14, weight='bold', transform=ax9.transAxes)
guidelines = [
'Single Dominant Eigenvalue:',
'• Power Method - Simple, memory efficient',
'• Suitable for large sparse matrices',
'',
'Specific Eigenvalue:',
'• Inverse Power Method - Fast convergence',
'• Requires an approximate value',
'',
'All Eigenvalues:',
'• QR Algorithm - Numerically stable',
'• Optimized versions for symmetric matrices',
'',
'Large-Scale Problems:',
'• Lanczos/Arnoldi Methods',
'• Krylov Subspace Methods'
]
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()
# Run all demonstrations
condition_nums = numerical_stability_demo()
matrix_factorizations_numerical()
iter_results = iterative_methods()
eigenvalue_algorithms()
# Final summary
def numerical_linear_algebra_summary():
"""Summary of Numerical Linear Algebra"""
print(f"\nSummary of Numerical Linear Algebra")
print("=" * 50)
summary = {
"Error Analysis": {
"Round-off Error": "Limitations of floating-point precision",
"Condition Number": "Sensitivity of the system to input perturbations",
"Numerical Stability": "Robustness of the algorithm to round-off errors",
"Backward Error": "Perturbed problem corresponding to the actual solution"
},
"Matrix Factorization": {
"LU Decomposition": "A = PLU, efficient for solving",
"QR Decomposition": "A = QR, numerically stable",
"SVD Decomposition": "A = UΣVᵀ, most stable and general",
"Cholesky": "A = LLᵀ, specific for positive definite matrices"
},
"Linear Equations": {
"Direct Methods": "Factorization for solving, exact but expensive",
"Jacobi": "Simple, parallelizable, slow convergence",
"Gauss-Seidel": "Faster convergence, sequential",
"SOR": "Accelerated convergence, requires parameter tuning"
},
"Eigenvalue Problems": {
"Power Method": "Dominant eigenvalue, simple and efficient",
"Inverse Power Method": "Specific eigenvalue, fast convergence",
"QR Algorithm": "All eigenvalues, industry standard",
"Krylov Methods": "Large sparse matrices"
}
}
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"\nPractical Recommendations:")
print("-" * 20)
recommendations = [
"• Prioritize numerical stability over speed",
"• Exploit special matrix structures (symmetric, sparse, etc.)",
"• Monitor condition numbers, be wary of ill-conditioned problems",
"• Choose appropriate factorization methods",
"• Use mature numerical libraries (LAPACK, BLAS)",
"• Validate the numerical accuracy of results"
]
for rec in recommendations:
print(rec)
numerical_linear_algebra_summary()
# Complete all chapters
print(f"\n" + "="*60)
print("Linear Algebra Complete Course Learning Notes Generation Completed!")
print("="*60)
print(f"Chapters generated:")
for i in range(8, 16):
print(f"• Chapter {i}: la-{i:03d}.md")
# Final update to todo list
from datetime import datetime
print(f"\nGeneration completed at: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")