Merge pull request #10 from kalmarek/mk/perf_tweaks

Various perf tweaks including moving to sparse matrices where possible
2024-11-25 16:50:27 +01:00 · 2022-12-02 14:21:24 +01:00 · 2022-12-02 14:21:24 +01:00 · 0127d05594
commit 0127d05594
parent 1072cc7074 2cc9444667
25 changed files with 604 additions and 409 deletions
--- a/.codecov.yml
+++ b/.codecov.yml
@ -1 +0,0 @@
-comment: false
--- a/.gitignore
+++ b/.gitignore
@ -2,7 +2,3 @@
 *.jl.*.cov
 *.jl.mem
 Manifest.toml
-test/SL*
-test/oSL*
-test/SAut*
-test/oSAut*
--- a/.travis.yml
+++ b/.travis.yml
@ -1,28 +0,0 @@
-# Documentation: http://docs.travis-ci.com/user/languages/julia/
-language: julia
-os:
-  - linux
-  - osx
-julia:
-  - 1.1
-  - 1.2
-  - 1.3
-  - nightly
-notifications:
-  email: true
-matrix:
-  fast_finish: true
-  allow_failures:
-    - julia: nightly
-    - os: osx
-
-addons:
-  apt:
-    packages:
-    - hdf5-tools
-
-## uncomment the following lines to override the default test
-# script:
-  # - julia -e 'using Pkg; Pkg.build(); Pkg.test(coverage=true);'
-
-codecov: true
--- a/Project.toml
+++ b/Project.toml
@ -8,6 +8,7 @@ Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
@ -17,9 +18,10 @@ COSMO = "0.8"
 Groups = "0.7"
 IntervalArithmetic = "0.20"
 JuMP = "1.3"
+ProgressMeter = "1.7"
 SCS = "1.1.0"
 StaticArrays = "1"
-SymbolicWedderburn = "0.3.1"
+SymbolicWedderburn = "0.3.2"
 julia = "1.6"

 [extras]
--- a/README.md
+++ b/README.md
@ -83,7 +83,10 @@ julia> include("test/optimizers.jl");

 Now we have everything what we need to solve the problem!
 ```julia
-julia> status, warmstart = PropertyT.solve(opt_problem, scs_optimizer(max_iters=5_000, accel=50, alpha=1.9));
+julia> status, warmstart = PropertyT.solve(
+    opt_problem,
+    scs_optimizer(max_iters=5_000, accel=50, alpha=1.9),
+);
 ------------------------------------------------------------------
               SCS v3.2.1 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012
@ -119,8 +122,12 @@ julia> status
 ALMOST_OPTIMAL::TerminationStatusCode = 7
 ```
 The solver didn't manage to solve the problem but it got quite close! (duality gap is ~`1.63e-6`). Let's try once again this time warmstarting the solver:
-```
-julia> status, warmstart = PropertyT.solve(opt_problem, scs_optimizer(max_iters=10_000, accel=50, alpha=1.9), warmstart);
+```julia
+julia> status, warmstart = PropertyT.solve(
+    opt_problem,
+    scs_optimizer(max_iters=10_000, accel=50, alpha=1.9),
+    warmstart,
+);
 ------------------------------------------------------------------
               SCS v3.2.1 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012
