From 971e07b81964902f75c73cff5d5d26e382557118 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:50:09 +0100 Subject: [PATCH] move to using sparse matrices in symmetrized sdp dense are faster for small sizes only --- src/sos_sdps.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 8f086c7..bbcba67 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -216,15 +216,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)) + M = spzeros.(T, size.(P)) M_orb = zeros(T, size(parent(elt).mstructure)...) - tmps = SymbolicWedderburn._tmps(wedderburn) end X = convert(Vector{T}, StarAlgebras.coeffs(elt)) @@ -235,16 +234,27 @@ function sos_problem_primal( @info "Adding $(length(invariant_vectors(wedderburn))) constraints" - for iv in invariant_vectors(wedderburn) + ds = SymbolicWedderburn.direct_summands(wedderburn) + Uπs = SymbolicWedderburn.image_basis.(ds) + T = eltype(first(Uπs)) + degrees = SymbolicWedderburn.degree.(ds) + + for (i, iv) in enumerate(invariant_vectors(wedderburn)) + # @debug i + i % 100 == 0 && print('.') + i % 10000 === 0 && print('\n') + x = dot(X, iv) u = dot(U, iv) M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv) + M = SymbolicWedderburn.diagonalize!(M, M_orb, Uπs, degrees) + SparseArrays.droptol!.(M, 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[π])) + # @debug [nnz(m) / length(m) for m in M] + # spM = sparse.(M) + # @time M_dot_P = sum(dot(spM[π], P[π]) for π in eachindex(spM) if !iszero(spM[π])) + # @info density = [count(!iszero, m) / sum(length, m) for m in M] if feasibility_problem JuMP.@constraint(model, x == __fast_recursive_dot!(JuMP.AffExpr(), P, M))