API Reference¶

def solve_qp(
    G: np.ndarray,
    a: np.ndarray,
    C: np.ndarray | None = None,
    b: np.ndarray | None = None,
    meq: int = 0,
    factorized: bool = False,
) -> Solution:
    r"""Solve a strictly convex quadratic program.

    .. math::
        \min_x \tfrac{1}{2} x^T G x - a^T x \quad\text{subject to}\quad C^T x \ge b

    The first ``meq`` constraints are treated as equalities.

    Args:
        G: ``(n, n)`` symmetric positive definite matrix of the quadratic term.
            If ``factorized`` is True, pass :math:`R^{-1}` instead, where
            :math:`G = R^T R` with ``R`` upper triangular.
        a: ``(n,)`` vector of the linear term.
        C: ``(n, m)`` constraint matrix, one column per constraint. Defaults to
            a single inactive constraint, giving the unconstrained problem.
        b: ``(m,)`` right-hand side of the constraints.
        meq: Number of leading constraints to treat as equalities.
        factorized: Whether ``G`` holds :math:`R^{-1}` rather than :math:`G`.

    Returns:
        A :class:`Solution` with the minimiser, the objective value, the
        unconstrained minimiser, the iteration counts, the Lagrange multipliers
        and the active set.

    Raises:
        ValueError: If the shapes are inconsistent, if ``meq`` is out of range,
            if ``G`` is not positive definite, or if the constraints admit no
            solution.
    """
    G = np.asarray(G, dtype=np.float64)
    a = np.asarray(a, dtype=np.float64)

    if C is None and b is None:
        C, b, meq = np.zeros((len(G), 1)), -np.ones(1), 0
    elif C is None or b is None:
        raise ValueError("C and b must be given together")

    C = np.asarray(C, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64)

    n, q = _validate(G, a, C, b, meq)
    r = min(n, q)

    # Initialisation. We want xv to hold G^-1 a, the unconstrained minimum, and
    # J to hold R^-1, so that J J^T = G^-1.
    J, xv = _factorize(G, a, factorized)

    # The objective at the unconstrained minimum. Kept as a running total: each
    # step updates it in closed form rather than re-evaluating the quadratic.
    obj = -float(a @ xv) / 2.0
    xu = xv.copy()

    # The norm of each column of C, used to scale the pivoting rule so that the
    # choice of constraint is invariant to how each one happens to be scaled.
    # A zero-norm column reads 0 >= b, which no x can influence; scoring it as
    # infinitely violated sends the solver to the infeasibility verdict instead
    # of dividing by zero below.
    nbv = np.sqrt(np.sum(C * C, axis=0))
    degenerate = nbv == 0.0
    nbv_safe = np.where(degenerate, 1.0, nbv)

    # Sparsity of C, detected once. Bound constraints make most columns a single
    # scaled unit vector, which turns three of the per-iteration operations from
    # O(n) or O(n*q) work into scalar indexing.
    single, srow, sval = _analyse_constraints(C)
    slack_of = _slack_evaluator(C, single)

    # Upper triangular, stored as packed columns -- see the note in _qr.
    R = np.zeros(r * (r + 1) // 2)
    uv = np.zeros(r)  # dual variables of the active constraints
    iact = np.zeros(q, dtype=np.int64)  # 1-based, first nact entries valid
    lagr = np.zeros(q)
    nact = 0
    iter_full, iter_partial = 0, 0

    while True:
        iter_full += 1

        # The slack of every constraint. Slacks of active constraints are forced
        # to exactly zero as a safeguard against rounding error.
        sv = slack_of(xv) - b
        sv[np.abs(sv) < VSMALL] = 0.0
        sv[iact[:nact] - 1] = 0.0

        iadd = _choose_constraint(sv, nbv_safe, degenerate, meq)

        if iadd == 0:
            # Every constraint is satisfied, so we are at the optimum.
            lagr[iact[:nact] - 1] = uv[:nact]
            iterations = np.array([iter_full, iter_partial], dtype=np.int64)
            return Solution(xv, obj, xu, iterations, lagr, iact[:nact])

        # An equality constraint may be violated from either side. When its
        # slack is positive we have to step in the opposite direction.
        slack = float(sv[iadd - 1])
        reverse_step = slack > 0.0
        u = 0.0

        # A column holding a single scaled unit vector e_row lets the three
        # products against it below be read off by index instead of computed.
        unit = bool(single[iadd - 1])
        if unit:
            row, val = int(srow[iadd - 1]), float(sval[iadd - 1])
        else:
            normal = C[:, iadd - 1]

        # Inner loop: walk towards the constraint boundary, dropping active
        # constraints whose multipliers would otherwise turn negative.
        while True:
            # dv = J^T n, split as (d_1, d_2) at the size of the active set.
            # For a unit column this is one scaled row of J, O(n) not O(n^2).
            dv = val * J[row, :] if unit else J.T @ normal

            # zv = J_2 d_2 is the step direction of the primal variable, the
            # component of the constraint normal orthogonal to the active set.
            zv = J[:, nact:] @ dv[nact:]

            # rv = R^-1 d_1 is the negated step direction of the dual variable.
            # Solved on a copy: dv is still needed intact for qr_insert below.
            rv = _TPSV(nact, R, dv[:nact].copy(), overwrite_x=True) if nact else _EMPTY

            # The largest step t1 that keeps the dual variables non-negative,
            # and the constraint idel that would be the first to bind at zero.
            t1, idel = _dual_step_limit(uv, rv, iact, nact, meq, reverse_step)
            t1inf = idel == 0

            # The step t2 that brings the slack of the entering constraint to
            # zero. ztn is the rate of change of that slack.
            t2inf = abs(float(zv @ zv)) <= VSMALL
            if not t2inf:
                ztn = val * float(zv[row]) if unit else float(zv @ normal)
                t2 = abs(slack) / ztn

            if t1inf and t2inf:
                # We can step infinitely far: the dual is unbounded, so the
                # primal is infeasible.
                raise ValueError("constraints are inconsistent, no solution")

            full_step = not t2inf and (t1inf or t1 >= t2)
            step_length = t2 if full_step else t1
            step = -step_length if reverse_step else step_length

            if not t2inf:
                xv += step * zv
                obj += step * ztn * (step / 2.0 + u)

            uv[:nact] -= step * rv
            u += step

            if full_step:
                break

            # Only a partial step: drop constraint idel from the active set.
            qr_delete(nact, idel, J, R)
            uv[idel - 1 : nact - 1] = uv[idel:nact].copy()
            iact[idel - 1 : nact - 1] = iact[idel:nact].copy()
            uv[nact - 1], iact[nact - 1] = 0.0, 0
            nact -= 1
            iter_partial += 1

            if not t2inf:
                # We moved in primal space, so the slack we are closing has
                # changed and must be recomputed.
                reached = val * float(xv[row]) if unit else float(xv @ normal)
                slack = reached - float(b[iadd - 1])

        # The entering constraint now holds with equality: add it.
        nact += 1
        uv[nact - 1], iact[nact - 1] = u, iadd
        qr_insert(nact, dv, J, R)