Coverage for src/cvx/linalg/operators/factor.py: 100%
75 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-03 18:56 +0000
1""":class:`FactorOperator`: diagonal-plus-low-rank ``A = diag(d) + U @ Delta @ U.T``."""
3from __future__ import annotations
5import numpy as np
7from ..core.exceptions import DimensionMismatchError, NonSquareMatrixError, NotAMatrixError
8from ..core.types import Matrix, Vector
9from ..decomposition.cholesky import cholesky_solve
10from .base import SymmetricOperator, _as_index
12_DIAGONAL_NDIM_MESSAGE = "diagonal must be a 1-D array"
13_DIAGONAL_POSITIVE_MESSAGE = "diagonal entries must be strictly positive"
16class FactorOperator(SymmetricOperator):
17 """Diagonal-plus-low-rank operator ``A = diag(d) + U @ Delta @ U.T``.
19 Free-block solves use the Woodbury identity, costing ``O(len(free) r**2 +
20 r**3)`` for a rank-``r`` factor rather than ``O(len(free)**3)``, and no
21 ``n x n`` matrix is formed (memory ``O(n r)``). With a strictly positive
22 diagonal *d* and positive-definite *Delta* every principal block is positive
23 definite, so :meth:`solve_free` is always well posed.
25 Args:
26 diagonal: The strictly positive diagonal ``d`` of length ``n``.
27 loadings: The ``n x r`` factor loadings ``U``.
28 inner: The ``r x r`` positive-definite inner matrix ``Delta``.
30 Example:
31 >>> import numpy as np
32 >>> from cvx.linalg import FactorOperator
33 >>> d = np.array([2.0, 3.0, 4.0])
34 >>> U = np.array([[1.0], [0.5], [-1.0]])
35 >>> Delta = np.array([[2.0]])
36 >>> op = FactorOperator(d, U, Delta)
37 >>> (op.n, op.k) # 3 assets, 1 factor
38 (3, 1)
39 >>> A = np.diag(d) + U @ Delta @ U.T
40 >>> free, rhs = np.array([0, 2]), np.array([1.0, 1.0])
41 >>> np.allclose(A[np.ix_(free, free)] @ op.solve_free(free, rhs), rhs)
42 True
43 """
45 def __init__(self, diagonal: Vector, loadings: Matrix, inner: Matrix) -> None:
46 """Store the diagonal, loadings, and inner block after shape checks."""
47 d = np.asarray(diagonal, dtype=np.float64)
48 u = np.asarray(loadings, dtype=np.float64)
49 delta = np.asarray(inner, dtype=np.float64)
50 if d.ndim != 1:
51 raise ValueError(_DIAGONAL_NDIM_MESSAGE)
52 if np.any(d <= 0.0):
53 raise ValueError(_DIAGONAL_POSITIVE_MESSAGE)
54 if u.ndim != 2:
55 raise NotAMatrixError(u.ndim, func="FactorOperator")
56 if u.shape[0] != d.shape[0]:
57 raise DimensionMismatchError(u.shape[0], d.shape[0])
58 if delta.ndim != 2:
59 raise NotAMatrixError(delta.ndim, func="FactorOperator")
60 if delta.shape[0] != delta.shape[1]:
61 raise NonSquareMatrixError(delta.shape[0], delta.shape[1])
62 if delta.shape[0] != u.shape[1]:
63 raise DimensionMismatchError(delta.shape[0], u.shape[1])
64 self._d = d
65 self._u = u
66 self._delta = delta
68 @property
69 def n(self) -> int:
70 """Dimension of the operator (length of the diagonal ``d``)."""
71 return int(self._d.shape[0])
73 @property
74 def k(self) -> int:
75 """Number of factors (rank ``r`` of the low-rank term; columns of ``U``)."""
76 return int(self._u.shape[1])
78 @property
79 def diag(self) -> Vector:
80 """The diagonal ``d_i + U[i] @ Delta @ U[i]``, at ``O(n r**2)`` without forming ``A``."""
81 result: Vector = self._d + np.einsum("ij,ij->i", self._u @ self._delta, self._u)
82 return result
84 def matvec(self, x: Vector | Matrix) -> Vector | Matrix:
85 """Return ``A @ x = d * x + U @ (Delta @ (U.T @ x))``."""
86 return (self._d * x.T).T + self._u @ (self._delta @ (self._u.T @ x))
88 def restricted(self, free: object) -> FactorOperator:
89 """Return ``FactorOperator(d[free], U[free], Delta)``: the free block, pre-sliced."""
90 free = _as_index(free)
91 return FactorOperator(self._d[free], np.ascontiguousarray(self._u[free, :]), self._delta)
93 def block_matvec(self, rows: object, cols: object, v: Vector | Matrix) -> Vector | Matrix:
94 """Return ``A[rows, cols] @ v`` from the low-rank term and the diagonal overlap."""
95 rows = _as_index(rows)
96 cols = _as_index(cols)
97 low_rank = self._u[rows] @ (self._delta @ (self._u[cols].T @ v))
98 # Diagonal couples only positions where a row index equals a column index.
99 common, r_idx, c_idx = np.intersect1d(rows, cols, return_indices=True)
100 diag = np.zeros_like(low_rank)
101 diag[r_idx] = (self._d[common] * np.asarray(v)[c_idx].T).T
102 result: Vector | Matrix = low_rank + diag
103 return result
105 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix:
106 """Solve the free block by the Woodbury identity on the ``r x r`` capacitance matrix."""
107 free = _as_index(free)
108 df = self._d[free]
109 uf = self._u[free]
110 # Woodbury: A_FF^{-1} = D^{-1} - D^{-1} U W^{-1} U.T D^{-1},
111 # with W = Delta^{-1} + U.T D^{-1} U.
112 dinv_rhs = (np.asarray(rhs, dtype=np.float64).T / df).T
113 delta_inv = np.linalg.solve(self._delta, np.eye(self._delta.shape[0]))
114 w = delta_inv + uf.T @ ((uf.T / df).T)
115 inner = cholesky_solve(w, uf.T @ dinv_rhs)
116 correction = (uf @ inner).T / df
117 result: Vector | Matrix = dinv_rhs - correction.T
118 return result
120 def rcond_free(self, free: object) -> float:
121 """Lower bound on the free block's reciprocal condition number, via Weyl's inequalities.
123 The free block ``diag(d_F) + U_F Delta U_F.T`` is positive definite (the
124 positive diagonal keeps it full rank). Rather than form it, bound
125 ``lambda_min >= min(d_F)`` and
126 ``lambda_max <= max(d_F) + ||U_F||_2^2 * lambda_max(Delta)``; their ratio is
127 a guaranteed lower bound on the true reciprocal condition number, at
128 ``O(len(free) r**2 + r**3)`` and without an ``n x n`` matrix.
129 """
130 free = _as_index(free)
131 if free.size == 0:
132 return 1.0
133 d_free = self._d[free]
134 u_free = self._u[free]
135 u_spectral_norm = float(np.linalg.svd(u_free, compute_uv=False)[0])
136 delta_max = float(np.linalg.eigvalsh(self._delta)[-1])
137 lam_max_upper = float(np.max(d_free)) + u_spectral_norm**2 * max(delta_max, 0.0)
138 return float(np.min(d_free)) / lam_max_upper