--- a/appveyor.yml
+++ b/appveyor.yml
@ -1,34 +0,0 @@
-environment:
-  matrix:
-  - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
-  - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe"
-
-branches:
-  only:
-    - master
-    - /release-.*/
-
-notifications:
-  - provider: Email
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: false
-
-install:
-# Download most recent Julia Windows binary
-  - ps: (new-object net.webclient).DownloadFile(
-        $("http://s3.amazonaws.com/"+$env:JULIAVERSION),
-        "C:\projects\julia-binary.exe")
-# Run installer silently, output to C:\projects\julia
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
-
-build_script:
-# Need to convert from shallow to complete for Pkg.clone to work
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"Property(T)\"); Pkg.build(\"Property(T)\")"
-
-test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"Property(T)\")"
--- a/scripts/SpN_Adj.jl
+++ b/scripts/SpN_Adj.jl
@ -0,0 +1,80 @@
+using LinearAlgebra
+BLAS.set_num_threads(8)
+
+ENV["OMP_NUM_THREADS"] = 1
+
+using Groups
+import Groups.MatrixGroups
+
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
+using PropertyT
+
+using PropertyT.SymbolicWedderburn
+using PropertyT.PermutationGroups
+using PropertyT.StarAlgebras
+
+include(joinpath(@__DIR__, "argparse.jl"))
+include(joinpath(@__DIR__, "utils.jl"))
+
+const N = parsed_args["N"]
+const HALFRADIUS = parsed_args["halfradius"]
+const UPPER_BOUND = parsed_args["upper_bound"]
+
+const GENUS = 2N
+
+G = MatrixGroups.SymplecticGroup{GENUS}(Int8)
+
+RG, S, sizes =
+    @time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
+
+wd = let RG = RG, N = N
+    G = StarAlgebras.object(RG)
+    P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
+    Σ = Groups.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
+    # Σ = P
+    act = PropertyT.action_by_conjugation(G, Σ)
+    @info "Computing WedderburnDecomposition"
+
+    wdfl = @time SymbolicWedderburn.WedderburnDecomposition(
+        Float64,
+        Σ,
+        act,
+        basis(RG),
+        StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
+    )
+end
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+Δs = PropertyT.laplacians(
+    RG,
+    S,
+    x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
+)
+
+# elt = Δ^2
+elt = PropertyT.Adj(Δs, :C₂)
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd,
+    upper_bound=UPPER_BOUND,
+    augmented=true,
+    show_progress=true
+)
+
+solve_in_loop(
+    model,
+    wd,
+    varP,
+    logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ",
+    optimizer=cosmo_optimizer(
+        eps=1e-10,
+        max_iters=20_000,
+        accel=50,
+        alpha=1.95,
+    ),
+    data=(elt=elt, unit=unit, halfradius=HALFRADIUS)
+)
+
--- a/scripts/argparse.jl
+++ b/scripts/argparse.jl
@ -0,0 +1,19 @@
+using ArgParse
+
+args_settings = ArgParseSettings()
+@add_arg_table! args_settings begin
+    "-N"
+    help = "the degree/genus/etc. parameter for a group"
+    arg_type = Int
+    default = 3
+    "--halfradius", "-R"
+    help = "the halfradius on which perform the sum of squares decomposition"
+    arg_type = Int
+    default = 2
+    "--upper_bound", "-u"
+    help = "set upper bound for the optimization problem to speed-up the convergence"
+    arg_type = Float64
+    default = Inf
+end
+
+parsed_args = parse_args(ARGS, args_settings)
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@ -0,0 +1,85 @@
+using Dates
+using Serialization
+using Logging
+
+import JuMP
+
+function get_solution(model)
+    λ = JuMP.value(model[:λ])
+    Q = real.(sqrt(JuMP.value.(model[:P])))
+    solution = Dict(:λ => λ, :Q => Q)
+    return solution
+end
+
+function get_solution(model, wd, varP; logdir)
+    λ = JuMP.value(model[:λ])
+
+    Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
+    Q = PropertyT.reconstruct(Qs, wd)
+    solution = Dict(:λ => λ, :Q => Q)
+    return solution
+end
+
+function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
+    @info "logging to $logdir"
+    status = JuMP.UNKNOWN_RESULT_STATUS
+    warm = try
+        solution = deserialize(joinpath(logdir, "solution.sjl"))
+        warm = solution[:warm]
+        @info "trying to warm-start model with λ=$(solution[:λ])..."
+        warm
+    catch
+        nothing
+    end
+    old_lambda = 0.0
+    while status != JuMP.OPTIMAL
+        date = now()
+        log_file = joinpath(logdir, "solver_$date.log")
+        @info "Current logfile is $log_file."
+        isdir(dirname(log_file)) || mkpath(dirname(log_file))
+
+        λ, flag, certified_λ = let
+            # logstream = current_logger().logger.stream
+            # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint
+            # @warn v
+            status, warm = @time PropertyT.solve(log_file, model, optimizer, warm)
+
+            solution = get_solution(model, args...; logdir=logdir)
+            solution[:warm] = warm
+
+            serialize(joinpath(logdir, "solution_$date.sjl"), solution)
+            serialize(joinpath(logdir, "solution.sjl"), solution)
+
+            flag, λ_cert = open(log_file, append=true) do io
+                with_logger(SimpleLogger(io)) do
+                    PropertyT.certify_solution(
+                        data.elt,
+                        data.unit,
+                        solution[:λ],
+                        solution[:Q],
+                        halfradius=data.halfradius,
+                    )
+                end
+            end
+
+            solution[:λ], flag, λ_cert
+        end
+
+        if flag == true && certified_λ ≥ 0
+            @info "Certification done with λ = $certified_λ"
+            return certified_λ
+        else
+            rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
+            @info "Certification failed with λ = $λ" certified_λ rel_change
+        end
+
+        old_lambda = certified_λ
+
+        if rel_change < 1e-9
+            @info "No progress detected, breaking"
+            break
+        end
+    end
+
+    return status == JuMP.OPTIMAL ? old_lambda : NaN
+end
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@ -1,4 +1,3 @@
-__precompile__()
 module PropertyT

 using LinearAlgebra
@ -15,6 +14,8 @@ import SymbolicWedderburn.PermutationGroups

 include("constraint_matrix.jl")
 include("sos_sdps.jl")
+include("solve.jl")
+include("reconstruct.jl")
 include("certify.jl")

 include("sqadjop.jl")
