Quartically regularized polynomials

Adapted from: Section 3.3.1 of (Cartis and Zhu, 2026)

Higher-order optimization methods minimize, at each iteration, a cubic Taylor model regularized by a quartic term:

\[m(s) = f_0 + g^\top s + \frac{1}{2} H[s]^2 + \frac{1}{6} T[s]^3 + \frac{\sigma}{4} \|s\|_2^4.\]

For any $\sigma > 0$, this model is bounded below so $m(s) - m^*$ is nonnegative where $m^*$ is its global minimum. If $m(s) - m^*$ is moreover a sum of squares then $m^*$ can be computed exactly with a single Sum-of-Squares program, i.e., the zeroth level of the Lasserre hierarchy. (Cartis and Zhu, 2026) shows that this is the case for $\sigma$ large enough and that, when $T = 0$, it is the case for all $\sigma > 0$. They also show that the geometry of the regularization matters: with the separable quartic regularization $\|s\|_4^4 = s_1^4 + \cdots + s_n^4$ instead of the Euclidean one $\|s\|_2^4 = (s_1^2 + \cdots + s_n^2)^2$, $m(s) - m^*$ may fail to be a sum of squares for every $\sigma > 0$. In this tutorial, we verify this numerically on their example.

using LinearAlgebra
using DynamicPolynomials
@polyvar s[1:3]

(DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[s₁, s₂, s₃],)

We consider the following nonconvex quadratic part whose Hessian has eigenvalues $-4$, $-4$ and $20$:

quad = 2 * sum(s .^ 2) + 8 * (s[1] * s[2] + s[1] * s[3] + s[2] * s[3])

\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} \]

We need to pick an SDP solver, see here for a list of the available choices.

using SumOfSquares
import Clarabel
solver = optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true)

MathOptInterface.OptimizerWithAttributes(Clarabel.MOIwrapper.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])

The following function computes the largest $\gamma$ such that multiplier * (p - γ) is a sum of squares. With the default multiplier, this is the zeroth level of the Lasserre hierarchy. It returns the lower bound $\gamma$ and the constraint reference, which we will use to extract minimizers from the moment matrix.

function sos_lower_bound(p; multiplier = 1)
    model = SOSModel(solver)
    @variable(model, γ)
    @objective(model, Max, γ)
    con_ref = @constraint(model, multiplier * (p - γ) >= 0)
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    return value(γ), con_ref
end

sos_lower_bound (generic function with 1 method)

Separable quartic regularization

We first consider the separable regularization with $\sigma = 4$. The resulting polynomial was shown in (Ahmadi et al., 2023; Theorem 3.3) to be such that $m_A(s) - m_A^*$ is not a sum of squares.

m_A = quad + sum(s .^ 4)

\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} + s_{3}^{4} + s_{2}^{4} + s_{1}^{4} \]

The zeroth level of the Lasserre hierarchy gives the following lower bound:

γ_A, _ = sos_lower_bound(m_A)
γ_A

-2.385952772480988

Is this lower bound tight ? Recall that a polynomial $p$ is nonnegative if $(s_1^2 + s_2^2 + s_3^2) p(s)$ is a sum of squares, see the Motzkin tutorial. Using this multiplier, we get a larger lower bound:

γ_A_mult, con_ref = sos_lower_bound(m_A, multiplier = sum(s .^ 2))
γ_A_mult

-2.1442784188346504

We can extract the minimizers from the moment matrix:

atoms = atomic_measure(moment_matrix(con_ref), 1e-4)

Atomic measure on the variables s[1], s[2], s[3] with 7 atoms:
 at [-0.667400995328262, 1.1497302403031486, -0.667398825519008] with weight 0.07528774640460814
 at [0.6674010075910534, 0.6673988228757881, -1.1497302568166359] with weight 0.07538269309709689
 at [-1.1497296867186602, 0.6673993968975332, 0.6673994305402426] with weight 0.0752951814451945
 at [1.1497296590300567, -0.6673993923853224, -0.6673993994935465] with weight 0.0752951888847822
 at [0.6674004416444095, -1.1497300395958414, 0.6673990095967173] with weight 0.0752878001957122
 at [-0.6674004262185959, -0.6673990280953056, 1.1497300416922287] with weight 0.07538275142485673
 at [0.0, 0.0, 0.0] with weight 0.0

We recover the six global minimizers $\pm(-1.1498, 0.6674, 0.6674)$ and their permutations. The moment matrix also contains an atom at the origin with zero weight. It comes from the multiplier which vanishes at the origin, so we ignore it.

minimizers = [atom.center for atom in atoms.atoms if atom.weight > 1e-4]

6-element Vector{Vector{Float64}}:
 [-0.667400995328262, 1.1497302403031486, -0.667398825519008]
 [0.6674010075910534, 0.6673988228757881, -1.1497302568166359]
 [-1.1497296867186602, 0.6673993968975332, 0.6673994305402426]
 [1.1497296590300567, -0.6673993923853224, -0.6673993994935465]
 [0.6674004416444095, -1.1497300395958414, 0.6673990095967173]
 [-0.6674004262185959, -0.6673990280953056, 1.1497300416922287]

Evaluating $m_A$ at the minimizers confirms that this is the global minimum and that the zeroth level of the hierarchy has a gap:

m_A_min = minimum(x -> m_A(s => x), minimizers)
m_A_min

-2.1442784126523144

