08. How does a solver reach the boundary?

Duality explained what primal and dual solutions must satisfy, but our code still jumped from optimize! to an answer. What happens between those calls? Production SDP solvers contain many safeguards and transformations, but one central idea can be seen on a problem small enough to solve by hand.

Why is a boundary solution awkward?

Consider

\[ \begin{aligned} \text{minimize}\quad &t\\ \text{subject to}\quad& \begin{bmatrix}1&t\\t&1\end{bmatrix}\succeq0. \end{aligned} \]

The feasible interval is \(-1\le t\le1\), so the solution is \(t=-1\). The optimal matrix is singular: one eigenvalue has reached zero. Algorithms based on inverses and Newton systems are better behaved at positive definite matrices, where every eigenvalue remains strictly positive.

A standard barrier for a positive definite matrix is \(-\mu\log\det X\). It grows without bound as an eigenvalue approaches zero, discouraging an iterate from leaving the interior. On this slice,

\[ \det \begin{bmatrix} 1&t\\ t&1 \end{bmatrix} =1-t^2, \]

so the matrix log-determinant barrier becomes the scalar problem

\[ \text{minimize}\quad \phi_\mu(t)=t-\mu\log(1-t^2),\qquad -1<t<1. \]

The logarithm tends to infinity near either boundary. The parameter \(\mu>0\) balances progress in the objective against distance from the boundary.

Can we approach without touching it?

For this scalar example,

\[ \phi_\mu'(t)=1+\frac{2\mu t}{1-t^2},\qquad \phi_\mu''(t)=\frac{2\mu(1+t^2)}{(1-t^2)^2}. \]

Newton’s method proposes \(t-\phi_\mu'(t)/\phi_\mu''(t)\). A backtracking step keeps the iterate inside the feasible interval and requires the barrier objective to decrease.

function barrier_value(t, μ)
    abs(t) < 1 || return Inf
    return t - μ * log(1 - t^2)
end

function center(μ; start = 0.0, tolerance = 1e-12)
    t = start
    history = Float64[t]
    for _ in 1:50
        gradient = 1 + 2μ * t / (1 - t^2)
        hessian = 2μ * (1 + t^2) / (1 - t^2)^2
        abs(gradient) <= tolerance && break
        direction = -gradient / hessian
        α = 1.0
        current = barrier_value(t, μ)
        accepted = false
        for _ in 1:60
            candidate = t + α * direction
            sufficient_decrease = current + 1e-4 * α * gradient * direction
            if abs(candidate) < 1 && barrier_value(candidate, μ) <= sufficient_decrease
                t = candidate
                accepted = true
                break
            end
            α /= 2
        end
        accepted || error("backtracking failed to find an interior descent step")
        push!(history, t)
    end
    return (; t, history)
end

path = [(μ = μ, center = center(μ).t) for μ in (1.0, 0.3, 0.1, 0.03, 0.01)]
path
5-element Vector{@NamedTuple{μ::Float64, center::Float64}}:
 (μ = 1.0, center = -0.41421356237309664)
 (μ = 0.3, center = -0.7440306508910554)
 (μ = 0.1, center = -0.9049875621120985)
 (μ = 0.03, center = -0.9704498987955369)
 (μ = 0.01, center = -0.9900499987500625)

As \(\mu\) shrinks, the centers approach -1 while remaining strictly inside the cone. This sequence is a simple central path. The method does not need to land on the singular boundary to approximate the optimal value accurately.

What must a production solver add?

The scalar barrier reveals the purpose of interior iterates, Newton directions, step lengths, and decreasing barrier parameters. A conic solver must also:

  • maintain primal and dual variables together;
  • solve a structured linear system at each iteration;
  • scale data and regularize nearly singular systems;
  • decide when residuals and the duality gap are small enough;
  • detect or report evidence of infeasibility and unboundedness.

Clarabel is a primal-dual interior-point conic solver. JuMP translates a model into MathOptInterface objects; bridges may reformulate those objects; the solver then receives cone data. These are separate layers. A model that prints nicely in JuMP is not necessarily stored by the solver in the same surface form.

NoteGo further: algorithms and software layers

The algorithm sections of Vandenberghe and Boyd’s SDP survey develop primal-dual interior-point methods. Clarabel’s official documentation describes its problem format, settings and tolerances, and the linear systems solved at each iteration. Use the former for mathematical depth and the latter to interpret this course’s particular solver—not as a substitute for the mathematics.

What can a solver log tell us?

using JuMP
using Clarabel
using LinearAlgebra

model = Model(Clarabel.Optimizer)
set_optimizer_attribute(model, "max_iter", 20)
@variable(model, X[1:2, 1:2], PSD)
@constraint(model, X[1, 1] == 1)
@constraint(model, X[2, 2] == 1)
@objective(model, Min, X[1, 2])
optimize!(model)

(status = termination_status(model),
 t = value(X[1, 2]),
 eigenvalues = eigvals(Symmetric(value.(X))))
-------------------------------------------------------------
           Clarabel.jl v0.11.1  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 3
  constraints   = 5
  nnz(P)        = 0
  nnz(A)        = 5
  cones (total) = 2
    : Zero        = 1,  numel = 2
    : PSDTriangle = 1,  numel = 3

settings:
  linear algebra: direct / qdldl, precision: 64 bit (1 thread)
  max iter = 20, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, 
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0   0.0000e+00  -0.0000e+00  0.00e+00  2.03e-02  4.47e-01  1.00e+00  1.35e+00   ------   
  1  -4.8817e-01  -4.8735e-01  8.20e-04  3.56e-03  1.05e-01  1.84e-01  2.55e-01  8.24e-01  
  2  -9.9456e-01  -9.9363e-01  9.28e-04  7.76e-05  2.40e-03  5.56e-03  6.43e-03  9.90e-01  
  3  -9.9995e-01  -9.9994e-01  9.28e-06  7.75e-07  2.39e-05  5.56e-05  6.43e-05  9.90e-01  
  4  -1.0000e+00  -1.0000e+00  9.28e-08  7.75e-09  2.39e-07  5.56e-07  6.43e-07  9.90e-01  
  5  -1.0000e+00  -1.0000e+00  9.28e-10  7.75e-11  2.39e-09  5.56e-09  6.43e-09  9.90e-01  
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  425ms
(status = OPTIMAL, t = -0.9999999945574471, eigenvalues = [5.442559125601321e-9, 1.9999999945574534])

The log is diagnostic evidence about an iterative process. The returned point still needs the independent feasibility and objective checks developed in Chapter 09.

Try it yourself

  1. Complete learner/08_how_a_solver_moves.jl.
  2. Run the barrier path with a starting point close to 1. Record which step sizes backtracking accepts.
  3. Derive the stationary point of \(\phi_\mu\) analytically and compare it with the Newton result.
  4. Explain why an optimum on the PSD boundary is normal rather than evidence that the model failed.