"""
MATLAB matrix group - Standard MATLAB test matrices
"""
import numpy as np
from scipy import sparse
from typing import Optional
def register_matlab_matrices(registry):
"""Register MATLAB-style test matrices"""
registry.register_matrix(
group='matlab',
name='magic',
generator=magic,
properties=[],
description='Magic square - rows, columns, diagonals sum to same value'
)
registry.register_matrix(
group='matlab',
name='pascal',
generator=pascal,
properties=['symmetric', 'positive definite', 'integer'],
description='Pascal matrix - binomial coefficients'
)
registry.register_matrix(
group='matlab',
name='rosser',
generator=rosser,
properties=['symmetric'],
description='Rosser matrix - classic eigenvalue test'
)
registry.register_matrix(
group='matlab',
name='wilkinson',
generator=wilkinson,
properties=['symmetric', 'tridiagonal'],
description='Wilkinson matrix W+ or W-'
)
registry.register_matrix(
group='matlab',
name='hilbert',
generator=hilbert,
properties=['symmetric', 'positive definite'],
description='Hilbert matrix - notoriously ill-conditioned'
)
registry.register_matrix(
group='matlab',
name='invhilb',
generator=invhilb,
properties=['symmetric', 'positive definite', 'integer'],
description='Inverse of Hilbert matrix'
)
registry.register_matrix(
group='matlab',
name='frank',
generator=frank,
properties=['upper hessenberg'],
description='Frank matrix - upper Hessenberg with determinant 1'
)
registry.register_matrix(
group='matlab',
name='companion',
generator=companion,
properties=[],
description='Companion matrix of a polynomial'
)
[docs]
def magic(n: int) -> np.ndarray:
"""
Magic square matrix
Args:
n: Size (n >= 3)
Returns:
n×n magic square
"""
if n < 3:
raise ValueError("Magic square size must be >= 3")
if n % 2 == 1:
# Odd order - De la Loubere method
M = np.zeros((n, n), dtype=int)
i, j = 0, n // 2
for num in range(1, n*n + 1):
M[i, j] = num
i_new, j_new = (i - 1) % n, (j + 1) % n
if M[i_new, j_new]:
i_new = (i + 1) % n
j_new = j
i, j = i_new, j_new
elif n % 4 == 0:
# Doubly even order
M = np.arange(1, n*n + 1).reshape(n, n)
# Create diagonal pattern
for i in range(n):
for j in range(n):
if (i % 4 == j % 4) or ((i % 4) + (j % 4) == 3):
M[i, j] = n*n + 1 - M[i, j]
else:
# Singly even order (n = 4k+2)
# Use Strachey method
k = n // 2
M_odd = magic(k)
M = np.zeros((n, n), dtype=int)
# Four copies with offsets
M[:k, :k] = M_odd * 4 - 3
M[:k, k:] = M_odd * 4 - 2
M[k:, :k] = M_odd * 4 - 1
M[k:, k:] = M_odd * 4
# Swap some elements for magic property
k2 = (n - 2) // 4
for i in range(k):
for j in range(k2):
if i != k2:
M[i, j], M[i + k, j] = M[i + k, j], M[i, j]
# Additional swap in middle
M[k2, k2], M[k2 + k, k2] = M[k2 + k, k2], M[k2, k2]
M[k2, k2 + k2], M[k2 + k, k2 + k2] = M[k2 + k, k2 + k2], M[k2, k2 + k2]
return M.astype(float)
[docs]
def pascal(n: int, k: int = 0) -> np.ndarray:
"""
Pascal matrix
Args:
n: Size
k: Type (0=symmetric, 1=lower, 2=upper)
Returns:
n×n Pascal matrix
"""
from scipy.special import comb
P = np.zeros((n, n))
if k == 0:
# Symmetric Pascal matrix
for i in range(n):
for j in range(n):
P[i, j] = comb(i + j, i, exact=True)
elif k == 1:
# Lower triangular
for i in range(n):
for j in range(i + 1):
P[i, j] = comb(i, j, exact=True)
else:
# Upper triangular
for i in range(n):
for j in range(i, n):
P[i, j] = comb(j, i, exact=True)
return P
[docs]
def rosser() -> np.ndarray:
"""
Rosser matrix - 8×8 symmetric test matrix
Classic eigenvalue test matrix with close eigenvalues
"""
R = np.array([
[611, 196, -192, 407, -8, -52, -49, 29],
[196, 899, 113, -192, -71, -43, -8, -44],
[-192, 113, 899, 196, 61, 49, 8, 52],
[407, -192, 196, 611, 8, 44, 59, -23],
[-8, -71, 61, 8, 411, -599, 208, 208],
[-52, -43, 49, 44, -599, 411, 208, 208],
[-49, -8, 8, 59, 208, 208, 99, -911],
[29, -44, 52, -23, 208, 208, -911, 99]
], dtype=float)
return R
[docs]
def wilkinson(n: int, mode: str = 'plus') -> np.ndarray:
"""
Wilkinson matrix
Args:
n: Size (must be odd)
mode: 'plus' for W+ or 'minus' for W-
Returns:
n×n symmetric tridiagonal Wilkinson matrix
"""
if n % 2 == 0:
raise ValueError("n must be odd for Wilkinson matrix")
m = (n - 1) // 2
if mode == 'plus':
# W+ : diagonal is |m, m-1, ..., 1, 0, 1, ..., m|
diag = np.concatenate([np.arange(m, 0, -1), [0], np.arange(1, m + 1)])
else:
# W- : diagonal is |m, m-1, ..., 1, 0, -1, ..., -m|
diag = np.concatenate([np.arange(m, 0, -1), [0], np.arange(-1, -m - 1, -1)])
# Superdiagonal and subdiagonal are all 1
W = np.diag(diag) + np.diag(np.ones(n - 1), 1) + np.diag(np.ones(n - 1), -1)
return W.astype(float)
[docs]
def hilbert(n: int) -> np.ndarray:
"""
Hilbert matrix - H(i,j) = 1/(i+j-1)
Notoriously ill-conditioned matrix
Args:
n: Size
Returns:
n×n Hilbert matrix
"""
i, j = np.meshgrid(range(1, n + 1), range(1, n + 1), indexing='ij')
H = 1.0 / (i + j - 1)
return H
[docs]
def invhilb(n: int) -> np.ndarray:
"""
Inverse of Hilbert matrix
Exact integer inverse of Hilbert matrix
Args:
n: Size
Returns:
n×n inverse Hilbert matrix
"""
from scipy.special import comb
InvH = np.zeros((n, n))
for i in range(n):
for j in range(n):
sign = (-1) ** (i + j)
term1 = (i + j + 1) * comb(n + i, n - j - 1, exact=True)
term2 = comb(n + j, n - i - 1, exact=True)
term3 = comb(i + j, i, exact=True) ** 2
InvH[i, j] = sign * term1 * term2 * term3
return InvH
[docs]
def frank(n: int, k: int = 0) -> np.ndarray:
"""
Frank matrix - upper Hessenberg with determinant 1
Args:
n: Size
k: Type (0=default, 1=reflected)
Returns:
n×n Frank matrix
"""
F = np.zeros((n, n))
if k == 0:
# Standard Frank matrix
for i in range(n):
for j in range(n):
if j >= i:
F[i, j] = n - max(i, j)
elif j == i - 1:
F[i, j] = n - i
else:
# Reflected Frank matrix
for i in range(n):
for j in range(n):
if j <= i:
F[i, j] = n - max(i, j)
elif j == i + 1:
F[i, j] = n - j
return F
[docs]
def companion(c: np.ndarray) -> np.ndarray:
"""
Companion matrix of polynomial
For polynomial p(x) = c[0] + c[1]*x + ... + c[n]*x^n,
creates companion matrix whose eigenvalues are roots of p
Args:
c: Polynomial coefficients [c0, c1, ..., cn]
Returns:
n×n companion matrix
"""
c = np.asarray(c)
n = len(c) - 1
if c[-1] == 0:
raise ValueError("Leading coefficient must be nonzero")
# Normalize
c = c / c[-1]
C = np.zeros((n, n))
C[:-1, 1:] = np.eye(n - 1)
C[-1, :] = -c[:-1]
return C