Coverage for src/cvx/linalg/operators/__init__.py: 100%

1"""Symmetric linear operators for active-set and Krylov solvers. 

2 

3Active-set and Krylov solvers never need a symmetric matrix ``A`` as an explicit 

4array. They touch it through only a handful of products: a full matrix-vector 

5product ``A x``, an arbitrary sub-block product ``A[rows, cols] v``, and a direct 

6solve ``A[free, free]^{-1} rhs`` against a principal ("free") sub-block. This 

7subpackage captures exactly that contract as :class:`SymmetricOperator` -- together 

8with :meth:`SymmetricOperator.rcond_free`, a conditioning check an active-set 

9solver uses to detect a rank-deficient free block -- and provides four backends 

10that implement it at very different cost: 

11 

12* :class:`DenseOperator` wraps an explicit ``n x n`` matrix. 

13* :class:`GramOperator` represents ``A = M.T @ M`` from a factor ``M`` and never 

14 forms the ``n x n`` Gram matrix. 

15* :class:`FactorOperator` represents a diagonal-plus-low-rank 

16 ``A = diag(d) + U @ Delta @ U.T`` and inverts free blocks by the Woodbury 

17 identity. 

18* :class:`IncrementalDenseOperator` is a :class:`DenseOperator` that maintains the 

19 free-block inverse across single-index changes for active-set sweeps. 

20 

21The protocol and the shared index/conditioning helpers live in :mod:`~.base`; 

22each backend lives in its own module (:mod:`~.dense`, :mod:`~.gram`, 

23:mod:`~.factor`). 

24 

25Two solving strategies compose with the same protocol. A direct solver consumes 

26:meth:`SymmetricOperator.solve_free`; a matrix-free conjugate-gradient solver 

27consumes :meth:`SymmetricOperator.apply_free` (the free-block forward product) 

28and does its own iterating. Because a solver reaches ``A`` only through the 

29protocol, switching backend requires no change to the solver. 

30 

31The two products a solver needs beyond the free-block solve are the free-block 

32forward product (:meth:`SymmetricOperator.apply_free`, i.e. ``block_matvec`` with 

33``rows == cols``) and a cross product ``A[bound, free] x_free`` -- the reduced 

34gradient at the bound indices -- obtained from :meth:`block_matvec` on the 

35disjoint index sets. ``block_matvec`` is defined for arbitrary ``rows`` and 

36``cols``, so both regimes and everything between are covered. 

37 

38Example: 

39 >>> import numpy as np 

40 >>> from cvx.linalg import GramOperator 

41 >>> rng = np.random.default_rng(0) 

42 >>> M = rng.standard_normal((20, 5)) 

43 >>> op = GramOperator(M) 

44 >>> A = M.T @ M 

45 >>> np.allclose(op.matvec(np.ones(5)), A @ np.ones(5)) 

46 True 

47 >>> free = np.array([0, 2, 4]) 

48 >>> rhs = np.array([1.0, 2.0, 3.0]) 

49 >>> x = op.solve_free(free, rhs) 

50 >>> np.allclose(A[np.ix_(free, free)] @ x, rhs) 

51 True 

52""" 

53 

54from .base import SymmetricOperator as SymmetricOperator 

55from .composite import SumOperator as SumOperator 

56from .dense import DenseOperator as DenseOperator 

57from .dense import IncrementalDenseOperator as IncrementalDenseOperator 

58from .factor import FactorOperator as FactorOperator 

59from .gram import GramOperator as GramOperator 

60 

61__all__ = [ 

62 "DenseOperator", 

63 "FactorOperator", 

64 "GramOperator", 

65 "IncrementalDenseOperator", 

66 "SumOperator", 

67 "SymmetricOperator", 

68]