Coverage for src / basanos / math / _linalg.py: 100%
57 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-19 05:23 +0000
1"""Linear algebra helpers used by the Basanos optimizer.
3This private module provides a small set of NumPy-based utilities for
4working with symmetric (correlation-like) matrices in a robust way:
5- valid(matrix): mask out rows/cols with non-finite diagonal entries and
6 return the corresponding sub-matrix.
7- is_positive_definite(matrix): explicitly test whether a matrix is
8 positive-definite using Cholesky decomposition.
9- inv_a_norm(vector, matrix=None, cond_threshold=...): compute the inverse
10 A-norm of a vector with respect to a (possibly masked) positive-definite
11 matrix; defaults to the Euclidean norm when no matrix is given. Emits
12 IllConditionedMatrixWarning when the condition number exceeds the threshold.
13- solve(matrix, rhs, cond_threshold=...): solve a linear system on the valid
14 subset indicated by finite diagonal entries, returning NaNs for invalid
15 positions. Emits IllConditionedMatrixWarning when the condition number
16 exceeds the threshold.
18Both solve and inv_a_norm use Cholesky decomposition as a numerically stable
19first attempt; they fall back to standard LU-based solving for matrices that
20are not positive-definite.
22These routines are intentionally lightweight and raise domain-specific
23exceptions to guard against shape mismatches and numerical errors.
24They are internal implementation details and may change without notice.
25"""
27from __future__ import annotations
29import warnings
31import numpy as np
32from scipy.linalg import cho_factor, cho_solve
34from basanos.exceptions import (
35 DimensionMismatchError,
36 IllConditionedMatrixWarning,
37 NonSquareMatrixError,
38 SingularMatrixError,
39)
41_DEFAULT_COND_THRESHOLD: float = 1e12
42"""Default condition-number threshold above which a warning is emitted."""
45def valid(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
46 """Validate a square matrix and return mask and valid sub-matrix.
48 Checks that the input is square and that diagonal entries are finite.
49 Returns a boolean mask for valid diagonal entries and the sub-matrix
50 consisting of valid rows/columns.
52 Args:
53 matrix (np.ndarray): Input square matrix.
55 Returns:
56 tuple[np.ndarray, np.ndarray]:
57 - mask: Boolean array indicating which diagonal elements are finite.
58 - submatrix: Matrix containing only rows/columns with finite diagonal.
60 Raises:
61 NonSquareMatrixError: If the input matrix is not square.
63 Examples:
64 >>> import numpy as np
65 >>> mask, sub = valid(np.eye(2))
66 >>> mask.tolist()
67 [True, True]
68 >>> sub.shape
69 (2, 2)
70 """
71 # make sure matrix is quadratic
72 if matrix.shape[0] != matrix.shape[1]:
73 raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])
75 v = np.isfinite(np.diag(matrix))
76 return v, matrix[:, v][v]
79def is_positive_definite(matrix: np.ndarray) -> bool:
80 """Return True if *matrix* is symmetric positive-definite, False otherwise.
82 The check is performed via an attempted Cholesky decomposition — the most
83 numerically reliable way to test positive-definiteness for symmetric
84 matrices. A matrix is positive-definite if and only if its Cholesky
85 decomposition exists and all diagonal entries of the Cholesky factor are
86 strictly positive.
88 This function is intentionally side-effect-free: it raises no exceptions
89 and emits no warnings. It is suitable for use as a guard before passing a
90 shrunk correlation matrix to a linear solver.
92 Args:
93 matrix (np.ndarray): Square matrix to test.
95 Returns:
96 bool: ``True`` if the matrix is positive-definite, ``False`` otherwise.
98 Examples:
99 >>> import numpy as np
100 >>> is_positive_definite(np.eye(3))
101 True
102 >>> is_positive_definite(np.array([[1.0, 2.0], [2.0, 1.0]]))
103 False
104 >>> is_positive_definite(np.array([[1.0, 0.5], [0.5, 1.0]]))
105 True
106 """
107 try:
108 cho_factor(matrix, check_finite=False)
109 except np.linalg.LinAlgError:
110 return False
111 else:
112 return True
115def _check_and_warn_condition(matrix: np.ndarray, threshold: float) -> None:
116 """Emit IllConditionedMatrixWarning when the condition number exceeds threshold.
118 Args:
119 matrix (np.ndarray): Square matrix whose condition number is checked.
120 threshold (float): Upper bound before a warning is issued.
121 """
122 cond = float(np.linalg.cond(matrix))
123 if cond > threshold:
124 warnings.warn(
125 f"Matrix condition number {cond:.3e} exceeds threshold {threshold:.3e}; "
126 "results may be numerically unreliable.",
127 IllConditionedMatrixWarning,
128 stacklevel=3,
129 )
132def _cholesky_solve(matrix: np.ndarray, rhs: np.ndarray) -> np.ndarray:
133 """Solve *matrix* x = *rhs* using Cholesky, falling back to LU if needed.
135 Cholesky decomposition is attempted first as it is numerically more stable
136 for positive-definite matrices. If the decomposition fails (the matrix is
137 not positive-definite), a standard LU-based solve is used instead.
139 Args:
140 matrix (np.ndarray): Square coefficient matrix.
141 rhs (np.ndarray): Right-hand side vector.
143 Returns:
144 np.ndarray: Solution vector.
146 Raises:
147 np.linalg.LinAlgError: If both Cholesky and LU-based solves fail.
148 """
149 try:
150 c, lower = cho_factor(matrix, check_finite=False)
151 return cho_solve((c, lower), rhs, check_finite=False)
152 except np.linalg.LinAlgError:
153 # Matrix is not positive-definite; fall back to general LU solve.
154 return np.linalg.solve(matrix, rhs)
157def inv_a_norm(
158 vector: np.ndarray,
159 matrix: np.ndarray | None = None,
160 cond_threshold: float = _DEFAULT_COND_THRESHOLD,
161) -> float:
162 """Compute inverse A-norm of a vector with an optional metric matrix.
164 If ``matrix`` is None, compute the Euclidean norm of finite entries.
165 Otherwise, validate the matrix is square and dimensionally compatible,
166 then return sqrt(v^T A^{-1} v) on the valid subset.
168 Cholesky decomposition is attempted first for numerical stability; the
169 solver falls back to LU decomposition for non-positive-definite matrices.
170 If the condition number of the valid sub-matrix exceeds *cond_threshold*,
171 an :class:`~basanos.exceptions.IllConditionedMatrixWarning` is emitted.
173 Args:
174 vector (np.ndarray): Input vector.
175 matrix (np.ndarray | None, optional): Positive-definite metric matrix.
176 cond_threshold (float, optional): Condition-number threshold above which
177 a warning is emitted. Defaults to ``1e12``.
179 Returns:
180 float: The computed norm. Returns np.nan if no valid entries exist.
182 Raises:
183 NonSquareMatrixError: If ``matrix`` is not square.
184 DimensionMismatchError: If ``vector`` length does not match the matrix dimension.
185 SingularMatrixError: If the valid sub-matrix is singular.
187 Examples:
188 >>> import numpy as np
189 >>> inv_a_norm(np.array([3.0, 4.0]))
190 5.0
191 """
192 if matrix is None:
193 return float(np.linalg.norm(vector[np.isfinite(vector)], 2))
195 # make sure matrix is quadratic
196 if matrix.shape[0] != matrix.shape[1]:
197 raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])
199 # make sure the vector has the right number of entries
200 if vector.size != matrix.shape[0]:
201 raise DimensionMismatchError(vector.size, matrix.shape[0])
203 v, mat = valid(matrix)
205 if v.any():
206 _check_and_warn_condition(mat, cond_threshold)
207 try:
208 return float(np.sqrt(np.dot(vector[v], _cholesky_solve(mat, vector[v]))))
209 except np.linalg.LinAlgError as exc:
210 raise SingularMatrixError(str(exc)) from exc
211 return float(np.nan)
214def solve(
215 matrix: np.ndarray,
216 rhs: np.ndarray,
217 cond_threshold: float = _DEFAULT_COND_THRESHOLD,
218) -> np.ndarray:
219 """Solve Ax=b on the valid subset indicated by finite diagonal entries.
221 Cholesky decomposition is attempted first for numerical stability; the
222 solver falls back to LU decomposition for non-positive-definite matrices.
223 If the condition number of the valid sub-matrix exceeds *cond_threshold*,
224 an :class:`~basanos.exceptions.IllConditionedMatrixWarning` is emitted.
226 Args:
227 matrix (np.ndarray): Square coefficient matrix A.
228 rhs (np.ndarray): Right-hand side vector b with matching size.
229 cond_threshold (float, optional): Condition-number threshold above which
230 a warning is emitted. Defaults to ``1e12``.
232 Returns:
233 np.ndarray: Solution vector with NaN for invalid positions.
235 Raises:
236 NonSquareMatrixError: If matrix is not square.
237 DimensionMismatchError: If rhs size does not match the matrix dimension.
238 SingularMatrixError: If the valid sub-matrix is singular.
240 Examples:
241 >>> import numpy as np
242 >>> solve(np.eye(2), np.array([1.0, 2.0])).tolist()
243 [1.0, 2.0]
244 """
245 # make sure matrix is quadratic
246 if matrix.shape[0] != matrix.shape[1]:
247 raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])
249 # make sure the vector rhs has the right number of entries
250 if rhs.size != matrix.shape[0]:
251 raise DimensionMismatchError(rhs.size, matrix.shape[0])
253 x = np.nan * np.ones(rhs.size)
254 v, mat = valid(matrix)
256 if v.any():
257 _check_and_warn_condition(mat, cond_threshold)
258 try:
259 x[v] = _cholesky_solve(mat, rhs[v])
260 except np.linalg.LinAlgError as exc:
261 raise SingularMatrixError(str(exc)) from exc
263 return x