kalmar/PropertyT.jl
1
0
Fork 0
mirror of https://github.com/kalmarek/PropertyT.jl.git synced 2024-12-27 02:45:30 +01:00
Code Issues Releases Wiki Activity

use sparse matrices for invariant constraint

Browse Source
This commit is contained in:
Marek Kaluba 2023-04-04 23:14:56 +02:00
parent 1a43a1b1be
commit 132802feeb
No known key found for this signature in database
GPG Key ID: 8BF1A3855328FC15
1 changed files with 30 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
src/sos_sdps.jl
View File

@ -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)
@ -190,7 +209,8 @@ 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())
@ -224,24 +244,23 @@ function sos_problem_primal(
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