src/constraint_matrix.jl

@@ -83,29 +83,26 @@ Base.@propagate_inbounds function Base.getindex(
     return pos - neg
 end
 
-struct NZPairsIter{T,I}
+struct NZPairsIter{T}
-    m::ConstraintMatrix{T,I}
+    m::ConstraintMatrix{T}
 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(
+function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int} = (1, 1))
-    itr::NZPairsIter,
-    state::Tuple{Int,Nothing} = (1, nothing),
-)
     k = iterate(itr.m.pos, state[1])
-    isnothing(k) && return iterate(itr, (nothing, 1))
+    isnothing(k) && return iterate(itr, state[2])
     idx, st = k
-    return idx => itr.m.val, (st, nothing)
+    return idx => itr.m.val, (st, 1)
 end
 
-function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
+function Base.iterate(itr::NZPairsIter, state::Int)
-    k = iterate(itr.m.neg, state[2])
+    k = iterate(itr.m.neg, state[1])
     isnothing(k) && return nothing
     idx, st = k
-    return idx => -itr.m.val, (nothing, st)
+    return idx => -itr.m.val, st
 end
 
 """

src/reconstruct.jl

@@ -12,10 +12,9 @@ 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)
+        res = _reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom)
         return res
     end
-    res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
     return res
 end
 
@@ -23,14 +22,16 @@ 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)
-        res = (U' * M * U) .* d
+    Θπ = (U' * M * U) .* d
-    end
+
-    return res
+    res .= zero(eltype(res))
+    Θπ = average!(res, Θπ, G, hom)
+    return Θπ
 end
 
 function __droptol!(M::AbstractMatrix, tol)
@@ -51,18 +52,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)
-                if !iszero(M[r, c])
+                res[r^p, c^p] += 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))
-    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
+    for i in SparseArrays.nonzeroinds(invariant_vec)
         g = basis[i]
         A = cnstrs[g]
         for (idx, v) in nzpairs(A)
@@ -109,25 +109,6 @@ 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)
@@ -194,23 +175,8 @@ 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, λ)
@@ -224,52 +190,58 @@ function sos_problem_primal(
         end
     end
 
-    # semidefinite constraints as described by wedderburn
+    P = map(direct_summands(wedderburn)) do ds
-    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 # Ms are preallocated for the constraints loop
+    begin # preallocating
         T = eltype(wedderburn)
-        Ms = [spzeros.(T, size(p)...) for p in Ps]
+        Ms = [spzeros.(T, size(p)...) for p in P]
-        _eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
+        M_orb = zeros(T, size(parent(elt).mstructure)...)
     end
 
-    X = StarAlgebras.coeffs(elt)
+    X = convert(Vector{T}, StarAlgebras.coeffs(elt))
-    U = StarAlgebras.coeffs(orderunit)
+    U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
 
     # defining constraints based on the multiplicative structure
     cnstrs = constraints(parent(elt); augmented = augmented)
 
-    # adding linear constraints: one per orbit
+    prog = ProgressMeter.Progress(
+        length(invariant_vectors(wedderburn));
+        dt = 1,
+        desc = "Adding constraints: ",
+        enabled = show_progress,
+        barlen = 60,
+        showspeed = true,
+    )
+
     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)
 
-        spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv)
+        M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
 
         Ms = SymbolicWedderburn.diagonalize!(
             Ms,
-            spM_orb,
+            M_orb,
             wedderburn;
             trace_preserving = true,
         )
-        for M in Ms
+        # SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
-            SparseArrays.droptol!(M, _eps)
+
-        end
+        # @info [nnz(m) / length(m) for m in Ms]
+
         if feasibility_problem
-            JuMP.@constraint(model, x == _dot(Ps, Ms))
+            JuMP.@constraint(model, x == _dot(P, Ms))
         else
-            JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
+            JuMP.@constraint(model, x - λ * u == _dot(P, Ms))
         end
     end
     ProgressMeter.finish!(prog)
-    return model, Ps
+    return model, P
 end

test/1812.03456.jl

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