Coverage for src/cvxcla/operators/factor.py: 100%
76 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-27 06:45 +0000
1"""Diagonal-plus-low-rank covariance backend.
3``FactorCovariance`` implements the :class:`~cvxcla.operators._core.QuadraticForm`
4protocol for factor risk models and RMT-cleaned covariances of the form
5``Sigma = diag(d) + U @ Delta @ U.T``, solving against the free block through the
6Woodbury identity so no ``n x n`` matrix is ever materialised.
7"""
9from __future__ import annotations
11from dataclasses import dataclass
12from functools import cached_property
13from typing import cast
15import numpy as np
16from numpy.typing import NDArray
19@dataclass(frozen=True)
20class FactorCovariance:
21 """Diagonal-plus-low-rank covariance ``Sigma = diag(d) + U @ Delta @ U.T``.
23 This backend implements the ``CovarianceOperator`` protocol for factor
24 risk models and RMT-cleaned covariances (constant-residual-eigenvalue /
25 eigenvalue clipping) without ever forming an ``n x n`` matrix. Solves
26 against the free block go through the Woodbury identity
28 ``Sigma_FF^{-1} = D_F^{-1} - D_F^{-1} U_F W^{-1} U_F^T D_F^{-1}``
30 with the ``k x k`` correction matrix ``W = Delta^{-1} + U_F^T D_F^{-1} U_F``,
31 so a solve costs ``O(n_F k^2 + k^3)`` instead of ``O(n_F^3)``. The
32 cross-product against blocked assets only involves the low-rank part,
33 since the diagonal contributes nothing off the diagonal.
35 Attributes:
36 d: Positive idiosyncratic variances of shape ``(n,)``.
37 u: Factor loadings of shape ``(n, k)``.
38 delta: Factor covariance, either as eigenvalues/variances of shape
39 ``(k,)`` (interpreted as a diagonal matrix) or as a symmetric
40 positive definite matrix of shape ``(k, k)``.
42 Examples:
43 >>> import numpy as np
44 >>> cov = FactorCovariance(d=np.array([1.0, 1.0]), u=np.eye(2), delta=np.array([1.0, 1.0]))
45 >>> cov.n
46 2
47 >>> cov.matvec(np.array([1.0, 2.0]))
48 array([2., 4.])
49 """
51 d: NDArray[np.float64]
52 u: NDArray[np.float64]
53 delta: NDArray[np.float64]
55 def __post_init__(self) -> None:
56 """Validate shapes, positivity of the diagonal, and symmetry of ``delta``.
58 Raises:
59 ValueError: If ``d`` is not a positive vector, ``u`` is not an
60 ``(n, k)`` matrix matching ``d``, or ``delta`` is neither a
61 ``(k,)`` vector nor a symmetric ``(k, k)`` matrix.
62 """
63 if self.d.ndim != 1 or not np.all(self.d > 0):
64 msg = "d must be a vector of positive idiosyncratic variances"
65 raise ValueError(msg)
66 if self.u.ndim != 2 or self.u.shape[0] != self.d.shape[0]:
67 msg = f"u must have shape (n, k) with n = {self.d.shape[0]}, got {self.u.shape}"
68 raise ValueError(msg)
69 k = self.u.shape[1]
70 if self.delta.ndim == 1:
71 if self.delta.shape[0] != k:
72 msg = f"delta must have {k} entries, got {self.delta.shape[0]}"
73 raise ValueError(msg)
74 if not np.all(self.delta > 0):
75 msg = "A diagonal delta must have positive entries"
76 raise ValueError(msg)
77 elif self.delta.ndim == 2:
78 if self.delta.shape != (k, k):
79 msg = f"delta must have shape ({k}, {k}), got {self.delta.shape}"
80 raise ValueError(msg)
81 if not np.allclose(self.delta, self.delta.T):
82 msg = "delta must be symmetric"
83 raise ValueError(msg)
84 else:
85 msg = f"delta must be a (k,) vector or (k, k) matrix, got ndim {self.delta.ndim}"
86 raise ValueError(msg)
88 @property
89 def n(self) -> int:
90 """Number of assets (the dimension of the covariance)."""
91 return int(self.d.shape[0])
93 @property
94 def k(self) -> int:
95 """Number of factors (the rank of the low-rank part)."""
96 return int(self.u.shape[1])
98 @cached_property
99 def _delta_matrix(self) -> NDArray[np.float64]:
100 """The factor covariance ``Delta`` as a ``(k, k)`` matrix."""
101 if self.delta.ndim == 1:
102 return np.diag(self.delta)
103 return self.delta
105 @cached_property
106 def _delta_inv(self) -> NDArray[np.float64]:
107 """The inverse of the ``(k, k)`` factor covariance."""
108 if self.delta.ndim == 1:
109 return np.diag(1.0 / self.delta)
110 return cast("NDArray[np.float64]", np.linalg.inv(self.delta))
112 def matvec(self, x: NDArray[np.float64]) -> NDArray[np.float64]:
113 """Compute the matrix-vector product ``Sigma @ x``.
115 Args:
116 x: Vector of shape ``(n,)`` or matrix of shape ``(n, r)``.
118 Returns:
119 ``Sigma @ x`` with the same trailing shape as ``x``, computed in
120 ``O(n k)`` per column.
121 """
122 low_rank = self.u @ (self._delta_matrix @ (self.u.T @ x))
123 diagonal = self.d * x if x.ndim == 1 else self.d[:, None] * x
124 return diagonal + low_rank
126 def solve_free(self, free: NDArray[np.bool_], rhs: NDArray[np.float64]) -> NDArray[np.float64]:
127 """Solve the free-block system ``Sigma[free][:, free] @ y = rhs`` via Woodbury.
129 Forms the ``k x k`` correction matrix
130 ``W = Delta^{-1} + U_F^T D_F^{-1} U_F`` and solves against it, so the
131 cost is ``O(n_F k^2 + k^3 + n_F k r)`` for ``r`` right-hand sides.
133 Args:
134 free: Boolean mask of shape ``(n,)`` selecting the free assets.
135 rhs: Right-hand side of shape ``(n_free,)`` or ``(n_free, r)``.
137 Returns:
138 The solution ``y`` with the same shape as ``rhs``.
139 """
140 d_free = self.d[free]
141 u_free = self.u[free]
143 rhs_2d = rhs if rhs.ndim == 2 else rhs[:, None]
144 dinv_rhs = rhs_2d / d_free[:, None]
145 dinv_u = u_free / d_free[:, None]
147 w = self._delta_inv + u_free.T @ dinv_u
148 solution = dinv_rhs - dinv_u @ np.linalg.solve(w, u_free.T @ dinv_rhs)
149 return solution if rhs.ndim == 2 else solution[:, 0]
151 def cross(self, free: NDArray[np.bool_], x: NDArray[np.float64]) -> NDArray[np.float64]:
152 """Compute the free-to-blocked cross product ``Sigma[free][:, ~free] @ x[~free]``.
154 Only the low-rank part contributes, since the diagonal of ``Sigma``
155 has no off-diagonal block: the result is
156 ``U_F @ Delta @ (U_out^T @ x_out)``, an ``O(n k)`` computation.
158 Args:
159 free: Boolean mask of shape ``(n,)`` selecting the free assets.
160 x: Full-length vector of shape ``(n,)``; only the blocked entries
161 ``x[~free]`` enter the product.
163 Returns:
164 Vector of shape ``(n_free,)``.
165 """
166 blocked = ~free
167 return self.u[free] @ (self._delta_matrix @ (self.u[blocked].T @ x[blocked]))
169 def rcond_free(self, free: NDArray[np.bool_]) -> float:
170 """Lower bound on the free block's reciprocal condition number.
172 The free block ``Sigma_FF = diag(d_F) + U_F @ Delta @ U_F^T`` is
173 positive definite by construction (the positive idiosyncratic floor
174 ``d`` keeps it full rank), which is exactly why this backend is the
175 documented remedy for rank deficiency. Rather than form the
176 ``n_F x n_F`` block, we bound the conditioning from the structure via
177 Weyl's inequalities:
179 * ``lambda_min(Sigma_FF) >= min(d_F)`` (the low-rank part is PSD), and
180 * ``lambda_max(Sigma_FF) <= max(d_F) + ||U_F||_2^2 * lambda_max(Delta)``.
182 Their ratio is a guaranteed *lower* bound on the true reciprocal
183 condition number, computed in ``O(n_F k^2 + k^3)`` without ever
184 materialising an ``n x n`` matrix. For a well-floored factor model it
185 sits far above the CLA's singularity threshold, so the guard never
186 trips; a pathologically tiny floor correctly drives it toward zero.
188 Args:
189 free: Boolean mask of shape ``(n,)`` selecting the free assets.
191 Returns:
192 A lower bound on the reciprocal 2-norm condition number, in ``[0, 1]``.
193 """
194 if not np.any(free):
195 return 1.0
196 d_free = self.d[free]
197 u_free = self.u[free]
198 u_spectral_norm = float(np.linalg.svd(u_free, compute_uv=False)[0])
199 delta_max = float(np.linalg.eigvalsh(self._delta_matrix)[-1])
200 lam_max_upper = float(np.max(d_free)) + u_spectral_norm**2 * max(delta_max, 0.0)
201 return float(np.min(d_free)) / lam_max_upper