Coverage for src/nncg/api.py: 100%
36 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-15 07:01 +0000
1"""One-call convenience entry points over the core solvers.
3:func:`solve_nnqp` and :func:`solve_nnqp_eq` compose the three pieces of the
4core API — wrap a plain SPD array in ``DenseOperator``, default-construct the
5inner solver from a bare string, bundle the outer-loop knobs into an
6:class:`~nncg.solver.ActiveSetConfig` — and delegate to
7:class:`~nncg.solver.ActiveSetSolver`. They hold no logic of their own: reach
8past them to ``ActiveSetSolver`` directly whenever you need to reuse a
9configured solver across problems, or an inner solver the string shortcut cannot
10express (``inner=Nystrom(nystrom=NystromConfig(rank=20))`` still works here,
11passed as an instance). :func:`solve_nnqp_mprgp` is the matching one-call wrapper
12over the projection-based :class:`~nncg.mprgp.MPRGP` solver for the same
13bound-constrained problem.
14"""
16from __future__ import annotations
18from collections.abc import Callable
19from typing import Literal
21import numpy as np
22from cvx.linalg import DenseOperator, Matrix, SymmetricOperator, Vector
23from numpy.typing import NDArray
25from .inner import CG, Exact, GlobalNystrom, Jacobi, Nystrom
26from .mprgp import MPRGP, MPRGPConfig, MPRGPResult
27from .solver import ActiveSetConfig, ActiveSetSolver, InnerSolver, Result
29#: Bare-string shortcuts mapping to a default-constructed inner solver.
30_INNER: dict[str, Callable[[], InnerSolver]] = {
31 "cg": CG,
32 "jacobi": Jacobi,
33 "nystrom": Nystrom,
34 "global_nystrom": GlobalNystrom,
35 "exact": Exact,
36}
38InnerKind = Literal["cg", "jacobi", "nystrom", "global_nystrom", "exact"]
39"""The bare-string shortcuts accepted for ``inner`` (keys of :data:`_INNER`)."""
42def _resolve_inner(inner: InnerSolver | InnerKind) -> InnerSolver:
43 """Return the inner solver, default-constructing it from a shortcut string.
45 Args:
46 inner: An :class:`~nncg.solver.InnerSolver` instance (returned as-is), or
47 one of the shortcut strings in :data:`_INNER`.
49 Returns:
50 An inner-solver instance.
52 Raises:
53 ValueError: When ``inner`` is a string outside the shortcut set.
54 """
55 if isinstance(inner, str):
56 try:
57 return _INNER[inner]()
58 except KeyError:
59 valid = ", ".join(map(repr, _INNER))
60 msg = f"unknown inner solver {inner!r}; pass an InnerSolver instance or one of {valid}"
61 raise ValueError(msg) from None
62 return inner
65def _as_operator(a: SymmetricOperator | NDArray[np.float64]) -> SymmetricOperator:
66 """Return ``a`` as a :class:`cvx.linalg.SymmetricOperator`, wrapping a plain array.
68 A :class:`~cvx.linalg.SymmetricOperator` is used unchanged; anything else is
69 treated as an explicit SPD array and wrapped in ``DenseOperator``. The
70 matrix-free ``A = M^T M + ridge I`` path is deliberately *not* inferred from
71 an ``M`` — pass ``GramOperator(M, ridge)`` explicitly for it.
73 Args:
74 a: A symmetric operator, or a 2-D SPD array.
76 Returns:
77 The operator form of ``a``.
78 """
79 if isinstance(a, SymmetricOperator):
80 return a
81 return DenseOperator(np.asarray(a, dtype=np.float64))
84def solve_nnqp(
85 a: SymmetricOperator | NDArray[np.float64],
86 b: Vector,
87 *,
88 inner: InnerSolver | InnerKind = "cg",
89 warm: tuple[NDArray[np.bool_], Vector] | None = None,
90 tol: float = 1e-8,
91 p_max: int = 3,
92 track: bool = False,
93 max_outer: int | None = None,
94) -> Result:
95 """Minimise ``1/2 x^T A x - b^T x`` over ``x >= 0`` — the one-call entry point.
97 A thin convenience wrapper that composes the three pieces of the layered API
98 for the common case: it wraps a plain SPD array in ``DenseOperator``, default-
99 constructs the inner solver from a bare string, and bundles the outer-loop
100 knobs into an :class:`ActiveSetConfig`, then delegates to
101 :meth:`ActiveSetSolver.solve`. It holds no logic of its own — reach past it to
102 :class:`ActiveSetSolver` directly whenever you need to reuse a configured
103 solver across problems, or an inner solver this shortcut cannot express.
105 Args:
106 a: The SPD quadratic term. A :class:`cvx.linalg.SymmetricOperator` is used
107 as-is; a plain 2-D array is wrapped in ``DenseOperator``. The matrix-
108 free ``A = M^T M + ridge I`` path is *not* inferred from an ``M`` —
109 pass ``GramOperator(M, ridge)`` explicitly for it.
110 b: The linear term ``b``.
111 inner: The inner solver for each free block, as an
112 :class:`nncg.inner.InnerSolver` instance (fully configurable — e.g.
113 ``Nystrom(nystrom=NystromConfig(rank=20))``), or one of the shortcut
114 strings ``"cg"``, ``"jacobi"``, ``"nystrom"``, ``"global_nystrom"``,
115 ``"exact"`` for its
116 default configuration.
117 warm: Optional ``(free_mask, x_prev)`` pair from a previous solve, forwarded
118 to :meth:`ActiveSetSolver.solve` — see there for the warm-start semantics.
119 tol: Threshold of the primal and dual KKT violator tests
120 (``ActiveSetConfig.tol``).
121 p_max: Patience budget before a least-index Bland fallback pivot
122 (``ActiveSetConfig.p_max``).
123 track: Record the visited free-set trajectory in ``Result.traj``.
124 max_outer: Optional cap on outer steps; when hit, the current iterate is
125 returned with ``converged=False``.
127 Returns:
128 A :class:`Result`; ``converged`` is True iff the KKT system was satisfied
129 to ``tol``, which certifies the unique global minimiser.
131 Raises:
132 TypeError: When ``a`` is neither a :class:`cvx.linalg.SymmetricOperator`
133 nor an array wrappable by ``DenseOperator``.
134 ValueError: When ``inner`` is a string outside the shortcut set, when the
135 operator dimension does not match ``len(b)``, or on the inner solver's
136 own conditions.
138 Examples:
139 The bound binds where the unconstrained minimiser would go negative. Here
140 ``A^-1 b = [1, -1]``, so the second coordinate is clamped to zero:
142 >>> import numpy as np
143 >>> a = np.array([[2.0, 0.0], [0.0, 2.0]])
144 >>> b = np.array([2.0, -2.0])
145 >>> result = solve_nnqp(a, b)
146 >>> result.converged
147 True
148 >>> result.x.round(6).tolist()
149 [1.0, 0.0]
150 """
151 config = ActiveSetConfig(tol=tol, p_max=p_max, track=track, max_outer=max_outer)
152 solver = ActiveSetSolver(inner=_resolve_inner(inner), config=config)
153 return solver.solve(_as_operator(a), b, warm=warm)
156def solve_nnqp_eq(
157 a: SymmetricOperator | NDArray[np.float64],
158 b: Vector,
159 b_eq: Matrix,
160 c_eq: Vector,
161 *,
162 inner: InnerSolver | InnerKind = "cg",
163 warm: tuple[NDArray[np.bool_], Vector] | None = None,
164 tol: float = 1e-8,
165 p_max: int = 3,
166 track: bool = False,
167 max_outer: int | None = None,
168) -> Result:
169 """Solve ``min 1/2 x^T A x - b^T x`` s.t. ``x >= 0`` and ``B x = c`` — one call.
171 The equality-augmented companion to :func:`solve_nnqp`, with identical
172 wrapping and configuration conventions; it delegates to
173 :meth:`ActiveSetSolver.solve_eq`, where the per-free-set saddle system and the
174 full-row-rank requirement on ``B`` are documented. The single normalisation
175 ``1^T x = beta`` is the ``p = 1`` case.
177 Args:
178 a: The SPD quadratic term — a :class:`cvx.linalg.SymmetricOperator`, or a
179 plain array wrapped in ``DenseOperator`` (see :func:`solve_nnqp`).
180 b: The linear term ``b``.
181 b_eq: Equality matrix ``B`` of shape ``(p, n)``, full row rank on the
182 visited free sets.
183 c_eq: Equality right-hand side ``c`` of shape ``(p,)``.
184 inner: The inner solver instance, or a shortcut string — see
185 :func:`solve_nnqp`.
186 warm: Optional ``(free_mask, x_prev)`` pair, forwarded to
187 :meth:`ActiveSetSolver.solve_eq`.
188 tol: KKT violator tolerance (``ActiveSetConfig.tol``).
189 p_max: Bland-fallback patience budget (``ActiveSetConfig.p_max``).
190 track: Record the visited free-set trajectory in ``Result.traj``.
191 max_outer: Optional outer-step cap; ``converged=False`` when hit.
193 Returns:
194 A :class:`Result` with the equality multipliers in ``lam``. The reduced
195 gradient underlying the dual test is ``s = A x - b - B^T lam``.
197 Raises:
198 TypeError: When ``a`` is neither a :class:`cvx.linalg.SymmetricOperator`
199 nor an array wrappable by ``DenseOperator``.
200 ValueError: When ``inner`` is a string outside the shortcut set, when the
201 operator dimension does not match ``len(b)``, or on the inner solver's
202 own conditions.
204 Examples:
205 The ``p = 1`` normalisation ``1^T x = 1`` — the minimum-norm point on the
206 simplex, here its centre:
208 >>> import numpy as np
209 >>> a = np.eye(2) * 2.0
210 >>> b = np.zeros(2)
211 >>> result = solve_nnqp_eq(a, b, np.array([[1.0, 1.0]]), np.array([1.0]))
212 >>> result.converged
213 True
214 >>> result.x.round(6).tolist()
215 [0.5, 0.5]
216 """
217 config = ActiveSetConfig(tol=tol, p_max=p_max, track=track, max_outer=max_outer)
218 solver = ActiveSetSolver(inner=_resolve_inner(inner), config=config)
219 return solver.solve_eq(_as_operator(a), b, b_eq, c_eq, warm=warm)
222def solve_nnqp_mprgp(
223 a: SymmetricOperator | NDArray[np.float64],
224 b: Vector,
225 *,
226 x0: Vector | None = None,
227 tol: float = 1e-8,
228 gamma: float = 1.0,
229 alpha_bar: float | None = None,
230 max_iter: int = 100_000,
231 seed: int = 0,
232) -> MPRGPResult:
233 """Minimise ``1/2 x^T A x - b^T x`` over ``x >= 0`` by MPRGP — one call.
235 The projection-based companion to :func:`solve_nnqp`: it solves the same
236 bound-constrained program with Dostál & Schöberl's MPRGP
237 (:class:`nncg.mprgp.MPRGP`) instead of the active-set loop — matrix-free and
238 factorisation-free, so it never forms or refactorises ``A``. Like
239 :func:`solve_nnqp` it wraps a plain SPD array in ``DenseOperator`` and bundles
240 the knobs into an :class:`nncg.mprgp.MPRGPConfig`, then delegates to
241 :meth:`nncg.mprgp.MPRGP.solve`. The equality-augmented variant is not covered
242 — use :func:`solve_nnqp_eq` for ``B x = c``.
244 Args:
245 a: The SPD quadratic term. A :class:`cvx.linalg.SymmetricOperator` is used
246 as-is; a plain 2-D array is wrapped in ``DenseOperator`` (the
247 matrix-free ``A = M^T M + ridge I`` path is *not* inferred — pass
248 ``GramOperator(M, ridge)`` explicitly for it).
249 b: The linear term ``b``.
250 x0: Optional feasible warm start, projected onto ``x >= 0``; ``None``
251 starts from the origin.
252 tol: Relative projected-gradient stopping tolerance
253 (``MPRGPConfig.tol``).
254 gamma: Proportioning constant ``Gamma > 0`` (``MPRGPConfig.gamma``).
255 alpha_bar: Fixed projected-gradient step in ``(0, 2/||A||]``; ``None``
256 estimates ``1/||A||`` matrix-free (``MPRGPConfig.alpha_bar``).
257 max_iter: Iteration cap; ``converged=False`` when hit
258 (``MPRGPConfig.max_iter``).
259 seed: Seed of the power-iteration ``||A||`` estimate
260 (``MPRGPConfig.seed``).
262 Returns:
263 An :class:`nncg.mprgp.MPRGPResult`; ``converged`` is True iff the
264 projected gradient fell below ``tol * ||b||``, which certifies the unique
265 global minimiser.
267 Raises:
268 TypeError: When ``a`` is neither a :class:`cvx.linalg.SymmetricOperator`
269 nor an array wrappable by ``DenseOperator``.
270 ValueError: When the operator dimension does not match ``len(b)``, when
271 ``gamma`` is not strictly positive, or when ``alpha_bar`` is set but
272 not strictly positive.
274 Examples:
275 The same program as the :func:`solve_nnqp` example, reached by projection
276 instead of the active-set loop — same unique minimiser:
278 >>> import numpy as np
279 >>> a = np.array([[2.0, 0.0], [0.0, 2.0]])
280 >>> b = np.array([2.0, -2.0])
281 >>> result = solve_nnqp_mprgp(a, b)
282 >>> result.converged
283 True
284 >>> result.x.round(6).tolist()
285 [1.0, 0.0]
286 """
287 config = MPRGPConfig(tol=tol, gamma=gamma, alpha_bar=alpha_bar, max_iter=max_iter, seed=seed)
288 return MPRGP(config=config).solve(_as_operator(a), b, x0=x0)