Ferkans — Interactive Telecom Tutor

ex-sp-ch06-01

Easy

Create a $4 \times 4$ random complex matrix $\mathbf{A}$ and a random complex vector $\mathbf{b}$ . Solve $\mathbf{A}\mathbf{x} = \mathbf{b}$ using both np.linalg.solve and np.linalg.inv. Verify that both give the same answer (up to numerical precision) and time both methods.

Show Hint

Use np.random.randn(4, 4) + 1j * np.random.randn(4, 4) for complex matrices.

Use np.allclose to compare results.

Solution

Implementation

import numpy as np
import time

np.random.seed(42)
A = np.random.randn(4, 4) + 1j * np.random.randn(4, 4)
b = np.random.randn(4) + 1j * np.random.randn(4)

t0 = time.perf_counter()
x_solve = np.linalg.solve(A, b)
t_solve = time.perf_counter() - t0

t0 = time.perf_counter()
x_inv = np.linalg.inv(A) @ b
t_inv = time.perf_counter() - t0

print(f"solve time: {t_solve:.6f}s")
print(f"inv time:   {t_inv:.6f}s")
print(f"Match: {np.allclose(x_solve, x_inv)}")
print(f"Residual (solve): {np.linalg.norm(A @ x_solve - b):.2e}")
print(f"Residual (inv):   {np.linalg.norm(A @ x_inv - b):.2e}")

ex-sp-ch06-02

Easy

Compute the eigenvalues and eigenvectors of the $3 \times 3$ correlation matrix $R_{ij} = 0.9^{|i-j|}$ . Verify that the eigenvalues are real, the eigenvectors are orthogonal, and $\mathbf{R} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^T$ .

Show Hint

Use np.linalg.eigh since R is symmetric.

Check orthogonality with np.allclose(V.T @ V, np.eye(3)).

Solution

Implementation

import numpy as np

rho = 0.9
R = np.array([[rho**abs(i-j) for j in range(3)] for i in range(3)])
eigenvalues, V = np.linalg.eigh(R)

print(f"Eigenvalues: {eigenvalues}")
print(f"All real: {np.all(np.isreal(eigenvalues))}")
print(f"Orthogonal: {np.allclose(V.T @ V, np.eye(3))}")

R_reconstructed = V @ np.diag(eigenvalues) @ V.T
print(f"Reconstruction error: {np.linalg.norm(R - R_reconstructed):.2e}")

ex-sp-ch06-03

Easy

Construct a $1000 \times 1000$ sparse tridiagonal matrix with $-2$ on the diagonal and $1$ on the sub/super-diagonals. Compare its memory usage (sparse vs dense) and compute a matrix-vector product timing for both.

Show Hint

Use scipy.sparse.diags([1, -2, 1], [-1, 0, 1], shape=(1000, 1000)).

Compare .nbytes for dense vs .data.nbytes + .indices.nbytes + .indptr.nbytes for sparse.

Solution

Implementation

import numpy as np
from scipy.sparse import diags
import time

n = 1000
A_sparse = diags([1, -2, 1], [-1, 0, 1], shape=(n, n), format='csr')
A_dense = A_sparse.toarray()
x = np.random.randn(n)

sparse_bytes = A_sparse.data.nbytes + A_sparse.indices.nbytes + A_sparse.indptr.nbytes
dense_bytes = A_dense.nbytes
print(f"Sparse: {sparse_bytes/1024:.1f} KB")
print(f"Dense:  {dense_bytes/1024:.1f} KB")
print(f"Ratio:  {dense_bytes/sparse_bytes:.0f}x")

t0 = time.perf_counter()
for _ in range(1000): y = A_sparse @ x
t_sparse = time.perf_counter() - t0

t0 = time.perf_counter()
for _ in range(1000): y = A_dense @ x
t_dense = time.perf_counter() - t0

print(f"Sparse matvec: {t_sparse:.4f}s")
print(f"Dense matvec:  {t_dense:.4f}s")

ex-sp-ch06-04

Easy

