Linear Algebra API#

Linear algebra routines operate on SpaceCore spaces and operators. Preconditions are mathematical: square, Hermitian, positive-definite, and residual statements refer to the domain and codomain inner products, not only coordinate arrays.

Linear solves#

`spacecore.linalg.cg`	Solve \(A x = b\) by conjugate gradients.
`spacecore.linalg.LinalgResult`	Result of every `spacecore.linalg` solver.

LinalgResult is returned by every solver below: x, converged, num_iters and residual_norm, plus eigenvalue or normal_residual_norm where the solver produces them.

Least squares#

spacecore.linalg.lsqr

Solve \(\min_x \|A x - b\|\) by LSQR.

Eigenvalue and spectral methods#

`spacecore.linalg.lanczos_smallest`	Approximate the smallest eigenpair of a Hermitian operator by Lanczos.
`spacecore.linalg.power_iteration`	Estimate the dominant eigenpair (largest in absolute value) of a self-adjoint action.

Matrix functions#

spacecore.linalg.expm_multiply

Compute \(\exp(t A) v\) by Lanczos projection, without forming A.

Native elements and JAX#

lanczos_smallest and expm_multiply retain native space elements and use zeros, add, scale, and inner together with the operator action. They do not require a coordinate shape or flatten/unflatten methods. Built-in dense spaces retain their vectorized basis operations. Custom elements are held directly during eager execution; on JAX they must be registered pytrees with fixed structure, leaf shapes, and dtypes throughout the recurrence, including space.zeros(). Basis capacity is fixed by the static max_iter argument. Mathematical linear combinations always use space algebra, even when the representation has multiple array leaves.

Lanczos requires a nonzero starting element. A zero start raises ValueError eagerly. Under JAX tracing, it returns a NaN eigenvalue, a zero eigenvector, num_iters=0, and converged=False. This replaces the previous coordinate vector fallback. expm_multiply accepts zero and returns zero with num_iters=0 and a zero residual_norm.

Autodoc#

spacecore.linalg.cg(A, b, *, x0=None, tol=1e-06, atol=0.0, maxiter=None, check_every=64, preconditioner=None)[source]#

Solve \(A x = b\) by conjugate gradients.

A must be square (A.domain == A.codomain), Hermitian and positive-definite in the domain’s inner product. Only A.apply and that inner product are used.

Parameters:

A (LinOp) – Hermitian positive-definite operator. One that is provably not Hermitian is rejected; definiteness is not checked.
b (array-like) – Right-hand side in A.codomain.
x0 (array-like or None, optional) – Initial guess; zero by default.
tol (float, optional) – Converged when the residual norm is below atol + tol * norm(b).
atol (float, optional) – Absolute residual tolerance.
maxiter (int or None, optional) – Iteration cap; prod(A.domain.shape) by default.
check_every (int, optional) – Test convergence every this many iterations, and on the last one.
preconditioner (LinOp or None, optional) – Operator K on A.domain approximating the inverse of A, self-adjoint and positive definite in A.domain; each residual r is preconditioned to K r. The residual is still measured in A.domain’s norm.

Returns:

(x, converged, num_iters, residual_norm).

Return type:

LinalgResult

Raises:

TypeError – If A is not a LinOp.
ValueError – If A is not square or an iteration parameter is invalid.

See also

lsqr: Least squares for rectangular operators.

Notes

JIT-compatible on JAX with maxiter and check_every static. Stops early when no useful update remains (residual at machine precision, or non-positive curvature). Trefethen & Bau, Numerical Linear Algebra, Algorithm 38.1. In an orthonormal basis the metric recurrence is ordinary CG; complex spaces use conjugate inner products, and a nonzero x0 starts from b - A x0.

Examples

Solve a small positive-definite system.

