ConvexOptimization: lower quadratic atoms through the rotated SOC

[ChrisRackauckas-Claude]

Summary

ConvexMOI could not solve least squares or QPs: SymbolicAnalysis certifies the quadratic atoms convex, but they were missing from LOWERABLE_ATOMS, so no epigraph was ever formed (and u'Pu surfaced as UnknownCurvature, sum(abs2.(w)) as an internal invalid index error). This PR adds abs2(w), w^2, sum(abs2.(w)) / sum(w .^ 2) / sum(abs2, w), SymbolicAnalysis.quad_form(u, P), u' * P * u, and self-products v' * v (u'*u, (A*u)'*(A*u), (A*u-b)'*(A*u-b)) as lowerable atoms through MOI.RotatedSecondOrderCone (‖w‖² <= τ is (τ, 1/2, w); u'Pu is factored as sym(P) = LᵀL once at build time). Unsound cases are rejected with clear errors: non-PSD or parameter-dependent P, non-affine atom arguments, -abs2 under MinSense, nested norm(w)^2, and non-square powers.

The quadratic-form matcher (_quad_form_parts) also handles how Symbolics flattens real input: a numeric scalar folded anywhere in the product (u'*(2u), (2u)'*u, u'*u*2, (A*u)'*(2*A*u), u'*A'*A*u standalone) is stripped and folded into P — so P = k·I for self-products and P = k·∏Mᵢ for quad forms, which routes both magnitude and sign through the PSD check (negative scale → rejected as non-PSD). A p-dependent scale correctly fails the numeric check. Products that cannot be traced to a scalar objective at all (u'*A*B*u - u1) are caught and re-thrown as a "could not be traced" error naming the supported spellings instead of leaking the internal shape ArgumentError.

A residual-objective pass (_expand_scalar_products) evaluates leftover scalar * terms over array arguments, so pure LPs such as c'*u, c'*M*u, (M*u)'*c, dot(c,u), sum(c .* u), and c'*(u .- 1) also solve (they previously failed with UnknownCurvature).

Parameter-only squares (p[1]^2, abs2(p[1]) with no optimization variable in the argument) are not epigraph atoms: each distinct one is lifted to an implicit parameter column evaluated at every θ by a compiled function (Symbolics.build_function over p), so non-polynomial arguments (exp(p1)^2, abs2(sin(p1)), (p1/p2)^2) work too. u[1] + p[1]^2 solves, u[1]^2 - p[1]^2 gets the right constant with the right sign (rather than a misleading epigraph-sign error), they work inside atom arguments and constraints, and reinit! recomputes them.

sym(P) is accepted only when its smallest eigenvalue is above -n*eps(Float64)*λmax; a strictly negative eigenvalue inside that band is clamped with a @debug log (rank-deficient BᵀB inputs hit the band routinely, so it must not warn), and anything more negative (e.g. [1e6 0; 0 -1e-6]) is rejected as non-convex.

Wrong-answer bug found in review, fixed

Round-3 review caught (A*u)' * (2*A*u) - c'*u reporting -2.1667 — the optimum of u'Mu - c'u, i.e. the factor 2 was dropped. Root cause: in _quad_form_parts, the self-product branch (rest multiplies back to v) returned mid = nothing, discarding the scalar scale that _strip_scalar_mul had extracted, while the v'*P*v branch folded it in. Fix: the self-product branch returns mid = [scale] (unless scale == 1), so P = scale·I and the PSD check sees the sign. N-argument products (adjoint(2*A*u) contains *(2, A, u) after Symbolics refactors 2A back out) are now stripped too. Every successful solve in core_tests.jl now goes through solve_checked/solve!_checked, which asserts f(sol.u, p) ≈ sol.objective — the guard that would have caught this (~120 assertions added).

Merge with #1367 (piecewise-linear atoms)

This branch merges origin/master at aaf82fe33 (ConvexOptimization v0.1.2, #1367 merged) — a merge commit, not a rebase, so no force-push. Both families coexist: LOWERABLE_ATOMS lists abs/max/min alongside abs2/quad_form; Mapreducer terms dispatch to the piecewise-linear check first, then the sum-of-squares check; _atom_lowering sends sum-of-squares reducers to RSOC and the rest to _reducer_lowering; _check_no_unsupported_reducer runs on the raw traced objective (before _scalar can drop dims/init). All #1367 solves route through the f(sol.u,p) ≈ sol.objective guard. One #1367 expected message changed: -abs2(u1) is now a recognized atom, so it is refused by the epigraph-sign check ("nondecreasing in it for MinSense") instead of the DCP gate ("not certified convex") — still a rejection, just from a later stage with a more specific reason. Mixed-family coverage added: sum(abs2.(u-b)) + sum(abs.(u)) → 1.75, maximum(u) + abs2(u1-1) → 0.75, u'P2u + |u1-0.5| - u1 → -0.125, all matching analytic optima.

Verification

Run in lib/ConvexOptimization/ with JULIA_PKG_SERVER="", Julia 1.12.4:

julia --project=. -e 'using Pkg; Pkg.test()'

Test Summary: | Pass  Total     Time
Core          |  529    529  4m45.1s
     Testing ConvexOptimization tests passed

julia --project=@runic -e 'using Runic; exit(Runic.main(ARGS))' -- --check src test

→ RUNIC_CHECK_EXIT=0. typos src test → TYPOS_EXIT=0.

Before/after on the same problems (parent commit 90a74c375):

sum(abs2.(A*u-b)): ERROR: ArgumentError: invalid index: _1 of type SymbolicUtils.BasicSymbolicImpl...
u'Pu: ERROR: Objective is not certified convex for MinSense: curvature = UnknownCurvature.

Both now solve to the analytic optimum within 1e-6.

Scalar-factor sweep (all observed this round)

A = [1 2; 0 1; 1 1], c = [1,-2], M = A'A, P2 = [2 1; 1 2], box [-10,10]². f(u) is the user objective evaluated at the returned point; "consistent" means f(u) == sol.objective.

k=0.5 k*(u'u) - c'u                obj=-2.5        f(u)=-2.5        consistent OK truth=-2.5
k=0.5 (k*u)'*u - c'u               obj=-2.5        f(u)=-2.5        consistent OK truth=-2.5
k=0.5 u'*(k*u) - c'u               obj=-2.5        f(u)=-2.5        consistent OK truth=-2.5
k=0.5 k*(u'P2u) - c'u              obj=-2.3333333  f(u)=-2.3333333  consistent OK truth=-2.3333333
k=0.5 u'*(k*P2)*u - c'u            obj=-2.3333333  f(u)=-2.3333333  consistent OK truth=-2.3333333
k=0.5 k*((Au)'(Au)) - c'u          obj=-4.3333333  f(u)=-4.3333333  consistent OK truth=-4.3333333
k=0.5 (k*(Au))'*(Au) - c'u         obj=-4.3333333  f(u)=-4.3333333  consistent OK truth=-4.3333333
k=0.5 (Au)'*(k*(Au)) - c'u         obj=-4.3333333  f(u)=-4.3333333  consistent OK truth=-4.3333333
k=2.0 k*(u'u) - c'u                obj=-0.625      f(u)=-0.625      consistent OK truth=-0.625
k=2.0 (k*u)'*u - c'u               obj=-0.625      f(u)=-0.625      consistent OK truth=-0.625
k=2.0 u'*(k*u) - c'u               obj=-0.625      f(u)=-0.625      consistent OK truth=-0.625
k=2.0 k*(u'P2u) - c'u              obj=-0.5833333  f(u)=-0.5833333  consistent OK truth=-0.5833333
k=2.0 u'*(k*P2)*u - c'u            obj=-0.5833333  f(u)=-0.5833333  consistent OK truth=-0.5833333
k=2.0 k*((Au)'(Au)) - c'u          obj=-1.0833333  f(u)=-1.0833333  consistent OK truth=-1.0833333
k=2.0 (k*(Au))'*(Au) - c'u         obj=-1.0833333  f(u)=-1.0833333  consistent OK truth=-1.0833333
k=2.0 (Au)'*(k*(Au)) - c'u         obj=-1.0833333  f(u)=-1.0833333  consistent OK truth=-1.0833333
k=-1.0 all 8 positions             REJECTED (PSD / non-convex-objective message)

p[1] scale, crossing zero (p=2 solves to analytic value; p=-2 must refuse):
p1*(u'u) - c'u                     p=[2.0] obj=-0.625 f(u)=-0.625 OK   p=[-2.0] REJECTED (epigraph sign)
p1*(u'P2u) - c'u                   p=[2.0] obj=-0.5833333 OK           p=[-2.0] REJECTED (epigraph sign)
p1*((Au)'(Au)) - c'u               p=[2.0] obj=-1.0833333 OK           p=[-2.0] REJECTED (epigraph sign)
(p1*u)'*u - c'u                    p=[2.0] REJECTED (affine-in-u)       p=[-2.0] REJECTED
u'*(p1*u) - c'u                    p=[2.0] REJECTED (affine-in-u)       p=[-2.0] REJECTED
u'*(p1*P2)*u - c'u                 p=[2.0] REJECTED (constant matrix)  p=[-2.0] REJECTED

Atom-level scales (all consistent): 2*sum(abs2.(u)) - c'u → -0.625, 0.5*sum(abs2.(u)) → -2.5, 2*sum(abs2,u) → -0.625, 2*quad_form(u,P2) → -0.5833, quad_form(u,2P2) → -0.5833, 2*(u1^2+u2^2) → -0.625, sum(abs2.(u))*2 → -0.625, sum(abs2.(u))/2 → -2.5, (u'u)*0.5 → -2.5, (u-b)'(u-b) - c'u → -1.25, 2*((u-b)'(u-b)) - c'u → -0.625, -2*sum(abs2.(u)) → rejected. Reviewer reproducer (A*u)'*(2*A*u) - c'*u → -1.0833 = -c'(M\c)/(4·2) ✓ (was -2.1667 before the fix). (2*A*u)'*(A*u) → -1.0833, u'*(3u) - c'u → -0.4167 = -‖c‖²/12, 2*((A*u)'*(A*u)) → -1.0833. (A*u)'*((2A)*u) (scalar pre-fused into a numeric matrix, a v'*M*w cross term) is refused with UnknownCurvature — sound.

New coverage in test/core_tests.jl: LP scalar products (-3, -11, -11, -2 on the unit box), v'*v self-products vs A\b / analytic optima, parameter-only squares in objectives/constraints/atom args and across reinit! including non-polynomial arguments (exp(p1)^2, abs2(sin(p1))), asymmetric P ([2 1; -1 2] uses only sym(P) = 2I), singular-PSD P, near-PSD eigenvalue clamped with a @debug log (@test_logs at Debug level, deterministic diagonal P), materially indefinite [1e6 0; 0 -1e-6] rejected, scalar-scaled and multi-matrix products (u'*(2u), (2u)'*P*u, u'*A'*A*u, indefinite u'*I*Bi*u, k∈{0.5,2} on both factors and outside, k inside P), trace-shape error messages, reinit! flipping the sign of p on abs2(u1)*p1 rejected while the cache still solves at p = 3, the rejection suite collapsed to a (message, f, u0, p, kwargs) loop with message matches, and a file-wide f(sol.u,p) ≈ sol.objective consistency guard. Flat-optimum primal checks assert the objective tightly and bound u by sqrt(dobj/λmin) instead of a magic tolerance.

Cross-check against Convex.jl 0.16.7 + Clarabel 0.11.1 in a scratch environment (same solver, tight tolerances):

least squares  min sum(abs2.(A*u-b)):  du=0.0  dobj=0.0
unconstrained QP min u'Pu + c'u:       du=8.5e-7   dobj=1.1e-10
equality-constrained QP:               du=2.9e-6   dobj=8.5e-11

The constrained-QP gap is Convex.jl's own solver error: against the independent KKT oracle, ConvexMOI is off by 3.1e-7 while Convex.jl is off by 2.6e-6.

Not verified

• Only the Core test group was run; the docs build and other groups were not exercised (this change adds no docs pages and no new dependencies).
• Convex.jl comparison was run in a scratch env, not committed; the committed tests pin solutions against A \ b and hand-computed KKT systems instead.
• u'*A*B*u inside a larger expression (e.g. - u[1]) still cannot be traced: Symbolics builds a (1,)-shaped term that cannot be added to scalar terms, inside the user's f. It now throws a clear "could not be traced" error naming u' * (A * B) * u and (A*u)'*(B*u) as supported spellings; standalone (1,)-shaped products do lower.
• (v)'*(w) with v ≠ w (e.g. (A*u)'*(B*u) with sym(A'B) PSD) is not decomposed into a quadratic form and is rejected as uncertified — a sound refusal, not a wrong answer. A scalar pre-fused into a matrix ((A*u)'*((2A)*u)) refuses the same way.
• Quadratic constraints remain out of scope — only objective epigraph atoms are lowered.

Reviewer pushback

• The PSD acceptance line is λmin >= -n*eps(Float64)*λmax on sym(P) — check that is the right boundary; anything below is rejected, anything negative inside is clamped with a @debug log.
• v' * v matching compares the factors after adjoint(v) to v structurally (isequal on materialized arrays); check for spellings that miss this and silently fall back to the generic rejection.
• Scalar factors are recognized via Symbolics.value(...) isa Number; a p-dependent factor correctly fails that check, but probe whether any other non-numeric symbolic slips through.
• Check epigraph sign bookkeeping for MaxSense and for atoms scaled by p that change sign between reinit! calls.

Closes #1365.

Please ignore this draft until reviewed by @ChrisRackauckas.

🤖 Generated with Devin — Devin CLI 3000.11.1, model swe-2-high; local session, transcript at /home/crackauc/.local/share/devin/cli/summaries/history_11b25520a0184738.md