SDP Relaxation
Quantum Density Matrix Optimization
Spectraplex(n) exactly represents the real density matrix feasible set
\[\{ X \in \mathbb{S}_+^n : \operatorname{tr}(X) = 1 \}.\]
This page uses that set directly: optimize the energy of a real symmetric Hamiltonian over density matrices. Problems such as the standard MaxCut SDP add extra constraints like $\operatorname{diag}(X) = 1$ and therefore require a richer oracle than Spectraplex alone.
For simplicity we work with real symmetric matrices. The package stores the matrix variable in vectorized form, so x = vec(X).
Random Hamiltonian
using Marguerite, LinearAlgebra, Random, UnicodePlots
Random.seed!(42)
n = 20
A = randn(n, n)
H = Symmetric((A + A') / 2)
H_vec = vec(Matrix(H))
println("Hilbert space dimension = $n")
println("Smallest eigenvalue of H = ", round(eigmin(H); sigdigits=6))Hilbert space dimension = 20
Smallest eigenvalue of H = -5.92417Linear Ground State SDP
The ground state problem is the linear SDP
\[\min_{X \succeq 0,\; \operatorname{tr}(X)=1} \operatorname{tr}(H X).\]
Its optimum is the minimum eigenvalue of $H$, attained by the rank-1 projector onto a ground state eigenvector. Frank-Wolfe finds that projector in one step.
m = n * n
lmo = Spectraplex(n)
X0 = Matrix{Float64}(I, n, n) / n
x0 = vec(X0)
f_linear(x) = dot(H_vec, x)
∇f_linear!(g, x) = (g .= H_vec)
x_lin, res_lin = solve(f_linear, lmo, x0;
grad=∇f_linear!, max_iters=100, tol=1e-10, verbose=false)
println("FW energy = ", round(res_lin.objective; sigdigits=6))
println("Ground state energy = ", round(eigmin(H); sigdigits=6))
println("FW gap = ", round(res_lin.gap; sigdigits=4))
println("Iterations = ", res_lin.iterations)FW energy = -5.92417
Ground state energy = -5.92417
FW gap = 2.665e-15
Iterations = 1Regularized Mixed State Problem
Adding a Frobenius penalty produces a nontrivial convex SDP whose solution is typically mixed rather than rank 1:
\[\min_{X \succeq 0,\; \operatorname{tr}(X)=1} \operatorname{tr}(H X) + \frac{\eta}{2}\|X\|_F^2.\]
eta = 1.5
f_reg(x) = dot(H_vec, x) + (eta / 2) * dot(x, x)
∇f_reg!(g, x) = (g .= H_vec .+ eta .* x)
x_reg, res_reg = solve(f_reg, lmo, x0;
grad=∇f_reg!, max_iters=1000, tol=1e-8, verbose=false)
println("Objective = ", round(res_reg.objective; sigdigits=6))
println("FW gap = ", round(res_reg.gap; sigdigits=4))
println("Iterations = ", res_reg.iterations)Objective = -5.24207
FW gap = 0.0004294
Iterations = 1000Convergence Trace
The hand-written Frank-Wolfe loop below records the duality gap for the regularized problem:
max_iters = 500
x = copy(x0)
g_buf = zeros(m)
v_buf = zeros(m)
step = MonotonicStepSize()
gaps = Vector{Float64}(undef, max_iters)
for t in 0:max_iters-1
∇f_reg!(g_buf, x)
lmo(v_buf, g_buf)
gaps[t+1] = dot(g_buf, x .- v_buf)
gamma = step(t)
x .= x .+ gamma .* (v_buf .- x)
end
scatterplot(1:max_iters, gaps;
yscale=:log10,
title="FW Gap — Regularized Density Matrix SDP (n=$n, eta=$eta)",
xlabel="iteration", ylabel="gap",
name="FW gap", width=60) ⠀⠀⠀⠀FW Gap — Regularized Density Matrix SDP (n=20, eta=1.5)⠀⠀⠀
┌────────────────────────────────────────────────────────────┐
10⁰⸱⁷⁶²⁸⁵ │⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ FW gap
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠅⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠤⢂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠥⠦⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠘⢤⠲⢩⠕⢂⣄⡀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
gap │⠀⠀⢀⠃⡘⡱⠡⠺⢪⠪⣒⢔⣔⢄⢄⣀⠀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠠⢀⢀⠩⠡⠁⢅⢂⠦⠡⡓⢑⠣⡪⡑⢝⠵⢕⢕⠶⡲⣢⡢⡢⢤⢤⢄⣄⢄⡠⣀⠄⡀⣀⢀⡀⡀⡀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
│⠀⠀⠀⠀⠁⠀⠀⡀⠔⠈⠀⠄⢁⠀⠆⠡⠐⡨⠂⡌⠤⠢⢑⠠⡑⡈⢌⠢⢁⠲⡑⢝⠊⢆⠪⠪⣑⠍⡐⢍⠺⡚⠝⡱⢝⢫⡪⠱⡣⢶⠵⡮⢢⡲⢖⠶⢢⠦⡢⢖│
│⠀⠀⠀⠁⠀⠄⠀⠀⠀⠀⠁⠠⠁⠀⠈⠀⠂⠀⠄⠠⠈⡀⠁⡀⢁⠈⡀⠊⡀⢂⠰⡀⢐⡀⢂⠑⡈⠐⡉⢊⠑⢐⠑⡌⢂⠑⢀⠑⠌⠢⠑⠢⠕⠌⠢⠁⠣⡑⢌⠢│
│⠀⠀⠄⠀⠀⠀⠂⠀⠀⠂⠀⠀⠀⠐⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⡀⡀⢀⠀⠄⠀⠄⠐⠀⠃⠠⠂⠐⠀⠈⡀⢁⠈⠅⠘⠄⠠⠀⠐│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⢀⠀⠀⠀⠀⠀⠀⡀⠀⠄⠠⠀⠂⠀⠁⠈⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⡀⠀⠄⠐⠁⠁⠀⠀⠀⠀⠀⠀⠠⠁⠂⠀│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠈⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠄⠀⠁⠀⠀⠀⠁│
10⁻⁴⸱⁷⁷⁰⁴² │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀│
└────────────────────────────────────────────────────────────┘
⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀500⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀iteration⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ Solution Structure
The regularized solution is still PSD with unit trace, but it can spread mass across several eigenvectors.
X_sol = reshape(x_reg, n, n)
X_sym = Symmetric((X_sol .+ X_sol') ./ 2)
lambda = eigvals(X_sym)
lambda_top = sort(lambda; rev=true)[1:min(10, n)]
println("Trace = ", round(tr(X_sym); sigdigits=6))
println("Min eig = ", round(minimum(lambda); sigdigits=6))
println("Purity = ", round(dot(x_reg, x_reg); sigdigits=6))
println("Numerical rank = ", count(>(1e-6), lambda), " / ", n)
barplot(["λ$i" for i in 1:length(lambda_top)], lambda_top;
title="Largest eigenvalues of X*",
xlabel="value", width=60) Largest eigenvalues of X*
┌ ┐
λ1 ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.787431
λ2 ┤■■■■■■■■■■■■■ 0.212569
λ3 ┤ 1.72361e⁻¹⁶
λ4 ┤ 6.1829e⁻¹⁷
λ5 ┤ 3.49199e⁻¹⁷
λ6 ┤ 2.98395e⁻¹⁷
λ7 ┤ 1.91492e⁻¹⁷
λ8 ┤ 1.09509e⁻¹⁷
λ9 ┤ 5.53531e⁻¹⁸
λ10 ┤ 2.3308e⁻¹⁸
└ ┘
value Spectraplex also supports implicit differentiation through these solves via a compact active_set representation; see Bilevel Optimization for the differentiated interface.
Why Frank-Wolfe Fits This SDP
Frank-Wolfe is effective for density matrix SDPs because:
- Cheap oracle: the linear minimization oracle only needs the minimum eigenvector of the current gradient matrix.
- Low-rank iterates: each FW step adds one rank-1 projector, so early iterates stay compact.
- No projections: projection onto the spectraplex needs a full eigendecomposition, while the FW step only solves a leading-eigenvector subproblem.
- Memory efficiency: the method stores a few dense matrices and vectors, not a full SDP interior-point system.