Coverage for src/cvx/quadprog/_structure.py: 100%
28 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-09 18:50 +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(C: np.ndarray, single: np.ndarray) -> Callable[[np.ndarray], np.ndarray]:
96 """Return the cheapest available way to evaluate ``C.T @ x``.
98 This runs once per outer iteration over all ``m`` constraints, so on a
99 box-constrained problem -- where ``m = 2n`` -- the dense product is the
100 single largest cost in the solver. Three strategies, in preference order:
102 * every column a single nonzero: one gather, ``O(m)``;
103 * big enough, and sparse enough, to pay for the bookkeeping: a CSR product,
104 ``O(nnz)``;
105 * otherwise the dense product, transposed once here rather than per call.
107 The CSR branch needs both tests. Density decides which product does fewer
108 flops, and it has to be genuinely low -- see :data:`_SPARSE_DENSITY_FACTOR`,
109 which is far stricter than the compiled matvec's speed alone would suggest.
110 But below :data:`_SPARSE_MIN_WORK` the flops are not what the product costs:
111 reaching scipy's matvec takes some twenty interpreter-level calls where a
112 dense product takes one, so a small enough problem is served better densely
113 however sparse it is.
115 Args:
116 C: ``(n, m)`` constraint matrix.
117 single: Mask of the columns holding exactly one nonzero.
119 Returns:
120 A callable mapping ``x`` to ``C.T @ x``.
121 """
122 n, m = C.shape
124 if m and single.all():
125 nonzero = C != 0.0
126 row = np.argmax(nonzero, axis=0)
127 val = C[row, np.arange(m)]
128 return lambda x: val * x[row]
130 # The size test is first because it is free, where count_nonzero is O(n * m).
131 if n * m >= _SPARSE_MIN_WORK and np.count_nonzero(C) * _SPARSE_DENSITY_FACTOR <= n * m:
132 ct = scipy.sparse.csr_matrix(C.T)
133 return lambda x: ct @ x
135 dense = np.ascontiguousarray(C.T)
136 return lambda x: dense @ x