Compute the SVD of a $6 \times 4$ random matrix. Print the singular values, verify $\mathbf{U}^H\mathbf{U} = \mathbf{I}$ and $\mathbf{V}^H\mathbf{V} = \mathbf{I}$ , and reconstruct the original matrix from the SVD factors.

Show Hint

Use np.linalg.svd(A, full_matrices=False) for the reduced SVD.

Reconstruct with (U * s) @ Vh.

Solution

Implementation

import numpy as np

A = np.random.randn(6, 4)
U, s, Vh = np.linalg.svd(A, full_matrices=False)

print(f"Singular values: {s}")
print(f"U^H U = I? {np.allclose(U.T @ U, np.eye(4))}")
print(f"V^H V = I? {np.allclose(Vh @ Vh.T, np.eye(4))}")

A_reconstructed = (U * s) @ Vh
print(f"Reconstruction error: {np.linalg.norm(A - A_reconstructed):.2e}")

ex-sp-ch06-05

Easy

Compute the Kronecker product of $\mathbf{A} = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}$ and $\mathbf{B} = \begin{bmatrix} 0 & 5 \\ 6 & 7 \end{bmatrix}$ . Verify the property $\det(\mathbf{A} \otimes \mathbf{B}) = (\det\mathbf{A})^2 (\det\mathbf{B})^2$ .

Show Hint

Use np.kron(A, B) and np.linalg.det.

Solution

Implementation

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 5], [6, 7]])
K = np.kron(A, B)

print(f"A kron B =\n{K}")
print(f"det(A kron B) = {np.linalg.det(K):.4f}")
print(f"det(A)^2 * det(B)^2 = {np.linalg.det(A)**2 * np.linalg.det(B)**2:.4f}")

ex-sp-ch06-06

Easy

Create a $4 \times 4$ Toeplitz matrix with first column $[2, 1, 0, 0]$ and first row $[2, -1, 0, 0]$ . Verify that it represents a specific filter by comparing T @ x with np.convolve(h, x).

Show Hint

Use scipy.linalg.toeplitz(c, r).

Solution

Implementation

from scipy.linalg import toeplitz
import numpy as np

c = [2, 1, 0, 0]
r = [2, -1, 0, 0]
T = toeplitz(c, r)
print(f"T =\n{T}")

x = np.array([1, 2, 3, 4])
print(f"T @ x = {T @ x}")

ex-sp-ch06-07

Medium

Implement the MMSE detector for a $4 \times 4$ MIMO system. Given $\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{n}$ , compute: $\hat{\mathbf{x}}_{\mathrm{MMSE}} = (\mathbf{H}^H\mathbf{H} + \sigma_n^2\mathbf{I})^{-1}\mathbf{H}^H\mathbf{y}$ Compare with ZF detection at SNR = 0, 10, 20 dB over 1000 trials. Plot the symbol error rate vs SNR.

Show Hint

Use np.linalg.solve for both ZF and MMSE.

QPSK symbols: (+/-1 +/- j) / sqrt(2).

Solution

Implementation

import numpy as np

N = 4
n_trials = 1000
snr_db_range = np.arange(0, 25, 2)

for snr_db in snr_db_range:
    snr_lin = 10 ** (snr_db / 10)
    sigma2 = 1.0 / snr_lin
    zf_errors, mmse_errors = 0, 0

    for _ in range(n_trials):
        H = (np.random.randn(N, N) + 1j * np.random.randn(N, N)) / np.sqrt(2)
        bits = np.random.choice([-1, 1], (N, 2))
        x = (bits[:, 0] + 1j * bits[:, 1]) / np.sqrt(2)
        n = np.sqrt(sigma2/2) * (np.random.randn(N) + 1j * np.random.randn(N))
        y = H @ x + n

        x_zf = np.linalg.solve(H.conj().T @ H, H.conj().T @ y)
        x_mmse = np.linalg.solve(H.conj().T @ H + sigma2 * np.eye(N),
                                  H.conj().T @ y)

        zf_errors += np.sum(np.sign(x_zf.real) != np.sign(x.real))
        mmse_errors += np.sum(np.sign(x_mmse.real) != np.sign(x.real))

    print(f"SNR={snr_db:2d} dB: ZF SER={zf_errors/(n_trials*N):.4f}, "
          f"MMSE SER={mmse_errors/(n_trials*N):.4f}")