@ -28,7 +29,7 @@ include("actions/actions.jl")
 function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool)
    S = union!(S, inv.(S))
    @info "generating wl-metric ball of radius $(2halfradius)"
-    @time E, sizes = Groups.wlmetric_ball_serial(S, radius=2halfradius)
+    @time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius)
    @info "sizes = $(sizes)"
    @info "computing the *-algebra structure for G"
    @time RG = StarAlgebras.StarAlgebra{twisted}(
--- a/src/actions/actions.jl
+++ b/src/actions/actions.jl
@ -4,6 +4,7 @@ StarAlgebras.star(g::Groups.GroupElement) = inv(g)
 include("alphabet_permutation.jl")

 include("sln_conjugation.jl")
+include("spn_conjugation.jl")
 include("autfn_conjugation.jl")

 function SymbolicWedderburn.action(
--- a/src/actions/spn_conjugation.jl
+++ b/src/actions/spn_conjugation.jl
@ -1,26 +1,43 @@
 ## Particular definitions for actions on Sp(n,ℤ)

 function _conj(
-    t::MatrixGroups.ElementarySymplectic{N,T},
+    s::MatrixGroups.ElementarySymplectic{N,T},
    σ::PermutationGroups.AbstractPerm,
 ) where {N,T}
    @assert iseven(N)
-    @assert degree(σ) == N ÷ 2 "Got degree = $(degree(σ)); N = $N"
-    i = mod1(t.i, N ÷ 2)
-    ib = i == t.i ? 0 : N ÷ 2
-    j = mod1(t.j, N ÷ 2)
-    jb = j == t.j ? 0 : N ÷ 2
-    return MatrixGroups.ElementarySymplectic{N}(t.symbol, i^inv(σ) + ib, j^inv(σ) + jb, t.val)
+    @assert PermutationGroups.degree(σ) == N ÷ 2 "Got degree = $(PermutationGroups.degree(σ)); N = $N"
+    n = N ÷ 2
+    @assert 1 ≤ s.i ≤ N
+    @assert 1 ≤ s.j ≤ N
+    if s.symbol == :A
+        @assert 1 ≤ s.i ≤ n
+        @assert 1 ≤ s.j ≤ n
+        i = s.i^inv(σ)
+        j = s.j^inv(σ)
+    else
+        @assert s.symbol == :B
+        @assert xor(s.i > n, s.j > n)
+        if s.i > n
+            i = (s.i - n)^inv(σ) + n
+            j = s.j^inv(σ)
+        elseif s.j > n
+            i = s.i^inv(σ)
+            j = (s.j - n)^inv(σ) + n
+        end
+    end
+    return MatrixGroups.ElementarySymplectic{N}(s.symbol, i, j, s.val)
 end

 function _conj(
-    t::MatrixGroups.ElementarySymplectic{N,T},
+    s::MatrixGroups.ElementarySymplectic{N,T},
    x::Groups.Constructions.DirectPowerElement,
 ) where {N,T}
    @assert Groups.order(Int, parent(x).group) == 2
    @assert iseven(N)
-    just_one_flips = xor(isone(x.elts[mod1(t.i, N ÷ 2)]), isone(x.elts[mod1(t.j, N ÷ 2)]))
-    return ifelse(just_one_flips, inv(t), t)
+    n = N ÷ 2
+    i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n)
+    just_one_flips = xor(isone(x.elts[i]), isone(x.elts[j]))
+    return ifelse(just_one_flips, inv(s), s)
 end

 action_by_conjugation(sln::Groups.MatrixGroups.SymplecticGroup, Σ::Groups.Group) =
--- a/src/certify.jl
+++ b/src/certify.jl
@ -5,80 +5,50 @@ function augment_columns!(Q::AbstractMatrix)
    return Q
 end

-function _fma_SOS_thr!(
-    result::AbstractVector{T},
-    mstructure::AbstractMatrix{<:Integer},
-    Q::AbstractMatrix{T},
-    acc_matrix=zeros(T, size(mstructure)...),
-) where {T}
-
-    s1, s2 = size(mstructure)
-
-    @inbounds for k = 1:s2
-        let k = k, s1 = s1, s2 = s2, Q = Q, acc_matrix = acc_matrix
-            Threads.@threads for j = 1:s2
-                for i = 1:s1
-                    @inbounds acc_matrix[i, j] =
-                        muladd(Q[i, k], Q[j, k], acc_matrix[i, j])
-                end
-            end
-        end
-    end
-
-    @inbounds for j = 1:s2
-        for i = 1:s1
-            result[mstructure[i, j]] += acc_matrix[i, j]
-        end
-    end
-
-    return result
-end
-
-function _cnstr_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix, cnstrs)
+function __sos_via_sqr!(
+    res::StarAlgebras.AlgebraElement,
+    P::AbstractMatrix;
+    augmented::Bool
+)
+    StarAlgebras.zero!(res)
+    A = parent(res)
+    b = basis(A)
+    @assert size(A.mstructure) == size(P)
+    e = b[one(b[1])]
+
+    for i in axes(A.mstructure, 1)
+        x = StarAlgebras._istwisted(A.mstructure) ? StarAlgebras.star(b[i]) : b[i]
+        for j in axes(A.mstructure, 2)
+            p = P[i, j]
+            xy = b[A.mstructure[i, j]]
+            # either result += P[x,y]*(x*y)
+            res[xy] += p
+            if augmented
+                # or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy)
+                y = b[j]
+                res[e] += p
+                res[x] -= p
+                res[y] -= p
+            end
+        end
+    end
+
+    return res
+end
+
+function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs)
    StarAlgebras.zero!(res)
