API reference¶

def kkt_violation(a: SymmetricOperator, b: Vector, x: Vector) -> float:
    """Maximum violation of the KKT system of ``min_{x>=0} 1/2 x'Ax - b'x``.

    Args:
        a: The SPD operator ``A`` (a :class:`cvx.linalg.SymmetricOperator`).
        b: The linear term ``b``.
        x: Candidate solution.

    Returns:
        ``max`` of the negativity violations of ``x`` and of the reduced
        gradient ``s = A x - b``, and of the complementarity products
        ``|x_i s_i|``. Zero certifies the unique global minimiser.
    """
    op = _require_operator(a)
    _check_dimension(op, b)
    s = op.matvec(x) - b
    return float(
        max(
            np.max(-x, initial=0.0),
            np.max(-s, initial=0.0),
            np.max(np.abs(x * s), initial=0.0),
        )
    )

def solve_nnqp(
    a: SymmetricOperator,
    b: Vector,
    tol: float = 1e-8,
    cg_tol: float = 1e-10,
    p_max: int = 3,
    inner: str = "cg",
    track: bool = False,
    cg_maxit: int = 100_000,
    max_outer: int | None = None,
    warm: tuple[NDArray[np.bool_], Vector] | None = None,
) -> Result:
    """Minimise ``1/2 x^T A x - b^T x`` over ``x >= 0`` by the active-set loop.

    Each free-block solve is matrix-free CG on ``v -> op.apply_free(F, v)``;
    the reduced matrix is never materialised and ``A`` is never refactorised.
    The batch block-pivot fast path is guarded by a least-index Bland
    fallback, so termination at the unique global minimiser is unconditional
    — no non-degeneracy assumption.

    Args:
        a: The SPD operator ``A`` (a :class:`cvx.linalg.SymmetricOperator`) —
            ``DenseOperator`` for an explicit array, ``GramOperator(M, ridge)``
            for ``A = M^T M + ridge I`` whose Gram matrix is never formed.
        b: The linear term ``b``.
        tol: Threshold of the primal and dual violator tests.
        cg_tol: Relative residual tolerance of the inner solves. Keep it a
            couple of orders below ``tol`` so the inexact loop makes the same
            sign decisions as the exact one (Lemma 5.1 of the paper).
        p_max: Patience budget — non-improving batch steps tolerated before a
            fallback pivot. Any value gives finite termination.
        inner: ``"cg"`` (matrix-free), ``"pcg"`` (Jacobi-preconditioned from
            ``op.diag``), or ``"exact"`` (direct solve of each free block via
            ``op.solve_free``). Match the inner solver to the backend: pick
            ``"exact"`` when ``solve_free`` is structured and cheap — e.g.
            ``FactorOperator``'s Woodbury solve at ``O(|F| r^2)`` — and CG
            when only products are cheap (large dense ``A``, Gram factors
            with many rows).
        track: Record the free-set trajectory in ``Result.traj``.
        cg_maxit: Iteration cap per inner solve.
        max_outer: Optional cap on outer steps; when hit, the current iterate
            is returned with ``converged=False``.
        warm: Optional ``(free_mask, x_prev)`` pair from a previous solve.
            Starts the loop from that free set and warm-starts every CG call
            from the newest iterate — across a support-stable parameter step
            the loop then terminates in a single outer step.

    Returns:
        A :class:`Result`; ``converged`` is True iff the KKT system was
        satisfied to ``tol``, which certifies the unique global minimiser.

    Raises:
        TypeError: When ``a`` is not a :class:`cvx.linalg.SymmetricOperator`.
        ValueError: When the operator dimension does not match ``len(b)``, or
            when ``inner="exact"`` meets a numerically singular free block
            (``op.rcond_free`` below 1e-12) — ``A`` is then not positive
            definite on that free set; add a ridge.
        NotImplementedError: When ``inner="pcg"`` meets a backend that does
            not expose ``diag`` (propagated from ``cvx.linalg``).
    """
    op = _require_operator(a)
    _check_dimension(op, b)
    dinv: Vector | None = None  # Jacobi preconditioner, read off op.diag on first use

    def sub_solve(idx: NDArray[np.int_], x0: Vector | None) -> tuple[Vector, Vector | None, int]:
        """Solve the reduced system ``A_F x_F = b_F`` with the chosen inner solver."""
        nonlocal dinv
        if inner == "exact":
            rcond = op.rcond_free(idx)
            if rcond < _RCOND_MIN:
                msg = f"free block of size {idx.size} is numerically singular (rcond={rcond:.2e})"
                raise ValueError(msg)
            return op.solve_free(idx, b[idx]), None, 1
        if inner == "pcg":
            if dinv is None:
                dinv = 1.0 / op.diag
            xf, k_step = pcg(_free_matvec(op, idx), b[idx], dinv[idx], tol=cg_tol, maxit=cg_maxit)
            return xf, None, k_step
        xf, k_step = cg(_free_matvec(op, idx), b[idx], tol=cg_tol, maxit=cg_maxit, x0=x0)
        return xf, None, k_step

    def reduced_gradient(x: Vector, lam: Vector | None) -> Vector:  # noqa: ARG001
        """Return the reduced gradient ``s = A x - b``."""
        return op.matvec(x) - b

    return _active_set_loop(
        len(b),
        sub_solve,
        reduced_gradient,
        tol=tol,
        p_max=p_max,
        track=track,
        max_outer=max_outer,
        warm=warm,
    )