ex-sp-ch06-08

Medium

Compute the condition number of a Hilbert matrix $H_{ij} = 1/(i+j-1)$ for sizes $n = 2, 4, 6, 8, 10, 12$ . Show how the condition number grows and explain why Hilbert matrices are notoriously ill-conditioned.

Show Hint

Use scipy.linalg.hilbert(n) and np.linalg.cond.

The condition number grows approximately as $e^{3.5n}$ .

Solution

Implementation

from scipy.linalg import hilbert
import numpy as np

for n in [2, 4, 6, 8, 10, 12]:
    H = hilbert(n)
    kappa = np.linalg.cond(H)
    print(f"n={n:2d}: cond = {kappa:.2e}")

ex-sp-ch06-09

Medium

Build a graph Laplacian for a $10$ -node ring graph (each node connected to its two neighbors) as a sparse matrix. Find the 3 smallest eigenvalues using eigsh. Verify that the smallest eigenvalue is 0 (connected graph).

Show Hint

Adjacency: node $i$ connects to $(i-1) \mod n$ and $(i+1) \mod n$ .

Laplacian: $\mathbf{L} = \mathbf{D} - \mathbf{A}$ where $\mathbf{D}$ is degree matrix.

Solution

Implementation

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import eigsh

n = 10
A = lil_matrix((n, n))
for i in range(n):
    A[i, (i+1) % n] = 1
    A[i, (i-1) % n] = 1
A = A.tocsr()
D = A.sum(axis=1).A1
L = -A + lil_matrix(np.diag(D))
L = L.tocsr()

eigenvalues, _ = eigsh(L.astype(float), k=3, which='SM')
print(f"3 smallest eigenvalues: {eigenvalues}")
print(f"Smallest ~0? {abs(eigenvalues[0]) < 1e-10}")

ex-sp-ch06-10

Medium

Implement kron_matvec(A, B, x) that computes $(\mathbf{A} \otimes \mathbf{B})\mathbf{x}$ without forming the full Kronecker product. Benchmark it against np.kron(A, B) @ x for $50 \times 50$ factors.

Show Hint

Reshape $\mathbf{x}$ to $(q, n)$ in Fortran order.

Compute $\mathbf{Y} = \mathbf{B} \mathbf{X} \mathbf{A}^T$ .

Solution

Implementation

import numpy as np
import time

def kron_matvec(A, B, x):
    m, n = A.shape
    p, q = B.shape
    X = x.reshape(q, n, order='F')
    Y = B @ X @ A.T
    return Y.ravel(order='F')

m = n = p = q = 50
A = np.random.randn(m, n)
B = np.random.randn(p, q)
x = np.random.randn(n * q)

t0 = time.perf_counter()
y_naive = np.kron(A, B) @ x
t_naive = time.perf_counter() - t0

t0 = time.perf_counter()
y_fast = kron_matvec(A, B, x)
t_fast = time.perf_counter() - t0

print(f"Naive: {t_naive:.4f}s, Efficient: {t_fast:.6f}s")
print(f"Speedup: {t_naive/t_fast:.0f}x")
print(f"Error: {np.linalg.norm(y_naive - y_fast):.2e}")

ex-sp-ch06-11

Medium

Verify the circulant diagonalization theorem: for a circulant matrix $\mathbf{C}$ with first column $\mathbf{c} = [1, 2, 3, 4, 5]$ , compute eigenvalues both via np.linalg.eig and via the DFT ( $\lambda = \sqrt{n}\,\mathbf{F}\mathbf{c}$ ). Show they match.

Show Hint

Use scipy.linalg.circulant and scipy.linalg.dft.

Sort eigenvalues before comparing (eig returns unsorted).

Solution

Implementation

import numpy as np
from scipy.linalg import circulant, dft

c = np.array([1, 2, 3, 4, 5], dtype=float)
n = len(c)
C = circulant(c)

