Coverage for src/cvx/quadprog/_structure.py: 100%
25 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 06:10 +0000
1"""Exploiting the shape of the constraint matrix.
3A bound constraint is a column holding a single nonzero, and a box is nothing
4else. Detecting that turns three per-iteration products into indexing, and
5choosing the slack product by size and density decides the rest.
6"""
8# G, C, R and J are the names used in Goldfarb & Idnani (1983) and in the
9# reference implementation's public signature `solve_qp(G, a, C, b, meq)`.
10# Lowercasing them would obscure the correspondence to the paper, so the
11# pep8-naming rules are waived here, as they are in _solve.py.
12# ruff: noqa: N803, TRY003
14from collections.abc import Callable
16import numpy as np
17import scipy.sparse
19# Size below which the slack product is dense even when C is sparse enough to
20# favour CSR on flops alone. scipy's CSR matvec is compiled, but reaching it
21# costs some twenty interpreter-level calls per product -- isinstance and abc
22# checks, sputils lookups, allocating the result -- against the single matmul a
23# dense product takes. That overhead is roughly fixed while the dense product
24# grows as n * m, so below some size it is not worth paying however sparse the
25# matrix is. Measured end to end on budget-plus-bounds problems, dense wins by
26# 1.34x at n = 10 and 1.10x at n = 100, ties at n = 400 (n * m = 320_400), and
27# loses from n = 450 (n * m = 405_450) on, reaching 0.62x by n = 1200.
28_SPARSE_MIN_WORK = 350_000
30# Reciprocal of the density above which the dense product wins outright, whatever
31# the size: CSR is used only when nnz * _SPARSE_DENSITY_FACTOR <= n * m. Measured
32# per product over n * m from 80_000 to 4_500_000, CSR wins at and below 2%
33# density at every size and loses at 5% for all but the largest, so the crossover
34# sits near 3-4% and moves little with size -- the two costs are both linear, in
35# nnz and in n * m respectively, so their ratio is what decides.
36_SPARSE_DENSITY_FACTOR = 25
39def _default_constraints(
40 G: np.ndarray, C: np.ndarray | None, b: np.ndarray | None, meq: int
41) -> tuple[np.ndarray, np.ndarray, int]:
42 """Fill in the unconstrained problem and coerce the constraint arrays.
44 Omitting both ``C`` and ``b`` asks for the unconstrained minimum. Rather
45 than branch on that everywhere below, it is expressed as a single constraint
46 ``0 >= -1``, which no ``x`` can violate: the solver then runs its ordinary
47 path and terminates on the first iteration. Supplying exactly one of the two
48 is an error rather than a shape crash further in.
50 Args:
51 G: ``(n, n)`` quadratic term, used only for its size.
52 C: ``(n, m)`` constraint matrix, or None.
53 b: ``(m,)`` right-hand side, or None.
54 meq: Number of leading constraints to treat as equalities.
56 Returns:
57 ``C``, ``b`` as float64 arrays and the ``meq`` that goes with them --
58 forced to 0 when the unconstrained placeholder is substituted, since
59 that constraint must not be read as an equality.
61 Raises:
62 ValueError: If exactly one of ``C`` and ``b`` is given.
63 """
64 if C is None and b is None:
65 return np.zeros((len(G), 1)), -np.ones(1), 0
66 if C is None or b is None:
67 raise ValueError("C and b must be given together")
68 return np.asarray(C, dtype=np.float64), np.asarray(b, dtype=np.float64), meq
71def _analyse_constraints(C: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
72 """Locate the columns of ``C`` that hold a single scaled unit vector.
74 Bound constraints -- the overwhelmingly common shape, and what
75 ``C = [I, -I]`` is -- make every column one nonzero. Recognising that lets
76 the products against a constraint normal become scalar indexing rather than
77 length-``n`` reductions.
79 Args:
80 C: ``(n, m)`` constraint matrix, one column per constraint.
82 Returns:
83 A boolean mask of the single-nonzero columns, the row index of that
84 nonzero per column, and its value. Entries of the latter two are
85 meaningless where the mask is False.
86 """
87 nonzero = C != 0.0
88 single = nonzero.sum(axis=0) == 1
89 # argmax gives the first nonzero row, which for a single-nonzero column is
90 # the only one. Columns failing the mask still index safely, just uselessly.
91 row = np.argmax(nonzero, axis=0)
92 return single, row, C[row, np.arange(C.shape[1])]
95def _slack_evaluator(
96 C: np.ndarray, single: np.ndarray, srow: np.ndarray, sval: np.ndarray
97) -> Callable[[np.ndarray], np.ndarray]:
98 """Return the cheapest available way to evaluate ``C.T @ x``.
100 This runs once per outer iteration over all ``m`` constraints, so on a
101 box-constrained problem -- where ``m = 2n`` -- the dense product is the
102 single largest cost in the solver. Three strategies, in preference order:
104 * every column a single nonzero: one gather, ``O(m)``;
105 * big enough, and sparse enough, to pay for the bookkeeping: a CSR product,
106 ``O(nnz)``;
107 * otherwise the dense product, transposed once here rather than per call.
109 The CSR branch needs both tests. Density decides which product does fewer
110 flops, and it has to be genuinely low -- see :data:`_SPARSE_DENSITY_FACTOR`,
111 which is far stricter than the compiled matvec's speed alone would suggest.
112 But below :data:`_SPARSE_MIN_WORK` the flops are not what the product costs:
113 reaching scipy's matvec takes some twenty interpreter-level calls where a
114 dense product takes one, so a small enough problem is served better densely
115 however sparse it is.
117 Args:
118 C: ``(n, m)`` constraint matrix.
119 single: Mask of the columns holding exactly one nonzero.
120 srow: Row index of that nonzero per column, from
121 :func:`_analyse_constraints`.
122 sval: Value of that nonzero per column, from the same place.
124 Returns:
125 A callable mapping ``x`` to ``C.T @ x``. Every branch returns a freshly
126 allocated array, which callers rely on: both this module's callers go on
127 to force entries of it to zero in place.
128 """
129 n, m = C.shape
131 if m and single.all():
132 # `srow` and `sval` are taken from the caller rather than rebuilt here.
133 # This used to recompute both, plus the n x m boolean `C != 0.0` they come
134 # from, on the same C the caller had just analysed -- two extra O(n * m)
135 # passes and a discarded n x m temporary per solve, measured at 0.70 ms of
136 # a 14.5 ms solve at n = 800 and 2.21 ms at n = 1400 (#108).
137 return lambda x: sval * x[srow]
139 # The size test is first because it is free, where count_nonzero is O(n * m).
140 if n * m >= _SPARSE_MIN_WORK and np.count_nonzero(C) * _SPARSE_DENSITY_FACTOR <= n * m:
141 ct = scipy.sparse.csr_matrix(C.T)
142 return lambda x: ct @ x
144 dense = np.ascontiguousarray(C.T)
145 return lambda x: dense @ x