2022-11-07 15:42:10 +01:00
|
|
|
"""
|
|
|
|
sos_problem_dual(X, [u = zero(X); upper_bound=Inf])
|
|
|
|
Formulate the dual to the sum of squares decomposition problem for `X - λ·u`.
|
|
|
|
|
|
|
|
See also [sos_problem_primal](@ref).
|
|
|
|
"""
|
|
|
|
function sos_problem_dual(
|
2022-11-07 17:01:06 +01:00
|
|
|
elt::StarAlgebras.AlgebraElement,
|
|
|
|
order_unit::StarAlgebras.AlgebraElement=zero(elt);
|
2022-11-07 15:42:10 +01:00
|
|
|
lower_bound=-Inf
|
|
|
|
)
|
|
|
|
@assert parent(elt) == parent(order_unit)
|
|
|
|
algebra = parent(elt)
|
|
|
|
mstructure = if StarAlgebras._istwisted(algebra.mstructure)
|
|
|
|
algebra.mstructure
|
|
|
|
else
|
|
|
|
StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
|
|
|
|
end
|
|
|
|
|
|
|
|
# 1 variable for every primal constraint
|
|
|
|
# 1 dual variable for every element of basis
|
|
|
|
# Symmetrized:
|
|
|
|
# 1 dual variable for every orbit of G acting on basis
|
|
|
|
model = Model()
|
|
|
|
@variable model y[1:length(basis(algebra))]
|
|
|
|
@constraint model λ_dual dot(order_unit, y) == 1
|
|
|
|
@constraint(model, psd, y[mstructure] in PSDCone())
|
|
|
|
|
|
|
|
if !isinf(lower_bound)
|
|
|
|
throw("Not Implemented yet")
|
|
|
|
@variable model λ_ub_dual
|
|
|
|
@objective model Min dot(elt, y) + lower_bound * λ_ub_dual
|
|
|
|
else
|
|
|
|
@objective model Min dot(elt, y)
|
|
|
|
end
|
|
|
|
|
|
|
|
return model
|
|
|
|
end
|
|
|
|
|
|
|
|
function constraints(
|
|
|
|
basis::StarAlgebras.AbstractBasis,
|
|
|
|
mstr::AbstractMatrix{<:Integer};
|
|
|
|
augmented::Bool=false,
|
|
|
|
table_size=size(mstr)
|
|
|
|
)
|
|
|
|
cnstrs = [signed(eltype(mstr))[] for _ in basis]
|
|
|
|
LI = LinearIndices(table_size)
|
|
|
|
|
|
|
|
for ci in CartesianIndices(table_size)
|
|
|
|
k = LI[ci]
|
|
|
|
a_star_b = basis[mstr[k]]
|
|
|
|
push!(cnstrs[basis[a_star_b]], k)
|
|
|
|
if augmented
|
|
|
|
# (1-a_star)(1-b) = 1 - a_star - b + a_star_b
|
|
|
|
|
|
|
|
i, j = Tuple(ci)
|
|
|
|
a, b = basis[i], basis[j]
|
|
|
|
|
|
|
|
push!(cnstrs[basis[one(a)]], k)
|
2022-11-07 17:01:06 +01:00
|
|
|
push!(cnstrs[basis[StarAlgebras.star(a)]], -k)
|
2022-11-07 15:42:10 +01:00
|
|
|
push!(cnstrs[basis[b]], -k)
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
|
|
|
end
|
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
return Dict(
|
|
|
|
basis[i] => ConstraintMatrix(c, table_size..., 1) for (i, c) in pairs(cnstrs)
|
|
|
|
)
|
|
|
|
end
|
2018-09-09 11:37:11 +02:00
|
|
|
|
2022-11-07 17:01:06 +01:00
|
|
|
function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool, twisted::Bool)
|
2022-11-07 15:42:10 +01:00
|
|
|
mstructure = if StarAlgebras._istwisted(A.mstructure) == twisted
|
|
|
|
A.mstructure
|
|
|
|
else
|
|
|
|
StarAlgebras.MTable{twisted}(basis(A), table_size=size(A.mstructure))
|
2018-09-09 11:37:11 +02:00
|
|
|
end
|
2022-11-07 15:42:10 +01:00
|
|
|
return constraints(basis(A), mstructure, augmented=augmented)
|
2018-09-09 11:37:11 +02:00
|
|
|
end
|
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
"""
|
|
|
|
sos_problem_primal(X, [u = zero(X); upper_bound=Inf])
|
|
|
|
Formulate sum of squares decomposition problem for `X - λ·u`.
|
|
|
|
|
|
|
|
Given element `X` of a star-algebra `A` and an order unit `u` of `Σ²A` the sum
|
|
|
|
of squares cone in `A`, fomulate sum of squares decomposition problem:
|
|
|
|
|
|
|
|
```
|
|
|
|
max: λ
|
|
|
|
subject to: X - λ·u ∈ Σ²A
|
|
|
|
```
|
|
|
|
|
|
|
|
If `upper_bound` is given a finite value, additionally `λ ≤ upper_bound` will
|
|
|
|
be added to the model. This may improve the accuracy of the solution if
|
|
|
|
`upper_bound` is less than the optimal `λ`.
|
|
|
|
|
|
|
|
The default `u = zero(X)` formulates a simple feasibility problem.
|
|
|
|
"""
|
|
|
|
function sos_problem_primal(
|
2022-11-07 17:01:06 +01:00
|
|
|
elt::StarAlgebras.AlgebraElement,
|
|
|
|
order_unit::StarAlgebras.AlgebraElement=zero(elt);
|
2022-11-07 15:42:10 +01:00
|
|
|
upper_bound=Inf,
|
2022-11-07 17:01:06 +01:00
|
|
|
augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit))
|
2022-11-07 15:42:10 +01:00
|
|
|
)
|
|
|
|
@assert parent(elt) === parent(order_unit)
|
|
|
|
|
|
|
|
N = LinearAlgebra.checksquare(parent(elt).mstructure)
|
|
|
|
model = JuMP.Model()
|
|
|
|
P = JuMP.@variable(model, P[1:N, 1:N], Symmetric)
|
|
|
|
JuMP.@constraint(model, psd, P in PSDCone())
|
|
|
|
|
|
|
|
if iszero(order_unit) && isfinite(upper_bound)
|
|
|
|
@warn "Setting `upper_bound` together with zero `order_unit` has no effect"
|
|
|
|
end
|
2017-03-13 14:49:55 +01:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
A = constraints(parent(elt), augmented=augmented, twisted=true)
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
if !iszero(order_unit)
|
|
|
|
λ = JuMP.@variable(model, λ)
|
|
|
|
if isfinite(upper_bound)
|
|
|
|
JuMP.@constraint model λ <= upper_bound
|
|
|
|
end
|
|
|
|
JuMP.@objective(model, Max, λ)
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
for b in basis(parent(elt))
|
|
|
|
JuMP.@constraint(model, elt(b) - λ * order_unit(b) == dot(A[b], P))
|
|
|
|
end
|
|
|
|
else
|
|
|
|
for b in basis(parent(elt))
|
|
|
|
JuMP.@constraint(model, elt(b) == dot(A[b], P))
|
|
|
|
end
|
|
|
|
end
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
return model
|
|
|
|
end
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
function invariant_constraint!(
|
|
|
|
result::AbstractMatrix,
|
|
|
|
basis::StarAlgebras.AbstractBasis,
|
|
|
|
cnstrs::AbstractDict{K,CM},
|
|
|
|
invariant_vec::SparseVector,
|
|
|
|
) where {K,CM<:ConstraintMatrix}
|
|
|
|
result .= zero(eltype(result))
|
|
|
|
for i in SparseArrays.nonzeroinds(invariant_vec)
|
|
|
|
g = basis[i]
|
|
|
|
A = cnstrs[g]
|
|
|
|
for (idx, v) in nzpairs(A)
|
|
|
|
result[idx] += invariant_vec[i] * v
|
|
|
|
end
|
2020-08-27 10:06:02 +02:00
|
|
|
end
|
2022-11-07 15:42:10 +01:00
|
|
|
return result
|
|
|
|
end
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
|
|
|
|
U = SymbolicWedderburn.image_basis(ds)
|
|
|
|
return isapprox(U * U', I)
|
2020-08-27 10:06:02 +02:00
|
|
|
end
|
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
sos_problem_primal(
|
2022-11-07 17:01:06 +01:00
|
|
|
elt::StarAlgebras.AlgebraElement,
|
2022-11-07 15:42:10 +01:00
|
|
|
wedderburn::WedderburnDecomposition;
|
|
|
|
kwargs...
|
|
|
|
) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
|
|
|
|
|
2022-11-14 19:47:38 +01:00
|
|
|
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
|
|
|
|
|
2022-11-15 00:27:09 +01:00
|
|
|
import ProgressMeter
|
|
|
|
__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
|
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
function sos_problem_primal(
|
2022-11-07 17:01:06 +01:00
|
|
|
elt::StarAlgebras.AlgebraElement,
|
|
|
|
orderunit::StarAlgebras.AlgebraElement,
|
2022-11-07 15:42:10 +01:00
|
|
|
wedderburn::WedderburnDecomposition;
|
|
|
|
upper_bound=Inf,
|
2022-11-07 17:01:06 +01:00
|
|
|
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
|
2022-11-15 00:27:09 +01:00
|
|
|
check_orthogonality=true,
|
|
|
|
show_progress=false
|
2022-11-07 15:42:10 +01:00
|
|
|
)
|
|
|
|
|
|
|
|
@assert parent(elt) === parent(orderunit)
|
2022-11-07 17:01:06 +01:00
|
|
|
if check_orthogonality
|
|
|
|
if any(!isorth_projection, direct_summands(wedderburn))
|
|
|
|
error("Wedderburn decomposition contains a non-orthogonal projection")
|
|
|
|
end
|
2022-11-07 15:42:10 +01:00
|
|
|
end
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
feasibility_problem = iszero(orderunit)
|
2018-08-20 03:46:44 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
model = JuMP.Model()
|
|
|
|
if !feasibility_problem # add λ or not?
|
|
|
|
λ = JuMP.@variable(model, λ)
|
|
|
|
JuMP.@objective(model, Max, λ)
|
2020-08-27 10:06:02 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
if isfinite(upper_bound)
|
|
|
|
JuMP.@constraint(model, λ <= upper_bound)
|
|
|
|
if feasibility_problem
|
|
|
|
@warn "setting `upper_bound` with zero `orderunit` has no effect"
|
|
|
|
end
|
|
|
|
end
|
2020-01-10 16:55:13 +01:00
|
|
|
end
|
2017-03-20 21:41:12 +01:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
P = map(direct_summands(wedderburn)) do ds
|
|
|
|
dim = size(ds, 1)
|
|
|
|
P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
|
2022-11-14 19:50:09 +01:00
|
|
|
JuMP.@constraint(model, P in PSDCone())
|
2022-11-07 15:42:10 +01:00
|
|
|
P
|
2017-03-13 14:49:55 +01:00
|
|
|
end
|
2019-04-12 23:18:48 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
begin # preallocating
|
|
|
|
T = eltype(wedderburn)
|
2022-11-15 00:27:09 +01:00
|
|
|
Ms = spzeros.(T, size.(P))
|
2022-11-07 15:42:10 +01:00
|
|
|
M_orb = zeros(T, size(parent(elt).mstructure)...)
|
|
|
|
end
|
2019-08-02 11:52:19 +02:00
|
|
|
|
2022-11-07 17:01:06 +01:00
|
|
|
X = convert(Vector{T}, StarAlgebras.coeffs(elt))
|
|
|
|
U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
|
2019-04-12 23:18:48 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
# defining constraints based on the multiplicative structure
|
|
|
|
cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)
|
2019-04-12 23:18:48 +02:00
|
|
|
|
2022-11-15 00:27:09 +01:00
|
|
|
prog = ProgressMeter.Progress(
|
|
|
|
length(invariant_vectors(wedderburn)),
|
|
|
|
dt=1,
|
|
|
|
desc="Adding constraints... ",
|
|
|
|
enabled=show_progress,
|
|
|
|
barlen=60,
|
|
|
|
showspeed=true
|
|
|
|
)
|
2022-11-14 19:50:09 +01:00
|
|
|
|
|
|
|
for (i, iv) in enumerate(invariant_vectors(wedderburn))
|
2022-11-15 00:27:09 +01:00
|
|
|
ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n))
|
2022-11-14 19:50:09 +01:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
x = dot(X, iv)
|
|
|
|
u = dot(U, iv)
|
2018-09-05 09:18:38 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
|
2022-11-15 00:27:09 +01:00
|
|
|
Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn)
|
|
|
|
SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
|
2018-09-05 09:18:38 +02:00
|
|
|
|
2022-11-15 00:27:09 +01:00
|
|
|
# @info [nnz(m) / length(m) for m in Ms]
|
2018-09-05 09:18:38 +02:00
|
|
|
|
2022-11-07 15:42:10 +01:00
|
|
|
if feasibility_problem
|
2022-11-15 00:27:09 +01:00
|
|
|
JuMP.@constraint(
|
|
|
|
model,
|
|
|
|
x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
|
|
|
|
)
|
2022-11-07 15:42:10 +01:00
|
|
|
else
|
2022-11-15 00:27:09 +01:00
|
|
|
JuMP.@constraint(
|
|
|
|
model,
|
|
|
|
x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
|
|
|
|
)
|
2022-11-07 15:42:10 +01:00
|
|
|
end
|
2018-09-05 09:18:38 +02:00
|
|
|
end
|
2022-11-15 00:27:09 +01:00
|
|
|
ProgressMeter.finish!(prog)
|
2022-11-07 15:42:10 +01:00
|
|
|
return model, P
|
|
|
|
end
|