Linear Algebra Review¶

Overview¶

Linear algebra plays a crucial role in numerical simulation. We will learn how to implement matrix operations, solving systems of linear equations, eigenvalue problems, and matrix decomposition using NumPy/SciPy.

1. Basic Matrix Operations¶

1.1 Matrix Creation and Operations¶

import numpy as np
from scipy import linalg

# Matrix creation
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 10]])

B = np.array([[1, 0, 0],
              [0, 1, 0],
              [0, 0, 1]])

# Basic operations
print("Matrix A:")
print(A)
print(f"\nTranspose: A.T =\n{A.T}")
print(f"\nMatrix multiplication: A @ B =\n{A @ B}")
print(f"\nElement-wise multiplication: A * B =\n{A * B}")
print(f"\nInverse: A⁻¹ =\n{np.linalg.inv(A)}")
print(f"\nDeterminant: det(A) = {np.linalg.det(A):.4f}")
print(f"\nTrace: tr(A) = {np.trace(A)}")

1.2 Special Matrix Creation¶

n = 4

# Identity matrix
I = np.eye(n)

# Zero matrix
Z = np.zeros((n, n))

# Diagonal matrix
D = np.diag([1, 2, 3, 4])

# Tridiagonal matrix
diag_main = 2 * np.ones(n)
diag_off = -1 * np.ones(n - 1)
T = np.diag(diag_main) + np.diag(diag_off, k=1) + np.diag(diag_off, k=-1)

print("Tridiagonal matrix:")
print(T)

# Random matrix
np.random.seed(42)
R = np.random.randn(3, 3)
print(f"\nRandom matrix:\n{R}")

1.3 Matrix Norms¶

A = np.array([[1, 2], [3, 4]])

# Various norms
print("Matrix norms:")
print(f"  1-norm (max column sum): {np.linalg.norm(A, 1)}")
print(f"  2-norm (spectral): {np.linalg.norm(A, 2)}")
print(f"  ∞-norm (max row sum): {np.linalg.norm(A, np.inf)}")
print(f"  Frobenius norm: {np.linalg.norm(A, 'fro')}")

# Vector norms
v = np.array([3, 4])
print(f"\nVector norms:")
print(f"  L1: {np.linalg.norm(v, 1)}")
print(f"  L2: {np.linalg.norm(v, 2)}")
print(f"  L∞: {np.linalg.norm(v, np.inf)}")

2. Solving Linear Systems¶

2.1 Direct Solution¶

# Solve Ax = b
A = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])

# numpy.linalg.solve (recommended)
x = np.linalg.solve(A, b)
print(f"Solution: x = {x}")
print(f"Verification: Ax = {A @ x}")

# Using inverse (not recommended - slow and unstable)
x_inv = np.linalg.inv(A) @ b
print(f"Inverse method: x = {x_inv}")

2.2 Overdetermined/Underdetermined Systems¶

# Overdetermined system (equations > unknowns) - least squares solution
A_over = np.array([[1, 1], [1, 2], [1, 3]])
b_over = np.array([1, 2, 2])

# Least squares solution
x_lstsq, residuals, rank, s = np.linalg.lstsq(A_over, b_over, rcond=None)
print(f"Least squares solution: {x_lstsq}")
print(f"Residuals: {residuals}")

# Underdetermined system (equations < unknowns) - minimum norm solution
A_under = np.array([[1, 2, 3]])
b_under = np.array([6])

# Using pseudoinverse
x_min_norm = np.linalg.pinv(A_under) @ b_under
print(f"\nMinimum norm solution: {x_min_norm}")
print(f"Norm: {np.linalg.norm(x_min_norm):.4f}")

3. Eigenvalues and Eigenvectors¶

3.1 Eigenvalue Decomposition¶

A = np.array([[4, -2],
              [1,  1]])

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

print("Eigenvalue decomposition:")
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# Verification: A @ v = λ @ v
for i in range(len(eigenvalues)):
    lam = eigenvalues[i]
    v = eigenvectors[:, i]
    print(f"\nλ_{i} = {lam:.4f}")
    print(f"  A @ v = {A @ v}")
    print(f"  λ * v = {lam * v}")

3.2 Eigenvalues of Symmetric Matrices¶

# Symmetric matrix - guarantees real eigenvalues
S = np.array([[2, 1, 0],
              [1, 3, 1],
              [0, 1, 2]])

# eigh is optimized for symmetric matrices (faster and more stable)
eigenvalues, eigenvectors = np.linalg.eigh(S)

