diff --git a/src/1712.07167.jl b/src/1712.07167.jl new file mode 100644 index 0000000..9f8ffc2 --- /dev/null +++ b/src/1712.07167.jl @@ -0,0 +1,248 @@ +############################################################################### +# +# Settings and filenames +# +############################################################################### + +struct Symmetrize end +struct Naive end + +abstract type PropertyTSettings end + +struct SolverSettings + sdpsolver::AbstractMathProgSolver + upper_bound::Float64 + warmstart::Bool + + SolverSettings(sol, ub, ws=true) = new(sol, upper_bound, ws) +end + +struct Naive <: PropertyTSettings + name::String + G::Group + S::Vector{GroupElem} + radius::Int + + solver::SolverSettings +end + +struct Symmetrized <: PropertyTSettings + name::String + G::Group + S::Vector{GroupElem} + autS::Group + radius::Int + + solver::SolverSettings +end + +function Settings(name::String, + G::Group, S::Vector{GEl}, r::Integer, + sol::Solver, ub, ws=true) where {GEl<:GroupElem, Solver<:AbstractMathProgSolver} + sol_sett = SolverSettings(sol, ub, ws) + return Naive(name, G, S, r, sol_sett) +end + +function Settings(name::String, + G::Group, S::Vector{GEl}, autS::Group, r::Integer, + sol::Solver, ub, ws=true) where {GEl<:GroupElem, Solver<:AbstractMathProgSolver} + sol_sett = SolverSettings(sol, ub, ws) + return Symmetrized(name, G, S, autS, r, sol_sett) +end + +prefix(s::Naive) = s.name +prefix(s::Symmetrized) = "o"*s.name +suffix(s::PropertyTSettings) = "$(s.upper_bound)" +prepath(s::PropertyTSettings) = prefix(s) +fullpath(s::PropertyTSettings) = joinpath(prefix(s), suffix(s)) + +filename(sett::PropertyTSettings, s::Symbol) = filename(sett, Val{s}) + +filename(sett::PropertyTSettings, ::Type{Val{:fulllog}}) = + joinpath(fullpath(sett), "full_$(string(now())).log") +filename(sett::PropertyTSettings, ::Type{Val{:solverlog}}) = + joinpath(fullpath(sett), "solver_$(string(now())).log") + +filename(sett::PropertyTSettings, ::Type{Val{:Δ}}) = + joinpath(prepath(sett), "delta.jld") +filename(sett::PropertyTSettings, ::Type{Val{:OrbitData}}) = + joinpath(prepath(sett), "OrbitData.jld") + +filename(sett::PropertyTSettings, ::Type{Val{:warmstart}}) = + joinpath(fullpath(sett), "warmstart.jld") +filename(sett::PropertyTSettings, ::Type{Val{:solution}}) = + joinpath(fullpath(sett), "solution.jld") + +############################################################################### +# +# λandP +# +############################################################################### + +function warmstart(sett::PropertyTSettings) + if sett.solver.warmstart && isfile(filename(sett, :warmstart)) + ws = load(filename(sett, :warmstart), "warmstart") + else + ws = nothing + end + return ws +end + +function computeλandP(sett::Naive, Δ::GroupRingElem; + solverlog=tempname()*".log") + + info("Creating SDP problem...") + SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.solver.upper_bound) + JuMP.setsolver(SDP_problem, sett.solver.sdpsolver) + info(Base.repr(SDP_problem)) + + ws = warmstart(sett) + @time status, (λ, P, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws) + @show status + save(filename(sett, :warmstart), "warmstart", ws, "P", P, "λ", λ) + + return λ, P +end + +function computeλandP(sett::Symmetrized, Δ::GroupRingElem; + solverlog=tempname()*".log") + + if !isfile(filename(sett, :OrbitData)) + isdefined(parent(Δ), :basis) || throw("You need to define basis of Group Ring to compute orbit decomposition!") + orbit_data = OrbitData(parent(Δ), sett.autS) + save(filename(sett, :OrbitData), "OrbitData", orbit_data) + end + orbit_data = load(filename(sett, :OrbitData), "OrbitData") + orbit_data = decimate(orbit_data) + + info("Creating SDP problem...") + + SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, orbit_data, upper_bound=sett.solver.upper_bound) + JuMP.setsolver(SDP_problem, sett.solver.sdpsolver) + info(Base.repr(SDP_problem)) + + ws = warmstart(sett) + @time status, (λ, Ps, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws) + @show status + save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ) + + info("Reconstructing P...") + @time P = reconstruct(Ps, orbit_data) + + return λ, P +end + +############################################################################### +# +# Checking solution +# +############################################################################### + +function distance_to_positive_cone(Δ::GroupRingElem, λ, Q; R::Int=2) + @info("------------------------------------------------------------") + @info("Checking in floating-point arithmetic...") + @info("λ = $λ") + eoi = Δ^2-λ*Δ + + @time residual = eoi - compute_SOS(parent(eoi), augIdproj(Q)) + @info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residual)))") + L1_norm = norm(residual,1) + @info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))") + + distance = λ - 2.0^(2ceil(log2(R)))*L1_norm + + @info("Floating point distance (to positive cone) ≈") + @info("$(@sprintf("%.10f", distance))") + @info("") + + if distance ≤ 0 + return distance + end + + @info("------------------------------------------------------------") + @info("Checking in interval arithmetic...") + @info("λ ∈ $λ") + + λ = @interval(λ) + eoi = Δ^2 - λ*Δ + + @info("Projecting columns of Q to the augmentation ideal...") + @time Q, check = augIdproj(Interval, Q) + @info("Checking that sum of every column contains 0.0... ") + @info((check? "They do." : "FAILED!")) + check || @warn("The following numbers are meaningless!") + + @time residual = eoi - compute_SOS(parent(eoi), Q) + @info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ∈ $(aug(residual))") + L1_norm = norm(residual,1) + @info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)") + + distance = λ - 2.0^(2ceil(log2(R)))*L1_norm + + @info("Interval distance (to positive cone) ∈") + @info("$(distance)") + @info("") + + return distance.lo +end + +############################################################################### +# +# Interpreting the numerical results +# +############################################################################### + +Kazhdan(λ::Number, N::Integer) = sqrt(2*λ/N) + +function interpret_results(sett::PropertyTSettings, sgap::Number) + + if sgap > 0 + Kazhdan_κ = Kazhdan(sgap, length(sett.S)) + if Kazhdan_κ > 0 + info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!") + return true + end + end + info("λ($(sett.name), S) ≥ $sgap < 0: Tells us nothing about property (T)") + return false +end + +function check_property_T(sett::PropertyTSettings) + fp = PropertyT.fullpath(sett) + isdir(fp) || mkpath(fp) + + if isfile(filename(sett,:Δ)) + # cached + Δ = loadLaplacian(filename(sett,:Δ), sett.G) + else + # compute + Δ = Laplacian(sett.S, sett.radius) + saveLaplacian(filename(sett, :Δ), Δ) + end + + if !sett.warmstart && isfile(filename(sett, :solution)) + λ, P = load(filename(sett, :solution), "λ", "P") + else + λ, P = computeλandP(sett, Δ, + solverlog=filename(sett, :solverlog)) + + save(filename(sett, :solution), "λ", λ, "P", P) + + if λ < 0 + warn("Solver did not produce a valid solution!") + end + end + + info("λ = $λ") + info("sum(P) = $(sum(P))") + info("maximum(P) = $(maximum(P))") + info("minimum(P) = $(minimum(P))") + + isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || + warn("The solution matrix doesn't seem to be positive definite!") + + @time Q = real(sqrtm((P+P')/2)) + sgap = distance_to_positive_cone(Δ, λ, Q, wlen=2*sett.radius) + + return interpret_results(sett, sgap) +end diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 89f5fd1..156dd4c 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -12,126 +12,10 @@ using JuMP import MathProgBase.SolverInterface.AbstractMathProgSolver -############################################################################### -# -# Settings and filenames -# -############################################################################### - -struct Symmetrize end -struct Naive end - -struct Settings{T, GEl<:GroupElem} - name::String - - G::Group - S::Vector{GEl} - radius::Int - - solver::AbstractMathProgSolver - upper_bound::Float64 - tol::Float64 - warmstart::Bool - - autS::Group - - function Settings(name::String, - G::Group, S::Vector{GEl}, r::Int, - sol::Sol, ub, tol, ws) where {GEl<:GroupElem, Sol<:AbstractMathProgSolver} - return new{Naive, GEl}(name, G, S, r, sol, ub, tol, ws) - end - - function Settings(name::String, - G::Group, S::Vector{GEl}, r::Int, - sol::Sol, ub, tol, ws, autS) where {GEl<:GroupElem, Sol<:AbstractMathProgSolver} - return new{Symmetrize, GEl}(name, G, S, r, sol, ub, tol, ws, autS) - end -end - - -prefix(s::Settings{Naive}) = s.name -prefix(s::Settings{Symmetrize}) = "o"*s.name -suffix(s::Settings) = "$(s.upper_bound)" -prepath(s::Settings) = prefix(s) -fullpath(s::Settings) = joinpath(prefix(s), suffix(s)) - -filename(sett::Settings, s::Symbol) = filename(sett, Val{s}) - -filename(sett::Settings, ::Type{Val{:fulllog}}) = - joinpath(fullpath(sett), "full_$(string(now())).log") -filename(sett::Settings, ::Type{Val{:solverlog}}) = - joinpath(fullpath(sett), "solver_$(string(now())).log") - -filename(sett::Settings, ::Type{Val{:Δ}}) = - joinpath(prepath(sett), "delta.jld") -filename(sett::Settings, ::Type{Val{:OrbitData}}) = - joinpath(prepath(sett), "OrbitData.jld") - -filename(sett::Settings, ::Type{Val{:warmstart}}) = - joinpath(fullpath(sett), "warmstart.jld") -filename(sett::Settings, ::Type{Val{:solution}}) = - joinpath(fullpath(sett), "solution.jld") - -function check_property_T(sett::Settings) - fp = PropertyT.fullpath(sett) - isdir(fp) || mkpath(fp) - - if isfile(filename(sett,:Δ)) - # cached - Δ = loadLaplacian(filename(sett,:Δ), sett.G) - else - # compute - Δ = Laplacian(sett.S, sett.radius) - saveLaplacian(filename(sett, :Δ), Δ) - end - - if !sett.warmstart && isfile(filename(sett, :solution)) - λ, P = load(filename(sett, :solution), "λ", "P") - else - λ, P = computeλandP(sett, Δ, - solverlog=filename(sett, :solverlog)) - - save(filename(sett, :solution), "λ", λ, "P", P) - - if λ < 0 - warn("Solver did not produce a valid solution!") - end - end - - info("λ = $λ") - info("sum(P) = $(sum(P))") - info("maximum(P) = $(maximum(P))") - info("minimum(P) = $(minimum(P))") - - isapprox(eigvals(P), abs.(eigvals(P)), atol=sett.tol) || - warn("The solution matrix doesn't seem to be positive definite!") - - @time Q = real(sqrtm((P+P')/2)) - sgap = distance_to_cone(Δ, λ, Q, wlen=2*sett.radius) - - return interpret_results(sett, sgap) -end - -Kazhdan(λ::Number, N::Integer) = sqrt(2*λ/N) - -function interpret_results(sett::Settings, sgap::Number) - - if sgap > 0 - Kazhdan_κ = Kazhdan(sgap, length(sett.S)) - if Kazhdan_κ > 0 - info("κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!") - return true - end - end - info("λ($(sett.name), S) ≥ $sgap < 0: Tells us nothing about property (T)") - return false -end - include("laplacians.jl") include("RGprojections.jl") include("orbitdata.jl") include("sos_sdps.jl") include("checksolution.jl") - end # module Property(T) diff --git a/src/checksolution.jl b/src/checksolution.jl index d61390e..91af92f 100644 --- a/src/checksolution.jl +++ b/src/checksolution.jl @@ -28,52 +28,5 @@ function augIdproj(Q::AbstractArray{T,2}) where {T<:Real} return result end -function distance_to_cone(Δ::GroupRingElem, λ, Q; wlen::Int=4) - info("------------------------------------------------------------") - info("Checking in floating-point arithmetic...") - info("λ = $λ") - @time sos = compute_SOS(parent(Δ), Q) - residue = Δ^2-λ*Δ - sos - info("ɛ(Δ² - λΔ - ∑ξᵢ*ξᵢ) ≈ $(@sprintf("%.10f", aug(residue)))") - L1_norm = norm(residue,1) - info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ≈ $(@sprintf("%.10f", L1_norm))") - - distance = λ - 2^(wlen-1)*L1_norm - - info("Floating point distance (to positive cone) ≈") - info("$(@sprintf("%.10f", distance))") - info("") - - if distance ≤ 0 - return distance end - - info("------------------------------------------------------------") - info("Checking in interval arithmetic...") - info("λ ∈ $λ") - - λ = @interval(λ) - eoi = Δ^2 - λ*Δ - info("Projecting columns of Q to the augmentation ideal...") - T = eltype(Q) - @time Q = augIdproj(Q) - - info("Checking that sum of every column contains 0.0... ") - check = all([zero(T) in sum(view(Q, :, i)) for i in 1:size(Q, 2)]) - info((check? "They do." : "FAILED!")) - - @assert check - - @time sos = compute_SOS(parent(Δ), Q) - residue = Δ^2-λ*Δ - sos - info("ɛ(∑ξᵢ*ξᵢ) ∈ $(aug(residue))") - L1_norm = norm(residue,1) - info("‖Δ² - λΔ - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)") - - distance = λ - 2^(wlen-1)*L1_norm - info("The Augmentation-projected distance (to positive cone) ∈") - info("$(distance)") - info("") - - return distance.lo end diff --git a/src/laplacians.jl b/src/laplacians.jl index 1e4010a..207d1ea 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -61,61 +61,3 @@ function loadGRElem(fname::String, G::Group) end return Δ end - -############################################################################### -# -# λandP -# -############################################################################### - -function computeλandP(sett::Settings{Naive}, Δ::GroupRingElem; - solverlog=tempname()*".log") - - info("Creating SDP problem...") - SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, upper_bound=sett.upper_bound) - JuMP.setsolver(SDP_problem, sett.solver) - info(Base.repr(SDP_problem)) - - ws = warmstart(sett) - @time status, (λ, P, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws) - save(filename(sett, :warmstart), "warmstart", ws) - - return λ, P -end - -function computeλandP(sett::Settings{Symmetrize}, Δ::GroupRingElem; - solverlog=tempname()*".log") - - if isfile(filename(sett, :OrbitData)) - orbit_data = load(filename(sett, :OrbitData), "OrbitData") - else - isdefined(parent(Δ), :basis) || throw("You need to define basis of Group Ring to compute orbit decomposition!") - orbit_data = OrbitData(parent(Δ), sett.autS) - save(filename(sett, :OrbitData), "OrbitData", orbit_data) - end - orbit_data = decimate(orbit_data) - - info("Creating SDP problem...") - - SDP_problem, varλ, varP = SOS_problem(Δ^2, Δ, orbit_data, upper_bound=sett.upper_bound) - JuMP.setsolver(SDP_problem, sett.solver) - info(Base.repr(SDP_problem)) - - ws = warmstart(sett) - @time status, (λ, Ps, ws) = PropertyT.solve(solverlog, SDP_problem, varλ, varP, ws) - save(filename(sett, :warmstart), "warmstart", ws, "Ps", Ps, "λ", λ) - - info("Reconstructing P...") - @time P = reconstruct(Ps, orbit_data) - - return λ, P -end - -function warmstart(sett::Settings) - if sett.warmstart && isfile(filename(sett, :warmstart)) - ws = load(filename(sett, :warmstart), "warmstart") - else - ws = nothing - end - return ws -end