-    Q² = Q' * Q
    for (g, A_g) in cnstrs
        res[g] = dot(A_g, Q²)
    end
    return res
 end

-function _augmented_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix)
-    A = parent(res)
-    StarAlgebras.zero!(res)
-    Q² = Q' * Q
-
-    N = LinearAlgebra.checksquare(A.mstructure)
-    augmented_basis = [A(1) - A(b) for b in @view basis(A)[1:N]]
-    tmp = zero(res)
-
-    for (j, y) in enumerate(augmented_basis)
-        for (i, x) in enumerate(augmented_basis)
-            # res += Q²[i, j] * x * y
-
-            StarAlgebras.mul!(tmp, x, y)
-            StarAlgebras.mul!(tmp, tmp, Q²[i, j])
-            StarAlgebras.add!(res, res, tmp)
-        end
-    end
-    return res
-end
-
 function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool)
-    if augmented
-        z = zeros(eltype(Q), length(basis(A)))
-        res = StarAlgebras.AlgebraElement(z, A)
-        return _augmented_sos!(res, Q)
-        cnstrs = constraints(basis(A), A.mstructure; augmented=true)
-        return _cnstr_sos!(res, Q, cnstrs)
-    else
-        @assert size(A.mstructure) == size(Q)
-        z = zeros(eltype(Q), length(basis(A)))
-
-        _fma_SOS_thr!(z, A.mstructure, Q)
-
-        return StarAlgebras.AlgebraElement(z, A)
-    end
+    Q² = Q' * Q
+    res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A)
+    res = __sos_via_sqr!(res, Q², augmented=augmented)
+    return res
 end

 function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
@ -159,7 +129,7 @@ function certify_solution(
        true, compute_sos(parent(elt), Q_int, augmented=augmented)
    end

-    @info "Checking in $(eltype(sos_int)) arithmetic with" λ
+    @info "Checking in $(eltype(sos_int)) arithmetic with" λ_int

    λ_certified =
        sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
--- a/src/gradings.jl
+++ b/src/gradings.jl
@ -6,7 +6,8 @@ Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
 function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
    if s.symbol === :A
        return Roots.𝕖(N ÷ 2, s.i) - Roots.𝕖(N ÷ 2, s.j)
-    else#if s.symbol === :B
+    else
+        @assert s.symbol === :B
        n = N ÷ 2
        i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n)
        return (-1)^(s.i > s.j) * (Roots.𝕖(n, i) + Roots.𝕖(n, j))
--- a/src/reconstruct.jl
+++ b/src/reconstruct.jl
@ -0,0 +1,69 @@
+__outer_dim(wd::WedderburnDecomposition) = size(first(direct_summands(wd)), 2)
+
+function __group_of(wd::WedderburnDecomposition)
+    # this is veeeery hacky... ;)
+    return parent(first(keys(wd.hom.cache)))
+end
+
+function reconstruct(
+    Ms::AbstractVector{<:AbstractMatrix},
+    wbdec::WedderburnDecomposition;
+    atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
+)
+    n = __outer_dim(wbdec)
+    res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
+        res = similar(M, n, n)
+        reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol)
+    end
+    return res
+end
+
+function reconstruct!(
+    res::AbstractMatrix,
+    M::AbstractMatrix,
+    ds::SymbolicWedderburn.DirectSummand,
+    G,
+    hom;
+    atol=eps(real(eltype(ds))) * 10max(size(res)...)
+)
+    res .= zero(eltype(res))
+    U = SymbolicWedderburn.image_basis(ds)
+    d = SymbolicWedderburn.degree(ds)
+    tmp = (U' * M * U) .* d
+
+    res = average!(res, tmp, G, hom)
+    if eltype(res) <: AbstractFloat
+        __droptol!(res, atol) # TODO: is this really necessary?!
+    end
+    return res
+end
+
+function __droptol!(M::AbstractMatrix, tol)
+    for i in eachindex(M)
+        if abs(M[i]) < tol
+            M[i] = zero(M[i])
+        end
+    end
+    return M
+end
+
+# implement average! for other actions when needed
+function average!(
+    res::AbstractMatrix,
+    M::AbstractMatrix,
+    G::Groups.Group,
+    hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations}
+)
+    @assert size(M) == size(res)
+    for g in G
+        gext = SymbolicWedderburn.induce(hom, g)
+        Threads.@threads for c in axes(res, 2)
+            for r in axes(res, 1)
+                res[r, c] += M[r^gext, c^gext]
+            end
+        end
+    end
+    o = Groups.order(Int, G)
+    res ./= o
+    return res
+end
--- a/src/roots.jl
+++ b/src/roots.jl
@ -48,10 +48,8 @@ function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
 end

 function _positive_direction(α::Root{N}) where {N}