# Method 1: direct eigendecomposition
eigs_direct = np.linalg.eigvals(C)

# Method 2: DFT of first column
F = dft(n, scale='sqrtn')
eigs_dft = np.sqrt(n) * (F @ c)

# Sort for comparison
idx1 = np.argsort(eigs_direct.real)
idx2 = np.argsort(eigs_dft.real)
print(f"Match: {np.allclose(eigs_direct[idx1], eigs_dft[idx2])}")

ex-sp-ch06-12

Medium

Solve an overdetermined system $\mathbf{A}\mathbf{x} \approx \mathbf{b}$ with $m = 100$ , $n = 5$ using three methods: normal equations, QR factorization, and np.linalg.lstsq. Compare the errors.

Show Hint

Use scipy.linalg.qr and scipy.linalg.solve_triangular.

Solution

Implementation

import numpy as np
from scipy.linalg import qr, solve_triangular

m, n = 100, 5
np.random.seed(42)
A = np.random.randn(m, n)
x_true = np.array([1.0, -2.0, 3.0, -4.0, 5.0])
b = A @ x_true + 0.01 * np.random.randn(m)

# Normal equations
x_ne = np.linalg.solve(A.T @ A, A.T @ b)

# QR
Q, R = qr(A, mode='economic')
x_qr = solve_triangular(R, Q.T @ b)

# lstsq
x_ls = np.linalg.lstsq(A, b, rcond=None)[0]

for name, x in [('Normal', x_ne), ('QR', x_qr), ('lstsq', x_ls)]:
    print(f"{name:8s}: error = {np.linalg.norm(x - x_true):.6f}")

ex-sp-ch06-13

Medium

Compute the matrix exponential $e^{\mathbf{A}t}$ for the $2 \times 2$ system $\mathbf{A} = \begin{bmatrix} -1 & 1 \\ 0 & -2 \end{bmatrix}$ at 20 time points from $t = 0$ to $t = 5$ . Verify that $e^{\mathbf{A} \cdot 0} = \mathbf{I}$ and that $\|e^{\mathbf{A}t}\| \to 0$ as $t \to \infty$ (stable system).

Show Hint

Use scipy.linalg.expm(A * t) in a loop.

Eigenvalues of A are $-1$ and $-2$ (both negative, so stable).

Solution

Implementation

import numpy as np
from scipy.linalg import expm

A = np.array([[-1, 1], [0, -2]])
times = np.linspace(0, 5, 20)

for t in times:
    Phi = expm(A * t)
    print(f"t={t:.2f}: ||expm(At)|| = {np.linalg.norm(Phi):.4f}")

assert np.allclose(expm(A * 0), np.eye(2))

ex-sp-ch06-14

Hard

Implement PCA (Principal Component Analysis) from scratch using eigendecomposition:

Generate 1000 points from a 2D Gaussian with covariance $\mathbf{R} = \begin{bmatrix} 3 & 1.5 \\ 1.5 & 1 \end{bmatrix}$ .
Center the data (subtract mean).
Compute the sample covariance matrix.
Find eigenvectors via eigh.
Project data onto the principal component and compute the fraction of variance explained.

Show Hint

Sample covariance: X.T @ X / (N - 1) after centering.

Variance explained: eigenvalue[i] / sum(eigenvalues).

Solution

Implementation

import numpy as np

R = np.array([[3, 1.5], [1.5, 1]])
L = np.linalg.cholesky(R)
X = np.random.randn(1000, 2) @ L.T

X_centered = X - X.mean(axis=0)
Sigma = X_centered.T @ X_centered / (len(X) - 1)
eigenvalues, V = np.linalg.eigh(Sigma)

# Sort descending
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
V = V[:, idx]

var_explained = eigenvalues / eigenvalues.sum()
print(f"Variance explained: {var_explained}")
print(f"PC1 direction: {V[:, 0]}")

# Project onto PC1
X_proj = X_centered @ V[:, :1] @ V[:, :1].T

ex-sp-ch06-15

Hard

