Coverage for src/cvx/linalg/operators/dense.py: 100%
92 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-03 18:56 +0000
1"""Dense-matrix backends: :class:`DenseOperator` and :class:`IncrementalDenseOperator`.
3:class:`DenseOperator` wraps an explicit ``n x n`` matrix and slices it directly.
4:class:`IncrementalDenseOperator` specialises it for active-set sweeps, maintaining
5the free-block inverse across single-index changes with rank-one bordered / deletion
6updates instead of refactorising each step.
7"""
9from __future__ import annotations
11import numpy as np
13from ..core.exceptions import NonSquareMatrixError, NotAMatrixError
14from ..core.types import Matrix, Vector
15from ..decomposition.cholesky import cholesky_solve
16from .base import SymmetricOperator, _as_index, _rcond_symmetric
19class DenseOperator(SymmetricOperator):
20 """Symmetric operator backed by an explicit dense matrix.
22 Args:
23 matrix: A symmetric ``n x n`` matrix. It is stored by reference, not
24 copied or symmetrised.
26 Example:
27 >>> import numpy as np
28 >>> from cvx.linalg import DenseOperator
29 >>> A = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
30 >>> op = DenseOperator(A)
31 >>> free = np.array([0, 2])
32 >>> np.allclose(op.apply_free(free, np.array([1.0, 1.0])), A[np.ix_(free, free)] @ np.ones(2))
33 True
34 """
36 def __init__(self, matrix: Matrix) -> None:
37 """Store the backing matrix after checking it is square."""
38 matrix = np.asarray(matrix, dtype=np.float64)
39 if matrix.ndim != 2:
40 raise NotAMatrixError(matrix.ndim, func="DenseOperator")
41 if matrix.shape[0] != matrix.shape[1]:
42 raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])
43 self._a = matrix
45 @property
46 def n(self) -> int:
47 """Dimension of the operator."""
48 return int(self._a.shape[0])
50 @property
51 def diag(self) -> Vector:
52 """The diagonal of the backing matrix (a read-only view)."""
53 result: Vector = np.diagonal(self._a)
54 return result
56 def matvec(self, x: Vector | Matrix) -> Vector | Matrix:
57 """Return ``A @ x`` by dense multiplication."""
58 return self._a @ x
60 def restricted(self, free: object) -> DenseOperator:
61 """Return ``DenseOperator(A[free, free])``: the free block, pre-sliced."""
62 free = _as_index(free)
63 return DenseOperator(self._a[np.ix_(free, free)])
65 def block_matvec(self, rows: object, cols: object, v: Vector | Matrix) -> Vector | Matrix:
66 """Return ``A[rows, cols] @ v`` by slicing the dense matrix."""
67 rows = _as_index(rows)
68 cols = _as_index(cols)
69 result: Vector | Matrix = self._a[np.ix_(rows, cols)] @ v
70 return result
72 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix:
73 """Solve the free block by Cholesky (LU fallback via :func:`cholesky_solve`)."""
74 free = _as_index(free)
75 return cholesky_solve(self._a[np.ix_(free, free)], rhs)
77 def rcond_free(self, free: object) -> float:
78 """Reciprocal condition number of the free block from its symmetric eigenvalues."""
79 free = _as_index(free)
80 return _rcond_symmetric(self._a[np.ix_(free, free)])
83class IncrementalDenseOperator(DenseOperator):
84 """Dense operator that maintains ``A[free, free]^{-1}`` across single-index flips.
86 A drop-in :class:`DenseOperator` whose :meth:`solve_free` reuses the previous
87 free-block inverse when the free set changed by exactly one index since the last
88 call, updating it with a rank-one bordered (index added) or deletion (index
89 removed) formula at ``O(len(free)**2)`` instead of refactorising at
90 ``O(len(free)**3)``. Any other change -- the first solve, a multi-index change,
91 or a non-positive/non-finite pivot -- recomputes the inverse from scratch.
93 This suits an active-set loop that changes its free set one index at a time. The
94 free index arrays must be **ascending** (as produced by ``np.flatnonzero`` of a
95 boolean mask) and *rhs* aligned to that order; the maintained inverse and the
96 returned solution follow the same order. :meth:`matvec`, :meth:`block_matvec`, and
97 :meth:`rcond_free` are the plain dense ones -- only :meth:`solve_free` differs.
99 A maintained inverse accumulates rounding over the ``O(n)`` updates of a sweep, so
100 on ill-conditioned problems the plain :class:`DenseOperator` (a clean solve each
101 step) is the safer choice.
103 Example:
104 >>> import numpy as np
105 >>> from cvx.linalg import IncrementalDenseOperator
106 >>> op = IncrementalDenseOperator(np.eye(3))
107 >>> np.allclose(op.solve_free(np.array([0, 1]), np.array([1.0, 2.0])), [1.0, 2.0])
108 True
109 """
111 def __init__(self, matrix: Matrix) -> None:
112 """Wrap ``matrix`` (validated as in :class:`DenseOperator`) and start with no cache."""
113 super().__init__(matrix)
114 self._free_idx: np.ndarray | None = None
115 self._inv: Matrix | None = None
117 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix:
118 """Solve ``A[free, free] @ y = rhs`` using the maintained (incrementally updated) inverse."""
119 cur = _as_index(free)
120 inv = self._inverse_for(cur)
121 self._free_idx, self._inv = cur, inv
122 return inv @ rhs
124 def _inverse_for(self, cur: np.ndarray) -> Matrix:
125 """Return ``A[cur, cur]^{-1}``, updating the cache incrementally when possible."""
126 prev, prev_inv = self._free_idx, self._inv
127 if prev is None or prev_inv is None:
128 return self._refactor(cur)
129 added = np.setdiff1d(cur, prev, assume_unique=True)
130 removed = np.setdiff1d(prev, cur, assume_unique=True)
131 if added.size == 1 and removed.size == 0:
132 updated = self._insert(prev, prev_inv, int(added[0]), cur)
133 elif removed.size == 1 and added.size == 0:
134 updated = self._delete(prev, prev_inv, int(removed[0]))
135 else:
136 updated = None # not a single-index flip; recompute
137 return updated if updated is not None else self._refactor(cur)
139 def _refactor(self, cur: np.ndarray) -> Matrix:
140 """Invert ``A[cur, cur]`` from scratch."""
141 if cur.size == 0:
142 return np.zeros((0, 0))
143 inv: Matrix = np.linalg.inv(self._a[np.ix_(cur, cur)])
144 return inv
146 def _insert(self, prev: np.ndarray, prev_inv: Matrix, asset: int, cur: np.ndarray) -> Matrix | None:
147 """Rank-one bordered update for one index entering the free set (``None`` if the pivot is bad)."""
148 c = self._a[prev, asset]
149 v = prev_inv @ c
150 schur = float(self._a[asset, asset] - c @ v)
151 if not np.isfinite(schur) or schur <= 0.0:
152 return None
153 k = prev.shape[0]
154 aug = np.empty((k + 1, k + 1))
155 aug[:k, :k] = prev_inv + np.outer(v, v) / schur
156 aug[:k, k] = -v / schur
157 aug[k, :k] = -v / schur
158 aug[k, k] = 1.0 / schur
159 # ``aug`` is ordered [prev..., asset]; permute to ascending ``cur`` order.
160 perm = np.empty(k + 1, dtype=np.intp)
161 is_new = cur == asset
162 perm[is_new] = k
163 perm[~is_new] = np.searchsorted(prev, cur[~is_new])
164 permuted: Matrix = aug[np.ix_(perm, perm)]
165 return permuted
167 def _delete(self, prev: np.ndarray, prev_inv: Matrix, asset: int) -> Matrix | None:
168 """Rank-one deletion update for one index leaving the free set (``None`` if the pivot is bad)."""
169 p = int(np.searchsorted(prev, asset))
170 pivot = float(prev_inv[p, p])
171 if not np.isfinite(pivot) or pivot <= 0.0:
172 return None
173 mask = np.ones(prev.shape[0], dtype=bool)
174 mask[p] = False
175 col = prev_inv[mask, p]
176 updated: Matrix = prev_inv[np.ix_(mask, mask)] - np.outer(col, col) / pivot
177 return updated