KKT Solvers

At each interior-point iteration, ConicIP solves a 3×3 block KKT system. The choice of solver for this system has a significant impact on performance. This guide covers the built-in solvers, when to use each, and how to write custom solvers.

The KKT System

The system solved at each iteration is:

┌             ┐ ┌   ┐   ┌   ┐
│ Q   G'  -A' │ │ a │   │ x │
│ G           │ │ b │ = │ y │
│ A       FᵀF │ │ c │   │ z │
└             ┘ └   ┘   └   ┘

where F is a Block diagonal matrix representing the Nesterov-Todd scaling. The block type depends on the cone:

Automatic Selection (default)

By default (kktsolver = default_kktsolver), ConicIP picks the solver per problem via choose_kktsolver, using cone mix, size, and structural sparsity — never the storage type of the inputs:

1. Any SDP cone → kktsolver_qr. The dense double-QR method is the numerically robust choice for the dense SDP scaling blocks.
2. Small problems (n + m + p < 1000) → kktsolver_qr. Dense factorization wins at small sizes; behavior matches the historical default exactly.
3. Dense-ish data (more than 10 structural nonzeros per column of [Q; A; G] on average) → kktsolver_qr. A sparse-typed matrix with dense columns (e.g. many small SOCs over a 10%-dense A) factors faster densely.
4. Otherwise → kktsolver_sparse. Large and genuinely sparse — the regime of issue #10, where the dense method is slower by orders of magnitude.

With verbose = true the solver prints the choice it made. Passing kktsolver explicitly always overrides the heuristic.

Built-in Solvers

kktsolver_qr

QR-based solver using the double QR method from CVXOPT: a thin QR of G' at setup, then per iteration a Cholesky factorization of the reduced Hessian Q₂'(Q + AᵀF⁻¹F⁻ᵀA)Q₂ on the null space of G.

When to use: SDP problems, and small or dense problems.

Trade-offs: Numerically robust, but the per-iteration factorization is dense — O((n-p)³) — so it is the wrong choice for large sparse problems.

kktsolver_sparse

Sparse LU solver that intelligently chooses between two internal strategies:

• Lifted formulation: Replaces large diagonal-plus-low-rank blocks with lifted variables, keeping the system sparse. Better for large second-order cones.
• Dense formulation: Converts all scaling blocks to dense matrices. Better when constraints are the product of many small cones.

The solver estimates the nonzero count for each strategy and picks the sparser one automatically. (SDP blocks cannot be lifted, so any "S" cone forces the dense formulation.) The sparse LU factorization reuses its symbolic analysis across iterations once the sparsity pattern stabilizes.

When to use: Large problems with sparse Q and A.

sol = conicIP(Q, c, A, b, cone_dims; kktsolver=kktsolver_sparse)

kktsolver_2x2 (with pivot)

A 2×2 sparse LU solver that works on the Schur complement system obtained by pivoting on the third block. Must be wrapped with pivot:

sol = conicIP(Q, c, A, b, cone_dims; kktsolver=pivot(kktsolver_2x2))

When to use: Problems where the Schur complement Q + Aᵀ(FᵀF)⁻¹A is sparser or better conditioned than the full 3×3 system.

Choosing a Solver Manually

The automatic default covers the common cases. Override it when you know better:

If a factorization breaks down (e.g. rank-deficient G passed directly to conicIP), the solver returns a Solution with status == :Error and the failure reason in sol.message rather than throwing; preprocess_conicIP removes redundant rows up front, and structurally degenerate inputs (an all-zero equality row, a variable appearing in no constraint) are detected exactly in O(nnz) and either deflated or answered with a certified :Infeasible/:Unbounded.

Writing a Custom Solver

A custom KKT solver is a three-level nested function:

function my_kktsolver(Q, A, G, cone_dims)
    # Level 1: One-time setup (symbolic factorization, preallocation)
    # Called once before the solve loop.

    function solve3x3gen(F, F⁻ᵀ)
        # Level 2: Per-iteration setup (F changes each iteration)
        # Compute numeric factorization using current scaling F.

        function solve3x3(x, y, z)
            # Level 3: Solve the 3×3 system given RHS (x, y, z).
            # Return (a, b, c) where:
            #   Qa + G'b - A'c = x
            #   Ga = y
            #   Aa + FᵀFc = z
            return (a, b, c)
        end
        return solve3x3
    end
    return solve3x3gen
end

Pass it to the solver via the kktsolver keyword:

sol = conicIP(Q, c, A, b, cone_dims; kktsolver=my_kktsolver)

Example: Diagonal QP

For minimize ½ xᵀQx - cᵀx subject to x ≥ 0 with diagonal Q, the KKT system simplifies (no equality constraints, G is empty) to:

┌            ┐ ┌   ┐   ┌   ┐
│ Q       -I │ │ a │   │ x │
│ I    FᵀF   │ │ c │ = │ z │
└            ┘ └   ┘   └   ┘

Since F is Diagonal for "R" cones, pivoting on the second block gives (Q + (FᵀF)⁻¹) a = x + (FᵀF)⁻¹ z, solvable by Cholesky:

using ConicIP, SparseArrays, LinearAlgebra, Random
Random.seed!(42)

n = 50
Q = sprandn(n, n, 0.3); Q = Q'Q + 0.1I  # make positive definite
c = ones(n)
A = sparse(1.0I, n, n)
b = zeros(n)
cone_dims = [("R", n)]

function my_kktsolver(Q, A, G, cone_dims)
    function solve3x3gen(F, F⁻ᵀ)
        invFᵀF = inv(F'F)
        QpD = cholesky(Q + spdiagm(0 => (F[1].diag).^(-2)))

        function solve3x3(x, y, z)
            a = QpD \ (x + A' * (invFᵀF * z))
            c = invFᵀF * (z - A * a)
            b = zeros(0)
            return (a, b, c)
        end
    end
end

sol = conicIP(Q, c, A, b, cone_dims; kktsolver=my_kktsolver, verbose=false)
sol.status

:Optimal

The pivot Wrapper

The pattern of reducing a 3×3 system to 2×2 by pivoting on the third block is common enough that ConicIP provides pivot to automate it.

A 2×2 solver has the signature:

function my_2x2_solver(Q, A, G, cone_dims)
    function solve2x2gen(F, F⁻ᵀ)
        # Build and factor the Schur complement: Q + Aᵀ(FᵀF)⁻¹A
        function solve2x2(y, w)
            # Solve for (Δy, Δw) and return them
            return (Δy, Δw)
        end
        return solve2x2
    end
    return solve2x2gen
end

Then pivot(my_2x2_solver) produces a valid 3×3 solver. Here's the same diagonal QP using pivot:

function my_2x2_solver(Q, A, G, cone_dims)
    function solve2x2gen(F, F⁻ᵀ)
        QpD = cholesky(Q + spdiagm(0 => (F[1].diag).^(-2)))
        return (y, w) -> (QpD \ y, zeros(0))
    end
end

sol2 = conicIP(Q, c, A, b, cone_dims;
               kktsolver=pivot(my_2x2_solver), verbose=false)
sol2.status

:Optimal

Performance Tips

1. Preallocate buffers in level 1 (the outer function) and reuse them in levels 2 and 3.
2. Reuse symbolic factorizations when the sparsity pattern doesn't change between iterations (only the numeric values of F change).
3. Avoid inv(F'F) for large blocks — compute the action of (FᵀF)⁻¹ on a vector instead.
4. For problems with a callback.ipynb example, see the examples/ directory in the repository.