Implement a function sparse_pagerank(A, alpha=0.85, tol=1e-6) that computes the PageRank vector for a sparse adjacency matrix using the power iteration method. Test on a random sparse graph with 10,000 nodes.

Show Hint

Normalize each column of $\mathbf{A}$ to get the transition matrix.

PageRank update: $\mathbf{r}^{(k+1)} = \alpha \mathbf{M}\mathbf{r}^{(k)} + (1-\alpha)/n \cdot \mathbf{1}$ .

Iterate until $\|\mathbf{r}^{(k+1)} - \mathbf{r}^{(k)}\| < \text{tol}$ .

Solution

Implementation

import numpy as np
from scipy.sparse import random as sparse_random

def sparse_pagerank(A, alpha=0.85, tol=1e-6, max_iter=100):
    n = A.shape[0]
    # Column-normalize
    col_sums = np.array(A.sum(axis=0)).flatten()
    col_sums[col_sums == 0] = 1
    M = A / col_sums

    r = np.ones(n) / n
    for iteration in range(max_iter):
        r_new = alpha * (M @ r) + (1 - alpha) / n
        if np.linalg.norm(r_new - r) < tol:
            print(f"Converged in {iteration + 1} iterations")
            return r_new
        r = r_new
    return r

A = sparse_random(10000, 10000, density=0.001, format='csr')
ranks = sparse_pagerank(A)
print(f"Top 5 node ranks: {np.argsort(ranks)[-5:]}")

ex-sp-ch06-16

Hard

Implement Tikhonov regularization with L-curve analysis:

Generate a noisy overdetermined system with a known solution.
Compute the Tikhonov solution for 50 values of $\alpha$ logarithmically spaced from $10^{-6}$ to $10^2$ .
Plot the L-curve: $\log\|\mathbf{Ax}_\alpha - \mathbf{b}\|$ vs $\log\|\mathbf{x}_\alpha\|$ .
Find the optimal $\alpha$ at the corner of the L-curve.

Show Hint

The corner is the point of maximum curvature.

A simple heuristic: find where the second derivative changes sign.

Solution

Implementation

import numpy as np

m, n = 100, 10
np.random.seed(42)
A = np.random.randn(m, n)
x_true = np.random.randn(n)
b = A @ x_true + 0.5 * np.random.randn(m)

alphas = np.logspace(-6, 2, 50)
residuals = []
norms = []

for alpha in alphas:
    x_tik = np.linalg.solve(A.T @ A + alpha * np.eye(n), A.T @ b)
    residuals.append(np.linalg.norm(A @ x_tik - b))
    norms.append(np.linalg.norm(x_tik))

# Find corner (maximum curvature heuristic)
log_res = np.log10(residuals)
log_nrm = np.log10(norms)
curvature = np.gradient(np.gradient(log_res)) + np.gradient(np.gradient(log_nrm))
best_idx = np.argmax(np.abs(curvature))
print(f"Optimal alpha: {alphas[best_idx]:.4e}")

ex-sp-ch06-17

Hard

Show that the 2D DFT of an $N \times N$ image can be computed as a Kronecker product: $\mathrm{vec}(\mathbf{F}\mathbf{X}\mathbf{F}^T) = (\mathbf{F} \otimes \mathbf{F})\,\mathrm{vec}(\mathbf{X})$ . Verify numerically for $N = 8$ and compare the timing of the direct Kronecker product vs np.fft.fft2.

Show Hint

Use scipy.linalg.dft(N, scale="sqrtn") for the DFT matrix.

np.fft.fft2 computes the 2D DFT in $O(N^2 \log N)$ vs $O(N^4)$ for Kronecker.

Solution

Implementation

import numpy as np
from scipy.linalg import dft

N = 8
F = dft(N, scale='sqrtn')
X = np.random.randn(N, N)

# Method 1: Direct 2D DFT matrix
Y1_vec = np.kron(F, F) @ X.ravel(order='F')
Y1 = Y1_vec.reshape(N, N, order='F')

# Method 2: Efficient F X F^T
Y2 = F @ X @ F.T

