replace checksoltuion by certify
This commit is contained in:
parent
ecea3dfbcb
commit
085f6bce3c
|
|
@ -5,6 +5,7 @@ using LinearAlgebra
|
|||
using SparseArrays
|
||||
using Dates
|
||||
|
||||
using IntervalArithmetic
|
||||
using JuMP
|
||||
|
||||
using Groups
|
||||
|
|
@ -14,7 +15,7 @@ using SymbolicWedderburn
|
|||
include("laplacians.jl")
|
||||
include("constraint_matrix.jl")
|
||||
include("sos_sdps.jl")
|
||||
include("checksolution.jl")
|
||||
include("certify.jl")
|
||||
|
||||
include("1712.07167.jl")
|
||||
include("1812.03456.jl")
|
||||
|
|
|
|||
|
|
@ -0,0 +1,168 @@
|
|||
function augment_columns!(Q::AbstractMatrix)
|
||||
for c in eachcol(Q)
|
||||
c .-= sum(c) ./ length(c)
|
||||
end
|
||||
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::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::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::StarAlgebra, Q::AbstractMatrix; augmented::Bool)
|
||||
if augmented
|
||||
z = zeros(eltype(Q), length(basis(A)))
|
||||
res = 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 AlgebraElement(z, A)
|
||||
end
|
||||
end
|
||||
|
||||
function sufficient_λ(residual::AlgebraElement, λ; halfradius)
|
||||
L1_norm = norm(residual, 1)
|
||||
suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
|
||||
|
||||
eq_sign = let T = eltype(residual)
|
||||
if T <: Interval
|
||||
"∈"
|
||||
elseif T <: Union{Rational,Integer}
|
||||
"="
|
||||
else # if T <: AbstractFloat
|
||||
"≈"
|
||||
end
|
||||
end
|
||||
|
||||
info_strs = [
|
||||
"Numerical metrics of the obtained SOS:",
|
||||
"ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(aug(residual))",
|
||||
"‖elt - λu - ∑ξᵢ*ξᵢ‖₁ $eq_sign $(L1_norm)",
|
||||
" λ $eq_sign $suff_λ",
|
||||
]
|
||||
@info join(info_strs, "\n")
|
||||
|
||||
return suff_λ
|
||||
end
|
||||
|
||||
function sufficient_λ(
|
||||
elt::AlgebraElement,
|
||||
order_unit::AlgebraElement,
|
||||
λ,
|
||||
sos::AlgebraElement;
|
||||
halfradius
|
||||
)
|
||||
|
||||
@assert parent(elt) === parent(order_unit) == parent(sos)
|
||||
residual = (elt - λ * order_unit) - sos
|
||||
|
||||
return sufficient_λ(residual, λ; halfradius=halfradius)
|
||||
end
|
||||
|
||||
function certify_solution(
|
||||
elt::AlgebraElement,
|
||||
orderunit::AlgebraElement,
|
||||
λ,
|
||||
Q::AbstractMatrix{<:AbstractFloat};
|
||||
halfradius,
|
||||
augmented=iszero(aug(elt)) && iszero(aug(orderunit))
|
||||
)
|
||||
|
||||
should_we_augment = !augmented && aug(elt) == aug(orderunit) == 0
|
||||
|
||||
Q = should_we_augment ? augment_columns!(Q) : Q
|
||||
@time sos = compute_sos(parent(elt), Q, augmented=augmented)
|
||||
|
||||
@info "Checking in $(eltype(sos)) arithmetic with" λ
|
||||
|
||||
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
|
||||
|
||||
if λ_flpoint ≤ 0
|
||||
return false, λ_flpoint
|
||||
end
|
||||
|
||||
λ_int = @interval(λ)
|
||||
Q_int = [@interval(q) for q in Q]
|
||||
|
||||
check, sos_int = @time if should_we_augment
|
||||
@info("Projecting columns of Q to the augmentation ideal...")
|
||||
Q_int = augment_columns!(Q_int)
|
||||
@info "Checking that sum of every column contains 0.0..."
|
||||
check_augmented = all(0 ∈ sum(c) for c in eachcol(Q_int))
|
||||
check_augmented || @error(
|
||||
"Augmentation failed! The following numbers are not certified!"
|
||||
)
|
||||
sos_int = compute_sos(parent(elt), Q_int; augmented=augmented)
|
||||
check_augmented, sos_int
|
||||
else
|
||||
true, compute_sos(parent(elt), Q_int, augmented=augmented)
|
||||
end
|
||||
|
||||
@info "Checking in $(eltype(sos_int)) arithmetic with" λ
|
||||
|
||||
λ_certified =
|
||||
sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
|
||||
|
||||
return check && inf(λ_certified) > 0.0, inf(λ_certified)
|
||||
end
|
||||
|
|
@ -1,77 +0,0 @@
|
|||
using IntervalArithmetic
|
||||
|
||||
IntervalArithmetic.setrounding(Interval, :tight)
|
||||
IntervalArithmetic.setformat(sigfigs=12)
|
||||
|
||||
function fma_SOS_thr!(result::AbstractVector{T}, pm::AbstractMatrix{<:Integer},
|
||||
Q::AbstractMatrix{T}, acc_matrix=zeros(T, size(pm)...)) where T
|
||||
|
||||
s1, s2 = size(pm)
|
||||
|
||||
@inbounds for k in 1:s2
|
||||
let k=k, s1=s1, s2=s2, Q=Q, acc_matrix=acc_matrix
|
||||
Threads.@threads for j in 1:s2
|
||||
for i in 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 in 1:s2
|
||||
for i in 1:s1
|
||||
result[pm[i,j]] += acc_matrix[i,j]
|
||||
end
|
||||
end
|
||||
|
||||
return result
|
||||
end
|
||||
|
||||
function compute_SOS(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
|
||||
result = zeros(eltype(Q), maximum(pm));
|
||||
return fma_SOS_thr!(result, pm, Q)
|
||||
end
|
||||
|
||||
function compute_SOS(RG::GroupRing, Q::AbstractMatrix{<:Real})
|
||||
result = compute_SOS(RG.pm, Q)
|
||||
return GroupRingElem(result, RG)
|
||||
end
|
||||
|
||||
function compute_SOS_square(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix)
|
||||
result = zeros(eltype(Q), maximum(pm));
|
||||
|
||||
for i in 1:size(Q,2)
|
||||
GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), pm)
|
||||
end
|
||||
|
||||
return result
|
||||
end
|
||||
|
||||
function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix)
|
||||
return GroupRingElem(compute_SOS_square(RG.pm, Q), RG)
|
||||
end
|
||||
|
||||
function augIdproj(Q::AbstractMatrix{T}) where T
|
||||
result = zeros(T, size(Q))
|
||||
l = size(Q, 2)
|
||||
Threads.@threads for j in 1:l
|
||||
col = sum(view(Q, :,j))/l
|
||||
for i in 1:size(Q, 1)
|
||||
result[i,j] = Q[i,j] - col
|
||||
end
|
||||
end
|
||||
return result
|
||||
end
|
||||
|
||||
function augIdproj(::Type{Interval}, Q::AbstractMatrix{T}) where {T<:Real}
|
||||
result = zeros(Interval{T}, size(Q))
|
||||
l = size(Q, 2)
|
||||
Threads.@threads for j in 1:l
|
||||
col = sum(view(Q, :,j))/l
|
||||
for i in 1:size(Q, 1)
|
||||
result[i,j] = @interval(Q[i,j] - col)
|
||||
end
|
||||
end
|
||||
check = all([zero(T) in sum(view(result, :, i)) for i in 1:size(result, 2)])
|
||||
return result, check
|
||||
end
|
||||
Loading…
Reference in New Issue