API Reference¶

class Solution(NamedTuple):
    """The outcome of a quadratic program.

    Iterating over an instance yields the same six values, in the same order, as
    the tuple returned by ``quadprog.solve_qp``, so it is a drop-in replacement.

    That ordering is what the example below pins down, on a problem small enough
    to check by hand: with ``G`` the identity the objective separates into
    ``x_i**2 / 2 - a_i x_i`` per coordinate, so subject to ``x >= 0`` the answer
    is ``a`` with its negative entries clipped to zero.

    >>> import numpy as np
    >>> from cvx.quadprog import solve_qp
    >>> solution = solve_qp(np.eye(2), np.array([1.0, -1.0]), np.eye(2), np.zeros(2))
    >>> x, f, xu, iterations, lagrangian, iact = solution
    >>> x.tolist()
    [1.0, 0.0]
    >>> xu.tolist()
    [1.0, -1.0]
    >>> round(f, 12)
    -0.5

    Only the second constraint binds, so ``iact`` names it -- 1-based -- and only
    its multiplier is non-zero. ``iterations`` is a pair, additions then removals,
    and its first entry counts *outer iterations*: it therefore exceeds the number
    of constraints that ended up active, one iteration having brought ``x_2`` onto
    its bound and a final one confirming there was nothing left to add.

    >>> iact.tolist()
    [2]
    >>> lagrangian.tolist()
    [0.0, 1.0]
    >>> iterations.tolist()
    [2, 0]

    Attributes:
        x: ``(n,)`` minimiser of the constrained problem.
        f: Value of the objective at ``x``.
        xu: ``(n,)`` minimiser of the unconstrained problem, ``G^-1 a``.
        iterations: ``(2,)`` count of constraints added to the active set (once
            per outer iteration) and of constraints removed from it.
        lagrangian: ``(m,)`` Lagrange multipliers, zero for inactive
            constraints.
        iact: 1-based indices of the constraints active at the solution.
    """

    x: np.ndarray
    f: float
    xu: np.ndarray
    iterations: np.ndarray
    lagrangian: np.ndarray
    iact: np.ndarray