def solve_nnqp_eq(
    a: SymmetricOperator,
    b: Vector,
    b_eq: Matrix,
    c_eq: Vector,
    tol: float = 1e-8,
    cg_tol: float = 1e-10,
    p_max: int = 3,
    warm: tuple[NDArray[np.bool_], Vector] | None = None,
) -> Result:
    """Solve ``min 1/2 x^T A x - b^T x`` subject to ``x >= 0`` and ``B x = c``.

    On each free set the saddle system is solved by eliminating the multiplier
    ``lambda`` in R^p through the p-by-p Schur complement
    ``S = B_F A_F^{-1} B_F^T``: the ``p + 1`` right-hand sides share the
    operator ``A_F`` and are each one matrix-free CG solve, then
    ``S lambda = c - B_F v0`` fixes the multipliers in closed form. The single
    normalisation ``1^T x = beta`` is the ``p = 1`` case. ``B`` must have full
    row rank on the visited free sets (automatic for ``p = 1``).

    Args:
        a: The SPD operator ``A`` (a :class:`cvx.linalg.SymmetricOperator`).
        b: The linear term ``b``.
        b_eq: Equality matrix ``B`` of shape ``(p, n)``, full row rank.
        c_eq: Equality right-hand side ``c`` of shape ``(p,)``.
        tol: Threshold of the primal and dual violator tests.
        cg_tol: Relative residual tolerance of the inner CG solves.
        p_max: Patience budget of the batch fast path.
        warm: Optional ``(free_mask, x_prev)`` pair from a previous solve —
            the same tuple :func:`solve_nnqp` accepts. Starts the loop from
            that free set and seeds the ``v0`` solve of every saddle step
            from the newest iterate; the ``v1`` columns are re-solved cold.
            Across a support-stable parameter step the loop then terminates
            in a single outer step.

    Returns:
        A :class:`Result` with the multipliers in ``lam``. The reduced
        gradient underlying the dual test is ``s = A x - b - B^T lam``.

    Raises:
        TypeError: When ``a`` is not a :class:`cvx.linalg.SymmetricOperator`.
        ValueError: When the operator dimension does not match ``len(b)``.
    """
    op = _require_operator(a)
    _check_dimension(op, b)
    p = b_eq.shape[0]

    def sub_solve(idx: NDArray[np.int_], x0: Vector | None) -> tuple[Vector, Vector | None, int]:
        """Solve the saddle system on the free set via the p-by-p Schur complement."""
        matvec_f = _free_matvec(op, idx)
        b_f = b_eq[:, idx]
        v0, k0 = cg(matvec_f, b[idx], tol=cg_tol, x0=x0)
        v1 = np.zeros((idx.size, p))
        k_cols = 0
        for j in range(p):
            v1[:, j], kj = cg(matvec_f, b_f[j], tol=cg_tol)
            k_cols += kj
        schur = b_f @ v1  # p-by-p Schur complement, SPD
        lam = cholesky_solve(schur, c_eq - b_f @ v0)
        xf = v0 + v1 @ lam  # x_F = A_F^{-1}(b_F + B_F^T lambda)
        return xf, lam, k0 + k_cols

    def reduced_gradient(x: Vector, lam: Vector | None) -> Vector:
        """Return the constrained reduced gradient ``s = A x - b - B^T lam``."""
        correction = b_eq.T @ lam if lam is not None else np.zeros_like(b)
        return op.matvec(x) - b - correction

    return _active_set_loop(len(b), sub_solve, reduced_gradient, tol=tol, p_max=p_max, warm=warm)