-    last = -1 / √2^(N - 1)
-    return Root{N,Float64}(
-        SVector(ntuple(i -> ifelse(i == N, last, (√2)^-i), N)),
-    )
+    v = α.coord + 1 / (N * 100) * rand(N)
+    return Root{N,Float64}(v / norm(v, 2))
 end

 function positive(roots::AbstractVector{<:Root{N}}) where {N}
--- a/src/solve.jl
+++ b/src/solve.jl
@ -0,0 +1,63 @@
+## Low-level solve
+
+setwarmstart!(model::JuMP.Model, ::Nothing) = model
+
+function setwarmstart!(model::JuMP.Model, warmstart)
+    constraint_map = Dict(
+        ct => JuMP.all_constraints(model, ct...) for
+        ct in JuMP.list_of_constraint_types(model)
+    )
+
+    JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
+
+    for (ct, idx) in pairs(constraint_map)
+        JuMP.set_start_value.(idx, warmstart.slack[ct])
+        JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
+    end
+    return model
+end
+
+function getwarmstart(model::JuMP.Model)
+    constraint_map = Dict(
+        ct => JuMP.all_constraints(model, ct...) for
+        ct in JuMP.list_of_constraint_types(model)
+    )
+
+    primal = value.(JuMP.all_variables(model))
+
+    slack = Dict(k => value.(v) for (k, v) in constraint_map)
+    duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map)
+
+    return (primal=primal, dual=duals, slack=slack)
+end
+
+function solve(m::JuMP.Model, optimizer, warmstart=nothing)
+
+    JuMP.set_optimizer(m, optimizer)
+    MOIU.attach_optimizer(m)
+
+    m = setwarmstart!(m, warmstart)
+
+    JuMP.optimize!(m)
+
+    status = JuMP.termination_status(m)
+
+    return status, getwarmstart(m)
+end
+
+function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
+
+    isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
+
+    Base.flush(Base.stdout)
+    Base.Libc.flush_cstdio()
+    status, warmstart = open(solverlog, "a+") do logfile
+        redirect_stdout(logfile) do
+            status, warmstart = solve(m, optimizer, warmstart)
+            Base.Libc.flush_cstdio()
+            status, warmstart
+        end
+    end
+
+    return status, warmstart
+end
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@ -160,13 +160,39 @@ sos_problem_primal(
    kwargs...
 ) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)

+function __fast_recursive_dot!(
+    res::JuMP.AffExpr,
+    Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
+    Ms::AbstractArray{<:AbstractSparseMatrix};
+)
+    @assert length(Ps) == length(Ms)
+
+    for (A, P) in zip(Ms, Ps)
+        iszero(Ms) && continue
+        rows = rowvals(A)
+        vals = nonzeros(A)
+        for cidx in axes(A, 2)
+            for i in nzrange(A, cidx)
+                ridx = rows[i]
+                val = vals[i]
+                JuMP.add_to_expression!(res, P[ridx, cidx], val)
+            end
+        end
+    end
+    return res
+end
+
+import ProgressMeter
+__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
+
 function sos_problem_primal(
    elt::StarAlgebras.AlgebraElement,
    orderunit::StarAlgebras.AlgebraElement,
    wedderburn::WedderburnDecomposition;
    upper_bound=Inf,
    augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
-    check_orthogonality=true
+    check_orthogonality=true,
+    show_progress=false
 )

    @assert parent(elt) === parent(orderunit)
@ -194,15 +220,14 @@ function sos_problem_primal(
    P = map(direct_summands(wedderburn)) do ds
        dim = size(ds, 1)
        P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
-        @constraint(model, P in PSDCone())
+        JuMP.@constraint(model, P in PSDCone())
        P
    end

    begin # preallocating
        T = eltype(wedderburn)
-        M = zeros.(T, size.(P))
+        Ms = [spzeros.(T, size(p)...) for p in P]
        M_orb = zeros(T, size(parent(elt).mstructure)...)
-        tmps = SymbolicWedderburn._tmps(wedderburn)
    end

    X = convert(Vector{T}, StarAlgebras.coeffs(elt))
@ -211,142 +236,39 @@ function sos_problem_primal(
    # defining constraints based on the multiplicative structure
    cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)

-    @info "Adding $(length(invariant_vectors(wedderburn))) constraints"
+    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))

