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.
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,
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.
functionbarrier_value(t, μ)abs(t) <1||returnInfreturn t - μ *log(1- t^2)endfunctioncenter(μ; start =0.0, tolerance =1e-12) t = start history =Float64[t]for _ in1:50 gradient =1+2μ * t / (1- t^2) hessian =2μ * (1+ t^2) / (1- t^2)^2abs(gradient) <= tolerance &&break direction =-gradient / hessian α =1.0 current =barrier_value(t, μ) accepted =falsefor _ in1:60 candidate = t + α * direction sufficient_decrease = current +1e-4* α * gradient * directionifabs(candidate) <1&&barrier_value(candidate, μ) <= sufficient_decrease t = candidate accepted =truebreakend α /=2end accepted ||error("backtracking failed to find an interior descent step")push!(history, t)endreturn (; t, history)endpath = [(μ = μ, 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.
-------------------------------------------------------------
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
Complete learner/08_how_a_solver_moves.jl.
Run the barrier path with a starting point close to 1. Record which step sizes backtracking accepts.
Derive the stationary point of \(\phi_\mu\) analytically and compare it with the Newton result.
Explain why an optimum on the PSD boundary is normal rather than evidence that the model failed.
---title: "08. How does a solver reach the boundary?"engine: juliajulia: exeflags: ["--project=@."]---Duality explained what primal and dual solutions must satisfy, but our codestill jumped from `optimize!` to an answer. What happens between those calls?Production SDP solvers contain many safeguards and transformations, but onecentral 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$. Theoptimal matrix is singular: one eigenvalue has reached zero. Algorithms basedon inverses and Newton systems are better behaved at positive definitematrices, 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 stepkeeps the iterate inside the feasible interval and requires the barrierobjective to decrease.```{julia}#| label: miniature-barrier-methodfunctionbarrier_value(t, μ)abs(t) <1||returnInfreturn t - μ *log(1- t^2)endfunctioncenter(μ; start =0.0, tolerance =1e-12) t = start history =Float64[t]for _ in1:50 gradient =1+2μ * t / (1- t^2) hessian =2μ * (1+ t^2) / (1- t^2)^2abs(gradient) <= tolerance &&break direction =-gradient / hessian α =1.0 current =barrier_value(t, μ) accepted =falsefor _ in1:60 candidate = t + α * direction sufficient_decrease = current +1e-4* α * gradient * directionifabs(candidate) <1&&barrier_value(candidate, μ) <= sufficient_decrease t = candidate accepted =truebreakend α /=2end accepted ||error("backtracking failed to find an interior descent step")push!(history, t)endreturn (; t, history)endpath = [(μ = μ, center =center(μ).t) for μ in (1.0, 0.3, 0.1, 0.03, 0.01)]path```As $\mu$ shrinks, the centers approach `-1` while remaining strictly insidethe cone. This sequence is a simple **central path**. The method does not needto 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 modelinto MathOptInterface objects; bridges may reformulate those objects; the solverthen receives cone data. These are separate layers. A model that prints nicelyin JuMP is not necessarily stored by the solver in the same surface form.::: {.callout-note title="Go further: algorithms and software layers"}The algorithm sections of[Vandenberghe and Boyd's SDP survey](https://web.stanford.edu/~boyd/papers/sdp.html)develop primal-dual interior-point methods. Clarabel's official documentationdescribes its [problem format](https://clarabel.org/stable/examples/jl/example_sdp/),[settings and tolerances](https://clarabel.org/stable/api_settings/), and the[linear systems solved at each iteration](https://clarabel.org/stable/user_guide_linsolvers/).Use the former for mathematical depth and the latter to interpret this course'sparticular solver—not as a substitute for the mathematics.:::## What can a solver log tell us?```{julia}#| label: visible-solver-logusingJuMPusingClarabelusingLinearAlgebramodel =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))))```The log is diagnostic evidence about an iterative process. The returned pointstill needs the independent feasibility and objective checks developed inChapter 09.## Try it yourself1. 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.