def cg(
    matvec: MatVec,
    rhs: Vector,
    tol: float = 1e-8,
    maxit: int = 100_000,
    x0: Vector | None = None,
) -> tuple[Vector, int]:
    """Solve an SPD system by conjugate gradients.

    Args:
        matvec: The action ``v -> A v`` of an SPD operator.
        rhs: Right-hand side ``b``.
        tol: Relative residual stopping tolerance ``||b - A x|| / ||b||``.
        maxit: Iteration cap; the current iterate is returned when it is hit.
        x0: Optional warm start. The initial residual is ``b - A x0``, so a
            good guess cuts the iteration count by the log of the initial
            error.

    Returns:
        The approximate solution and the number of iterations taken.
    """
    if x0 is None:
        x = np.zeros_like(rhs)
        r = rhs.copy()
    else:
        x = x0.astype(np.float64, copy=True)
        r = rhs - matvec(x)
    p = r.copy()
    rs = float(r @ r)
    bnorm = float(np.linalg.norm(rhs))
    if bnorm == 0.0:
        return np.zeros_like(rhs), 0
    for it in range(1, maxit + 1):
        ap = matvec(p)
        alpha = rs / float(p @ ap)
        x += alpha * p
        r -= alpha * ap
        rs_new = float(r @ r)
        if np.sqrt(rs_new) / bnorm <= tol:
            return x, it
        p = r + (rs_new / rs) * p
        rs = rs_new
    return x, maxit

def pcg(
    matvec: MatVec,
    rhs: Vector,
    dinv: Vector,
    tol: float = 1e-8,
    maxit: int = 100_000,
) -> tuple[Vector, int]:
    """Solve an SPD system by Jacobi-preconditioned conjugate gradients.

    The preconditioner is ``M^{-1} = diag(dinv)``; for operators that are a
    well-conditioned core under a bad diagonal scaling, PCG runs at the core's
    condition number regardless of the scaling.

    Args:
        matvec: The action ``v -> A v`` of an SPD operator.
        rhs: Right-hand side ``b``.
        dinv: Elementwise inverse of the operator's diagonal.
        tol: Relative residual stopping tolerance.
        maxit: Iteration cap; the current iterate is returned when it is hit.

    Returns:
        The approximate solution and the number of iterations taken.
    """
    x = np.zeros_like(rhs)
    r = rhs.copy()
    z = dinv * r
    p = z.copy()
    rz = float(r @ z)
    bnorm = float(np.linalg.norm(rhs))
    if bnorm == 0.0:
        return x, 0
    for it in range(1, maxit + 1):
        ap = matvec(p)
        alpha = rz / float(p @ ap)
        x += alpha * p
        r -= alpha * ap
        if float(np.linalg.norm(r)) / bnorm <= tol:
            return x, it
        z = dinv * r
        rz_new = float(r @ z)
        p = z + (rz_new / rz) * p
        rz = rz_new
    return x, maxit