mirror of
https://github.com/kalmarek/PropertyT.jl.git
synced 2024-09-07 07:11:42 +02:00
Merge remote-tracking branch 'origin/enh/rename_sgap' into 1703.09680v1
This commit is contained in:
commit
c5741f0ce3
143
src/PropertyT.jl
143
src/PropertyT.jl
@ -20,14 +20,14 @@ function pmΔfilenames(name::String)
|
|||||||
return pm_filename, Δ_coeff_filename
|
return pm_filename, Δ_coeff_filename
|
||||||
end
|
end
|
||||||
|
|
||||||
function κSDPfilenames(name::String)
|
function λSDPfilenames(name::String)
|
||||||
if !isdir(name)
|
if !isdir(name)
|
||||||
mkdir(name)
|
mkdir(name)
|
||||||
end
|
end
|
||||||
prefix = name
|
prefix = name
|
||||||
κ_filename = joinpath(prefix, "kappa.jld")
|
λ_filename = joinpath(prefix, "lambda.jld")
|
||||||
SDP_filename = joinpath(prefix, "SDPmatrixA.jld")
|
SDP_filename = joinpath(prefix, "SDPmatrix.jld")
|
||||||
return κ_filename, SDP_filename
|
return λ_filename, SDP_filename
|
||||||
end
|
end
|
||||||
|
|
||||||
function ΔandSDPconstraints(name::String)
|
function ΔandSDPconstraints(name::String)
|
||||||
@ -46,27 +46,36 @@ function ΔandSDPconstraints(name::String)
|
|||||||
return Δ, sdp_constraints
|
return Δ, sdp_constraints
|
||||||
end
|
end
|
||||||
|
|
||||||
function ΔandSDPconstraints(name::String, generating_set::Function)
|
function ΔandSDPconstraints(name::String, generating_set::Function, radius::Int)
|
||||||
|
try
|
||||||
|
return ΔandSDPconstraints(name)
|
||||||
|
catch err
|
||||||
|
if isa(err, ArgumentError)
|
||||||
pm_fname, Δ_fname = pmΔfilenames(name)
|
pm_fname, Δ_fname = pmΔfilenames(name)
|
||||||
S, ID = generating_set()
|
S, Id = generating_set()
|
||||||
Δ, sdp_constraints = Main.ΔandSDPconstraints(ID, S)
|
Δ, sdp_constraints = Main.ΔandSDPconstraints(Id, S, radius)
|
||||||
save(pm_fname, "pm", Δ.product_matrix)
|
save(pm_fname, "pm", Δ.product_matrix)
|
||||||
save(Δ_fname, "Δ", Δ.coefficients)
|
save(Δ_fname, "Δ", Δ.coefficients)
|
||||||
return Δ, sdp_constraints
|
return Δ, sdp_constraints
|
||||||
|
else
|
||||||
|
error(logger, err)
|
||||||
|
end
|
||||||
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
function κandA(name::String)
|
function λandP(name::String)
|
||||||
κ_fname, SDP_fname = κSDPfilenames(name)
|
λ_fname, SDP_fname = λSDPfilenames(name)
|
||||||
f₁ = isfile(κ_fname)
|
f₁ = isfile(λ_fname)
|
||||||
f₂ = isfile(SDP_fname)
|
f₂ = isfile(SDP_fname)
|
||||||
|
|
||||||
if f₁ && f₂
|
if f₁ && f₂
|
||||||
info(logger, "Loading precomputed κ, A...")
|
info(logger, "Loading precomputed λ, P...")
|
||||||
κ = load(κ_fname, "κ")
|
λ = load(λ_fname, "λ")
|
||||||
A = load(SDP_fname, "A")
|
P = load(SDP_fname, "P")
|
||||||
else
|
else
|
||||||
throw(ArgumentError("You need to precompute κ and SDP matrix A to load it!"))
|
throw(ArgumentError("You need to precompute λ and SDP matrix P to load it!"))
|
||||||
end
|
end
|
||||||
return κ, A
|
return λ, P
|
||||||
end
|
end
|
||||||
|
|
||||||
function timed_msg(t)
|
function timed_msg(t)
|
||||||
@ -78,41 +87,62 @@ function timed_msg(t)
|
|||||||
return "took: $elapsed s, allocated: $bytes_alloc bytes ($(gc_diff.poolalloc) allocations)."
|
return "took: $elapsed s, allocated: $bytes_alloc bytes ($(gc_diff.poolalloc) allocations)."
|
||||||
end
|
end
|
||||||
|
|
||||||
function κandA(name::String, sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver; upper_bound=Inf)
|
|
||||||
if isfile("$name/solver.log")
|
function λandP(name::String, opts...)
|
||||||
rm("$name/solver.log")
|
try
|
||||||
|
return λandP(name)
|
||||||
|
catch err
|
||||||
|
if isa(err, ArgumentError)
|
||||||
|
if isfile(joinpath(name, "solver.log"))
|
||||||
|
rm(joinpath(name, "solver.log"))
|
||||||
end
|
end
|
||||||
|
|
||||||
add_handler(solver_logger, DefaultHandler("./$name/solver.log", DefaultFormatter("{date}| {msg}")), "solver_log")
|
add_handler(solver_logger, DefaultHandler(joinpath(name, "solver.log"), DefaultFormatter("{date}| {msg}")), "solver_log")
|
||||||
|
|
||||||
info(logger, "Creating SDP problem...")
|
info(logger, "Creating SDP problem...")
|
||||||
|
|
||||||
|
λ, P = compute_λandP(opts...)
|
||||||
|
|
||||||
|
remove_handler(solver_logger, "solver_log")
|
||||||
|
|
||||||
|
λ_fname, P_fname = λSDPfilenames(name)
|
||||||
|
|
||||||
|
if λ > 0
|
||||||
|
save(λ_fname, "λ", λ)
|
||||||
|
save(P_fname, "P", P)
|
||||||
|
else
|
||||||
|
throw(ErrorException("Solver $solver did not produce a valid solution!: λ = $λ"))
|
||||||
|
end
|
||||||
|
return λ, P
|
||||||
|
|
||||||
|
else
|
||||||
|
# throw(err)
|
||||||
|
error(logger, err)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function compute_λandP(sdp_constraints, Δ::GroupAlgebraElement, solver::AbstractMathProgSolver, upper_bound=Inf)
|
||||||
|
|
||||||
t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound)
|
t = @timed SDP_problem = create_SDP_problem(sdp_constraints, Δ; upper_bound=upper_bound)
|
||||||
info(logger, timed_msg(t))
|
info(logger, timed_msg(t))
|
||||||
|
|
||||||
κ = 0.0
|
λ = 0.0
|
||||||
A = nothing
|
P = nothing
|
||||||
while κ == 0.0
|
while λ == 0.0
|
||||||
try
|
try
|
||||||
κ, A = solve_SDP(SDP_problem, solver)
|
λ, P = solve_SDP(SDP_problem, solver)
|
||||||
catch y
|
catch y
|
||||||
warn(solver_logger, y)
|
warn(solver_logger, y)
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
return λ, P
|
||||||
remove_handler(solver_logger, "solver_log")
|
|
||||||
|
|
||||||
κ_fname, A_fname = κSDPfilenames(name)
|
|
||||||
if κ > 0
|
|
||||||
save(κ_fname, "κ", κ)
|
|
||||||
save(A_fname, "A", A)
|
|
||||||
else
|
|
||||||
throw(ErrorException("Solver $solver did not produce a valid solution!: κ = $κ"))
|
|
||||||
end
|
|
||||||
return κ, A
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
Kazhdan_from_sgap(λ,N) = sqrt(2*λ/N)
|
||||||
|
|
||||||
function check_property_T(name::String, generating_set::Function,
|
function check_property_T(name::String, generating_set::Function,
|
||||||
solver, upper_bound, tol=1e-6)
|
solver, upper_bound, tol, radius)
|
||||||
|
|
||||||
if !isdir(name)
|
if !isdir(name)
|
||||||
mkdir(name)
|
mkdir(name)
|
||||||
@ -124,47 +154,36 @@ function check_property_T(name::String, generating_set::Function,
|
|||||||
info(logger, "Precision: $tol")
|
info(logger, "Precision: $tol")
|
||||||
info(logger, "Upper bound: $upper_bound")
|
info(logger, "Upper bound: $upper_bound")
|
||||||
|
|
||||||
Δ, sdp_constraints = try
|
Δ, sdp_constraints = ΔandSDPconstraints(name, generating_set, radius)
|
||||||
ΔandSDPconstraints(name)
|
|
||||||
catch err
|
|
||||||
if isa(err, ArgumentError)
|
|
||||||
ΔandSDPconstraints(name, generating_set)
|
|
||||||
else
|
|
||||||
error(logger, err)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
S = countnz(Δ.coefficients) - 1
|
S = countnz(Δ.coefficients) - 1
|
||||||
info(logger, "|S| = $S")
|
info(logger, "|S| = $S")
|
||||||
info(logger, "length(Δ) = $(length(Δ))")
|
info(logger, "length(Δ) = $(length(Δ))")
|
||||||
info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))")
|
info(logger, "size(Δ.product_matrix) = $(size(Δ.product_matrix))")
|
||||||
|
|
||||||
κ, A = try
|
λ, P = λandP(name, sdp_constraints, Δ, solver, upper_bound)
|
||||||
κandA(name)
|
|
||||||
catch err
|
|
||||||
if isa(err, ArgumentError)
|
|
||||||
κandA(name, sdp_constraints, Δ, solver; upper_bound=upper_bound)
|
|
||||||
else
|
|
||||||
# throw(err)
|
|
||||||
error(logger, err)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
info(logger, "κ = $κ")
|
info(logger, "λ = $λ")
|
||||||
info(logger, "sum(A) = $(sum(A))")
|
info(logger, "sum(P) = $(sum(P))")
|
||||||
info(logger, "maximum(A) = $(maximum(A))")
|
info(logger, "maximum(P) = $(maximum(P))")
|
||||||
info(logger, "minimum(A) = $(minimum(A))")
|
info(logger, "minimum(P) = $(minimum(P))")
|
||||||
|
|
||||||
if κ > 0
|
if λ > 0
|
||||||
spectral_gap = check_distance_to_positive_cone(Δ, κ, A, tol=tol, rational=false)
|
spectral_gap = check_distance_to_positive_cone(Δ, λ, P, tol=tol, rational=false)
|
||||||
|
if isa(spectral_gap, Interval)
|
||||||
|
spectral_gap = spectral_gap.lo
|
||||||
|
end
|
||||||
if spectral_gap > 0
|
if spectral_gap > 0
|
||||||
Kazhdan_κ = sqrt(2*spectral_gap/S)
|
@show spectral_gap
|
||||||
|
Kazhdan_κ = Kazhdan_from_sgap(spectral_gap, S)
|
||||||
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
|
||||||
info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
info(logger, "κ($name, S) ≥ $Kazhdan_κ: Group HAS property (T)!")
|
||||||
else
|
else
|
||||||
|
spectral_gap = Float64(trunc(spectral_gap, 12))
|
||||||
info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!")
|
info(logger, "λ($name, S) ≥ $spectral_gap: Group may NOT HAVE property (T)!")
|
||||||
end
|
end
|
||||||
else
|
else
|
||||||
info(logger, "κ($name, S) ≥ $κ < 0: Tells us nothing about property (T)")
|
info(logger, "κ($name, S) ≥ $λ < 0: Tells us nothing about property (T)")
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -2,10 +2,12 @@ import Base: rationalize
|
|||||||
|
|
||||||
using ValidatedNumerics
|
using ValidatedNumerics
|
||||||
ValidatedNumerics.setrounding(Interval, :correct)
|
ValidatedNumerics.setrounding(Interval, :correct)
|
||||||
|
# ValidatedNumerics.setrounding(Interval, :fast) #which is slower??
|
||||||
ValidatedNumerics.setformat(:standard)
|
ValidatedNumerics.setformat(:standard)
|
||||||
|
# setprecision(Interval, 53) # slightly faster than 256
|
||||||
|
|
||||||
function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T)
|
function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, λ::T)
|
||||||
return Δ*Δ - κ*Δ
|
return Δ*Δ - λ*Δ
|
||||||
end
|
end
|
||||||
|
|
||||||
function algebra_square(vector, elt)
|
function algebra_square(vector, elt)
|
||||||
@ -46,7 +48,7 @@ end
|
|||||||
import ValidatedNumerics.±
|
import ValidatedNumerics.±
|
||||||
|
|
||||||
function (±){T<:Number}(X::AbstractArray{T}, tol::Real)
|
function (±){T<:Number}(X::AbstractArray{T}, tol::Real)
|
||||||
r{T}(x::T) = (x == zero(T)? @biginterval(0) : x ± tol)
|
r{T}(x::T) = (x == zero(T)? @interval(0) : x ± tol)
|
||||||
return r.(X)
|
return r.(X)
|
||||||
end
|
end
|
||||||
|
|
||||||
@ -60,91 +62,97 @@ end
|
|||||||
|
|
||||||
ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol)
|
ℚ(x, tol::Real) = rationalize(BigInt, x, tol=tol)
|
||||||
|
|
||||||
function distance_to_cone{T<:Rational}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
|
function distance_to_cone{T<:Rational}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
|
||||||
SOS = compute_SOS(sqrt_matrix, Δ)
|
SOS = compute_SOS(sqrt_matrix, Δ)
|
||||||
|
|
||||||
SOS_diff = EOI(Δ, κ) - SOS
|
SOS_diff = EOI(Δ, λ) - SOS
|
||||||
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
||||||
|
|
||||||
info(logger, "κ = $κ (≈$(@sprintf("%.10f", float(κ)))")
|
info(logger, "λ = $λ (≈$(@sprintf("%.10f", float(λ)))")
|
||||||
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
||||||
if ɛ_dist ≠ 0//1
|
if ɛ_dist ≠ 0//1
|
||||||
warn(logger, "The SOS is not in the augmentation ideal, number below are meaningless!")
|
warn(logger, "The SOS is not in the augmentation ideal, number below are meaningless!")
|
||||||
end
|
end
|
||||||
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist")
|
info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) = $ɛ_dist")
|
||||||
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L₁_dist)))")
|
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ = $(@sprintf("%.10f", float(eoi_SOS_L₁_dist)))")
|
||||||
|
|
||||||
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist
|
distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
|
||||||
return distance_to_cone
|
return distance_to_cone
|
||||||
end
|
end
|
||||||
|
|
||||||
function distance_to_cone{T<:Rational, S<:Interval}(κ::T, sqrt_matrix::Array{S,2}, Δ::GroupAlgebraElement{T})
|
function distance_to_cone{T<:Rational, S<:Interval}(λ::T, sqrt_matrix::Array{S,2}, Δ::GroupAlgebraElement{T})
|
||||||
SOS = compute_SOS(sqrt_matrix, Δ)
|
SOS = compute_SOS(sqrt_matrix, Δ)
|
||||||
info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))")
|
info(logger, "ɛ(∑ξᵢ*ξᵢ) ∈ $(GroupAlgebras.ɛ(SOS))")
|
||||||
|
λⁱⁿᵗ = @interval(λ)
|
||||||
SOS_diff = EOI(Δ, κ) - SOS
|
Δⁱⁿᵗ = GroupAlgebraElement([@interval(c) for c in Δ.coefficients], Δ.product_matrix)
|
||||||
|
SOS_diff = EOI(Δⁱⁿᵗ, λⁱⁿᵗ) - SOS
|
||||||
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
||||||
|
|
||||||
info(logger, "κ = $κ (≈$(@sprintf("%.10f",float(κ))))")
|
info(logger, "λ = $λ (≈≥$(@sprintf("%.10f",float(λ))))")
|
||||||
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
||||||
|
|
||||||
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
|
info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(ɛ_dist)")
|
||||||
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L₁_dist)")
|
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(eoi_SOS_L₁_dist)")
|
||||||
|
|
||||||
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist
|
distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
|
||||||
return distance_to_cone
|
return distance_to_cone
|
||||||
end
|
end
|
||||||
|
|
||||||
function distance_to_cone{T<:AbstractFloat}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
|
function distance_to_cone{T<:AbstractFloat}(λ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T})
|
||||||
SOS = compute_SOS(sqrt_matrix, Δ)
|
SOS = compute_SOS(sqrt_matrix, Δ)
|
||||||
|
|
||||||
SOS_diff = EOI(Δ, κ) - SOS
|
SOS_diff = EOI(Δ, λ) - SOS
|
||||||
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
eoi_SOS_L₁_dist = norm(SOS_diff,1)
|
||||||
|
|
||||||
info(logger, "κ = $κ")
|
info(logger, "λ = $λ")
|
||||||
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
ɛ_dist = GroupAlgebras.ɛ(SOS_diff)
|
||||||
info(logger, "ɛ(Δ² - κΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
|
info(logger, "ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", ɛ_dist))")
|
||||||
info(logger, "‖Δ² - κΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L₁_dist))")
|
info(logger, "‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", eoi_SOS_L₁_dist))")
|
||||||
|
|
||||||
distance_to_cone = κ - 2^3*eoi_SOS_L₁_dist
|
distance_to_cone = λ - 2^3*eoi_SOS_L₁_dist
|
||||||
return distance_to_cone
|
return distance_to_cone
|
||||||
end
|
end
|
||||||
|
|
||||||
function check_distance_to_positive_cone(Δ::GroupAlgebraElement, κ, A;
|
function check_distance_to_positive_cone(Δ::GroupAlgebraElement, λ, P;
|
||||||
tol=1e-7, rational=false)
|
tol=1e-7, rational=false)
|
||||||
|
|
||||||
isapprox(eigvals(A), abs(eigvals(A)), atol=tol) ||
|
isapprox(eigvals(P), abs(eigvals(P)), atol=tol) ||
|
||||||
warn("The solution matrix doesn't seem to be positive definite!")
|
warn("The solution matrix doesn't seem to be positive definite!")
|
||||||
@assert A == Symmetric(A)
|
@assert P == Symmetric(P)
|
||||||
A_sqrt = real(sqrtm(A))
|
Q = real(sqrtm(P))
|
||||||
|
|
||||||
info(logger, "------------------------------------------------------------")
|
info(logger, "------------------------------------------------------------")
|
||||||
info(logger, "")
|
info(logger, "")
|
||||||
info(logger, "Checking in floating-point arithmetic...")
|
info(logger, "Checking in floating-point arithmetic...")
|
||||||
t = @timed fp_distance = distance_to_cone(κ, A_sqrt, Δ)
|
t = @timed fp_distance = distance_to_cone(λ, Q, Δ)
|
||||||
info(logger, timed_msg(t))
|
info(logger, timed_msg(t))
|
||||||
info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))")
|
info(logger, "Floating point distance (to positive cone) ≈ $(@sprintf("%.10f", fp_distance))")
|
||||||
info(logger, "------------------------------------------------------------")
|
info(logger, "------------------------------------------------------------")
|
||||||
|
|
||||||
info(logger, "Projecting columns of rationalized A_sqrt to the augmentation ideal...")
|
# if fp_distance ≤ 0
|
||||||
δ = eps(κ)
|
# return fpdistance
|
||||||
A_sqrt_ℚ = ℚ(A_sqrt, δ)
|
# end
|
||||||
A_sqrt_ℚ_aug = correct_to_augmentation_ideal(A_sqrt_ℚ)
|
|
||||||
κ_ℚ = ℚ(κ, δ)
|
info(logger, "Projecting columns of rationalized Q to the augmentation ideal...")
|
||||||
|
δ = eps(λ)
|
||||||
|
Q_ℚ = ℚ(Q, δ)
|
||||||
|
t = @timed Q_ℚω = correct_to_augmentation_ideal(Q_ℚ)
|
||||||
|
info(logger, timed_msg(t))
|
||||||
|
λ_ℚ = ℚ(λ, δ)
|
||||||
Δ_ℚ = ℚ(Δ, δ)
|
Δ_ℚ = ℚ(Δ, δ)
|
||||||
|
|
||||||
info(logger, "Checking in interval arithmetic")
|
info(logger, "Checking in interval arithmetic")
|
||||||
A_sqrt_ℚ_augᴵ = A_sqrt_ℚ_aug ± δ
|
Q_ℚωⁱⁿᵗ = Float64.(Q_ℚω) ± δ
|
||||||
t = @timed Interval_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_augᴵ, Δ_ℚ)
|
t = @timed Interval_dist_to_Σ² = distance_to_cone(λ_ℚ, Q_ℚωⁱⁿᵗ, Δ_ℚ)
|
||||||
info(logger, timed_msg(t))
|
info(logger, timed_msg(t))
|
||||||
info(logger, "The Augmentation-projected actual distance (to positive cone) ≥ $(@sprintf("%.10f", Interval_dist_to_Σ².lo))")
|
info(logger, "The Augmentation-projected actual distance (to positive cone) ∈ $(Interval_dist_to_Σ²)")
|
||||||
info(logger, "------------------------------------------------------------")
|
info(logger, "------------------------------------------------------------")
|
||||||
|
|
||||||
if Interval_dist_to_Σ².lo ≤ 0 || !rational
|
if Interval_dist_to_Σ².lo ≤ 0 || !rational
|
||||||
return Interval_dist_to_Σ².lo
|
return Interval_dist_to_Σ²
|
||||||
else
|
else
|
||||||
info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...")
|
info(logger, "Checking Projected SOS decomposition in exact rational arithmetic...")
|
||||||
t = @timed ℚ_dist_to_Σ² = distance_to_cone(κ_ℚ, A_sqrt_ℚ_aug, Δ_ℚ)
|
t = @timed ℚ_dist_to_Σ² = distance_to_cone(λ_ℚ, Q_ℚω, Δ_ℚ)
|
||||||
info(logger, timed_msg(t))
|
info(logger, timed_msg(t))
|
||||||
@assert isa(ℚ_dist_to_Σ², Rational)
|
@assert isa(ℚ_dist_to_Σ², Rational)
|
||||||
info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))")
|
info(logger, "Augmentation-projected rational distance (to positive cone) ≥ $(Float64(trunc(ℚ_dist_to_Σ²,8)))")
|
||||||
|
20
src/sdps.jl
20
src/sdps.jl
@ -52,21 +52,21 @@ function create_SDP_problem(matrix_constraints, Δ::GroupAlgebraElement; upper_b
|
|||||||
Δ² = Δ*Δ
|
Δ² = Δ*Δ
|
||||||
@assert length(Δ) == length(matrix_constraints)
|
@assert length(Δ) == length(matrix_constraints)
|
||||||
m = JuMP.Model();
|
m = JuMP.Model();
|
||||||
JuMP.@variable(m, A[1:N, 1:N], SDP)
|
JuMP.@variable(m, P[1:N, 1:N], SDP)
|
||||||
JuMP.@SDconstraint(m, A >= 0)
|
JuMP.@SDconstraint(m, P >= 0)
|
||||||
JuMP.@constraint(m, sum(A[i] for i in eachindex(A)) == 0)
|
JuMP.@constraint(m, sum(P[i] for i in eachindex(P)) == 0)
|
||||||
|
|
||||||
if upper_bound < Inf
|
if upper_bound < Inf
|
||||||
JuMP.@variable(m, 0.0 <= κ <= upper_bound)
|
JuMP.@variable(m, 0.0 <= λ <= upper_bound)
|
||||||
else
|
else
|
||||||
JuMP.@variable(m, κ >= 0)
|
JuMP.@variable(m, λ >= 0)
|
||||||
end
|
end
|
||||||
|
|
||||||
for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients)
|
for (pairs, δ², δ) in zip(matrix_constraints, Δ².coefficients, Δ.coefficients)
|
||||||
JuMP.@constraint(m, sum(A[i,j] for (i,j) in pairs) == δ² - κ*δ)
|
JuMP.@constraint(m, sum(P[i,j] for (i,j) in pairs) == δ² - λ*δ)
|
||||||
end
|
end
|
||||||
|
|
||||||
JuMP.@objective(m, Max, κ)
|
JuMP.@objective(m, Max, λ)
|
||||||
|
|
||||||
return m
|
return m
|
||||||
end
|
end
|
||||||
@ -95,7 +95,7 @@ function solve_SDP(SDP_problem, solver)
|
|||||||
end
|
end
|
||||||
info(logger, solution_status)
|
info(logger, solution_status)
|
||||||
|
|
||||||
κ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :κ))
|
λ = JuMP.getvalue(JuMP.getvariable(SDP_problem, :λ))
|
||||||
A = JuMP.getvalue(JuMP.getvariable(SDP_problem, :A))
|
P = JuMP.getvalue(JuMP.getvariable(SDP_problem, :P))
|
||||||
return κ, A
|
return λ, P
|
||||||
end
|
end
|
||||||
|
Loading…
Reference in New Issue
Block a user