cosmo fails at high precision with zero condition

src/constraint_matrix.jl

@@ -83,26 +83,29 @@ Base.@propagate_inbounds function Base.getindex(
     return pos - neg
 end
 
-struct NZPairsIter{T}
-    m::ConstraintMatrix{T}
+struct NZPairsIter{T,I}
+    m::ConstraintMatrix{T,I}
 end
 
 Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
 Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
 
 # TODO: iterate over (idx=>val) pairs combining vals
-function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int} = (1, 1))
+function Base.iterate(
+    itr::NZPairsIter,
+    state::Tuple{Int,Nothing} = (1, nothing),
+)
     k = iterate(itr.m.pos, state[1])
-    isnothing(k) && return iterate(itr, state[2])
+    isnothing(k) && return iterate(itr, (nothing, 1))
     idx, st = k
-    return idx => itr.m.val, (st, 1)
+    return idx => itr.m.val, (st, nothing)
 end
 
-function Base.iterate(itr::NZPairsIter, state::Int)
-    k = iterate(itr.m.neg, state[1])
+function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
+    k = iterate(itr.m.neg, state[2])
     isnothing(k) && return nothing
     idx, st = k
-    return idx => -itr.m.val, st
+    return idx => -itr.m.val, (nothing, st)
 end
 
 """

src/reconstruct.jl

@@ -12,9 +12,10 @@ function reconstruct(
     n = __outer_dim(wbdec)
     res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
         res = similar(M, n, n)
-        res = _reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom)
+        res = _reconstruct!(res, M, ds)
         return res
     end
+    res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
     return res
 end
 
@@ -22,16 +23,14 @@ function _reconstruct!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     ds::SymbolicWedderburn.DirectSummand,
-    G,
-    hom,
 )
+    res .= zero(eltype(res))
+    if !iszero(M)
         U = SymbolicWedderburn.image_basis(ds)
         d = SymbolicWedderburn.degree(ds)
-    Θπ = (U' * M * U) .* d
-
-    res .= zero(eltype(res))
-    Θπ = average!(res, Θπ, G, hom)
-    return Θπ
+        res = (U' * M * U) .* d
+    end
+    return res
 end
 
 function __droptol!(M::AbstractMatrix, tol)
@@ -52,18 +51,18 @@ function average!(
         <:SymbolicWedderburn.ByPermutations,
     },
 )
+    res .= zero(eltype(res))
     @assert size(M) == size(res)
+    o = Groups.order(Int, G)
     for g in G
         p = SymbolicWedderburn.induce(hom, g)
         Threads.@threads for c in axes(res, 2)
-            # note: p is a permutation,
-            # so accesses below are guaranteed to be disjoint
             for r in axes(res, 1)
-                res[r^p, c^p] += M[r, c]
+                if !iszero(M[r, c])
+                    res[r^p, c^p] += M[r, c] / o
+                end
             end
         end
     end
-    o = Groups.order(Int, G)
-    res ./= o
     return res
 end

src/sos_sdps.jl

@@ -99,7 +99,7 @@ function invariant_constraint!(
     invariant_vec::SparseVector,
 ) where {K}
     result .= zero(eltype(result))
-    for i in SparseArrays.nonzeroinds(invariant_vec)
+    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
         g = basis[i]
         A = cnstrs[g]
         for (idx, v) in nzpairs(A)
@@ -109,6 +109,25 @@ function invariant_constraint!(
     return result
 end
 
+function invariant_constraint(basis, cnstrs, invariant_vec)
+    I = UInt32[]
+    J = UInt32[]
+    V = Float64[]
+    _M = first(values(cnstrs))
+    CI = CartesianIndices(_M)
+    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
+        g = basis[i]
+        A = cnstrs[g]
+        for (idx, v) in nzpairs(A)
+            ci = CI[idx]
+            push!(I, ci[1])
+            push!(J, ci[2])
+            push!(V, invariant_vec[i] * v)
+        end
+    end
+    return sparse(I, J, V, size(_M)...)
+end
+
 function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
     U = SymbolicWedderburn.image_basis(ds)
     return isapprox(U * U', I)
@@ -175,8 +194,23 @@ function sos_problem_primal(
         end
     end
 
+    id_one = findfirst(invariant_vectors(wedderburn)) do v
+        b = basis(parent(elt))
+        return sparsevec([b[one(first(b))]], [1 // 1], length(v)) == v
+    end
+
+    prog = ProgressMeter.Progress(
+        length(invariant_vectors(wedderburn));
+        dt = 1,
+        desc = "Adding constraints: ",
+        enabled = show_progress,
+        barlen = 60,
+        showspeed = true,
+    )
+
     feasibility_problem = iszero(orderunit)
 
+    # problem creation starts here
     model = JuMP.Model()
     if !feasibility_problem # add λ or not?
         λ = JuMP.@variable(model, λ)
@@ -190,58 +224,52 @@ function sos_problem_primal(
         end
     end
 
-    P = map(direct_summands(wedderburn)) do ds
+    # semidefinite constraints as described by wedderburn
+    Ps = map(direct_summands(wedderburn)) do ds
         dim = size(ds, 1)
         P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
         JuMP.@constraint(model, P in PSDCone())
         return P
     end
 
-    begin # preallocating
+    begin # Ms are preallocated for the constraints loop
         T = eltype(wedderburn)
-        Ms = [spzeros.(T, size(p)...) for p in P]
-        M_orb = zeros(T, size(parent(elt).mstructure)...)
+        Ms = [spzeros.(T, size(p)...) for p in Ps]
+        _eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
     end
 
-    X = convert(Vector{T}, StarAlgebras.coeffs(elt))
-    U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
+    X = StarAlgebras.coeffs(elt)
+    U = StarAlgebras.coeffs(orderunit)
 
     # defining constraints based on the multiplicative structure
     cnstrs = constraints(parent(elt); augmented = augmented)
 
-    prog = ProgressMeter.Progress(
-        length(invariant_vectors(wedderburn));
-        dt = 1,
-        desc = "Adding constraints: ",
-        enabled = show_progress,
-        barlen = 60,
-        showspeed = true,
-    )
-
+    # adding linear constraints: one per orbit
     for (i, iv) in enumerate(invariant_vectors(wedderburn))
         ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n))
+        augmented && i == id_one && continue
+        # i == 500 && break
 
         x = dot(X, iv)
         u = dot(U, iv)
 
-        M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
+        spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv)
 
         Ms = SymbolicWedderburn.diagonalize!(
             Ms,
-            M_orb,
+            spM_orb,
             wedderburn;
             trace_preserving = true,
         )
-        # SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
-
-        # @info [nnz(m) / length(m) for m in Ms]
-
+        for M in Ms
+            SparseArrays.droptol!(M, _eps)
+        end
         if feasibility_problem
-            JuMP.@constraint(model, x == _dot(P, Ms))
+            JuMP.@constraint(model, x == _dot(Ps, Ms))
         else
-            JuMP.@constraint(model, x - λ * u == _dot(P, Ms))
+            JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
         end
     end
     ProgressMeter.finish!(prog)
-    return model, P
+    return model, Ps
 end

test/1812.03456.jl

@@ -176,7 +176,7 @@ end
                 wd;
                 upper_bound = UB,
                 halfradius = 2,
-                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
+                optimizer = scs_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test !certified