"""
Gallery matrix group - Matrices from MATLAB gallery
Python implementations of common MATLAB gallery matrices.
"""
import numpy as np
from scipy import sparse
from typing import Optional
def register_gallery_matrices(registry):
"""Register all gallery group matrices"""
registry.register_matrix(
group='gallery',
name='lehmer',
generator=lehmer,
properties=['symmetric', 'positive definite'],
description='Lehmer matrix - symmetric positive definite'
)
registry.register_matrix(
group='gallery',
name='minij',
generator=minij,
properties=['symmetric', 'positive definite'],
description='MIN(i,j) matrix'
)
registry.register_matrix(
group='gallery',
name='moler',
generator=moler,
properties=['symmetric', 'positive definite'],
description='Moler matrix'
)
registry.register_matrix(
group='gallery',
name='pei',
generator=pei,
properties=['symmetric'],
description='Pei matrix'
)
registry.register_matrix(
group='gallery',
name='clement',
generator=clement,
properties=['tridiagonal', 'symmetric'],
description='Clement matrix - tridiagonal with zero diagonal'
)
registry.register_matrix(
group='gallery',
name='kms',
generator=kms,
properties=['toeplitz', 'symmetric'],
description='Kac-Murdock-Szego Toeplitz matrix'
)
[docs]
def lehmer(n: int) -> np.ndarray:
"""
Lehmer matrix - A(i,j) = min(i,j)/max(i,j)
Symmetric positive definite matrix.
Args:
n: Matrix size
Returns:
n×n Lehmer matrix
"""
i, j = np.meshgrid(range(1, n+1), range(1, n+1), indexing='ij')
A = np.minimum(i, j) / np.maximum(i, j)
return A
[docs]
def minij(n: int) -> np.ndarray:
"""
MIN(i,j) matrix - A(i,j) = min(i,j)
Symmetric positive definite matrix.
Args:
n: Matrix size
Returns:
n×n minij matrix
"""
i, j = np.meshgrid(range(1, n+1), range(1, n+1), indexing='ij')
A = np.minimum(i, j).astype(float)
return A
[docs]
def moler(n: int, alpha: float = -1.0) -> np.ndarray:
"""
Moler matrix - a symmetric positive definite matrix
Args:
n: Matrix size
alpha: Parameter (default -1)
Returns:
n×n Moler matrix
"""
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
if i == j:
A[i, j] = i + 1
else:
A[i, j] = min(i, j) + alpha
return A
[docs]
def pei(n: int, alpha: float = 1.0) -> np.ndarray:
"""
Pei matrix - alpha*I + ones(n)
Args:
n: Matrix size
alpha: Scalar parameter
Returns:
n×n Pei matrix
"""
return alpha * np.eye(n) + np.ones((n, n))
[docs]
def clement(n: int, kind: int = 0) -> np.ndarray:
"""
Clement matrix - tridiagonal with zero diagonal
Args:
n: Matrix size
kind: 0 for symmetric, 1 for nonsymmetric
Returns:
n×n Clement matrix
"""
if kind == 0:
# Symmetric
A = np.diag(np.arange(n-1, 0, -1), 1) + np.diag(np.arange(1, n), -1)
else:
# Nonsymmetric
A = np.diag(np.arange(1, n), 1) + np.diag(np.arange(n-1, 0, -1), -1)
return A.astype(float)
[docs]
def kms(n: int, rho: float = 0.5) -> np.ndarray:
"""
Kac-Murdock-Szego Toeplitz matrix
A symmetric Toeplitz matrix with A(i,j) = rho^|i-j|
Args:
n: Matrix size
rho: Parameter (0 < rho < 1)
Returns:
n×n KMS matrix
"""
i, j = np.meshgrid(range(n), range(n), indexing='ij')
A = rho ** np.abs(i - j)
return A