-    for iv in invariant_vectors(wedderburn)
        x = dot(X, iv)
        u = dot(U, iv)

        M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
+        Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn)
+        SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))

-        M = SymbolicWedderburn.diagonalize!(M, M_orb, wedderburn, tmps)
-        SymbolicWedderburn.zerotol!.(M, atol=1e-12)
-
-        M_dot_P = sum(dot(M[π], P[π]) for π in eachindex(M) if !iszero(M[π]))
+        # @info [nnz(m) / length(m) for m in Ms]

        if feasibility_problem
-            JuMP.@constraint(model, x == M_dot_P)
+            JuMP.@constraint(
+                model,
+                x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
+            )
        else
-            JuMP.@constraint(model, x - λ * u == M_dot_P)
+            JuMP.@constraint(
+                model,
+                x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
+            )
        end
    end
+    ProgressMeter.finish!(prog)
    return model, P
 end
-
-function reconstruct(Ps, wd::WedderburnDecomposition)
-    N = size(first(direct_summands(wd)), 2)
-    P = zeros(eltype(wd), N, N)
-    return reconstruct!(P, Ps, wd)
-end
-
-function group_of(wd::WedderburnDecomposition)
-    # this is veeeery hacky... ;)
-    return parent(first(keys(wd.hom.cache)))
-end
-
-# TODO: move to SymbolicWedderburn
-SymbolicWedderburn.action(wd::WedderburnDecomposition) =
-    SymbolicWedderburn.action(wd.hom)
-
-function reconstruct!(
-    res::AbstractMatrix,
-    Ps,
-    wedderburn::WedderburnDecomposition,
-)
-    G = group_of(wedderburn)
-
-    act = SymbolicWedderburn.action(wedderburn)
-
-    @assert act isa SymbolicWedderburn.ByPermutations
-
-    for (π, ds) in pairs(direct_summands(wedderburn))
-        Uπ = SymbolicWedderburn.image_basis(ds)
-
-        # LinearAlgebra.mul!(tmp, Uπ', P[π])
-        # LinearAlgebra.mul!(tmp2, tmp, Uπ)
-        tmp2 = Uπ' * Ps[π] * Uπ
-        if eltype(res) <: AbstractFloat
-            SymbolicWedderburn.zerotol!(tmp2, atol=1e-12)
-        end
-        tmp2 .*= SymbolicWedderburn.degree(ds)
-
-        @assert size(tmp2) == size(res)
-
-        for g in G
-            p = SymbolicWedderburn.induce(wedderburn.hom, g)
-            for c in axes(res, 2)
-                for r in axes(res, 1)
-                    res[r, c] += tmp2[r^p, c^p]
-                end
-            end
-        end
-    end
-    res ./= Groups.order(Int, G)
-
-    return res
-end
-
-##
-# Low-level solve
-
-setwarmstart!(model::JuMP.Model, ::Nothing) = model
-
-function setwarmstart!(model::JuMP.Model, warmstart)
-    constraint_map = Dict(
-        ct => JuMP.all_constraints(model, ct...) for
-        ct in JuMP.list_of_constraint_types(model)
-    )
-
-    JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
-
-    for (ct, idx) in pairs(constraint_map)
-        JuMP.set_start_value.(idx, warmstart.slack[ct])
-        JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
-    end
-    return model
-end
-
-function getwarmstart(model::JuMP.Model)
-    constraint_map = Dict(
-        ct => JuMP.all_constraints(model, ct...) for
-        ct in JuMP.list_of_constraint_types(model)
-    )
-
-    primal = value.(JuMP.all_variables(model))
-
-    slack = Dict(k => value.(v) for (k, v) in constraint_map)
-    duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map)
-
-    return (primal=primal, dual=duals, slack=slack)
-end
-
-function solve(m::JuMP.Model, optimizer, warmstart=nothing)
-
-    JuMP.set_optimizer(m, optimizer)
-    MOIU.attach_optimizer(m)
-
-    m = setwarmstart!(m, warmstart)
-
-    JuMP.optimize!(m)
-    Base.Libc.flush_cstdio()
-
-    status = JuMP.termination_status(m)
-
-    return status, getwarmstart(m)
-end
-
-function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
-
-    isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
-
-    Base.flush(Base.stdout)
-    Base.Libc.flush_cstdio()
-    status, warmstart = open(solverlog, "a+") do logfile
-        redirect_stdout(logfile) do
-            status, warmstart = solve(m, optimizer, warmstart)
-            status, warmstart
-        end
-    end
-
-    return status, warmstart
-end
--- a/test/1703.09680.jl
+++ b/test/1703.09680.jl
@ -1,20 +1,3 @@
-function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
-    @time sos_problem =
-        PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound)
-
-    status, _ = PropertyT.solve(sos_problem, optimizer)
-    P = JuMP.value.(sos_problem[:P])
-    Q = real.(sqrt(P))
-    certified, λ_cert = PropertyT.certify_solution(
-        elt,
-        unit,
-        JuMP.objective_value(sos_problem),
-        Q,
-        halfradius=halfradius,
-    )
-    return status, certified, λ_cert
-end
-
@testset "1703.09680 Examples" begin

    @testset "SL(2,Z)" begin
