"""
Regtools matrix group - Matrices from regularization problems
Based on Per Christian Hansen's Regularization Tools
"""
import numpy as np
from scipy import sparse
from typing import Optional, Tuple
def register_regtools_matrices(registry):
"""Register regularization problem matrices"""
registry.register_matrix(
group='regtools',
name='deriv2',
generator=deriv2,
properties=['sparse', 'symmetric'],
description='Second derivative operator (discrete Laplacian)'
)
registry.register_matrix(
group='regtools',
name='shaw',
generator=shaw,
properties=[],
description='Shaw test problem - 1D image restoration'
)
registry.register_matrix(
group='regtools',
name='phillips',
generator=phillips,
properties=[],
description='Phillips test problem - image deblurring'
)
registry.register_matrix(
group='regtools',
name='blur',
generator=blur,
properties=['sparse'],
description='Blurring matrix for image deconvolution'
)
registry.register_matrix(
group='regtools',
name='gravity',
generator=gravity,
properties=[],
description='Gravity surveying test problem'
)
[docs]
def deriv2(n: int, example: int = 1) -> sparse.spmatrix:
"""
Second derivative operator
Discrete approximation to second derivative (1D Laplacian)
Args:
n: Size
example: 1 for standard, 2 for periodic BC
Returns:
n×n sparse matrix
"""
h = 1.0 / (n + 1)
if example == 1:
# Standard second derivative with zero BC
diags = [
-2 * np.ones(n),
np.ones(n - 1),
np.ones(n - 1)
]
L = sparse.diags(diags, [0, 1, -1], format='csr') / (h ** 2)
else:
# Periodic boundary conditions
diags = [
-2 * np.ones(n),
np.ones(n - 1),
np.ones(n - 1)
]
L = sparse.diags(diags, [0, 1, -1], format='lil')
L[0, -1] = 1 / (h ** 2)
L[-1, 0] = 1 / (h ** 2)
L = L.tocsr()
return L
[docs]
def shaw(n: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Shaw test problem - 1D image restoration
Args:
n: Size
Returns:
Tuple of (A, b, x_true) where A@x_true ≈ b
"""
h = np.pi / n
# Discretization points
t = -np.pi/2 + h * np.arange(n)
s = t.copy()
# Create matrix
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
if abs(s[i] - t[j]) < np.pi / 2:
A[i, j] = (np.cos(s[i]) + np.cos(t[j])) * \
(np.sin((s[i] + t[j]) / 2)) ** 2
else:
A[i, j] = 0
A = A * h
# True solution (smooth function)
x_true = 2 * np.exp(-6 * (t - 0.8) ** 2) + np.exp(-2 * (t + 0.5) ** 2)
# Right-hand side
b = A @ x_true
return A, b, x_true
[docs]
def phillips(n: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Phillips test problem - image deblurring
Args:
n: Size
Returns:
Tuple of (A, b, x_true)
"""
h = 12.0 / n
# Discretization
t = -6 + h * (np.arange(n) + 0.5)
s = t.copy()
# Create matrix (convolution-type)
A = np.zeros((n, n))
c = np.cos(t)
d = np.sin(t)
for i in range(n):
for j in range(n):
A[i, j] = h * c[i] * c[j] + h * d[i] * d[j]
# True solution
x_true = np.zeros(n)
x_true[t > -0.5] = 1
x_true[t > 0.5] = 0
# Right-hand side
b = A @ x_true
return A, b, x_true
[docs]
def blur(n: int, band: int = 3, sigma: float = 0.7) -> sparse.spmatrix:
"""
Blurring matrix for image deconvolution
Creates a sparse matrix representing Gaussian blur
Args:
n: Image size (total size will be n×n)
band: Bandwidth of blur kernel
sigma: Standard deviation of Gaussian
Returns:
n²×n² sparse blurring matrix
"""
# 1D Gaussian kernel
x = np.arange(-band, band + 1)
kernel = np.exp(-x**2 / (2 * sigma**2))
kernel = kernel / kernel.sum()
# Create 1D blurring matrix
diags = []
offsets = []
for i, k in enumerate(kernel):
offset = i - band
diags.append(np.ones(n - abs(offset)) * k)
offsets.append(offset)
B1D = sparse.diags(diags, offsets, shape=(n, n), format='csr')
# 2D blurring via Kronecker product
I = sparse.eye(n, format='csr')
B2D = sparse.kron(B1D, I) + sparse.kron(I, B1D)
B2D = B2D / 2 # Average of row and column blur
return B2D
[docs]
def gravity(n: int, example: int = 1, a: float = 0, b: float = 1,
d: float = 0.25) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Gravity surveying test problem
Vertical gravity surveying over a 2D domain
Args:
n: Number of discretization points
example: Problem variant (1 or 2)
a: Left endpoint
b: Right endpoint
d: Depth
Returns:
Tuple of (A, b, x_true)
"""
h = (b - a) / n
t = a + h * (np.arange(n) + 0.5) # Measurement points
s = t.copy() # Source points
# Kernel matrix
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
A[i, j] = h * d / ((t[i] - s[j])**2 + d**2)**1.5
# True solution (density distribution)
if example == 1:
x_true = np.sin(np.pi * s) + 0.5 * np.sin(2 * np.pi * s)
else:
x_true = np.zeros(n)
mask = (s >= 0.3) & (s <= 0.7)
x_true[mask] = 1
# Measurements
b = A @ x_true
return A, b, x_true