>>> import numpy as np
>>> import spacecore as sc
>>> ctx = sc.Context(sc.NumpyOps(), dtype=np.float64)
>>> X = sc.DenseCoordinateSpace((3,), ctx)
>>> M = ctx.asarray([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
>>> A = sc.DenseLinOp(M, X, X, ctx)
>>> b = ctx.asarray([1.0, 2.0, 3.0])
>>> result = sc.cg(A, b, tol=1e-10)
>>> bool(result.converged)
True
>>> np.allclose(A.apply(result.x), b)
True

class spacecore.linalg.LinalgResult(x, converged, num_iters, residual_norm, eigenvalue=None, normal_residual_norm=None)[source]#

Bases: NamedTuple

Result of every spacecore.linalg solver.

Fields a solver does not produce are None.

Parameters:

x (array-like) – The computed element: the solution (cg(), lsqr()), the eigenvector estimate (power_iteration(), lanczos_smallest()), or the action exp(t * A) @ v (expm_multiply()).
converged (bool-like) – Whether the solver’s stopping test was met.
num_iters (int-like) – Iterations executed; for the Krylov methods lanczos_smallest() and expm_multiply(), the Krylov dimension reached.
residual_norm (scalar) – The solver’s residual measure: norm(b - A x) for cg() and lsqr() (exact or estimated per residual_mode), norm(A x - eigenvalue x) for power_iteration(), the Ritz residual estimate for lanczos_smallest(), and for expm_multiply() the projection indicator abs(beta[m] phi[m-1]), which omits norm(v) and does not bound the error.
eigenvalue (scalar or None) – Eigenvalue estimate (power_iteration(), lanczos_smallest()).
normal_residual_norm (scalar or None) – norm(A.H (A x - b)) (lsqr()).

x: Any#: Alias for field number 0

converged: Any#: Alias for field number 1

num_iters: Any#: Alias for field number 2

residual_norm: Any#: Alias for field number 3

eigenvalue: Any#: Alias for field number 4

normal_residual_norm: Any#: Alias for field number 5

count(value, /)#: Return number of occurrences of value.

index(value, start=0, stop=sys.maxsize, /)#

Return first index of value.

Raises ValueError if the value is not present.

spacecore.linalg.lsqr(A, b, *, x0=None, tol=1e-06, atol=0.0, maxiter=None, check_every=64, residual_mode='exact')[source]#

Solve \(\min_x \|A x - b\|\) by LSQR.

A may map between different spaces. Only A.apply and A.H.apply are used; the normal equations are never formed.

Parameters:

A (LinOp) – Operator from A.domain to A.codomain.
b (array-like) – Right-hand side in A.codomain.
x0 (array-like or None, optional) – Initial guess; zero by default.
tol (float, optional) – Converged when norm(A.H @ (A @ x - b)) is below atol + tol * norm(b).
atol (float, optional) – Absolute tolerance on that normal-equation residual.
maxiter (int or None, optional) – Iteration cap; prod(A.domain.shape) by default.
check_every (int, optional) – Refresh residual diagnostics every this many iterations, and on the last one.
residual_mode ({"exact", "recurrence"}, optional) – "exact" recomputes both residuals at each check, costing one extra forward and adjoint apply; "recurrence" uses LSQR’s scalar estimates for free.

Returns:

With normal_residual_norm set.

Return type:

LinalgResult

Raises:

TypeError – If A is not a LinOp.
ValueError – If an iteration parameter or residual_mode is invalid.

See also

cg: Square Hermitian positive-definite systems.

Notes

JIT-compatible on JAX with maxiter, check_every and residual_mode static. A small normal-equation residual does not imply a small error when A is ill-conditioned: the error solves A.H A e = A.H (A x - b). Golub & Van Loan, Matrix Computations, 4th ed., Sec. 11.4.2. Orthonormal bases in both spaces give ordinary LSQR; complex spaces replace transposes by conjugate transposes in the bidiagonalization.

Examples

Solve a small overdetermined least-squares problem.

>>> import numpy as np
>>> import spacecore as sc
>>> ctx = sc.Context(sc.NumpyOps(), dtype=np.float64)
>>> X = sc.DenseCoordinateSpace((2,), ctx)
>>> Y = sc.DenseCoordinateSpace((3,), ctx)
>>> M = ctx.asarray([[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
>>> A = sc.DenseLinOp(M, X, Y, ctx)
>>> b = ctx.asarray([1.0, 2.0, 3.0])
>>> result = sc.lsqr(A, b, tol=1e-10)
>>> np.allclose(result.x, [1.0, 2.0])
True

spacecore.linalg.lanczos_smallest(A, initial_vector, *, max_iter=100, tol=1e-06, check_every=64)[source]#

Approximate the smallest eigenpair of a Hermitian operator by Lanczos.

The eigenvalue returned is the Rayleigh quotient of the Ritz vector, reconstructed in A.domain.

Parameters:

A (LinOp) – Square operator, Hermitian in the domain’s inner product. Known non-Hermitian operators are rejected; unknown ones are trusted.
initial_vector (array-like) – Nonzero starting element in A.domain.
max_iter (int, optional) – Maximum Krylov dimension; a static argument under jax.jit.
tol (float, optional) – Breakdown threshold, and the residual below which converged is set.
check_every (int, optional) – Test for breakdown every this many iterations, and on the last one.

Returns:

x is the Ritz vector, eigenvalue the Ritz value, num_iters the Krylov dimension reached.

Return type:

LinalgResult

Raises:

TypeError – If A is not a LinOp.
ValueError – If A is not square, is known non-Hermitian, max_iter is invalid, or the starting element is zero (eager execution).

See also

power_iteration: The dominant eigenpair.

Notes

residual_norm is the recurrence estimate \(\beta_m |y_{m-1}|\). The smallest Ritz value can recover the smallest eigenvalue only if initial_vector overlaps its eigenspace. Trefethen & Bau, Numerical Linear Algebra, Algorithm 36.1 and Thm 36.1; multiplying the Lanczos relation by a unit Ritz vector gives the residual. Uses the space algebra and native elements. Custom JAX elements must be registered pytrees whose structure, leaf shapes, and dtypes stay fixed. max_iter and check_every must be static under jax.jit. A zero start raises eagerly; under JAX tracing it returns a NaN eigenvalue, a zero eigenvector, and converged=False.

Examples

Approximate the smallest eigenpair of a diagonal operator.

>>> import numpy as np
>>> import spacecore as sc
>>> ctx = sc.Context(sc.NumpyOps(), dtype=np.float64)
>>> X = sc.DenseCoordinateSpace((3,), ctx)
>>> A = sc.DiagonalLinOp(ctx.asarray([1.0, 2.0, 4.0]), X, ctx)
>>> result = sc.lanczos_smallest(A, ctx.asarray([1.0, 1.0, 1.0]), max_iter=3)
>>> np.allclose(result.eigenvalue, 1.0)
True

spacecore.linalg.power_iteration(A, *, x0=None, tol=1e-06, maxiter=None)[source]#

Estimate the dominant eigenpair (largest in absolute value) of a self-adjoint action.

Parameters:

A (LinOp or QuadraticForm) – Square operator, or quadratic form with hess_apply. A form may also supply hess_quad(x, Hx=None) and hess_residual_norm(x, Hx, eigenvalue) for the diagnostics. For a spectral norm of a rectangular B, pass B.H @ B.
x0 (array-like or None, optional) – Initial vector; normalized ones by default.
tol (float, optional) – Converged when norm(A @ x - eigenvalue * x) < tol.
maxiter (int or None, optional) – Iteration cap; prod(A.domain.shape) by default.

Returns:

x is the eigenvector estimate; eigenvalue its Rayleigh quotient.

Return type:

LinalgResult

Raises:

TypeError – If A is neither a LinOp nor a QuadraticForm.
ValueError – If a linear operator is not square or known non-Hermitian, or an iteration parameter is invalid.

See also

lanczos_smallest: The smallest eigenpair.

Notes

Convergence needs a separated dominant modulus and nonzero initial overlap with its eigenspace. Trefethen & Bau, Numerical Linear Algebra, Thm 27.1; expansion in a metric-orthonormal eigenbasis gives the same eigenvalue-ratio argument. With both lambda and -lambda present at the dominant modulus the iterate can oscillate. JIT-compatible on JAX with maxiter static.

Examples

Estimate the largest eigenvalue of a diagonal operator.

>>> import numpy as np
>>> import spacecore as sc
>>> ctx = sc.Context(sc.NumpyOps(), dtype=np.float64)
>>> X = sc.DenseCoordinateSpace((3,), ctx)
>>> A = sc.DiagonalLinOp(ctx.asarray([1.0, 3.0, 2.0]), X, ctx)
>>> result = sc.power_iteration(A, maxiter=20, tol=1e-10)
>>> np.allclose(result.eigenvalue, 3.0)
True

spacecore.linalg.expm_multiply(A, v, t=1.0, *, max_iter=30, tol=1e-10)[source]#

Compute \(\exp(t A) v\) by Lanczos projection, without forming A.

Parameters:

A (LinOp) – Square operator, Hermitian in the domain’s inner product. Known non-Hermitian operators are rejected; unknown ones are trusted.
v (array-like) – Element in A.domain; zero returns zero without applying A.
t (float or complex, optional) – Scalar multiplying A; a complex t needs a complex ctx.dtype.
max_iter (int, optional) – Maximum Krylov dimension. Static under jax.jit.
tol (float, optional) – Breakdown threshold, and the residual below which converged is set.

Returns:

x is the result vector, num_iters the Krylov dimension reached, and residual_norm the normalized projection indicator \(|\beta_m \phi_{m-1}|\). The approximation is norm(v) V_m exp(t T_m) e_1 in a metric-orthonormal basis. Higham, Functions of Matrices, Sec. 13.2.2, Eq. (13.7).

Return type:

LinalgResult

Raises:

TypeError – If A is not a LinOp.
ValueError – If A is not square, is known non-Hermitian, or max_iter is invalid.

See also

lanczos_smallest: The same Krylov projection, for eigenpairs.

Examples

Apply the exponential of a diagonal operator.

>>> import numpy as np
>>> import spacecore as sc
>>> ctx = sc.Context(sc.NumpyOps(), dtype=np.float64)
>>> X = sc.DenseCoordinateSpace((2,), ctx)
>>> A = sc.DiagonalLinOp(ctx.asarray([0.0, 1.0]), X, ctx)
>>> v = ctx.asarray([2.0, 3.0])
>>> result = sc.expm_multiply(A, v, t=0.5, max_iter=5)
>>> np.allclose(result.x, [2.0, 3.0 * np.exp(0.5)], atol=1e-10)
True