@ -79,15 +62,15 @@ end
        @test λ > 1

        m = PropertyT.sos_problem_dual(elt, unit)
-        PropertyT.solve(m, scs_optimizer(
-            eps=1e-10,
+        PropertyT.solve(m, cosmo_optimizer(
+            eps=1e-6,
            max_iters=5_000,
            accel=50,
            alpha=1.9,
        ))

        @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
-        @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-3
+        @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2
    end

    @testset "SAut(F₂)" begin
--- a/test/1712.07167.jl
+++ b/test/1712.07167.jl
@ -1,44 +1,3 @@
-function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
-    @assert aug(elt) == aug(unit) == 0
-    @time sos_problem, Ps =
-        PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
-
-    @time status, _ = PropertyT.solve(sos_problem, optimizer)
-
-    Q = let Ps = Ps
-        flPs = [real.(sqrt(JuMP.value.(P))) for P in Ps]
-        PropertyT.reconstruct(flPs, wd)
-    end
-
-    λ = JuMP.value(sos_problem[:λ])
-
-    sos = let RG = parent(elt), Q = Q
-        z = zeros(eltype(Q), length(basis(RG)))
-        res = AlgebraElement(z, RG)
-        cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true)
-        PropertyT._cnstr_sos!(res, Q, cnstrs)
-    end
-
-    residual = elt - λ * unit - sos
-    λ_fl = PropertyT.sufficient_λ(residual, λ, halfradius=2)
-
-    λ_fl < 0 && return status, false, λ_fl
-
-    sos = let RG = parent(elt), Q = [PropertyT.IntervalArithmetic.@interval(q) for q in Q]
-        z = zeros(eltype(Q), length(basis(RG)))
-        res = AlgebraElement(z, RG)
-        cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true)
-        PropertyT._cnstr_sos!(res, Q, cnstrs)
-    end
-
-    λ_int = PropertyT.IntervalArithmetic.@interval(λ)
-
-    residual_int = elt - λ_int * unit - sos
-    λ_int = PropertyT.sufficient_λ(residual_int, λ_int, halfradius=2)
-
-    return status, λ_int > 0, PropertyT.IntervalArithmetic.inf(λ_int)
-end
-
@testset "1712.07167 Examples" begin

    @testset "SAut(F₃)" begin
--- a/test/SOS_correctness.jl
+++ b/test/SOS_correctness.jl
@ -1,59 +0,0 @@
-@testset "Correctness of HPC SOS computation" begin
-
-    function prepare(G_name, λ, S_size)
-        pm = load("$G_name/delta.jld", "pm")
-        P = load("$G_name/$λ/solution.jld", "P")
-        @time Q = real(sqrt(P))
-
-        Δ_coeff = SparseVector(maximum(pm), collect(1:1+S_size), [S_size; ((-1.0) for i in 1:S_size)...])
-
-        Δ²_coeff = GroupRings.GRmul!(spzeros(length(Δ_coeff)), Δ_coeff, Δ_coeff, pm)
-
-        eoi = Δ²_coeff - λ*Δ_coeff
-
-        Q = PropertyT.augIdproj(Q)
-
-        return eoi, pm, Q
-    end
-
-    #########################################################
-    NAME = "SL(3,Z)"
-    eoi, pm, Q = prepare(NAME, 0.1, 3*2*2)
-
-    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
-    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
-
-    @test norm(sos_sqr - sos_hpc, 1) < 3e-12
-    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
-
-    #########################################################
-    NAME = "SL(3,Z)_orbit"
-    eoi, pm, Q = prepare(NAME, 0.27, 3*2*2)
-
-    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
-    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
-
-    @test norm(sos_sqr - sos_hpc, 1) < 5e-12
-    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
-
-    #########################################################
-    NAME = "SL(4,Z)_orbit"
-    eoi, pm, Q = prepare(NAME, 1.3, 4*3*2)
-
-    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
-    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
-
-    @test norm(sos_sqr - sos_hpc, 1) < 2e-10
-    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
-
-    #########################################################
-    NAME = "SAut(F3)_orbit"
-    eoi, pm, Q = prepare(NAME, 0.15, 4*3*2*2)
-
-    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
-    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
-
-    @test norm(sos_sqr - sos_hpc, 1) < 6e-11
-    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
-
-end
--- a/test/check_positivity.jl
+++ b/test/check_positivity.jl
@ -0,0 +1,41 @@
+function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
+    @time sos_problem =
+        PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound)
+
+    status, _ = PropertyT.solve(sos_problem, optimizer)
+    P = JuMP.value.(sos_problem[:P])
+    Q = real.(sqrt(P))
+    certified, λ_cert = PropertyT.certify_solution(
+        elt,
+        unit,
+        JuMP.objective_value(sos_problem),
+        Q,
+        halfradius=halfradius,
+    )
+    return status, certified, λ_cert
+end
+
+function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
+    @assert aug(elt) == aug(unit) == 0
+    @time sos_problem, Ps =
+        PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
+
+    @time status, _ = PropertyT.solve(sos_problem, optimizer)
+
+    Q = let Ps = Ps
+        Qs = [real.(sqrt(JuMP.value.(P))) for P in Ps]
+        PropertyT.reconstruct(Qs, wd)
+    end
+
+    λ = JuMP.value(sos_problem[:λ])
+
+    certified, λ_cert = PropertyT.certify_solution(
+        elt,
+        unit,
+        λ,
+        Q,
+        halfradius=halfradius
+    )
+    return status, certified, λ_cert
+end
+
--- a/test/graded_adj.jl
+++ b/test/graded_adj.jl
@ -46,6 +46,28 @@


    @testset "Symplectic group" begin