class Sweep:
    """Solve a family of QPs sharing ``G``, ``C``, ``b`` and ``meq``.

    Only the linear term changes between calls. The first call solves from
    scratch; later ones reuse the factorisation when the active set still holds.

    >>> import numpy as np
    >>> from cvx.quadprog import Sweep, solve_qp
    >>> G = np.eye(3)
    >>> C = np.array([[-4.0, 2.0, 0.0], [-3.0, 1.0, -2.0], [0.0, 0.0, 1.0]])
    >>> b = np.array([-8.0, 2.0, 0.0])
    >>> sweep = Sweep(G, C, b)
    >>> a = np.array([0.0, 5.0, 0.0])
    >>> bool(np.allclose(sweep.solve(a).x, solve_qp(G, a, C, b).x))
    True
    >>> bool(np.allclose(sweep.solve(1.01 * a).x, solve_qp(G, 1.01 * a, C, b).x))
    True
    """

    def __init__(
        self,
        G: np.ndarray,
        C: np.ndarray | None = None,
        b: np.ndarray | None = None,
        meq: int = 0,
        check_finite: bool = False,
        blas_threads: int | None = None,
    ) -> None:
        """Fix the part of the problem that does not vary.

        Args:
            G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
            C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
                the unconstrained problem.
            b: ``(m,)`` right-hand side of the constraints.
            meq: Number of leading constraints to treat as equalities.
            check_finite: Whether to reject NaN and infinity in ``G``, ``C`` and
                ``b``, and in each ``a`` passed to :meth:`solve`. Off by default,
                matching :func:`~cvx.quadprog.solve_qp`.
            blas_threads: Cap the BLAS thread count for the expensive parts of this
                sweep, as :func:`~cvx.quadprog.solve_qp`'s argument of the same name
                does for one solve: the factorisation below, and every
                :meth:`solve` that misses the cache. A hit is deliberately left
                outside the context, which costs ~100 microseconds against an
                ``O(n^2)`` recovery -- and against the handful of microseconds a
                hit actually takes at the small ``n`` this class is most worth
                using at, where that arithmetic is still below the dispatch
                overhead.

                Decided once here rather than per call, because ``n`` is fixed for
                this object's lifetime and so the automatic gate's answer is too.
                Left unset, that gate is consulted exactly as it is for
                ``solve_qp`` -- see there for what it does and does not change, and
                :func:`~cvx.quadprog._threads.auto_cap_threads` for the conditions.

        Raises:
            ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
                if ``G`` is not positive definite, or if ``blas_threads`` is not at
                least 1.
            ImportError: If ``blas_threads`` is given and ``threadpoolctl`` is not
                installed.
        """
        G = np.asarray(G, dtype=np.float64)
        self.C, self.b, self.meq = _default_constraints(G, C, b, meq)
        self.n, self._q = _validate(G, np.zeros(len(G)), self.C, self.b, self.meq, check_finite)
        self._check_finite = check_finite
        self.G = G

        # C, b and meq are fixed for this object's lifetime, so the shape analysis
        # runs once here and is amortised over every call -- where solve_qp has to
        # pay it per solve. Before this, the hit path re-derived the slacks with a
        # dense `C.T @ x` and never reached the bound-constraint gather at all,
        # which on a box family is 13% of a hit at n = 800 (#109).
        self._single, self._srow, self._sval = _analyse_constraints(self.C)
        self._slack_of = _slack_evaluator(self.C, self._single, self._srow, self._sval)

        # An explicit count is used as given; otherwise the automatic gate decides,
        # and it is asked once because `n` cannot change under it. None means "leave
        # the process alone", which is what `scoped_limit` turns into a no-op.
        #
        # Sweep is the API most exposed to the OpenBLAS collapse -- large problems,
        # solved repeatedly -- and until #107 it was the one path with no guard,
        # because it calls `_solve_with_factors` below the level solve_qp installs
        # the context at.
        self._blas_threads = blas_threads if blas_threads is not None else _threads.auto_cap_threads(self.n, fast=False)

        # The Cholesky is a property of G alone, so it is done once here and every
        # later cold solve is handed the factor instead, via `factorized=True`.
        # That is the same reuse the reference package offers, and it is the part
        # of the saving that applies even when the active set does change.
        #
        # It is also the single largest BLAS call this object ever makes, at
        # O(n^3), so it is inside the cap. A bad `blas_threads` therefore raises
        # here, at construction, rather than at the first solve.
        with _threads.scoped_limit(self._blas_threads):
            self._Rinv, _xu = _factorize(G, np.zeros(self.n), False)
        self._cache: _Cache | None = None
        self.hits = 0
        self.misses = 0

    def solve(self, a: np.ndarray) -> Solution:
        """Solve for a new linear term.

        Args:
            a: ``(n,)`` vector of the linear term.

        Returns:
            The same :class:`~cvx.quadprog.Solution` that
            :func:`~cvx.quadprog.solve_qp` would return for this problem, except
            that ``iterations`` is ``(0, 0)`` when the cached factorisation was
            reused outright -- no active-set iteration was performed.

        Raises:
            ValueError: If ``a`` has the wrong shape, if the constraints admit no
                solution, or if ``check_finite`` is set and ``a`` holds a
                non-finite value.
        """
        a = np.asarray(a, dtype=np.float64)
        warm = None
        cache = self._cache
        if cache is not None and self._usable(a):
            hit = self._reuse(a, cache)
            if hit is not None:
                self.hits += 1
                return hit
            # The cached set is stale, but it is still a far better place to start
            # than the unconstrained minimum: repairing it into a dual-feasible
            # state costs a few drops, where a cold solve re-walks the whole set.
            warm = self._repair(a, cache)

        self.misses += 1
        # Only the miss is wrapped. A hit is an O(n^2) recovery plus a KKT check,
        # and entering a threadpoolctl context costs ~100 microseconds, which would
        # be a tax on exactly the path this class exists to make cheap.
        with _threads.scoped_limit(self._blas_threads):
            solution, J, R = _solve_with_factors(
                self._Rinv, a, self.C, self.b, self.meq, True, self._check_finite, warm
            )
        self._cache = _Cache(J, R, solution.iact)
        return solution

    def _usable(self, a: np.ndarray) -> bool:
        """Return whether the cache may be consulted at all for this ``a``.

        Args:
            a: ``(n,)`` vector of the linear term.

        Returns:
            False when ``a`` is the wrong length, or when ``check_finite`` is set
            and it is not finite -- every KKT comparison against NaN is False, so
            without this the fast path would *accept* a non-finite point instead
            of rejecting it. Falling back lets the full solve raise, which is what
            the caller asked for.
        """
        if len(a) != self.n:
            return False
        return not (self._check_finite and not np.isfinite(a).all())

    def _recover(
        self, a: np.ndarray, J: np.ndarray, R: np.ndarray, iact: np.ndarray, nact: int
    ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Return the minimiser over a given active set, and its multipliers.

        Costs ``O(n^2)``, all of it in ``x_u = J J^T a`` below: the factors already
        encode everything about ``G`` and the active constraints, so only the
        right-hand side has changed, but ``J`` is a dense ``(n, n)`` and those two
        products do not shrink with the active set. The rest is ``O(nk + k^2)``.

        Args:
            a: ``(n,)`` linear term.
            J: Inverse Cholesky factor for this active set.
            R: Packed triangular factor for this active set.
            iact: 1-based active set, first ``nact`` entries valid.
            nact: Size of the active set.

        Returns:
            ``(x, lam, xu)`` -- the minimiser subject to the active set held as
            equalities, its multipliers, and the unconstrained minimiser.
        """
        xu = J @ (J.T @ a)
        if nact == 0:
            # Distinct arrays even though the values coincide: a resumed solve
            # updates the iterate in place, and would otherwise corrupt ``xu``
            # along with it. The cold path copies here for the same reason.
            return xu.copy(), _EMPTY, xu
        active = iact[:nact] - 1
        y = dtpsv(nact, R, self.b[active] - self._active_product(active, xu), lower=0, trans=1, overwrite_x=True)
        x = xu + J[:, :nact] @ y
        lam = dtpsv(nact, R, y.copy(), lower=0, trans=0, overwrite_x=True)
        return x, lam, xu

    def _active_product(self, active: np.ndarray, xu: np.ndarray) -> np.ndarray:
        """Return ``C_A^T xu`` for the active columns, by gather where it can.

        Where every active column holds a single nonzero the product is ``k``
        multiplications (#109). The test is on the *active* columns rather than on
        all of ``C``, so a mixed matrix -- a budget row plus bounds -- still takes
        that path whenever the set happens to be all bounds. It is ``O(k)``
        against what it guards.

        What it guards is no longer a block of ``C``. Fancy-indexing an ``(n, k)``
        block out and multiplying against it costs ``O(nk)`` in flops but a copy
        of the block in bandwidth, and the copy is what dominates: at
        ``n = 800``, ``m = 400``, ``k = 50`` it measured 0.024 ms against 0.005 ms
        for evaluating all ``m`` products and keeping ``k`` of them, a product the
        evaluator of :mod:`._structure` has already chosen the cheapest form for.
        Doing the arithmetic for constraints whose answers are then discarded is
        the faster route by a factor of five, and on a reused solve of a
        dense-``C`` family it was 54% of the whole cost.

        Args:
            active: 0-based indices of the active constraints.
            xu: ``(n,)`` unconstrained minimiser.

        Returns:
            The length-``k`` vector of active constraint values at ``xu``.
        """
        # Annotated on the way out for the reason given in _threads.limit: indexing
        # an ndarray by an ndarray is typed as Any, so returning either expression
        # directly is an untyped escape under --strict.
        if self._single[active].all():
            gathered: np.ndarray = self._sval[active] * xu[self._srow[active]]
            return gathered
        product: np.ndarray = self._slack_of(xu)[active]
        return product

    def _reuse(self, a: np.ndarray, cache: "_Cache") -> Solution | None:
        """Return the solution from the cached factorisation, or None if it is stale.

        Args:
            a: ``(n,)`` vector of the linear term.
            cache: The factorisation a previous solve ended on.

        Returns:
            A :class:`~cvx.quadprog.Solution` when the cached active set still
            satisfies the KKT conditions for this ``a``, otherwise None.
        """
        J, R, iact = cache
        x, lam, xu = self._recover(a, J, R, iact, len(iact))
        lagr = np.zeros(self._q)
        lagr[iact - 1] = lam
        return self._verified(a, x, xu, lagr, lam, iact)

    def _repair(self, a: np.ndarray, cache: "_Cache") -> _WarmEntry:
        """Turn a stale active set into a dual-feasible state to resume from.

        A multiplier that has gone negative marks a constraint that no longer
        belongs in the active set. Dropping it and recomputing is exactly the
        step the solver's own inner loop takes, and repeating until none is
        negative restores the invariant the iteration requires. Whatever is left
        may still be primally infeasible -- constraints outside the set may be
        violated -- and driving that to zero is what the resumed loop is for.

        Terminates because each pass either stops or shrinks the active set; in
        the worst case everything is dropped and the resumed loop starts from the
        unconstrained minimum, which is the cold start.

        Args:
            a: ``(n,)`` linear term.
            cache: The stale factorisation.

        Returns:
            A :class:`~cvx.quadprog._solve._WarmEntry` satisfying that invariant.
        """
        # Copied because a repair mutates them, and the cache must survive intact
        # if the resumed solve then fails.
        J, R = cache.J.copy(), cache.R.copy()
        nact = len(cache.iact)
        iact = np.zeros(self._q, dtype=np.int64)
        iact[:nact] = cache.iact
        uv = np.zeros(min(self.n, self._q))

        while True:
            x, lam, xu = self._recover(a, J, R, iact, nact)
            uv[:nact] = lam
            if nact == 0:
                break
            # Equalities carry unrestricted multipliers, so only inequalities can
            # mark themselves as no longer belonging.
            candidates = np.where(iact[:nact] > self.meq, lam, np.inf)
            worst = int(np.argmin(candidates))
            if candidates[worst] >= 0.0:
                break
            nact = _drop_constraint(worst + 1, nact, uv, iact, J, R)

        obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)
        return _WarmEntry(J, R, iact, nact, x, uv, obj, xu)

    def _verified(
        self,
        a: np.ndarray,
        x: np.ndarray,
        xu: np.ndarray,
        lagr: np.ndarray,
        lam: np.ndarray,
        iact: np.ndarray,
    ) -> Solution | None:
        """Return a Solution if ``x`` satisfies the KKT conditions, else None.

        For a strictly convex QP the KKT conditions are sufficient, so this is a
        proof rather than a heuristic: dual feasibility on the inequalities, and
        primal feasibility of everything not held active.

        Args:
            a: ``(n,)`` linear term.
            x: Candidate minimiser.
            xu: Unconstrained minimiser.
            lagr: Full-length multiplier vector.
            lam: Multipliers of the active constraints only.
            iact: 1-based active set.

        Returns:
            The :class:`~cvx.quadprog.Solution`, or None if the cache is stale.
        """
        scale = _STALE_MARGIN * VSMALL * max(1.0, float(np.max(np.abs(x))))

        if lam.size and np.any(lam[iact > self.meq] < -scale):
            return None

        # Fresh array from every branch of the evaluator, which matters because the
        # active entries are forced to zero in place on the next line.
        sv = self._slack_of(x) - self.b
        if iact.size:
            sv[iact - 1] = 0.0
        if np.any(sv[self.meq :] < -scale) or np.any(np.abs(sv[: self.meq]) > scale):
            return None

        obj = 0.5 * float(x @ (self.G @ x)) - float(a @ x)
        return Solution(x, obj, xu, np.zeros(2, dtype=np.int64), lagr, iact)