Compare commits

..

No commits in common. "005ffc29cb8278648d85015b46743695c177ed03" and "f9f852439f8a113e429336a0e73f0ae15824ed1b" have entirely different histories.

4 changed files with 48 additions and 78 deletions

View File

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

View File

@ -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,
)
U = SymbolicWedderburn.image_basis(ds)
d = SymbolicWedderburn.degree(ds)
Θπ = (U' * M * U) .* d
res .= zero(eltype(res))
if !iszero(M)
U = SymbolicWedderburn.image_basis(ds)
d = SymbolicWedderburn.degree(ds)
res = (U' * M * U) .* d
end
return 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] / o
end
res[r^p, c^p] += M[r, c]
end
end
end
o = Groups.order(Int, G)
res ./= o
return res
end

View File

@ -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
Ps = map(direct_summands(wedderburn)) do ds
P = 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]
_eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
Ms = [spzeros.(T, size(p)...) for p in P]
M_orb = zeros(T, size(parent(elt).mstructure)...)
end
X = StarAlgebras.coeffs(elt)
U = StarAlgebras.coeffs(orderunit)
X = convert(Vector{T}, StarAlgebras.coeffs(elt))
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!(M, _eps)
end
# SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
# @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

View File

@ -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