+        @testset "Sp2(ℤ)" begin
+            genus = 2
+            halfradius = 1
+
+            SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
+
+            RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
+
+            Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
+                Δ = RG(length(S)) - sum(RG(s) for s in S)
+                Δs = PropertyT.laplacians(
+                    RG,
+                    S,
+                    x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])),
+                )
+                Δ, Δs
+            end
+
+            sq = sum(Δᵢ^2 for Δᵢ in values(Δs))
+            @test PropertyT.Adj(Δs, :C₂) + sq == Δ^2
+        end
+
        genus = 3
        halfradius = 1

@ -63,7 +85,7 @@
            Δ, Δs
        end

-        @testset "Adj correctness: genus=$genus" begin
+        @testset "Adj numerics for genus=$genus" begin

            all_subtypes = (
                :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
--- a/test/quick_tests.jl
+++ b/test/quick_tests.jl
@ -0,0 +1,80 @@
+@testset "Quick tests" begin
+
+    @testset "SL(2,F₇)" begin
+        N = 2
+        p = 7
+        halfradius = 3
+        G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p})
+        RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true)
+
+        Δ = let RG = RG, S = S
+            RG(length(S)) - sum(RG(s) for s in S)
+        end
+
+        elt = Δ^2
+        unit = Δ
+        ub = 0.58578# Inf# 1.5
+
+        @testset "standard formulation" begin
+            status, certified, λ_cert = check_positivity(
+                elt,
+                unit,
+                upper_bound=ub,
+                halfradius=2,
+                optimizer=cosmo_optimizer(
+                    eps=1e-7,
+                    max_iters=5_000,
+                    accel=50,
+                    alpha=1.95,
+                )
+            )
+
+            @test status == JuMP.OPTIMAL
+            @test certified
+            @test λ_cert > 5857 // 10000
+
+            m = PropertyT.sos_problem_dual(elt, unit)
+            PropertyT.solve(m, cosmo_optimizer(
+                eps=1e-7,
+                max_iters=10_000,
+                accel=50,
+                alpha=1.95,
+            ))
+
+            @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
+            @test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2
+        end
+
+        @testset "Wedderburn decomposition" begin
+            P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
+            Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
+            act = PropertyT.action_by_conjugation(G, Σ)
+
+            wd = WedderburnDecomposition(
+                Float64,
+                Σ,
+                act,
+                basis(RG),
+                StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]),
+            )
+
+            status, certified, λ_cert = check_positivity(
+                elt,
+                unit,
+                wd,
+                upper_bound=ub,
+                halfradius=2,
+                optimizer=cosmo_optimizer(
+                    eps=1e-7,
+                    max_iters=10_000,
+                    accel=50,
+                    alpha=1.9,
+                ),
+            )
+
+            @test status == JuMP.OPTIMAL
+            @test certified
+            @test λ_cert > 5857 // 10000
+        end
+    end
+end
--- a/test/runtests.jl
+++ b/test/runtests.jl
@ -1,8 +1,6 @@
 using Test
 using LinearAlgebra
 using SparseArrays
-BLAS.set_num_threads(1)
-ENV["OMP_NUM_THREADS"] = 4

 using Groups
 using Groups.GroupsCore
@ -14,9 +12,11 @@ using SymbolicWedderburn.StarAlgebras
 using SymbolicWedderburn.PermutationGroups

 include("optimizers.jl")
+include("check_positivity.jl")
+include("quick_tests.jl")

+if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
    @testset "PropertyT" begin
-
        include("constratint_matrices.jl")
        include("actions.jl")

@ -26,3 +26,4 @@ include("optimizers.jl")

        include("graded_adj.jl")
    end
+end