API Reference

conicIP(Q, c, A, b, conedims, G, d; kktsolver = defaultkktsolver, optTol = 1e-6, DTB = 0.01, verbose = true, maxRefinementSteps = 3, maxIters = 100, cache_nestodd = false, infeasTol = 1e-7, infeasAbsTol = 1e-9, staticReg = 0.0, certFallback = true, certFallbackIters = 50, refineRelTol = 1e-13, refineAbsTol = 1e-12, timeLimit = Inf, centralityCorrectors = 0, timing = nothing)

Interior point solver for the system

minimize    ½yᵀQy - cᵀy
s.t         Ay >= b
            Gy  = d

c, b, d are vectors (or any AbstractVector)

cone_dims is an array of tuples (Cone Type, Dimension)

e.g. [("R",2),("Q",4)] means
(y₁, y₂)          in  R+
(y₃, y₄, y₅, y₆)  in  Q

A semidefinite block is ("S", k) with k = n(n+1)/2 for an n × n symmetric matrix. Those k rows carry the matrix in the vecm form: the upper triangle read row by row, with the off-diagonal entries scaled by √2, so that dot(vecm(X), vecm(Y)) == tr(X*Y). See vecm and mat.

Returns a Solution whose status is one of

• :Optimal — max(rDu, rPr, rCp, rEq, rGap) < optTol, and every cone and equality row also passes the row-wise test |rᵢ| / (1 + |bᵢ| + (|A||y|)ᵢ + |sᵢ|) < optTol (|rᵢ| / (1 + |dᵢ| + (|G||y|)ᵢ) for equalities) in the original coordinates; the aggregate 2-norm test alone can accept a point that violates a small-scale row when the data span many orders of magnitude.
• :Infeasible / :DualInfeasible — a ray passed a screen and was accepted by the corresponding validator; has_certificate is then true.
• :AlmostInfeasible / :AlmostDualInfeasible — set only at loop exhaustion, when the best iterate carries a ray that validates at 100*infeasTol but not at infeasTol. The best iterate is retained.
• :Abandoned — iteration limit reached with no verdict.
• :TimeLimit — timeLimit seconds elapsed; the best iterate is retained.
• :Error — nonfinite residuals, a nonfinite search direction or iterate, or a KKT factorization failure (the reason is recorded in sol.message). Rank-deficient G handed directly to conicIP typically lands here; use preprocess_conicIP to trim redundant rows first. staticReg regularizes only the Q block of the KKT system and cannot repair a rank-deficient G.

