"""
Core matrix group - Miscellaneous test matrices
This module contains Python implementations of matrices from the
anymatrix 'core' group.
"""
import numpy as np
from scipy import sparse
from typing import Optional, Union
def register_core_matrices(registry):
"""Register all core group matrices with the registry"""
# Beta matrix
registry.register_matrix(
group='core',
name='beta',
generator=beta,
properties=['symmetric', 'positive definite'],
description='Beta matrix - a symmetric positive definite matrix'
)
# Fourier matrix
registry.register_matrix(
group='core',
name='fourier',
generator=fourier,
properties=['unitary', 'symmetric'],
description='Fourier matrix - discrete Fourier transform matrix'
)
# Nilpotent triangular
registry.register_matrix(
group='core',
name='nilpot_triang',
generator=nilpot_triang,
properties=['upper triangular', 'nilpotent'],
description='Nilpotent upper triangular matrix'
)
# Nilpotent tridiagonal
registry.register_matrix(
group='core',
name='nilpot_tridiag',
generator=nilpot_tridiag,
properties=['tridiagonal', 'nilpotent'],
description='Nilpotent tridiagonal matrix'
)
# Augmented matrix
registry.register_matrix(
group='core',
name='augment',
generator=augment,
properties=['symmetric'],
description='Augmented system matrix'
)
# Vandermonde matrix
registry.register_matrix(
group='core',
name='vand',
generator=vand,
properties=[],
description='Vandermonde matrix'
)
# Wilson matrix
registry.register_matrix(
group='core',
name='wilson',
generator=wilson,
properties=['symmetric', 'positive definite'],
description='Wilson matrix'
)
# Circulant binomial
registry.register_matrix(
group='core',
name='circul_binom',
generator=circul_binom,
properties=['circulant', 'toeplitz'],
description='Circulant matrix with binomial coefficients'
)
# Stochastic Cesaro
registry.register_matrix(
group='core',
name='stoch_cesaro',
generator=stoch_cesaro,
properties=['stochastic', 'lower triangular'],
description='Stochastic Cesaro matrix'
)
# Tournament matrix
registry.register_matrix(
group='core',
name='tournament',
generator=tournament,
properties=['binary'],
description='Random tournament matrix'
)
# Perfect shuffle
registry.register_matrix(
group='core',
name='perfect_shuffle',
generator=perfect_shuffle,
properties=['orthogonal', 'sparse'],
description='Perfect shuffle permutation matrix'
)
# Block circulant
registry.register_matrix(
group='core',
name='blockcirc',
generator=blockcirc,
properties=['circulant'],
description='Block circulant matrix'
)
# Collatz matrix
registry.register_matrix(
group='core',
name='collatz',
generator=collatz,
properties=['sparse', 'integer'],
description='Collatz conjecture matrix'
)
[docs]
def beta(n: int, beta_val: float = 0.5) -> np.ndarray:
"""
Beta matrix - symmetric positive definite
Args:
n: Matrix size
beta_val: Beta parameter (default 0.5)
Returns:
n×n beta matrix
"""
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
A[i, j] = 1.0 / (i + j + beta_val)
return A
[docs]
def fourier(n: int) -> np.ndarray:
"""
Discrete Fourier transform matrix
Args:
n: Matrix size
Returns:
n×n Fourier matrix
"""
omega = np.exp(-2j * np.pi / n)
i, j = np.meshgrid(range(n), range(n), indexing='ij')
F = omega ** (i * j) / np.sqrt(n)
return F
[docs]
def nilpot_triang(n: int) -> np.ndarray:
"""
Nilpotent upper triangular matrix
Args:
n: Matrix size
Returns:
n×n nilpotent upper triangular matrix
"""
A = np.zeros((n, n))
for i in range(n - 1):
A[i, i + 1] = 1
return A
[docs]
def nilpot_tridiag(n: int) -> np.ndarray:
"""
Nilpotent tridiagonal matrix
A tridiagonal matrix that is nilpotent. The matrix has the form:
[ 0 1 0 0 ... 0 ]
[-1 0 2 0 ... 0 ]
[ 0 -2 0 3 ... 0 ]
[ 0 0 -3 0 ... 0 ]
...
[ 0 0 0 ... 0 n-1]
[ 0 0 0 ... -(n-1) 0 ]
This is NOT nilpotent! Let me check the original implementation...
Actually, looking at the original anymatrix code, this matrix is
skew-symmetric tridiagonal, which means A^T = -A, so A^2 is symmetric.
For a truly nilpotent tridiagonal, we need a different construction.
Args:
n: Matrix size
Returns:
n×n nilpotent tridiagonal matrix
"""
# The original version from my previous code was:
# Superdiagonal: 1, 2, 3, ..., n-1
# Subdiagonal: -1, -2, -3, ..., -(n-1)
# But this is NOT nilpotent!
# A proper nilpotent tridiagonal matrix:
# Option 1: Just superdiagonal (like Jordan block with eigenvalue 0)
A = np.diag(np.ones(n - 1), 1)
# This is nilpotent with index n: A^n = 0
return A
[docs]
def augment(A: np.ndarray, b: Optional[np.ndarray] = None) -> np.ndarray:
"""
Augmented system matrix [A b; b' 0]
Args:
A: Original matrix
b: Right-hand side vector (default: ones)
Returns:
Augmented matrix
"""
n = A.shape[0]
if b is None:
b = np.ones((n, 1))
else:
b = b.reshape(-1, 1)
# Create augmented matrix
top = np.hstack([A, b])
bottom = np.hstack([b.T, np.zeros((1, 1))])
return np.vstack([top, bottom])
[docs]
def vand(c: Optional[np.ndarray] = None, n: Optional[int] = None) -> np.ndarray:
"""
Vandermonde matrix
Args:
c: Vector of values (default: 0, 1, ..., n-1)
n: Matrix size (used if c is None)
Returns:
Vandermonde matrix
"""
if c is None:
if n is None:
raise ValueError("Either c or n must be provided")
c = np.arange(n)
else:
c = np.asarray(c)
n = len(c)
V = np.vander(c, n, increasing=True)
return V
[docs]
def wilson(n: int = 4) -> np.ndarray:
"""
Wilson matrix - symmetric positive definite
Args:
n: Matrix size (default 4)
Returns:
n×n Wilson matrix
"""
W = np.zeros((n, n))
for i in range(n):
for j in range(n):
W[i, j] = (i + 1) ** (n - j)
return W
[docs]
def circul_binom(n: int) -> np.ndarray:
"""
Circulant matrix with binomial coefficients
Args:
n: Matrix size
Returns:
n×n circulant matrix
"""
from scipy.special import comb
# First row contains binomial coefficients
first_row = np.array([comb(n-1, k, exact=True) for k in range(n)])
# Create circulant matrix
C = np.zeros((n, n))
for i in range(n):
C[i, :] = np.roll(first_row, i)
return C
[docs]
def stoch_cesaro(n: int) -> np.ndarray:
"""
Stochastic Cesaro matrix (lower triangular, row stochastic)
Args:
n: Matrix size
Returns:
n×n stochastic Cesaro matrix
"""
C = np.zeros((n, n))
for i in range(n):
for j in range(i + 1):
C[i, j] = 1.0 / (i + 1)
return C
def tournament(n: int, seed: Optional[int] = None) -> np.ndarray:
"""
Random tournament matrix (random binary skew-symmetric + I)
Args:
n: Matrix size
seed: Random seed
Returns:
n×n tournament matrix
"""
if seed is not None:
np.random.seed(seed)
# Random binary matrix
A = np.random.randint(0, 2, size=(n, n))
# Make skew-symmetric
A = np.triu(A, 1)
A = A - A.T
# Add identity (self-loops)
A = A + np.eye(n)
# Make binary
A = (A > 0).astype(int)
return A
[docs]
def perfect_shuffle(n: int) -> sparse.spmatrix:
"""
Perfect shuffle permutation matrix
The perfect shuffle permutation interleaves the first and second
halves of a vector.
Args:
n: Matrix size (should be even)
Returns:
n×n sparse permutation matrix
"""
if n % 2 != 0:
raise ValueError("n must be even for perfect shuffle")
perm = np.zeros(n, dtype=int)
m = n // 2
# Interleave indices
perm[::2] = np.arange(m)
perm[1::2] = np.arange(m, n)
# Create permutation matrix
P = sparse.eye(n, format='csr')
P = P[perm, :]
return P
[docs]
def blockcirc(blocks: list, k: int = 1) -> np.ndarray:
"""
Block circulant matrix
Args:
blocks: List of square matrices
k: Shift parameter
Returns:
Block circulant matrix
"""
num_blocks = len(blocks)
block_size = blocks[0].shape[0]
n = num_blocks * block_size
C = np.zeros((n, n))
for i in range(num_blocks):
for j in range(num_blocks):
idx = (i - j * k) % num_blocks
row_start = i * block_size
col_start = j * block_size
C[row_start:row_start+block_size, col_start:col_start+block_size] = blocks[idx]
return C
def collatz(n: int) -> sparse.spmatrix:
"""
Collatz conjecture matrix
Matrix representation of the Collatz (3n+1) problem.
Args:
n: Matrix size
Returns:
n×n sparse integer matrix
"""
row = []
col = []
data = []
for i in range(1, n + 1):
# Even: i -> i/2
if i % 2 == 0:
target = i // 2
if target <= n:
row.append(i - 1)
col.append(target - 1)
data.append(1)
# Odd: i -> 3i+1
else:
target = 3 * i + 1
if target <= n:
row.append(i - 1)
col.append(target - 1)
data.append(1)
return sparse.csr_matrix((data, (row, col)), shape=(n, n))