Coverage for src/nncg/mprgp.py: 100%
117 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-15 07:01 +0000
1"""MPRGP: modified proportioning with reduced gradient projections.
3A matrix-free, projection-based alternative outer solver for the same
4strictly convex non-negative quadratic program the active-set loop targets,
6 min_{x >= 0} 1/2 x^T A x - b^T x, A symmetric positive definite.
8Where :class:`nncg.solver.ActiveSetSolver` toggles a working set and solves an
9unconstrained system on each free block, MPRGP (Dostál & Schöberl, 2005) never
10factorises anything: it interleaves three cheap first-order moves, each costing
11one or two Hessian products,
13* a **conjugate-gradient step** that minimises within the current face while it
14 stays feasible (the free set unchanged),
15* an **expansion step** that walks to the nearest bound and takes one fixed-step
16 projected-gradient move to *add* constraints to the active set, and
17* a **proportioning step** along the chopped gradient that *removes* constraints
18 from the active set,
20switched by the proportioning test ``||beta(x)||^2 <= gamma^2 phi~(x)^T phi(x)``.
21With the projected-gradient step bounded by ``alpha_bar in (0, 2/||A||]`` the
22iteration converges for any feasible start, and — because it identifies the
23active set of the minimiser in finitely many steps and then reduces to plain CG
24on the optimal face — it terminates finitely in exact arithmetic. ``A`` enters
25only through :meth:`cvx.linalg.SymmetricOperator.matvec`, so the ``n x n`` matrix
26is never formed; ``||A||`` for the step bound is estimated matrix-free by power
27iteration.
29This is the bound-constrained solver; the equality-augmented variant ``B x = c``
30is out of scope here (it needs an augmented-Lagrangian outer wrap, SMALBE/SMALSE
31around MPRGP) — use :meth:`nncg.solver.ActiveSetSolver.solve_eq` for that.
33Reference: Z. Dostál and J. Schöberl, "Minimizing quadratic functions subject to
34bound constraints with the rate of convergence and finite termination",
35Comput. Optim. Appl. 30 (2005), 23-43.
36"""
38from __future__ import annotations
40from collections.abc import Callable
41from dataclasses import dataclass
43import numpy as np
44from cvx.linalg import SymmetricOperator, Vector, power_iteration
45from numpy.typing import NDArray
47from .certificate import _require_operator
49MatVec = Callable[[Vector], Vector]
50"""The action ``v -> A v`` of the SPD operator — MPRGP's only access to ``A``."""
52Iterate = tuple[Vector, Vector, Vector]
53"""The MPRGP iteration state ``(x, g, p)``: iterate, gradient ``A x - b``, CG direction.
55Every move consumes one state and returns the next, so the loop in :func:`_mprgp`
56carries no other mutable numerics — only the counters the result reports.
57"""
60def _free_gradient(x: Vector, g: Vector) -> Vector:
61 """Free gradient ``phi``: the gradient on the free set, zero on the active set.
63 ``phi_i = g_i`` where ``x_i > 0`` (free) and ``0`` where ``x_i = 0`` (active).
64 It drives the conjugate-gradient and expansion steps, which move only the
65 free variables.
66 """
67 return np.where(x > 0.0, g, 0.0)
70def _chopped_gradient(x: Vector, g: Vector) -> Vector:
71 """Chopped gradient ``beta``: the releasing part of the gradient on the active set.
73 ``beta_i = min(g_i, 0)`` where ``x_i = 0`` (active) and ``0`` where
74 ``x_i > 0`` (free). A negative gradient at a bound means the objective still
75 decreases along the feasible ``+e_i`` direction, so the proportioning step
76 follows ``beta`` to release such constraints.
77 """
78 return np.where(x > 0.0, 0.0, np.minimum(g, 0.0))
81def _reduced_free_gradient(x: Vector, g: Vector, alpha_bar: float) -> Vector:
82 """Reduced free gradient ``phi~``: the free gradient capped by the feasible step.
84 ``phi~_i = min(g_i, x_i / alpha_bar)`` on the free set and ``0`` on the active
85 set. It measures the decrease a single ``alpha_bar`` projected-gradient step
86 can realise on each free variable (a downhill move is limited by the distance
87 ``x_i`` to the bound), and enters only the proportioning test.
88 """
89 return np.where(x > 0.0, np.minimum(g, x / alpha_bar), 0.0)
92def _max_feasible_step(x: Vector, p: Vector) -> float:
93 """Largest ``alpha`` with ``x - alpha p >= 0``: ``min_{p_i > 0} x_i / p_i``.
95 Only components that decrease (``p_i > 0``) can reach the bound; when none do,
96 the whole ray is feasible and the step is unbounded (``+inf``).
97 """
98 decreasing = p > 0.0
99 if not decreasing.any():
100 return np.inf
101 return float(np.min(x[decreasing] / p[decreasing]))
104def _initial_iterate(matvec: MatVec, b: Vector, x0: Vector | None) -> Iterate:
105 """Project the warm start onto ``x >= 0`` and seed the gradient and CG direction.
107 Args:
108 matvec: The action ``v -> A v`` of the SPD operator.
109 b: The linear term ``b``.
110 x0: Optional warm start, projected onto the feasible set; ``None`` starts
111 at the origin.
113 Returns:
114 The initial ``(x, g, p)``. Costs the one Hessian product of ``g = A x - b``.
115 """
116 x = np.zeros_like(b) if x0 is None else np.maximum(np.asarray(x0, dtype=np.float64), 0.0)
117 g = matvec(x) - b
118 return x, g, _free_gradient(x, g)
121def _stopping_threshold(b: Vector, tol: float) -> float:
122 """Absolute exit threshold ``tol * ||b||`` of the projected-gradient stopping test.
124 ``||b||`` is replaced by ``1`` when ``b = 0``, so the test stays meaningful on
125 the degenerate problem whose minimiser is the origin.
126 """
127 return tol * (float(np.linalg.norm(b)) or 1.0)
130def _cg_step(x: Vector, g: Vector, p: Vector, ap: Vector, p_ap: float, alpha: float) -> Iterate:
131 """Conjugate-gradient step: minimise within the current face, leaving it unchanged.
133 The caller has already checked ``alpha`` feasible, so no projection is needed
134 and the gradient updates linearly. The next direction is the new free gradient
135 made ``A``-conjugate to ``p`` by one Gram-Schmidt sweep. Costs no Hessian
136 product of its own — ``ap`` is the one the caller computed to form ``alpha``.
137 """
138 x = x - alpha * p
139 g = g - alpha * ap
140 phi = _free_gradient(x, g)
141 beta_gs = float(phi @ ap) / p_ap
142 return x, g, phi - beta_gs * p
145def _expansion_step(
146 matvec: MatVec,
147 b: Vector,
148 x: Vector,
149 g: Vector,
150 p: Vector,
151 ap: Vector,
152 alpha_f: float,
153 alpha_bar: float,
154) -> Iterate:
155 """Expansion step: walk to the nearest bound, then one projected-gradient move.
157 Taken when the conjugate-gradient step would leave the feasible set. The
158 iterate advances only as far as the bound (``alpha_f``) and then takes one
159 fixed-step ``alpha_bar`` projected-gradient move, which *adds* the newly
160 active constraints. The projection is non-linear, so the gradient is
161 recomputed rather than updated — the extra Hessian product this step costs —
162 and CG restarts from the new free gradient.
163 """
164 x_half = x - alpha_f * p
165 g = g - alpha_f * ap
166 phi_half = _free_gradient(x_half, g)
167 x = np.maximum(x_half - alpha_bar * phi_half, 0.0)
168 g = matvec(x) - b
169 return x, g, _free_gradient(x, g)
172def _proportioning_step(matvec: MatVec, x: Vector, g: Vector, beta: Vector) -> Iterate:
173 """Proportioning step: release active constraints along the chopped gradient.
175 Taken on a disproportional iterate. ``beta`` is non-positive on the active set
176 and zero on the free set, and the exact line minimiser along it has
177 ``alpha >= 0``, so the move points inward and stays feasible without a
178 projection — the gradient therefore updates linearly. Costs one Hessian
179 product; CG restarts afterwards.
180 """
181 ad = matvec(beta)
182 alpha = float(g @ beta) / float(beta @ ad)
183 x = x - alpha * beta
184 g = g - alpha * ad
185 return x, g, _free_gradient(x, g)
188def _proportional_step(matvec: MatVec, b: Vector, iterate: Iterate, alpha_bar: float) -> tuple[Iterate, int, str]:
189 """Take the conjugate-gradient or expansion move on a proportional iterate.
191 Forms the trial conjugate-gradient step and compares it with the largest
192 feasible step along ``p``: within the face it is a plain :func:`_cg_step`,
193 otherwise the face must change and :func:`_expansion_step` walks to the bound
194 and projects.
196 Args:
197 matvec: The action ``v -> A v`` of the SPD operator.
198 b: The linear term ``b``.
199 iterate: The current ``(x, g, p)``.
200 alpha_bar: The fixed projected-gradient step of the expansion move.
202 Returns:
203 The next ``(x, g, p)``, the Hessian products consumed (one for the
204 conjugate-gradient step, two for the expansion step), and which move was
205 taken — ``"cg"`` or ``"expansion"`` — for the result's step counters.
206 """
207 x, g, p = iterate
208 ap = matvec(p)
209 p_ap = float(p @ ap)
210 alpha_cg = float(g @ p) / p_ap
211 alpha_f = _max_feasible_step(x, p)
212 if alpha_cg <= alpha_f:
213 return _cg_step(x, g, p, ap, p_ap, alpha_cg), 1, "cg"
214 return _expansion_step(matvec, b, x, g, p, ap, alpha_f, alpha_bar), 2, "expansion"
217def _resolve_alpha_bar(a: SymmetricOperator, alpha_bar: float | None, seed: int) -> float:
218 """Return the fixed projected-gradient step, estimating ``1/||A||`` when unset.
220 The convergence proof requires ``alpha_bar in (0, 2/||A||]``. Power iteration
221 approaches ``||A|| = lambda_max`` from below, so ``2/lambda_est`` could exceed
222 the bound; the conservative default ``1/lambda_est`` stays safely inside it for
223 any estimate with ``lambda_est >= lambda_max / 2``.
225 Args:
226 a: The SPD operator ``A``.
227 alpha_bar: An explicit step (returned as-is after validation), or ``None``
228 for the ``1/lambda_max`` estimate.
229 seed: Seed of the power-iteration test vector, for reproducibility.
231 Returns:
232 The positive fixed step ``alpha_bar``.
234 Raises:
235 ValueError: When an explicit ``alpha_bar`` is not strictly positive.
236 """
237 if alpha_bar is not None:
238 if alpha_bar <= 0.0:
239 msg = f"alpha_bar must be strictly positive; got {alpha_bar:.2e}"
240 raise ValueError(msg)
241 return float(alpha_bar)
242 lam_max, _ = power_iteration(a, seed=seed)
243 return 1.0 / float(lam_max)
246@dataclass(frozen=True)
247class MPRGPConfig:
248 """Configuration of the MPRGP solver (:class:`MPRGP`).
250 Attributes:
251 tol: Relative stopping tolerance on the projected gradient — the loop
252 exits when ``||phi(x) + beta(x)|| <= tol * ||b||`` (``||b||`` replaced
253 by ``1`` when ``b = 0``), which certifies the KKT conditions.
254 gamma: Proportioning constant ``Gamma > 0`` balancing expansion against
255 proportioning. ``1.0`` is the standard, near-optimal choice; larger
256 values expand more eagerly, smaller ones release more eagerly.
257 alpha_bar: Fixed projected-gradient step, which must satisfy
258 ``0 < alpha_bar <= 2/||A||`` for convergence. ``None`` estimates the
259 safe default ``1/||A||`` by matrix-free power iteration.
260 max_iter: Iteration cap; the current iterate is returned with
261 ``converged=False`` when it is hit.
262 seed: Seed of the power-iteration ``||A||`` estimate (only used when
263 ``alpha_bar is None``), fixed so a solve is reproducible.
265 Raises:
266 ValueError: When ``gamma`` is not strictly positive.
267 """
269 tol: float = 1e-8
270 gamma: float = 1.0
271 alpha_bar: float | None = None
272 max_iter: int = 100_000
273 seed: int = 0
275 def __post_init__(self) -> None:
276 """Validate that the proportioning constant is strictly positive."""
277 if self.gamma <= 0.0:
278 msg = f"gamma must be strictly positive; got {self.gamma:.2e}"
279 raise ValueError(msg)
282#: Shared default so ``MPRGPConfig()`` is not called in argument defaults (ruff B008).
283_DEFAULT_MPRGP = MPRGPConfig()
286@dataclass(frozen=True)
287class MPRGPResult:
288 """Outcome of an MPRGP solve.
290 The iteration counts are broken out by move because they carry the algorithm's
291 signature: expansion and proportioning steps are the ones that change the
292 active set, while a run of conjugate-gradient steps is plain CG on a fixed
293 face. ``hessian_products`` is the honest cost of a matrix-free method — one
294 product per CG or proportioning step, two per expansion step, plus one for the
295 initial gradient.
297 Attributes:
298 x: The minimiser (or the final iterate if ``converged`` is False).
299 iterations: Total MPRGP steps taken (the sum of the three move counts).
300 hessian_products: Number of operator matrix-vector products consumed.
301 cg_steps: Conjugate-gradient (minimisation-within-the-face) steps.
302 expansion_steps: Expansion (bound-hitting projected-gradient) steps.
303 proportioning_steps: Proportioning (constraint-releasing) steps.
304 converged: True when the projected-gradient stopping test was met; False
305 when ``max_iter`` stopped the loop first.
306 free: Boolean mask of the final free set (``x > 0``).
307 """
309 x: Vector
310 iterations: int
311 hessian_products: int
312 cg_steps: int
313 expansion_steps: int
314 proportioning_steps: int
315 converged: bool
316 free: NDArray[np.bool_]
319@dataclass(frozen=True)
320class MPRGP:
321 """The MPRGP solver for the non-negative quadratic program.
323 A matrix-free, factorisation-free alternative to
324 :class:`nncg.solver.ActiveSetSolver` on the bound-constrained problem
325 ``min_{x>=0} 1/2 x^T A x - b^T x``: it interleaves conjugate-gradient,
326 expansion and proportioning steps under the proportioning test, never forming
327 or factorising ``A``. Holds only its :class:`MPRGPConfig`; the operator and
328 right-hand side are passed to :meth:`solve`.
330 Attributes:
331 config: Solver configuration (tolerance, proportioning constant,
332 projected-gradient step, iteration cap, seed).
334 Examples:
335 The same operator interface as the active-set loop:
337 >>> import numpy as np
338 >>> from cvx.linalg import DenseOperator
339 >>> from nncg import MPRGP, MPRGPConfig, kkt_violation
340 >>> a = DenseOperator(np.array([[2.0, 0.0], [0.0, 2.0]]))
341 >>> b = np.array([2.0, -2.0])
342 >>> res = MPRGP().solve(a, b)
343 >>> res.converged
344 True
345 >>> bool(np.allclose(res.x, [1.0, 0.0]))
346 True
347 >>> round(kkt_violation(a, b, res.x), 12)
348 0.0
350 The counts break the run down by move, and always sum to
351 ``iterations``; ``hessian_products`` is the honest matrix-free cost —
352 at least one product per step, plus the initial gradient:
354 >>> res.iterations == res.cg_steps + res.expansion_steps + res.proportioning_steps
355 True
356 >>> res.hessian_products > res.iterations
357 True
359 Nothing is factorised, so the whole configuration is a handful of
360 scalars — here a tighter tolerance and an explicit projected-gradient
361 step, which skips the power-iteration estimate of ``1/||A||``:
363 >>> tight = MPRGP(config=MPRGPConfig(tol=1e-12, alpha_bar=0.5)).solve(a, b)
364 >>> tight.converged
365 True
366 >>> bool(np.allclose(tight.x, [1.0, 0.0]))
367 True
368 """
370 config: MPRGPConfig = _DEFAULT_MPRGP
372 def solve(self, a: SymmetricOperator, b: Vector, x0: Vector | None = None) -> MPRGPResult:
373 """Minimise ``1/2 x^T A x - b^T x`` over ``x >= 0`` by MPRGP.
375 Args:
376 a: The SPD operator ``A`` (a :class:`cvx.linalg.SymmetricOperator`) —
377 ``DenseOperator`` for an explicit array, ``GramOperator(M, ridge)``
378 for ``A = M^T M + ridge I`` whose Gram matrix is never formed.
379 b: The linear term ``b``.
380 x0: Optional feasible warm start; it is projected onto ``x >= 0`` and
381 the iteration begins there. ``None`` starts from the origin.
383 Returns:
384 An :class:`MPRGPResult`; ``converged`` is True iff the projected
385 gradient fell below ``config.tol * ||b||``, which certifies the unique
386 global minimiser.
388 Raises:
389 TypeError: When ``a`` is not a :class:`cvx.linalg.SymmetricOperator`.
390 ValueError: When the operator dimension does not match ``len(b)``, or
391 when ``config.alpha_bar`` is set but not strictly positive.
392 """
393 _require_operator(a, b)
394 cfg = self.config
395 alpha_bar = _resolve_alpha_bar(a, cfg.alpha_bar, cfg.seed)
396 return _mprgp(a.matvec, b, x0, alpha_bar, cfg.gamma, cfg.tol, cfg.max_iter)
399def _mprgp(
400 matvec: MatVec,
401 b: Vector,
402 x0: Vector | None,
403 alpha_bar: float,
404 gamma: float,
405 tol: float,
406 max_iter: int,
407) -> MPRGPResult:
408 """Run the MPRGP iteration and assemble its :class:`MPRGPResult`.
410 The pure algorithm behind :meth:`MPRGP.solve`: it takes the resolved operator
411 action and step bound and owns the loop — the projected-gradient stopping
412 test and the proportioning switch — while the three moves themselves live in
413 :func:`_cg_step`, :func:`_expansion_step` and :func:`_proportioning_step`,
414 with the first two selected by :func:`_proportional_step`. Kept as a plain
415 function (operator access reduced to ``matvec``) so the numerics can be
416 exercised directly.
418 Args:
419 matvec: The action ``v -> A v`` of the SPD operator.
420 b: The linear term ``b``.
421 x0: Optional warm start (projected onto ``x >= 0``); ``None`` starts at 0.
422 alpha_bar: The fixed projected-gradient step, ``0 < alpha_bar <= 2/||A||``.
423 gamma: The proportioning constant ``Gamma > 0``.
424 tol: Relative projected-gradient stopping tolerance (against ``||b||``).
425 max_iter: Iteration cap.
427 Returns:
428 The completed :class:`MPRGPResult`.
429 """
430 x, g, p = _initial_iterate(matvec, b, x0)
431 products = 1 # the gradient A x - b of the initial iterate
432 stop = _stopping_threshold(b, tol)
433 counts = {"cg": 0, "expansion": 0, "proportioning": 0}
434 iterations = 0
435 converged = False
437 while iterations < max_iter:
438 phi = _free_gradient(x, g)
439 beta = _chopped_gradient(x, g)
440 if float(np.linalg.norm(phi + beta)) <= stop:
441 converged = True
442 break
443 iterations += 1
445 phi_tilde = _reduced_free_gradient(x, g, alpha_bar)
446 if float(beta @ beta) <= gamma * gamma * float(phi_tilde @ phi):
447 # Proportional iterate: minimise within the face, or expand it.
448 (x, g, p), used, move = _proportional_step(matvec, b, (x, g, p), alpha_bar)
449 products += used
450 counts[move] += 1
451 else:
452 # Disproportional iterate: release constraints along the chopped gradient.
453 x, g, p = _proportioning_step(matvec, x, g, beta)
454 products += 1
455 counts["proportioning"] += 1
457 return MPRGPResult(
458 x=x,
459 iterations=iterations,
460 hessian_products=products,
461 cg_steps=counts["cg"],
462 expansion_steps=counts["expansion"],
463 proportioning_steps=counts["proportioning"],
464 converged=converged,
465 free=x > 0.0,
466 )