print("Symmetric matrix eigenvalue decomposition:")
print(f"Eigenvalues: {eigenvalues}")
print(f"\nEigenvector orthogonality verification:")
print(f"V^T @ V =\n{eigenvectors.T @ eigenvectors}")  # Identity matrix

3.3 Power Method¶

def power_method(A, max_iter=100, tol=1e-10):
    """Calculate largest eigenvalue and eigenvector using power method"""
    n = A.shape[0]
    v = np.random.randn(n)
    v = v / np.linalg.norm(v)

    for _ in range(max_iter):
        v_new = A @ v
        eigenvalue = np.dot(v, v_new)  # Rayleigh quotient
        v_new = v_new / np.linalg.norm(v_new)

        if np.linalg.norm(v_new - v) < tol:
            break
        v = v_new

    return eigenvalue, v_new

A = np.array([[4, 1], [2, 3]])
lam, v = power_method(A)
print(f"Power method result:")
print(f"  Largest eigenvalue: {lam:.6f}")
print(f"  Eigenvector: {v}")

# Compare with numpy result
lam_np, v_np = np.linalg.eig(A)
print(f"\nnumpy result:")
print(f"  Eigenvalues: {lam_np}")

4. Matrix Decomposition¶

4.1 LU Decomposition¶

from scipy.linalg import lu, lu_factor, lu_solve

A = np.array([[2, 1, 1],
              [4, 3, 3],
              [8, 7, 9]], dtype=float)

# LU decomposition
P, L, U = lu(A)

print("LU decomposition:")
print(f"P (permutation matrix):\n{P}")
print(f"L (lower triangular):\n{L}")
print(f"U (upper triangular):\n{U}")
print(f"\nVerification: P @ L @ U =\n{P @ L @ U}")

# Using for solving linear systems
b = np.array([4, 10, 24])
lu_piv = lu_factor(A)
x = lu_solve(lu_piv, b)
print(f"\nSolution: {x}")

4.2 Cholesky Decomposition¶

# Only possible for positive definite symmetric matrices
A_spd = np.array([[4, 2, 2],
                  [2, 5, 1],
                  [2, 1, 6]], dtype=float)

# Cholesky decomposition: A = L @ L.T
L = np.linalg.cholesky(A_spd)

print("Cholesky decomposition:")
print(f"L:\n{L}")
print(f"\nVerification: L @ L.T =\n{L @ L.T}")

# Solving linear systems (2x faster than LU)
b = np.array([8, 8, 9])
# Solve L @ y = b, then L.T @ x = y
y = linalg.solve_triangular(L, b, lower=True)
x = linalg.solve_triangular(L.T, y)
print(f"\nSolution: {x}")

4.3 QR Decomposition¶

A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 10],
              [10, 11, 12]], dtype=float)

# QR decomposition
Q, R = np.linalg.qr(A)

print("QR decomposition:")
print(f"Q (orthogonal matrix):\n{Q}")
print(f"\nR (upper triangular):\n{R}")
print(f"\nVerification: Q @ R =\n{Q @ R}")
print(f"\nOrthogonality of Q: Q.T @ Q =\n{Q.T @ Q}")

# Using for least squares problem
b = np.array([1, 2, 3, 4])
x = linalg.solve_triangular(R, Q.T @ b)
print(f"\nLeast squares solution: {x}")

4.4 SVD (Singular Value Decomposition)¶

A = np.array([[1, 2],
              [3, 4],
              [5, 6]])

# SVD: A = U @ Σ @ V.T
U, s, Vt = np.linalg.svd(A)

print("SVD decomposition:")
print(f"U (m×m orthogonal):\n{U}")
print(f"\nSingular values: {s}")
print(f"\nV.T (n×n orthogonal):\n{Vt}")

# Construct Σ matrix
Sigma = np.zeros_like(A, dtype=float)
np.fill_diagonal(Sigma, s)

print(f"\nVerification: U @ Σ @ V.T =\n{U @ Sigma @ Vt}")

# Condition number calculation
cond = s[0] / s[-1]
print(f"\nCondition number: {cond:.4f}")

5. Sparse Matrices¶

5.1 Sparse Matrix Formats¶

from scipy import sparse

# COO format (coordinate format)
row = [0, 1, 2, 2]
col = [0, 1, 0, 2]
data = [1, 2, 3, 4]
A_coo = sparse.coo_matrix((data, (row, col)), shape=(3, 3))

print("COO format:")
print(A_coo)
print(f"\nDense form:\n{A_coo.toarray()}")

# CSR format (compressed sparse row, efficient for matrix-vector multiplication)
A_csr = A_coo.tocsr()
print(f"\nCSR format: {A_csr}")

# CSC format (compressed sparse column, efficient for column slicing)
A_csc = A_coo.tocsc()