# Method 3: FFT
Y3 = np.fft.fft2(X) / N  # scale to match unitary DFT

print(f"Kron vs FXF^T: {np.allclose(Y1, Y2)}")

ex-sp-ch06-18

Hard

Implement the Lanczos algorithm for finding the $k$ largest eigenvalues of a symmetric matrix. Compare your implementation with scipy.sparse.linalg.eigsh on a $1000 \times 1000$ sparse matrix.

Show Hint

Build the tridiagonal matrix $\mathbf{T}$ iteratively.

At each step: $\mathbf{w} = \mathbf{A}\mathbf{q}_k - \beta_{k-1}\mathbf{q}_{k-1}$ , $\alpha_k = \mathbf{q}_k^T\mathbf{w}$ , $\mathbf{w} = \mathbf{w} - \alpha_k\mathbf{q}_k$ , $\beta_k = \|\mathbf{w}\|$ .

Solution

Implementation

import numpy as np
from scipy.sparse import random as sparse_random
from scipy.sparse.linalg import eigsh

def lanczos(A, k, m=None):
    n = A.shape[0]
    if m is None:
        m = min(2 * k + 10, n)

    Q = np.zeros((n, m))
    alpha = np.zeros(m)
    beta = np.zeros(m)

    q = np.random.randn(n)
    q /= np.linalg.norm(q)
    Q[:, 0] = q

    for j in range(m):
        w = A @ Q[:, j]
        if j > 0:
            w -= beta[j-1] * Q[:, j-1]
        alpha[j] = Q[:, j] @ w
        w -= alpha[j] * Q[:, j]
        beta[j] = np.linalg.norm(w)
        if j < m - 1 and beta[j] > 1e-12:
            Q[:, j+1] = w / beta[j]

    T = np.diag(alpha[:m]) + np.diag(beta[:m-1], 1) + np.diag(beta[:m-1], -1)
    evals = np.sort(np.linalg.eigvalsh(T))[-k:]
    return evals

A = sparse_random(1000, 1000, density=0.01, format='csr')
A = (A + A.T) / 2  # make symmetric

evals_lanczos = lanczos(A, k=5)
evals_scipy, _ = eigsh(A, k=5, which='LM')

print(f"Lanczos: {np.sort(evals_lanczos)}")
print(f"eigsh:   {np.sort(evals_scipy)}")

ex-sp-ch06-19

Challenge

Implement a complete SVD-based MIMO system:

Generate a $4 \times 4$ Rayleigh channel $\mathbf{H}$ .
Compute $\mathbf{H} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^H$ .
Implement SVD precoding: transmit $\mathbf{x} = \mathbf{V}\tilde{\mathbf{x}}$ , receive $\tilde{\mathbf{y}} = \mathbf{U}^H\mathbf{y}$ .
Implement water-filling power allocation across the parallel channels $\tilde{y}_i = \sigma_i \tilde{x}_i + \tilde{n}_i$ .
Compute and plot the capacity vs SNR, comparing uniform power allocation with water-filling.

Show Hint

Water-filling: $P_i = \max(0, \mu - \sigma_n^2/\sigma_i^2)$ where $\mu$ satisfies $\sum P_i = P_{\mathrm{total}}$ .

Use bisection to find the water level $\mu$ .

Solution

Implementation

import numpy as np

def water_filling(gains, noise_var, total_power):
    """Water-filling power allocation."""
    n = len(gains)
    mu_low, mu_high = 0, total_power + noise_var / min(gains)**2
    for _ in range(100):
        mu = (mu_low + mu_high) / 2
        powers = np.maximum(0, mu - noise_var / gains**2)
        if np.sum(powers) > total_power:
            mu_high = mu
        else:
            mu_low = mu
    return powers

Nr, Nt = 4, 4
snr_db_range = np.arange(0, 31, 2)
n_trials = 500
cap_uniform = np.zeros(len(snr_db_range))
cap_wf = np.zeros(len(snr_db_range))