Every factorization or cone line search that can fail on a boundary iterate is guarded, so such a failure is reported as an :Error status with a reason in sol.message and never as an escaped exception (issue #10). A guarded failure returns immediately with the current iterate (the initial point, with Iter = 0 and pobj = Inf, if no iteration completed); it does not go through the certificate-fallback path, which runs only when the iteration limit is exhausted.

Structurally degenerate inputs are handled exactly before any factorization: an all-zero row of G is deflated (dᵢ = 0) or answered with a certified :Infeasible (dᵢ ≠ 0), and a variable absent from Q, A, and G is deflated (cⱼ = 0) or answered with a certified :DualInfeasible (cⱼ ≠ 0).

Selected keyword arguments:

• infeasTol — infeasibility-certificate tolerance, decoupled from optTol.
• infeasAbsTol — absolute tolerance for certificate validation.
• staticReg — static KKT regularization scale; 0 (default) disables it. preprocess_conicIP enables it when it detects rank deficiency.
• certFallback — enable fallback certificate solve on stall.
• maxRefinementSteps, refineRelTol, refineAbsTol — the predictor and corrector steps are refined against the 4×4 KKT system until ‖r − KΔz‖ ≤ refineAbsTol + refineRelTol·‖r‖ or the step budget is spent. The residual is always evaluated: the KKT solver's own residual bound (LDLDiagnostics.last_bound) is reported for diagnosis but does not stand in for it.
• timeLimit — wall-clock budget in seconds, checked once per iteration (a single factorization can overrun it). On expiry the status is :TimeLimit and the solution holds the best iterate so far; the certificate-fallback solves are skipped.
• objective_offset — a constant added to the objective for the purpose of the relative gap test only (⟨v,s⟩/(1 + |pobj + objective_offset|)). preprocess_conicIP passes the constant carried by fixed variables, so a reduced problem terminates by the same criterion as the full one.
• centralityCorrectors — number of Gondzio multiple centrality correctors tried per iteration (default 0, off). Each corrector re-solves the current KKT factorization once for a direction that pushes the trial complementarity toward the box [0.1·σμ, 10·σμ] (in the Jordan frame of each cone) and is kept only if it lengthens the step; the loop stops at the first rejected corrector, and nothing is tried when the step is already full. Extra solves are counted in kkt_solves; the verbose cc column shows accepted/tried.
• timing — a PhaseTimes object to accumulate per-phase wall times, allocation bytes, and counts into (see src/timing.jl for the contract); nothing (default) leaves the solver uninstrumented.

The parameter solve3x3gen allows the passing of a custom solver for the KKT System, as follows

julia> L = solve3x3gen(F,F⁻ᵀ,Q,A,G)

Then this

julia> (a,b,c) = L(y,w,v)

solves the system
┌             ┐ ┌   ┐   ┌   ┐
│ Q   G'  -A' │ │ a │ = │ y │
│ G           │ │ b │   │ w │
│ A       FᵀF │ │ c │   │ v │
└             ┘ └   ┘   └   ┘

a, b and c may be fresh vectors or views into the solver's own workspace, valid only until the next call to L — conicIP copies them into its direction buffers before calling again. kktsolver_ldl returns views; the other built-in solvers return fresh vectors.

We can also wrap a 2x2 solver using pivot3gen(solve2x2gen) The 2x2 solves the system

julia> L = solve2x2gen(F,F⁻ᵀ,Q,A,G)

Then this

julia> (a,b) = L(y,w)

solves the system

┌                     ┐ ┌   ┐   ┌   ┐
│ Q + Aᵀinv(FᵀF)A  G' │ │ a │ = │ y │
│ G                   │ │ b │   │ w │
└                     ┘ └   ┘   └   ┘

equilibrate = true (default) applies Ruiz equilibration to the data before solving and maps the solution, rays, and residuals back to the original coordinates; see equilibrate_conicIP. A custom kktsolver then receives the scaled data.

ConicIP with preprocessing to ensure the following rank constraints

Primal equailty constraints : Gx = d Rank condition : rank(G) = size(G,1)

Dual equality constraints : [ Q A' G'] = c Rank condition : rank([Q A' G']) = size(Q,1)

Inconsistent data is reported with a certificate whenever one can be constructed and verified against the original problem data:

• Gy = d inconsistent → :Infeasible with a Farkas ray (w,v)
• c ∉ range([Q Aᵀ Gᵀ]) → :DualInfeasible with a recession ray y

Rank deficiency that is not an inconsistency is handled by opting into conicIP's static KKT regularization rather than by perturbing Q.

rank_check controls the sparse-QR rank detection, which is the expensive part: :always runs it; :never skips it and relies on the KKT solver's regularization (what kktsolver_ldl provides); :auto (default) runs it only when the KKT solver that will be used needs a full-rank G — the dense QR solver, chosen for small and semidefinite problems, and any explicit solver other than kktsolver_ldl.

Before either, singleton equality rows (gᵢⱼ yⱼ = dᵢ) fix their variable and are removed together with the column; the fixed values, their equality duals, and the objective values are restored on the way out (fix_singletons = false disables this step).

timing = pt with a PhaseTimes charges the singleton and rank work to t_presolve, the restoration to t_postsolve, and forwards pt to the nested solve, which fills the remaining phases.

Solution

Return type of conicIP and preprocess_conicIP.

Fields

• y::Vector{Float64} – primal variables
• w::Vector{Float64} – dual variables for equality constraints (Gy = d)
• v::Vector{Float64} – dual variables for inequality constraints (Ay ≥_K b)
• s::Vector{Float64} – cone slack variables (Ay - s = b, s ∈ K)
• status::Symbol – :Optimal, :Infeasible, :DualInfeasible, :AlmostInfeasible, :AlmostDualInfeasible, :Abandoned, :TimeLimit, or :Error
• Iter::Integer – number of interior-point iterations
• Mu::Real – final complementarity gap parameter
• prFeas::Real – primal feasibility residual
• duFeas::Real – dual feasibility residual
• muFeas::Real – complementarity residual
• pobj::Real – primal objective value
• dobj::Real – dual objective value
• has_certificate::Bool – the returned vectors carry a verified ray certifying infeasibility or unboundedness (see the table below)

Field conventions by status

One exception: when a ray found on equilibrated or presolved data fails revalidation against the original data (sol.message says so), the :Almost* or :Abandoned solution holds that ray in the ray fields, not the best iterate.

The infeasibility ray is normalized so that dᵀw̄ - bᵀv̄ = -1 with Gᵀw̄ - Aᵀv̄ ≈ 0; the unboundedness ray is normalized so that cᵀȳ = +1 with Qȳ ≈ 0, Gȳ ≈ 0 and Aȳ ∈ K. See validate_infeasibility_certificate and validate_unboundedness_certificate.

kkt_solves counts the KKT back-solves the main loop performed (initial point, predictor, corrector, and refinement corrections); solves made by the certificate-fallback auxiliary problems are not included.

Residual and diagnostic tail:

• rEq::Real – relative equality residual ‖Gy − d‖ / (1 + max(‖d‖, ‖|G||y|‖)) of the returned point (prFeas = max(rPr, rEq)); NaN when no iterate was evaluated or when a certificate is returned (the iterate is discarded)
• rGap::Real – relative duality gap |vᵀs| / (1 + |pobj + objective_offset|), the quantity the termination test compares with optTol; NaN when unset or when a certificate is returned
• kkt_repaired::Int – pivots the KKT solver dynamically regularized, summed over every factorization of the solve (0 unless the solver reports diagnostics; see kkt_diagnostics)
• kkt_refactors::Int – refactorizations after a regularization bump (kktsolver_ldl with retry_max > 0), summed over the solve

Constructors with 12, 13, 14, or 15 positional arguments default the trailing fields to has_certificate = false, message = "", kkt_solves = 0, rEq = rGap = NaN, kkt_repaired = kkt_refactors = 0.

CertificateCheck

Verdict returned by the certificate validators.

Fields

• valid::Bool – all checks passed
• farkas_residual::Float64 – ‖·‖_∞ of the linear residual of the ray
• separation::Float64 – the (pre-normalization) separation value; must be > 0
• cone_margin::Float64 – blockwise minimum cone margin of the ray (≥ 0 if in K)
• finite::Bool – the candidate ray had only finite entries

cone_margin(x, cone_dims)

Blockwise minimum margin of x with respect to the cone K described by cone_dims. Nonnegative iff x ∈ K; the magnitude measures the distance to the boundary (violation depth when negative).

Per block:

"R"  minimum(x[I])
"Q"  x[I][1] - ‖x[I][2:end]‖
"S"  eigmin(Symmetric(mat(x[I])))

Returns +Inf when there are no blocks.

validate_infeasibility_certificate(Q, c, A, b, cone_dims, G, d, w, v;
                                   abstol, reltol)

Validate (w,v) as a Farkas ray proving primal infeasibility of {y : Ay ≥_K b, Gy = d}. A valid ray satisfies

Gᵀw̄ - Aᵀv̄ ≈ 0,    v̄ ∈ K,    dᵀw̄ - bᵀv̄ = -1

Returns (check::CertificateCheck, w̄, v̄) with the ray normalized so that dᵀw̄ - bᵀv̄ = -1. If the separation -(dᵀw - bᵀv) is nonpositive, or the candidate has nonfinite entries, the verdict is invalid and the candidates are returned unchanged (no normalization).

validate_unboundedness_certificate(Q, c, A, b, cone_dims, G, d, y;
                                   abstol, reltol)

Validate y as a recession ray proving ½yᵀQy - cᵀy is unbounded below over the feasible set. A valid ray satisfies

Qȳ ≈ 0,    Gȳ ≈ 0,    Aȳ ∈ K,    cᵀȳ = +1

so that the objective decreases without bound along ȳ. Returns (check::CertificateCheck, ȳ) with the ray normalized so that cᵀȳ = +1. If the separation cᵀy is nonpositive, or the candidate has nonfinite entries, the verdict is invalid and the candidate is returned unchanged.

fallback_infeasibility_ray(Q, c, A, b, cone_dims, G, d;
                           kktsolver = default_kktsolver, maxIters = 50)

Solve the minimum-norm Farkas auxiliary QP and return a candidate infeasibility ray (w, v), or nothing if the auxiliary solve does not reach :Optimal (which includes the common case that the original problem is feasible, making the auxiliary problem infeasible).

The returned pair is not validated and not normalized — pass it to validate_infeasibility_certificate.

fallback_unbounded_ray(Q, c, A, b, cone_dims, G, d;
                       kktsolver = default_kktsolver, maxIters = 50)

Solve the minimum-norm recession auxiliary QP and return a candidate unboundedness ray y, or nothing if the auxiliary solve does not reach :Optimal (in particular when no recession ray exists).

The Qy = 0 rows are omitted when Q is identically zero, which keeps the auxiliary system at its smallest for the LP case. The returned ray is not validated and not normalized — pass it to validate_unboundedness_certificate.

Optimizer(; kwargs...)

MathOptInterface optimizer wrapping the ConicIP interior-point solver. Use as a JuMP solver via Model(ConicIP.Optimizer).

Options

Settable as constructor keywords or through MOI.RawOptimizerAttribute / JuMP's set_attribute:

• verbose::Bool – print solver iterations (default: false)
• optTol::Float64 – optimality tolerance (default: 1e-6)
• maxIters::Int – maximum iterations (default: 100)
• infeasTol::Float64 – infeasibility/unboundedness certificate tolerance (default: 1e-7)
• kktsolver – "auto" (default; picks by cone mix and predicted factorization cost via choose_kktsolver), "ldl", "qr", "sparse", "2x2", or any callable solver object (a function or an instance such as cached_kktsolver_ldl())
• preprocess::Bool – remove redundant equality rows via preprocess_conicIP before solving (default: true)
• equilibrate::Bool – Ruiz-scale the data before solving (default: true)
• timeLimit::Float64 – wall-clock budget in seconds (also MOI.TimeLimitSec)
• assemble_only::Bool – stop optimize! once the solver's matrices (Q_int, c_int, ineq_A, ineq_b, cone_dims, eq_G, eq_d) are assembled, without calling the solver; the model then reports OPTIMIZE_NOT_CALLED and ResultCount == 0 (default: false)
• centralityCorrectors::Int – Gondzio centrality correctors tried per iteration, each one extra back-solve of the current KKT factorization (default: 0, off)
• timing – nothing (default) or a ConicIP.PhaseTimes that optimize! adds its per-phase wall times to (model assembly is t_frontend; the solver fills the rest); the caller owns and resets it
• plus infeasAbsTol, DTB, maxRefinementSteps, refineRelTol, refineAbsTol, staticReg, certFallback, certFallbackIters, cache_nestodd — forwarded to conicIP

MOI.Silent is supported and overrides verbose.

Supported Objectives

Affine and convex quadratic (ScalarQuadraticFunction), both handled natively: a quadratic objective becomes the solver's Q rather than a second-order-cone reformulation, so positive semidefinite but singular Hessians are fine. A Hessian that is not positive semidefinite for the given sense (a nonconvex QP) is rejected before the solve with TerminationStatus == INVALID_MODEL and no result.

Supported Constraints

• Vector: Zeros, Nonnegatives, Nonpositives, SecondOrderCone, PositiveSemidefiniteConeTriangle
• Scalar: EqualTo, GreaterThan, LessThan

equilibrate_conicIP(Q, c, A, b, cone_dims, G, d;
                    iters = 10, tol = 1e-2, bound = 1e6,
                    σ_range = (1e-3, 1e3))

Ruiz-equilibrate the problem data. Returns a named tuple with the scaled data Q, c, A, b, G, d and the scalings Dc (variables), Dr (cone rows, constant on each second-order and semidefinite block), De (equality rows) and σ (objective), such that

Q = Dc⁻¹ Q̃ Dc⁻¹ / σ,  c = Dc⁻¹ c̃ / σ,  A = Dr⁻¹ Ã Dc⁻¹,  b = Dr⁻¹ b̃,
G = De⁻¹ G̃ Dc⁻¹,  d = De⁻¹ d̃.

At most iters sweeps are made, stopping early when every row and column ∞-norm of [Q Aᵀ Gᵀ; A 0 0; G 0 0] is within tol of one. Individual scale factors are kept within [1/bound, bound]; a zero row or column is left unscaled. The objective is rescaled to unit ∞-norm (σ = 1/‖c̃‖∞) only when ‖c̃‖∞ falls outside σ_range: the iteration is not invariant to the objective scale, and on well-scaled problems a σ ≠ 1 costs iterations. The matrices keep their storage type (dense stays dense, sparse stays sparse).

default_kktsolver(Q, A, G, cone_dims)

The default kktsolver for conicIP: dispatches to the solver picked by choose_kktsolver, reusing the KKT pattern the choice analysed when the answer is kktsolver_ldl. Satisfies the standard kktsolver interface, so it can be passed anywhere a concrete solver can.

choose_kktsolver(Q, A, G, cone_dims;
                 size_min = 200, dense_bytes_max = 4 * 2^30,
                 ldl_flop_weight = 10.0)

Pick a KKT solver from the problem's cone mix, size, and predicted factorization cost (issue #10). Returns one of the solver constructors, chosen by:

1. the dense solver's storage estimate dense_kkt_bytes above dense_bytes_max ⇒ kktsolver_ldl, whatever the rules below would say — dense QR at that size is an out-of-memory error, not a slow solve. If the problem also has a semidefinite cone an ArgumentError is thrown instead: the sparse solver materializes each k×k SDP scaling block as a dense k²×k² block, so it would trade one huge allocation for another. Such problems currently need the dense path; raise the budget explicitly with kktsolver = (Q, A, G, cd) -> choose_kktsolver(Q, A, G, cd; dense_bytes_max = …)(Q, A, G, cd) or pass kktsolver = kktsolver_qr directly;
2. any SDP cone ⇒ kktsolver_qr — the dense double-QR method is the numerically robust choice for the dense SDP scaling blocks, and the sparse solver's SDP path is dense in k(k+1)/2;
3. n + m + p < size_min ⇒ kktsolver_qr — dense factorization wins at small sizes (the symbolic analysis below would cost more than it saves there);
4. otherwise the two per-iteration flop estimates decide: kktsolver_qr if dense_kkt_flops is below ldl_flop_weight times the LDLᵀ estimate Σⱼ nnz(L₍:,ⱼ₎)² taken from a symbolic analysis of the quasi-definite KKT pattern, else kktsolver_ldl. The weight accounts for the dense path running in BLAS and the LDLᵀ in scalar code; 10 reproduces the measured crossover on the benchmark set (many small SOCs over a 10%-dense A go dense, banded and block-structured problems go sparse).

kktsolver_sparse (UMFPACK LU) is no longer selected automatically; it remains available explicitly.

dense_kkt_bytes(n, m, p)

Routing estimate of the persistent dense arrays kktsolver_qr holds at once, used by choose_kktsolver to keep hopeless sizes off the dense path: the n×n orthogonal factor and the n×p dense copy of Gᵀ at setup, the m×(n−p) A*Q2 and the (n−p)² reduced Hessian, plus one more of each per iteration (W and Lmat). Float64 throughout.

This is not a peak-memory bound: LAPACK workspaces, the temporaries of F.Q * I and Q2' * (Q * Q2), the Cholesky (or fallback QR) copy of the reduced Hessian, and any fill in the caller's own matrices are all excluded. Peak usage can be a small multiple of this figure.

dense_kkt_flops(n, m, p; nnzA = 0, nnzQ = 0, iters = 25)

Per-iteration flop estimate for kktsolver_qr, comparable with the per-factorization estimate Σⱼ nnz(L₍:,ⱼ₎)² used for kktsolver_ldl. With r = n − p it counts

• per iteration: W = F⁻ᵀ(AQ2) (2mr), the reduced Hessian S22 + WᵀW (mr²) and its Cholesky (r³/3), plus three solves (predictor, corrector, one refinement) that each apply Q1, Q1ᵀ, Q2, Q2ᵀ (4np + 4nr), the Cholesky and R1 triangles (2r² + 2p²) and two Hmul products (2(nnzQ + 2nnzA));
• setup, amortized over iters iterations: the Householder QR of Gᵀ (2np² − 2p³/3), materializing the n×n Q0 = F.Q * I (4n²p), A*Q2 (2·nnzA·r) and Q2' * (Q * Q2) (2·nnzQ·r + 2nr²).

The setup terms matter when equalities nearly determine the variables: at p = n the per-iteration terms vanish but the n×n QR and the dense Q0 are still paid, so the old model (mr² + r³/3 alone) reported zero cost and routed n = m = p problems to the dense solver.

nnzA and nnzQ are structural nonzero counts; when omitted those terms are dropped, which only makes the estimate optimistic for the dense path.

kktsolver_ldl(Q, A, G, cone_dims;
              static_reg = 1e-8, cone_reg = 0.0,
              dynamic_eps = 1e-13, dynamic_delta = 2e-7,
              lift_min = 6, refine_steps = 2, refine_tol = 1e-13,
              retry_max = 0, retry_factor = 10.0,
              shift_floor = 1e-8, shift_max = 1e-4)

Sparse LDLᵀ KKT solver for the quasi-definite form of the 3×3 system. The upper triangle of

K_δ = [ Q + δ_p I   Gᵀ   −Aᵀ ;  G   −δ_e I   0 ;  −A   0   −(FᵀF + δ_c I) ]

(second-order cones of dimension at least lift_min lifted to diagonal + two columns + two pivots, smaller SOC and SDP blocks dense) is assembled once with every entry stored structurally, ordered by AMD, and analysed symbolically once; each iteration rewrites the scaling entries and refactorizes numerically in place (QDLDL.jl).

Regularization: static_reg is δ_p = δ_e, cone_reg is δ_c; the factorization also applies QDLDL's dynamic regularization, replacing any pivot whose sign disagrees with the quasi-definite pattern or whose magnitude is below dynamic_eps by ±dynamic_delta. Each solve is then refined against the unregularized matrix (δ = 0) for up to refine_steps corrections or until the residual is below refine_tol · (1 + ‖rhs‖), so the perturbation acts as a preconditioner rather than a change of problem. The residual is evaluated after every correction; a correction that does not reduce it is discarded and ends the refinement, so the returned solution is the best one seen and never worse than the unrefined solve. The tolerance is a target, not a guarantee. conicIP's own refinement against the 4×4 system runs on top of this.

Retry (off by default, retry_max = 0): when a solve is still above the refinement tolerance and either the factorization repaired a pivot or the first correction failed to reduce the residual, the static shifts are bumped, δ ← min(max(δ·retry_factor, shift_floor), shift_max) for δ_p and δ_e (a zero base shift starts at shift_floor), the matrix is refactorized, and the same right-hand side is solved again; the better of the two results by unregularized residual is returned. At most retry_max bumps per factorization. Bumped shifts last for the remaining solves with that factorization only: the next solve3x3gen call restores the base shifts. δ_c and dynamic_delta are never bumped.

Diagnostics: the object solve3x3gen returns is callable as before and carries an LDLDiagnostics record, reachable through ConicIP.kkt_diagnostics(solve3x3); conicIP prints its counts in the verbose kkt column and sums them into Solution.kkt_repaired and Solution.kkt_refactors. ConicIP.kkt_attach_timing!(solve3x3, pt) makes the backend add the time and count of every numeric refactorization, triangular solve, and refinement residual to the t_ldl_* / n_ldl_* fields of pt::PhaseTimes (conicIP does this when called with timing = pt); the record is shared across the factorizations of one kktsolver_ldl instance, so the attachment persists.

No rank assumption on G: dependent equality rows are handled by δ_e and the refinement. Semidefinite blocks are supported through a dense k×k block, which is O(k⁴) in memory and not meant for large k; choose_kktsolver keeps SDP problems on kktsolver_qr.

Returns solve3x3gen(F, F⁻ᵀ) per the conicIP KKT solver interface.

cached_kktsolver_ldl(; kwargs...)

A kktsolver_ldl that remembers the fill-reducing ordering of the last problem it saw and reuses it when the next problem has the same structure (dimensions, cone list, and sparsity patterns of Q, A, G), which is the situation in model-predictive control, sequential convex programming, or any loop that re-solves with new data. The AMD ordering is the only super-linear part of the symbolic setup; the rest is rebuilt in O(nnz) each solve. On banded and block-structured problems AMD is a few percent of a solve, so the saving is modest; it grows with the fill of the pattern. kwargs are forwarded to kktsolver_ldl.

ks = ConicIP.cached_kktsolver_ldl()
for t in 1:T
    sol = conicIP(Q, c[t], A, b[t], cone_dims; kktsolver = ks)
end
ks.hits   # number of solves that reused the ordering

soc_uv(Blk::SymWoodbury) -> (d², u, v)

Decompose the square of a diagonal-plus-rank-one scaling F = D + c·w·wᵀ (the second-order-cone Nesterov–Todd scaling from nestod_soc, with D = Diagonal([−β; β; …; β]) and c = 1) as

FᵀF = D² + u uᵀ − v vᵀ,     uᵀv = 0,

returning the vector d² = diag(D)² and the two vectors. Writing w̃ = √c·w, FᵀF − D² = B M Bᵀ with B = [D w̃ w̃] and M = [0 1; 1 w̃ᵀw̃]; det M = −1, so the rank-two term has exactly one positive and one negative eigenvalue when B has full column rank.

The construction goes through a thin Householder QR B = Qb R (Qb k×2 with orthonormal columns, R 2×2 upper triangular), so that B M Bᵀ = Qb (R M Rᵀ) Qbᵀ. The 2×2 symmetric matrix R M Rᵀ = V Λ Vᵀ is eigendecomposed and u = √λ₊·Qb V[:, +], v = √(−λ₋)·Qb V[:, −]. Because Qb V has orthonormal columns, uᵀv = 0 to rounding. Nothing is inverted: when Dw and w̃ are nearly (or exactly) parallel the second row of R is ~0 and the formula degrades gracefully to the rank-one term, with reconstruction error O(ε‖B‖²‖M‖) regardless of the conditioning of B. (The earlier route through T = BᵀB and T^{±1/2} squared the condition number and needed a scale-sensitive degeneracy branch that could discard part of the rank-two term.)

For the NT scaling the negative part satisfies vᵀv < β² (its eigenvalues on the cone's two-dimensional subspace multiply to β²), which is what keeps the lifted system quasi-definite.

Solves the 3x3 system

┌             ┐ ┌    ┐   ┌   ┐
│ Q   G'  -A' │ │ y' │ = │ y │
│ G           │ │ w' │   │ w │
│ A       FᵀF │ │ v' │   │ v │
└             ┘ └    ┘   └   ┘

by the double QR method described in CVXOPT http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf section 10.2

Solves the 3x3 system

┌             ┐ ┌    ┐   ┌   ┐
│ Q   G'  -A' │ │ y' │ = │ y │
│ G           │ │ w' │   │ w │
│ A       FᵀF │ │ v' │   │ v │
└             ┘ └    ┘   └   ┘

By lifting the large diagonal plus rank 3 blocks of FᵀF

Intelligently chooses between solve3x3gensparselift and solve3x3gensparsedense by approximating the number of non-zeros in both and choosing the form with more sparsity. The former is better for large second order cones, while the latter is better if the constraints are the product of many small cones.

Solves the 2x2 system

┌                   ┐ ┌    ┐   ┌   ┐
│ Q + A'F⁻¹F⁻ᵀA  G' │ │ y' │ = │ y │
│ G                 │ │ w' │   │ w │
└                   ┘ └    ┘   └   ┘

pivot(kktsolver_2x2)

Wrap a 2-by-2 KKT solver into a 3-by-3 solver by pivoting on the third component. The inner solver handles the Schur complement system; pivot reconstructs the full solution.

See also conicIP for the KKT solver interface specification.

Block(size::Int)
Block(Blk::Vector)

Block diagonal matrix type. Each diagonal block can be a different matrix type (Diagonal, SymWoodbury, VecCongurance, or dense Matrix).

Used internally to represent the Nesterov-Todd scaling matrix, where each block corresponds to a cone in the cone specification.

Supports arithmetic (*, +, -, inv, adjoint, ^), conversion to sparse and Matrix, and block-wise function application via broadcastf.

Indexing

• B[i] returns the i-th diagonal block
• B[i] = M sets the i-th diagonal block

block_idx(A::Block)

Return a vector of UnitRange{Int} giving the row/column index ranges for each diagonal block of A.

broadcastf(op, A::Block)
broadcastf(op, A::Block, B::Block)
broadcastf(op, A::Block, x::Union{Vector,Matrix})

Apply function op block-wise to the diagonal blocks of A (and optionally B or the corresponding segments of x).

Id(n)

Create an n-by-n identity matrix as Diagonal(ones(n)).

VecCongurance(R)

Linear operator representing a congruence transform in vectorized form. The action W * x computes vecm(R' * mat(x) * R).

Used internally as the Nesterov-Todd scaling matrix for semidefinite cones.

mat(x)

Convert a vectorized symmetric matrix (scaled lower-triangular form) back to a full symmetric matrix. Inverse of vecm.

mat!(Z, x)

In-place mat: fill the symmetric matrix Z from the vectorized form x. Z must be ord(x) square.

vecm(Z)

Vectorize a symmetric matrix Z into scaled lower-triangular form. Off-diagonal entries are scaled by √2 so that dot(vecm(X), vecm(Y)) == tr(X*Y). Inverse of mat.

vecm!(x, Z)

In-place vecm: write the vectorized form of the symmetric matrix Z into x, which must have length n(n+1)/2.

imcols(A, b, ϵ = 1e-8)

Removes redundant inequalities in a system of equations

Ax = b

and checks if the equations are consistent. Returns (R, consistent) where R are the indices of a maximal independent row set.

PhaseTimes()

Accumulator for the opt-in per-phase instrumentation of conicIP: pass one as timing = pt (or the MOI option "timing") and the solver adds wall time in nanoseconds (t_* fields), allocation bytes (b_*) and counts (n_*) for each phase; see the header of src/timing.jl for the field contract. Use one object per measured call, reset! to reuse it, and phase_table to print it. nothing (the default) leaves the solver uninstrumented.

reset!(pt::PhaseTimes)

Zero every field in place.

phase_table(pt::PhaseTimes) -> Vector{Pair{Symbol,Float64}}

Seconds per whole-call phase and per loop child, in the contract's order, for printing and CSV export. Diagnostics and counts are not included.

mul_adjoint!(y, A::Block, x)

Write A'x into y without forming the adjoint Block. Equivalent to mul!(y, A', x), and provided because A' on a Block is an eager broadcastf that copies the block vector.

inv_adjoint_block!(dest::Block, src::Block, i)

One block of inv_adjoint!. A Diagonal destination of the right length is overwritten in place — inv(::Diagonal) would otherwise allocate a vector per cone per iteration — with inv's singularity check kept. Any other block type is rebuilt.

SOCScratch(n)

Per-cone buffers for the in-place second-order-cone NT scaling: the SymWoodbury factors of F (j, w) and of F⁻ᵀ (ij, iw), the normalized iterates, and the internal temporaries both wrappers need.

nestod_soc!(sc::SOCScratch, z, s)

In-place nestod_soc: the same scaling matrix, with its factors written into sc instead of freshly allocated vectors.

soc_inv_adjoint!(sc::SOCScratch, W)

In-place adjoint(inv(W)) for the second-order-cone scaling block W built by nestod_soc! (a SymWoodbury of real type is its own adjoint, so this is inv(W) written into sc).

kkt_attach_timing!(solve3x3, pt::PhaseTimes)

Ask a KKT solver object to report its internal timings and counts into pt (the t_ldl_* / n_ldl_* fields). The default does nothing; kktsolver_ldl implements it for both the generator returned by kktsolver(Q, A, G, cone_dims) and the per-factorization object returned by solve3x3gen(F, F⁻ᵀ). The main loop calls it on the generator right after construction (so the initial factorization is counted) and on every factorization, only when timing is on.

@phase timing t_field b_field expr
@phase timing t_field expr

Evaluate expr, and when timing !== nothing add its wall time to timing.t_field (and its allocated bytes to timing.b_field when given). Timestamps are local to the expansion, so nested @phase blocks are safe.

Restrictions. expr is expanded twice (once per branch), so it must not define a named local function (the timed branch would see the method overwritten and throw UndefVarError); use an anonymous function or hoist the definition. A return, break or continue inside expr keeps its control flow but leaves the TIMED span unaccounted, because the update runs only when expr completes normally; where a span can exit early, use @phase_start / @phase_stop and close the span on each exit path.

(t0, b0) = @phase_start timing
@phase_stop timing t_field b_field t0 b0
@phase_stop timing t_field t0

Manual form of @phase for spans that can exit early: open once, close on every exit path. The stamps are (UInt64(0), Int64(0)) when timing === nothing, and no clock is read.

gc_start(timing) -> UInt64
gc_stop!(timing, gc0, x) -> x

t_gc bookkeeping for an entry point (conicIP, preprocess_conicIP, _preprocess_core, the MOI optimize!): read Base.gc_time_ns() on entry and ASSIGN the delta on every exit. Nested entry points each assign, and the outermost assigns last, so the recorded value is the outermost call's.

inv_adjoint!(dest::Block, src::Block)

Compute adjoint(inv(src)) block-wise, reusing the dest Block shell. Avoids allocating two intermediate Blocks for inv(F)'.

Wrapper around solve2xegen to solve 3x3 systems by pivoting on the third component.

Creates a matrix with the same sparsity structure as F

Checks if two sparse matrices have the same sparse structure

Estimates for the number of nonzeros of lift(F)

Estimates the number of nonzeros of F

Return (perm, is_offdiag) where perm[moi_k] is the vecm position for MOI triangle position moi_k, and is_offdiag[moi_k] is true when position moi_k corresponds to an off-diagonal entry.

kkt_diagnostics(solve3x3) -> diagnostics or nothing

Hook through which a KKT solver reports per-factorization diagnostics to conicIP. solve3x3 is the object a solver's solve3x3gen(F, F⁻ᵀ) returned; the default answers nothing (no diagnostics). A solver that returns an object with integer fields repaired and refactors (for the current factorization) and repaired_total and refactors_total (summed over the solve) has the former printed in the verbose kkt column as repaired/refactors and the latter stored in Solution.kkt_repaired and Solution.kkt_refactors. kktsolver_ldl implements it with LDLDiagnostics. Not exported.

Its companion kkt_attach_timing!(solve3x3, pt::PhaseTimes) (timing.jl) is the hook through which a KKT solver reports per-call timings: when conicIP runs with timing = pt it calls it once per factorization, right after solve3x3gen, and the solver may from then on add to pt.t_ldl_factor / pt.n_ldl_factor (numeric refactorizations, wall nanoseconds and count), pt.t_ldl_solve / pt.n_ldl_solve (triangular solves with the factorization) and pt.t_ldl_resid / pt.n_ldl_resid (residual evaluations of its internal refinement). These are inclusive diagnostics, never summed with the loop phases that contain them. Both hooks are optional: the defaults answer nothing and do nothing, so a custom kktsolver need not know about either; a solver that implements kkt_attach_timing! must keep its untimed path free of the timing work when nothing is attached.

LDLDiagnostics

Per-factorization and per-solve counters of kktsolver_ldl, reachable through kkt_diagnostics(solve3x3) on the object its solve3x3gen returns.

• repaired – pivots QDLDL's dynamic regularization replaced in the current factorization; repaired_total sums them over the solve
• pos_inertia – positive pivots of the current factorization. Recorded only: with Dsigns QDLDL forces every pivot's sign, so the count always matches the quasi-definite pattern and detects nothing
• δp, δe, δc – static shifts in effect (δp, δe grow under the retry policy; δc is fixed)
• refactors – shift bumps applied to the current factorization; refactors_total sums them over the solve
• last_residual – unregularized residual norm ‖rhs − K₀x‖ of the last solve3x3 return, over the LIFTED system (auxiliary rows included)
• last_rtol – the tolerance refine_tol·(1 + ‖rhs‖) that the internal refinement of that solve aimed at, so a caller can tell a solve that met its own target from one that gave up
• last_bound – upper bound on the residual norm of the UNLIFTED 3×3 system ‖(bx, by, −bz) − K₃ₓ₃ x₃ₓ₃‖ for the same solution. It is ‖r_u‖ + lift_gain·‖r_a‖, where r_u/r_a split the lifted residual into its 3×3 rows and its auxiliary rows (eliminating the auxiliaries a = uᵀz, b = vᵀz from a residual (r_u, r_a) leaves r_u + [u v]·r_a on the 3×3 rows, and [u v] per lifted block has spectral norm at most √(‖u‖² + ‖v‖²)). Each norm is measured on its own rows, so last_bound equals last_residual when nothing was lifted and stays faithful when one part is orders of magnitude below the other
• lift_gain – maxᵦ √(‖uᵦ‖² + ‖vᵦ‖²) over the lifted second-order-cone blocks of the current factorization; 0 when nothing is lifted
• timing – the PhaseTimes the backend reports into, or nothing (the default). Set by kkt_attach_timing!; while attached, every numeric refactorization, triangular solve, and residual evaluation adds its wall time and a count to the t_ldl_* / n_ldl_* fields. With nothing each of those sites costs one pointer comparison.

spectral_map!(out, w, f, cone_dims)

Apply the scalar map f to the Jordan-algebraic eigenvalues of w, block by block over the cone product cone_dims, keeping the Jordan frame of w. Writes the result into out (which may alias w) and returns it.

clip_spectral!(out, w, lo, hi, cone_dims)

Project the Jordan-algebraic eigenvalues of w onto [lo, hi] in the frame of w, block by block over cone_dims ("R": entrywise clamp; "Q": both eigenvalues w₁ ± ‖w̄‖; "S": the eigenvalues of mat(w)). The result is written into out and returned; out may alias w.

centrality_correction!(Δw, w, lo, hi, cap, cone_dims)

Gondzio's corrector target for the trial complementarity w: the spectral difference Π_[lo,hi](w) − w, with every component below −cap raised to −cap (in the frame of w). Returns Δw; the correction is identically zero when every eigenvalue of w already lies in the box.

Convert a vector from vecm order (solver convention) to MOI triangle order, dividing off-diagonal entries by √2.