Coverage for src/nncg/api.py: 100%

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.

104

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``.

126

127 Returns:

128 A :class:`Result`; ``converged`` is True iff the KKT system was satisfied

129 to ``tol``, which certifies the unique global minimiser.

130

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.

137

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:

141

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)

154

155

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.

170

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.

176

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.

192

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``.

196

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)``, when ``b_eq`` is not of

202 shape ``(p, n)`` or ``c_eq`` not of shape ``(p,)``, or on the inner

203 solver's own conditions.

204

205 Examples:

206 The ``p = 1`` normalisation ``1^T x = 1`` — the minimum-norm point on the

207 simplex, here its centre:

208

209 >>> import numpy as np

210 >>> a = np.eye(2) * 2.0

211 >>> b = np.zeros(2)

212 >>> result = solve_nnqp_eq(a, b, np.array([[1.0, 1.0]]), np.array([1.0]))

213 >>> result.converged

214 True

215 >>> result.x.round(6).tolist()

216 [0.5, 0.5]

217 """

218 config = ActiveSetConfig(tol=tol, p_max=p_max, track=track, max_outer=max_outer)

219 solver = ActiveSetSolver(inner=_resolve_inner(inner), config=config)

220 return solver.solve_eq(_as_operator(a), b, b_eq, c_eq, warm=warm)

221

222

223def solve_nnqp_mprgp(

224 a: SymmetricOperator | NDArray[np.float64],

225 b: Vector,

226 *,

227 x0: Vector | None = None,

228 tol: float = 1e-8,

229 gamma: float = 1.0,

230 alpha_bar: float | None = None,

231 max_iter: int = 100_000,

232 seed: int = 0,

233) -> MPRGPResult:

234 """Minimise ``1/2 x^T A x - b^T x`` over ``x >= 0`` by MPRGP — one call.

235

236 The projection-based companion to :func:`solve_nnqp`: it solves the same

237 bound-constrained program with Dostál & Schöberl's MPRGP

238 (:class:`nncg.mprgp.MPRGP`) instead of the active-set loop — matrix-free and

239 factorisation-free, so it never forms or refactorises ``A``. Like

240 :func:`solve_nnqp` it wraps a plain SPD array in ``DenseOperator`` and bundles

241 the knobs into an :class:`nncg.mprgp.MPRGPConfig`, then delegates to

242 :meth:`nncg.mprgp.MPRGP.solve`. The equality-augmented variant is not covered

243 — use :func:`solve_nnqp_eq` for ``B x = c``.

244

245 Args:

246 a: The SPD quadratic term. A :class:`cvx.linalg.SymmetricOperator` is used

247 as-is; a plain 2-D array is wrapped in ``DenseOperator`` (the

248 matrix-free ``A = M^T M + ridge I`` path is *not* inferred — pass

249 ``GramOperator(M, ridge)`` explicitly for it).

250 b: The linear term ``b``.

251 x0: Optional feasible warm start, projected onto ``x >= 0``; ``None``

252 starts from the origin.

253 tol: Relative projected-gradient stopping tolerance

254 (``MPRGPConfig.tol``).

255 gamma: Proportioning constant ``Gamma > 0`` (``MPRGPConfig.gamma``).

256 alpha_bar: Fixed projected-gradient step in ``(0, 2/||A||]``; ``None``

257 estimates ``1/||A||`` matrix-free (``MPRGPConfig.alpha_bar``).

258 max_iter: Iteration cap; ``converged=False`` when hit

259 (``MPRGPConfig.max_iter``).

260 seed: Seed of the power-iteration ``||A||`` estimate

261 (``MPRGPConfig.seed``).

262

263 Returns:

264 An :class:`nncg.mprgp.MPRGPResult`; ``converged`` is True iff the

265 projected gradient fell below ``tol * ||b||``, which certifies the unique

266 global minimiser.

267

268 Raises:

269 TypeError: When ``a`` is neither a :class:`cvx.linalg.SymmetricOperator`

270 nor an array wrappable by ``DenseOperator``.

271 ValueError: When the operator dimension does not match ``len(b)``, when

272 ``gamma`` is not strictly positive, or when ``alpha_bar`` is set but

273 not strictly positive.

274

275 Examples:

276 The same program as the :func:`solve_nnqp` example, reached by projection

277 instead of the active-set loop — same unique minimiser:

278

279 >>> import numpy as np

280 >>> a = np.array([[2.0, 0.0], [0.0, 2.0]])

281 >>> b = np.array([2.0, -2.0])

282 >>> result = solve_nnqp_mprgp(a, b)

283 >>> result.converged

284 True

285 >>> result.x.round(6).tolist()

286 [1.0, 0.0]

287 """

288 config = MPRGPConfig(tol=tol, gamma=gamma, alpha_bar=alpha_bar, max_iter=max_iter, seed=seed)

289 return MPRGP(config=config).solve(_as_operator(a), b, x0=x0)