for trial in range(n_trials):
    H = (np.random.randn(Nr, Nt) + 1j * np.random.randn(Nr, Nt)) / np.sqrt(2)
    _, s, _ = np.linalg.svd(H)

    for i, snr_db in enumerate(snr_db_range):
        snr_lin = 10 ** (snr_db / 10)
        sigma2 = 1.0 / snr_lin

        # Uniform
        p_uni = snr_lin / Nt * np.ones(Nt)
        cap_uniform[i] += np.sum(np.log2(1 + p_uni * s**2))

        # Water-filling
        p_wf = water_filling(s**2, sigma2, snr_lin)
        cap_wf[i] += np.sum(np.log2(1 + p_wf * s**2))

cap_uniform /= n_trials
cap_wf /= n_trials

for i, db in enumerate(snr_db_range):
    print(f"SNR={db:2d} dB: Uniform={cap_uniform[i]:.2f}, "
          f"WF={cap_wf[i]:.2f} bits/s/Hz")

ex-sp-ch06-20

Challenge

Build a preconditioned conjugate gradient (PCG) solver for a large sparse system:

Construct the 2D Laplacian on a $100 \times 100$ grid (10,000 unknowns) as a sparse matrix.
Implement CG from scratch.
Add Jacobi preconditioning ( $\mathbf{M}^{-1} = \mathrm{diag}(\mathbf{A})^{-1}$ ).
Compare convergence (iterations to tolerance $10^{-8}$ ) of CG, PCG, and scipy.sparse.linalg.cg.

Show Hint

The 2D Laplacian on an $n \times n$ grid is $\mathbf{I} \otimes \mathbf{L}_n + \mathbf{L}_n \otimes \mathbf{I}$ where $\mathbf{L}_n$ is the 1D Laplacian.

CG algorithm: initialize $\mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0$ , $\mathbf{p}_0 = \mathbf{r}_0$ , then iterate.

Solution

Implementation

import numpy as np
from scipy.sparse import diags, eye, kron
from scipy.sparse.linalg import cg as scipy_cg

def cg_solve(A, b, tol=1e-8, maxiter=1000, M_inv=None):
    x = np.zeros_like(b)
    r = b - A @ x
    if M_inv is not None:
        z = M_inv @ r
    else:
        z = r.copy()
    p = z.copy()
    rz = r @ z

    for k in range(maxiter):
        Ap = A @ p
        alpha = rz / (p @ Ap)
        x += alpha * p
        r -= alpha * Ap
        if np.linalg.norm(r) < tol:
            return x, k + 1
        if M_inv is not None:
            z = M_inv @ r
        else:
            z = r.copy()
        rz_new = r @ z
        beta = rz_new / rz
        p = z + beta * p
        rz = rz_new
    return x, maxiter

n = 100
L1d = diags([1, -2, 1], [-1, 0, 1], shape=(n, n))
L2d = kron(eye(n), L1d) + kron(L1d, eye(n))
L2d = -L2d.tocsr()

b = np.ones(n * n)
diag_inv = 1.0 / L2d.diagonal()

x_cg, iter_cg = cg_solve(L2d, b)
x_pcg, iter_pcg = cg_solve(L2d, b, M_inv=lambda r: diag_inv * r)

print(f"CG:  {iter_cg} iterations")
print(f"PCG: {iter_pcg} iterations")

Exercises

ex-sp-ch06-01

Implementation

ex-sp-ch06-02

Implementation

ex-sp-ch06-03

Implementation

ex-sp-ch06-04

Implementation

ex-sp-ch06-05

Implementation

ex-sp-ch06-06

Implementation

ex-sp-ch06-07

Implementation

ex-sp-ch06-08

Implementation

ex-sp-ch06-09

Implementation

ex-sp-ch06-10

Implementation

ex-sp-ch06-11

Implementation

ex-sp-ch06-12

Implementation

ex-sp-ch06-13

Implementation

ex-sp-ch06-14

Implementation

ex-sp-ch06-15

Implementation

ex-sp-ch06-16

Implementation

ex-sp-ch06-17

Implementation

ex-sp-ch06-18

Implementation

ex-sp-ch06-19

Implementation

ex-sp-ch06-20

Implementation