As shown in (Cartis and Zhu, 2026; Lemma 3.3), the gap does not close as $\sigma$ increases. Indeed, with $\tilde{s} = \sqrt{\sigma}s/2$, we have $\mathrm{quad}(s) + \sigma \|s\|_4^4/4 = 4 m_A(\tilde{s})/\sigma$ so the ratio between the lower bound and the minimum is independent of $\sigma$:

for σ in [1, 10, 100]
    m = quad + σ / 4 * sum(s .^ 4)
    γ, _ = sos_lower_bound(m)
    println("σ = $σ: lower bound = $γ, minimum = $(4m_A_min / σ)")
end

σ = 1: lower bound = -9.54381117121182, minimum = -8.577113650609258
σ = 10: lower bound = -0.9543811237701892, minimum = -0.8577113650609258
σ = 100: lower bound = -0.09543811863221864, minimum = -0.08577113650609258

Euclidean quartic regularization

Let us now replace the separable regularization with the Euclidean one:

σ = 4
m_E = quad + σ / 4 * sum(s .^ 2)^2

\[ 2.0s_{3}^{2} + 8.0s_{2}s_{3} + 2.0s_{2}^{2} + 8.0s_{1}s_{3} + 8.0s_{1}s_{2} + 2.0s_{1}^{2} + s_{3}^{4} + 2.0s_{2}^{2}s_{3}^{2} + s_{2}^{4} + 2.0s_{1}^{2}s_{3}^{2} + 2.0s_{1}^{2}s_{2}^{2} + s_{1}^{4} \]

This is a quadratic model with Euclidean quartic regularization ($T = 0$) so $m_E(s) - m_E^*$ is a sum of squares for every $\sigma > 0$. Indeed, it can be checked that

\[m_E(s) + \frac{4}{\sigma} = \frac{\sigma}{4}\left(\|s\|_2^2 - \frac{4}{\sigma}\right)^2 + 4(s_1 + s_2 + s_3)^2.\]

The zeroth level of the hierarchy finds this lower bound $-4/\sigma = -1$:

γ_E, con_ref = sos_lower_bound(m_E)
γ_E

-0.9999999989968844

The global minimizers form the circle $s_1 + s_2 + s_3 = 0$, $\|s\|_2 = 2/\sqrt{\sigma}$. Even if there are infinitely many of them, the low rank decomposition of the moment matrix still finds a few points of this circle:

atoms = atomic_measure(moment_matrix(con_ref), 1e-4)

Atomic measure on the variables s[1], s[2], s[3] with 4 atoms:
 at [0.5773381513376452, 0.21135177661655416, -0.7886899275789001] with weight 0.24602530944083023
 at [-0.5773381504170676, 0.7886899270637537, -0.21135177675460198] with weight 0.24600154846920752
 at [-0.5773381513376424, -0.2113517766165557, 0.7886899275788989] with weight 0.24602530944083129
 at [0.5773381504170648, -0.7886899270637523, 0.21135177675460348] with weight 0.24600154846920813

We can verify that they belong to the circle:

for atom in atoms.atoms
    x = atom.center
    println("sum = $(sum(x)), norm = $(norm(x))")
end

sum = 3.752992361327756e-10, norm = 1.0000103581130317
sum = -1.0791584292846323e-10, norm = 1.0000103572044425
sum = -3.752992361327756e-10, norm = 1.0000103581130295
sum = 1.0791600946191693e-10, norm = 1.00001035720444

Adding a cubic term

We now add a small cubic term $0.1(s_1^3 + s_2^3 + s_3^3)$ to both models. With the separable regularization, there is still a gap between the lower bound of the zeroth level and the global minimum:

cubic = 0.1 * sum(s .^ 3)
m_A2 = m_A + cubic
γ_A2, _ = sos_lower_bound(m_A2)
γ_A2

-2.39899054004181

The multiplier closes the gap again and gives the global minimum:

γ_A2_mult, con_ref = sos_lower_bound(m_A2, multiplier = sum(s .^ 2))
γ_A2_mult

-2.2409484856506636

With the Euclidean regularization, the model is now a cubic model so (Cartis and Zhu, 2026; Theorem 2.3) only guarantees that $m(s) - m^*$ is a sum of squares for $\sigma$ large enough. We compare the lower bound found by the zeroth level of the hierarchy with the value of $m$ at the minimizers extracted from the moment matrix. If they match, then the lower bound is the global minimum.

for σ in [0.1, 1, 4, 10, 40]
    m = quad + cubic + σ / 4 * sum(s .^ 2)^2
    γ, c = sos_lower_bound(m)
    μ = atomic_measure(moment_matrix(c), 1e-4)
    m_min = minimum(atom -> m(s => atom.center), μ.atoms)
    println("σ = $σ: lower bound = $γ, value at extracted minimizers = $m_min")
end

σ = 0.1: lower bound = -53.43765171392503, value at extracted minimizers = -53.437653996164805
σ = 1.0: lower bound = -4.353277438110028, value at extracted minimizers = -4.35327756450587
σ = 4.0: lower bound = -1.042438226705649, value at extracted minimizers = -1.0424382384182889
σ = 10.0: lower bound = -0.4105830614302021, value at extracted minimizers = -0.410583051107534
σ = 40.0: lower bound = -0.1013068122685663, value at extracted minimizers = -0.1013067761618109

For this example, the zeroth level is exact even for small values of $\sigma$. The theoretical threshold of (Cartis and Zhu, 2026; Theorem 2.3) is sufficient but not necessary.

This page was generated using Literate.jl.