5.2 Sparse Matrix Operations¶

# Create large sparse matrix
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"Sparse matrix size: {A_sparse.shape}")
print(f"Number of non-zero elements: {A_sparse.nnz}")
print(f"Sparsity: {1 - A_sparse.nnz / (n*n):.4%}")

# Sparse linear system
b = np.random.randn(n)
x = sparse.linalg.spsolve(A_sparse, b)
print(f"\nNorm of sparse linear system solution: {np.linalg.norm(x):.4f}")

5.3 Iterative Solvers¶

from scipy.sparse.linalg import cg, gmres, bicgstab

# Create symmetric positive definite matrix
n = 100
A = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(n, n), format='csr')
A = A.T @ A + 0.1 * sparse.eye(n)  # Make positive definite
b = np.random.randn(n)

# CG (Conjugate Gradient) - optimal for symmetric positive definite
x_cg, info_cg = cg(A, b, tol=1e-10)
print(f"CG: convergence status = {info_cg}, residual = {np.linalg.norm(A @ x_cg - b):.2e}")

# GMRES - general matrices
x_gmres, info_gmres = gmres(A, b, tol=1e-10)
print(f"GMRES: convergence status = {info_gmres}, residual = {np.linalg.norm(A @ x_gmres - b):.2e}")

6. Application Examples¶

6.1 Discretization of Heat Conduction Problem¶

def heat_equation_matrix(n, alpha, dx, dt):
    """Implicit discretization matrix for 1D heat equation"""
    r = alpha * dt / dx**2

    # Create tridiagonal matrix
    main_diag = (1 + 2*r) * np.ones(n)
    off_diag = -r * np.ones(n - 1)

    A = np.diag(main_diag) + np.diag(off_diag, 1) + np.diag(off_diag, -1)
    return A

n = 10
A = heat_equation_matrix(n, alpha=0.1, dx=0.1, dt=0.01)
print("Heat equation discretization matrix:")
print(A[:5, :5])  # Print only part

6.2 Image Compression (SVD)¶

def svd_compression(image, k):
    """Compress image using SVD (using top k singular values)"""
    U, s, Vt = np.linalg.svd(image, full_matrices=False)

    # Use only top k components
    compressed = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

    # Calculate compression ratio
    original_size = image.shape[0] * image.shape[1]
    compressed_size = k * (image.shape[0] + image.shape[1] + 1)
    compression_ratio = compressed_size / original_size

    return compressed, compression_ratio

# Example (random image)
image = np.random.randn(100, 100)
for k in [5, 10, 20, 50]:
    comp, ratio = svd_compression(image, k)
    error = np.linalg.norm(image - comp, 'fro') / np.linalg.norm(image, 'fro')
    print(f"k={k:2d}: compression ratio={ratio:.2%}, relative error={error:.4f}")

Practice Problems¶

Problem 1¶

Find the eigenvalues of a 3x3 tridiagonal matrix and compare with the power method result.

# Solution
T = np.array([[2, -1, 0],
              [-1, 2, -1],
              [0, -1, 2]], dtype=float)

# numpy
eig_np, _ = np.linalg.eig(T)
print(f"numpy eigenvalues: {sorted(eig_np)}")

# Power method
lam, v = power_method(T)
print(f"Power method largest eigenvalue: {lam:.6f}")

Problem 2¶

Solve a 100x100 sparse matrix linear system using direct and iterative methods and compare execution times.

import time

n = 1000
A = sparse.diags([1, -2, 1], [-1, 0, 1], shape=(n, n), format='csr')
A = A + 3 * sparse.eye(n)  # Make diagonally dominant
b = np.random.randn(n)

# Direct method
start = time.time()
x_direct = sparse.linalg.spsolve(A, b)
time_direct = time.time() - start

# Iterative method (CG)
A_spd = A.T @ A
b_spd = A.T @ b
start = time.time()
x_cg, _ = cg(A_spd, b_spd, tol=1e-10)
time_cg = time.time() - start

print(f"Direct method: {time_direct:.4f}s")
print(f"CG: {time_cg:.4f}s")

Summary¶

Decomposition Method	Purpose	Requirements
LU	Solving linear systems	Square matrix
Cholesky	Solving linear systems (fast)	Symmetric positive definite
QR	Least squares, eigenvalues	Any matrix
SVD	Compression, pseudoinverse	Any matrix

Solver	Matrix Type	Characteristics
spsolve	Sparse	Direct method
CG	Symmetric positive definite	Iterative method
GMRES	General	Iterative method
BiCGSTAB	Non-symmetric	Iterative method