From d096d06b0da3be8bc08bc91f7ac3c3f6da6be862 Mon Sep 17 00:00:00 2001 From: Joaquim Garcia Date: Tue, 22 Sep 2020 13:13:29 -0300 Subject: [PATCH 1/4] Organize eigsolvers and options --- src/eigsolver.jl | 207 +++++++++++++++++++++++++----------------- src/pdhg.jl | 2 +- src/prox_operators.jl | 58 ++++++++---- src/structs.jl | 20 ++-- test/moitest.jl | 8 ++ 5 files changed, 189 insertions(+), 106 deletions(-) diff --git a/src/eigsolver.jl b/src/eigsolver.jl index 7d1e773..c9825a0 100644 --- a/src/eigsolver.jl +++ b/src/eigsolver.jl @@ -333,12 +333,24 @@ c WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1). c c----------------------------------------------------------------------------- =# -mutable struct ARPACKAlloc{T} +mutable struct EigSolverAlloc{T} + converged::Bool n::Int nev::Int ncv::Int maxiter::Int + tol::Float64 + resid::Vector{T} + converged_eigs::Int + + rng::Random.MersenneTwister + + # KrylovKit data + vals::Vector{Float64} + vecs::Vector{Vector{Float64}} + + # Arpack data bmat::String which::String mode::Int @@ -348,7 +360,6 @@ mutable struct ARPACKAlloc{T} v::Matrix{T} workd::Vector{T} workl::Vector{T} - resid::Vector{T} info::Vector{BlasInt} iparam::Vector{BlasInt} ipntr::Vector{BlasInt} @@ -359,35 +370,63 @@ mutable struct ARPACKAlloc{T} info_e::Vector{BlasInt} d::Vector{T} sigmar::Vector{T} - converged::Bool arpackerror::Bool x::Vector{T} y::Vector{T} - rng::Random.MersenneTwister - function ARPACKAlloc{T}() where T + function EigSolverAlloc{T}() where T new{T}() end end -hasconverged(arc::ARPACKAlloc) = arc.converged +hasconverged(arc::EigSolverAlloc) = arc.converged -function unsafe_getvalues(arc::ARPACKAlloc) - return arc.d +function EigSolverAlloc(T::DataType, n::Int64, opt::Options)::EigSolverAlloc + arc = EigSolverAlloc{T}() + @timeit "init_arc" _init_eigsolver!(arc, Symmetric(Matrix{T}(I, n, n), :U), 1, n, opt) + return arc end -function unsafe_getvectors(arc::ARPACKAlloc) - return arc.v +function eigsolver_init_resid!(arc, opt, n) + arc.resid = zeros(n) + arc.rng = Random.MersenneTwister(opt.eigsolver_resid_seed) + return nothing +end +function eigsolver_update_resid!(arc, opt, init) + if init == 2 + Random.seed!(arc.rng, opt.eigsolver_resid_seed) + Random.rand!(arc.rng, arc.resid) + elseif init == 3 + Random.seed!(arc.rng, opt.eigsolver_resid_seed) + Random.randn!(arc.rng, arc.resid) + LinearAlgebra.normalize!(arc.resid) + elseif init == 1 + fill!(arc.resid, 1.0) + else + fill!(arc.resid, 0.0) + end + return nothing end -function ARPACKAlloc(T::DataType, n::Int64, opt::Options)::ARPACKAlloc - arc = ARPACKAlloc{T}() - @timeit "init_arc" _init_arc!(arc, Symmetric(Matrix{T}(I, n, n), :U), 1, n, opt) - return arc +function _init_eigsolver!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T + if opt.eigsolver == 1 + arpack_init!(arc, A, nev, n, opt) + else + krylovkit_init!(arc, A, nev, n, opt) + end + return nothing +end + +function arpack_getvalues(arc::EigSolverAlloc) + return arc.d +end + +function arpack_getvectors(arc::EigSolverAlloc) + return arc.v end -function _init_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T +function arpack_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T # IDO - Integer. (INPUT/OUTPUT) (in julia its is a integer vector) # Reverse communication flag. @@ -425,18 +464,8 @@ function _init_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, # On OUTPUT: # RESID contains the final residual vector. # reset the curent mersenne twister to keep determinism - if opt.arpack_resid_init == 2 - arc.rng = Random.MersenneTwister(opt.arpack_resid_seed) - arc.resid = Random.rand(arc.rng, arc.n) - elseif opt.arpack_resid_init == 3 - arc.rng = Random.MersenneTwister(opt.arpack_resid_seed) - arc.resid = Random.randn(arc.rng, arc.n) - LinearAlgebra.normalize!(arc.resid) - elseif opt.arpack_resid_init == 1 - arc.resid = ones(arc.n) - else - arc.resid = zeros(arc.n) - end + eigsolver_init_resid!(arc, opt, n) + eigsolver_update_resid!(arc, opt, opt.arpack_resid_init) # NCV - Integer. (INPUT) # How many Lanczos vectors are generated @@ -450,7 +479,7 @@ function _init_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, # also increases the work and storage required to maintain the orthogonal # basis vectors. The optimal "cross-over" with respect to CPU time # is problem dependent and must be determined empirically. - arc.ncv = max(2 * arc.nev + 1, opt.arpack_min_lanczos) + arc.ncv = max(2 * arc.nev + 1, opt.eigsolver_min_lanczos) # TODO: 10 might be a way to tight bound # why not 20 @@ -539,14 +568,14 @@ function _init_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, return nothing end -function _update_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, opt::Options, up_ncv::Bool) where T +function arpack_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, opt::Options, up_ncv::Bool) where T need_resize = up_ncv if !need_resize need_resize = nev != arc.nev end if !need_resize - need_resize = arc.ncv != max(2 * arc.nev + 1, opt.arpack_min_lanczos) + need_resize = arc.ncv != max(2 * arc.nev + 1, opt.eigsolver_min_lanczos) end # IDO must be zero on the first call to dsaupd @@ -560,18 +589,7 @@ function _update_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64 # RESID - Double precision array of length N. (INPUT/OUTPUT) if opt.arpack_reset_resid - if opt.arpack_resid_init == 2 - Random.seed!(arc.rng, opt.arpack_resid_seed) - Random.rand!(arc.rng, arc.resid) - elseif opt.arpack_resid_init == 3 - Random.seed!(arc.rng, opt.arpack_resid_seed) - Random.randn!(arc.rng, arc.resid) - LinearAlgebra.normalize!(arc.resid) - elseif opt.arpack_resid_init == 1 - fill!(arc.resid, 1.0) - else - fill!(arc.resid, 0.0) - end + eigsolver_update_resid!(arc, opt, opt.arpack_resid_init) end # How many Lanczos vectors are generated. Needs to be NCV > NEV. It is recommended that NCV >= 2*NEV. @@ -581,7 +599,7 @@ function _update_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64 if up_ncv arc.ncv += max(arc.nev, 10) else - arc.ncv = max(2 * arc.nev + 1, opt.arpack_min_lanczos) + arc.ncv = max(2 * arc.nev + 1, opt.eigsolver_min_lanczos) end end @@ -646,7 +664,7 @@ function _update_arc!(arc::ARPACKAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64 return nothing end -function _saupd!(arc::ARPACKAlloc{T})::Nothing where T +function _saupd!(arc::EigSolverAlloc{T})::Nothing where T while true Arpack.saupd(arc.ido, arc.bmat, arc.n, arc.which, arc.nev, arc.TOL, arc.resid, arc.ncv, arc.v, arc.n, @@ -693,7 +711,7 @@ function _saupd!(arc::ARPACKAlloc{T})::Nothing where T return nothing end -function _seupd!(arc::ARPACKAlloc{T})::Nothing where T +function _seupd!(arc::EigSolverAlloc{T})::Nothing where T # Check if _AUPD has converged properly if !arc.converged @@ -726,50 +744,77 @@ function _seupd!(arc::ARPACKAlloc{T})::Nothing where T return nothing end -function eig!(arc::ARPACKAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T +function arpack_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T up_ncv = false - if opt.eigsolver == 1 - for i in 1:1 - # Initialize parameters and do memory allocation - @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing + for i in 1:1 + # Initialize parameters and do memory allocation + @timeit "update_arc" arpack_update!(arc, A, nev, opt, up_ncv)::Nothing - # Top level reverse communication interface to solve real double precision symmetric problems. - @timeit "saupd" _saupd!(arc)::Nothing + # Top level reverse communication interface to solve real double precision symmetric problems. + @timeit "saupd" _saupd!(arc)::Nothing - if arc.nev > arc.iparam[5] - up_ncv = true - else - break - end + if arc.nev > arc.iparam[5] + up_ncv = true + else + break end + end - # Post processing routine (eigenvalues and eigenvectors purification) - @timeit "seupd" _seupd!(arc)::Nothing - else - @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing - @timeit "krylovkit" begin - vals, vecs, info = eigsolve( - A, arc.resid, arc.nev, :LR,Lanczos(krylovdim = arc.ncv, tol = 1e-12, maxiter=100) - ) + # Post processing routine (eigenvalues and eigenvectors purification) + @timeit "seupd" _seupd!(arc)::Nothing + return nothing +end - if info.converged == 0 - @show arc.nev, info - end - fill!(arc.d, 0.0) - fill!(arc.v, 0.0) - for i in 1:min(arc.nev, info.converged) - arc.d[i] = vals[i] - for j in 1:arc.n - arc.v[j,i] = vecs[i][j] - end - end - arc.converged = true - arc.arpackerror = false - # if info.converged == 1 - # arc.converged = false - # end +function krylovkit_getvalues(arc::EigSolverAlloc) + return arc.vals +end + +function krylovkit_getvector(arc::EigSolverAlloc, i) + return arc.v[i] +end + +function krylovkit_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T + arc.n = n + arc.nev = nev + eigsolver_init_resid!(arc, opt, n) + eigsolver_update_resid!(arc, opt, opt.krylovkit_resid_init) + arc.ncv = max(2 * arc.nev + 1, opt.eigsolver_min_lanczos) + return nothing +end + +function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T + arc.n = n + arc.nev = nev + if opt.krylovkit_reset_resid + eigsolver_update_resid!(arc, opt, opt.krylovkit_resid_init) + end + arc.ncv = max(2 * arc.nev + 1, opt.eigsolver_min_lanczos) + return nothing +end + +function krylovkit_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T + up_ncv = false + arc.converged = true + @timeit "update_arc" krylovkit_update!(arc, A, nev, opt, up_ncv)::Nothing + @timeit "krylovkit" begin + vals, vecs, info = KrylovKit.eigsolve( + A, arc.resid, arc.nev, :LR, + KrylovKit.Lanczos( + KrylovKit.KrylovDefaults.orth, + arc.ncv, + opt.krylovkit_max_iter, + opt.krylovkit_tol, + opt.krylovkit_eager, + opt.krylovkit_verbose + ) + ) + arc.converged_eigs = info.converged + if info.converged == 0 + arc.converged = false end + arc.vals = vals + arc.vecs = vecs end return nothing end \ No newline at end of file diff --git a/src/pdhg.jl b/src/pdhg.jl index bc4085a..111c969 100644 --- a/src/pdhg.jl +++ b/src/pdhg.jl @@ -16,7 +16,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP p.target_rank = 2 * ones(length(conic_sets.sdpcone)) p.current_rank = 2 * ones(length(conic_sets.sdpcone)) p.min_eig = zeros(length(conic_sets.sdpcone)) - arc_list = [ARPACKAlloc(Float64, sdp.sq_side, opt) for (idx, sdp) in enumerate(conic_sets.sdpcone)] + arc_list = [EigSolverAlloc(Float64, sdp.sq_side, opt) for (idx, sdp) in enumerate(conic_sets.sdpcone)] ada_count = 0 # Print header diff --git a/src/prox_operators.jl b/src/prox_operators.jl index b35cb50..67dcab4 100644 --- a/src/prox_operators.jl +++ b/src/prox_operators.jl @@ -28,29 +28,19 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, p.min_eig[idx] = a.m[idx][1] elseif !opt.full_eig_decomp && - p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && - sdp.sq_side > opt.min_size_krylov_eigs && - mod(p.iter, opt.full_eig_freq) > opt.full_eig_len # for full from time to time - @timeit "eigs" begin - eig!(arc_list[idx], a.m[idx], p.target_rank[idx], opt) - - if hasconverged(arc_list[idx]) - fill!(a.m[idx].data, 0.) - for i in 1:p.target_rank[idx] - if unsafe_getvalues(arc_list[idx])[i] > 0. - p.current_rank[idx] += 1 - vec = view(unsafe_getvectors(arc_list[idx]), :, i) - LinearAlgebra.BLAS.gemm!('N', 'T', unsafe_getvalues(arc_list[idx])[i], vec, vec, 1., a.m[idx].data) - end - end - end + p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && + sdp.sq_side > opt.min_size_krylov_eigs && + mod(p.iter, opt.full_eig_freq) > opt.full_eig_len # for full from time to time + @timeit "eigs" if opt.eigsolver == 1 + arpack_eig!(arc_list[idx], a, idx, opt, p) + else + end if hasconverged(arc_list[idx]) - @timeit "get min eig" p.min_eig[idx] = minimum(unsafe_getvalues(arc_list[idx])) + @timeit "get min eig" p.min_eig[idx] = minimum(arpack_getvalues(arc_list[idx])) else @timeit "eigfact" full_eig!(a, idx, opt, p) end - else p.min_eig[idx] = 0. @timeit "eigfact" full_eig!(a, idx, opt, p) @@ -72,6 +62,38 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, return nothing end +function arpack_eig!(solver::EigSolverAlloc, a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing + arpack_eig!(solver, a.m[idx], p.target_rank[idx], opt) + if hasconverged(solver) + fill!(a.m[idx].data, 0.) + for i in 1:p.target_rank[idx] + val = arpack_getvalues(solver)[i] + if val > 0. + p.current_rank[idx] += 1 + vec = view(arpack_getvectors(solver), :, i) + LinearAlgebra.BLAS.gemm!('N', 'T', val, vec, vec, 1., a.m[idx].data) + end + end + end + return nothing +end + +function krylovkit_eig!(solver::EigSolverAlloc, a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing + krylovkit_eig!(solver, a.m[idx], p.target_rank[idx], opt) + if hasconverged(solver) + fill!(a.m[idx].data, 0.) + for i in 1:p.target_rank[idx] + val = krylovkit_getvalues(solver)[i] + if val > 0. + p.current_rank[idx] += 1 + vec = krylovkit_getvector(solver, i) + LinearAlgebra.BLAS.gemm!('N', 'T', val, vec, vec, 1., a.m[idx].data) + end + end + end + return nothing +end + function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing p.current_rank[idx] = 0 fact = eigen!(a.m[idx]) diff --git a/src/structs.jl b/src/structs.jl index ecb7234..f77d262 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -69,22 +69,30 @@ Base.@kwdef mutable struct Options 2: KrylovKit [eigsolve/Lanczos] (DEFAULT) =# eigsolver::Int = 2 + eigsolver_min_lanczos::Int = 25 + eigsolver_resid_seed::Int = 1234 # Arpack arpack_tol::Float64 = 1e-10 #= 0: arpack random [usually faster - NON-DETERMINISTIC - slightly] 1: all ones [???] - 2: julia random uniform (arpack_resid_seed) [medium for DETERMINISTIC] - 3: julia normalized random normal (arpack_resid_seed) [best for DETERMINISTIC] + 2: julia random uniform (eigsolver_resid_seed) [medium for DETERMINISTIC] + 3: julia normalized random normal (eigsolver_resid_seed) [best for DETERMINISTIC] =# arpack_resid_init::Int = 3 - arpack_resid_seed::Int = 1234 - arpack_reset_resid::Bool = false # true for determinism - arpack_max_iter::Int = 10_000 - arpack_min_lanczos::Int = 25 + arpack_reset_resid::Bool = true # true for determinism # larger is more stable to converge and more deterministic + arpack_max_iter::Int = 10_000 # see remark for of dsaupd + + # KrylovKit + krylovkit_reset_resid::Bool = false + krylovkit_resid_init::Int = 3 + krylovkit_tol::Float64 = 1e-12 + krylovkit_max_iter::Int = 100 + krylovkit_eager::Bool = false + krylovkit_verbose::Int = 0 # Reduce rank [warning: heuristics] reduce_rank::Bool = false diff --git a/test/moitest.jl b/test/moitest.jl index 07b5d43..0e98779 100644 --- a/test/moitest.jl +++ b/test/moitest.jl @@ -19,6 +19,10 @@ const optimizer_full = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(full_eig_decomp = true, tol_primal = 1e-4, tol_dual = 1e-4)) const optimizer_print = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(log_freq = 10, log_verbose = true, timer_verbose = true, tol_primal = 1e-4, tol_dual = 1e-4)) +const optimizer_lowacc_arpack = MOIU.CachingOptimizer(cache, + ProxSDP.Optimizer(eigsolver = 1, tol_primal = 1e-3, tol_dual = 1e-3, log_verbose = false)) +const optimizer_lowacc_krylovkit = MOIU.CachingOptimizer(cache, + ProxSDP.Optimizer(eigsolver = 2, tol_primal = 1e-3, tol_dual = 1e-3, log_verbose = false)) const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, infeas_certificates = false) const config_conic = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false, infeas_certificates = false) @@ -483,11 +487,15 @@ end # badly conditioned # moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-1.dat-s"), test = true) moi_sdplib(optimizer_low_acc, joinpath(datapath, "gpp124-2.dat-s"), test = true) + moi_sdplib(optimizer_lowacc_arpack, joinpath(datapath, "gpp124-2.dat-s"), test = true) + moi_sdplib(optimizer_lowacc_krylovkit, joinpath(datapath, "gpp124-2.dat-s"), test = true) # moi_sdplib(optimizer, joinpath(datapath, "gpp124-3.dat-s"), test = true) # moi_sdplib(optimizer, joinpath(datapath, "gpp124-4.dat-s"), test = true) end @testset "MAX CUT" begin moi_sdplib(optimizer_low_acc, joinpath(datapath, "mcp124-1.dat-s"), test = true) + moi_sdplib(optimizer_lowacc_arpack, joinpath(datapath, "mcp124-1.dat-s"), test = true) + moi_sdplib(optimizer_lowacc_krylovkit, joinpath(datapath, "mcp124-1.dat-s"), test = true) # moi_sdplib(optimizer, joinpath(datapath, "mcp124-2.dat-s"), test = true) # moi_sdplib(optimizer, joinpath(datapath, "mcp124-3.dat-s"), test = true) # moi_sdplib(optimizer, joinpath(datapath, "mcp124-4.dat-s"), test = true) From 27da113245a6b6d5216b88dd129211d1e41090fd Mon Sep 17 00:00:00 2001 From: Joaquim Garcia Date: Tue, 22 Sep 2020 13:14:50 -0300 Subject: [PATCH 2/4] up version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5fc119d..8fa846c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ProxSDP" uuid = "65e78d25-6039-50a4-9445-38022e3d2eb3" repo = "https://github.com/mariohsouto/ProxSDP.jl.git" -version = "1.4.1" +version = "1.5.0" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" From ab69b2ec5d3afc4c037863151083709c58e316bb Mon Sep 17 00:00:00 2001 From: Joaquim Garcia Date: Tue, 22 Sep 2020 16:26:34 -0300 Subject: [PATCH 3/4] fix krylov kit and disable logger --- Project.toml | 1 + src/MOI_wrapper.jl | 11 +++++++++++ src/ProxSDP.jl | 1 + src/eigsolver.jl | 12 +++++++----- src/pdhg.jl | 6 +++--- src/printing.jl | 2 +- src/prox_operators.jl | 39 +++++++++++++++++++++++++++------------ src/residuals.jl | 9 +++++++-- src/structs.jl | 3 +++ test/jump_mimo.jl | 2 -- 10 files changed, 61 insertions(+), 25 deletions(-) diff --git a/Project.toml b/Project.toml index 8fa846c..b0d73aa 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "1.5.0" Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 14d016a..9913394 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -561,8 +561,19 @@ function MOI.optimize!(optimizer::Optimizer) # warm = WarmStart() + if options.disable_julia_logger + # disable logger + global_log = Logging.current_logger() + Logging.global_logger(Logging.NullLogger()) + end + sol = @timeit "Main" chambolle_pock(aff, con, options) + if options.disable_julia_logger + # re-enable logger + Logging.global_logger(global_log) + end + #= Unload solution =# diff --git a/src/ProxSDP.jl b/src/ProxSDP.jl index bb3bf92..629a83d 100644 --- a/src/ProxSDP.jl +++ b/src/ProxSDP.jl @@ -9,6 +9,7 @@ module ProxSDP using SparseArrays using LinearAlgebra import Random + import Logging import LinearAlgebra: BlasInt diff --git a/src/eigsolver.jl b/src/eigsolver.jl index c9825a0..bfbd313 100644 --- a/src/eigsolver.jl +++ b/src/eigsolver.jl @@ -761,6 +761,8 @@ function arpack_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Intege end end + arc.converged_eigs = arc.iparam[5] + # Post processing routine (eigenvalues and eigenvectors purification) @timeit "seupd" _seupd!(arc)::Nothing return nothing @@ -771,7 +773,7 @@ function krylovkit_getvalues(arc::EigSolverAlloc) end function krylovkit_getvector(arc::EigSolverAlloc, i) - return arc.v[i] + return arc.vecs[i] end function krylovkit_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T @@ -783,8 +785,7 @@ function krylovkit_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev: return nothing end -function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T - arc.n = n +function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, opt::Options) where T arc.nev = nev if opt.krylovkit_reset_resid eigsolver_update_resid!(arc, opt, opt.krylovkit_resid_init) @@ -794,9 +795,8 @@ function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, ne end function krylovkit_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T - up_ncv = false arc.converged = true - @timeit "update_arc" krylovkit_update!(arc, A, nev, opt, up_ncv)::Nothing + @timeit "update_arc" krylovkit_update!(arc, A, nev, opt)::Nothing @timeit "krylovkit" begin vals, vecs, info = KrylovKit.eigsolve( A, arc.resid, arc.nev, :LR, @@ -809,6 +809,8 @@ function krylovkit_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Int opt.krylovkit_verbose ) ) + arc.vals = vals + arc.vecs = vecs arc.converged_eigs = info.converged if info.converged == 0 arc.converged = false diff --git a/src/pdhg.jl b/src/pdhg.jl index 111c969..025c573 100644 --- a/src/pdhg.jl +++ b/src/pdhg.jl @@ -146,7 +146,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP p.rank_update += 1 if residuals.dual_gap <= opt.tol_primal && residuals.equa_feasibility <= opt.tol_primal - if convergedrank(p, conic_sets, opt) && soc_convergence(a, conic_sets, pair, opt, p) && p.iter > opt.min_iter + if convergedrank(a, p, conic_sets, opt) && soc_convergence(a, conic_sets, pair, opt, p) && p.iter > opt.min_iter p.stop_reason = 1 # Optimal p.stop_reason_string = "Optimal solution found" if opt.log_verbose @@ -160,7 +160,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP if p.update_cont > 0 for (idx, sdp) in enumerate(conic_sets.sdpcone) if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx] - if p.min_eig[idx] > opt.tol_psd + if min_eig(a, idx, p) > opt.tol_psd if opt.rank_increment == 0 p.target_rank[idx] = min(opt.rank_increment_factor * p.target_rank[idx], sdp.sq_side) else @@ -183,7 +183,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP full_rank_flag = false end if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx] - if p.min_eig[idx] > opt.tol_psd + if min_eig(a, idx, p) > opt.tol_psd if opt.rank_increment == 0 p.target_rank[idx] = min(opt.rank_increment_factor * p.target_rank[idx], sdp.sq_side) else diff --git a/src/printing.jl b/src/printing.jl index df851da..24d6e63 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -3,7 +3,7 @@ function print_header_1() println("=======================================================================================") println(" ProxSDP : Proximal Semidefinite Programming Solver ") println(" (c) Mario Souto and Joaquim D. Garcia, 2020 ") - println(" v1.4.1 ") + println(" v1.5.0 ") println("---------------------------------------------------------------------------------------") end diff --git a/src/prox_operators.jl b/src/prox_operators.jl index 67dcab4..abb9d7c 100644 --- a/src/prox_operators.jl +++ b/src/prox_operators.jl @@ -1,5 +1,3 @@ - - function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, opt::Options, p::Params, arc_list, iter::Int64)::Nothing p.min_eig, current_rank, sqrt_2 = zeros(length(cones.sdpcone)), 0, sqrt(2.) @@ -23,10 +21,8 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, p.current_rank[idx] = 0 if sdp.sq_side == 1 - a.m[idx][1] = max(0., a.m[idx][1]) p.min_eig[idx] = a.m[idx][1] - elseif !opt.full_eig_decomp && p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && sdp.sq_side > opt.min_size_krylov_eigs && @@ -34,15 +30,12 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, @timeit "eigs" if opt.eigsolver == 1 arpack_eig!(arc_list[idx], a, idx, opt, p) else - + krylovkit_eig!(arc_list[idx], a, idx, opt, p) end - if hasconverged(arc_list[idx]) - @timeit "get min eig" p.min_eig[idx] = minimum(arpack_getvalues(arc_list[idx])) - else + if !hasconverged(arc_list[idx]) @timeit "eigfact" full_eig!(a, idx, opt, p) end else - p.min_eig[idx] = 0. @timeit "eigfact" full_eig!(a, idx, opt, p) end end @@ -66,9 +59,14 @@ function arpack_eig!(solver::EigSolverAlloc, a::AuxiliaryData, idx::Int, opt::Op arpack_eig!(solver, a.m[idx], p.target_rank[idx], opt) if hasconverged(solver) fill!(a.m[idx].data, 0.) + # TODO: how to measure this when convergen eigs is less than target? + p.min_eig[idx] = minimum(arpack_getvalues(solver)) + # if solver.converged_eigs < p.target_rank[idx] && p.min_eig[idx] > 0.0 + # p.min_eig[idx] = -Inf + # end for i in 1:p.target_rank[idx] val = arpack_getvalues(solver)[i] - if val > 0. + if val > 0. p.current_rank[idx] += 1 vec = view(arpack_getvectors(solver), :, i) LinearAlgebra.BLAS.gemm!('N', 'T', val, vec, vec, 1., a.m[idx].data) @@ -82,9 +80,15 @@ function krylovkit_eig!(solver::EigSolverAlloc, a::AuxiliaryData, idx::Int, opt: krylovkit_eig!(solver, a.m[idx], p.target_rank[idx], opt) if hasconverged(solver) fill!(a.m[idx].data, 0.) - for i in 1:p.target_rank[idx] + # TODO: how to measure this when convergen eigs is less than target? + # ? min_eig is just checking if the rank projection is going far enough to reach zeros and negatives + p.min_eig[idx] = minimum(krylovkit_getvalues(solver)) + # if solver.converged_eigs < p.target_rank[idx] #&& p.min_eig[idx] > 0.0 + # p.min_eig[idx] = -Inf + # end + for i in 1:min(p.target_rank[idx], solver.converged_eigs) val = krylovkit_getvalues(solver)[i] - if val > 0. + if val > 0. p.current_rank[idx] += 1 vec = krylovkit_getvector(solver, i) LinearAlgebra.BLAS.gemm!('N', 'T', val, vec, vec, 1., a.m[idx].data) @@ -97,6 +101,7 @@ end function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing p.current_rank[idx] = 0 fact = eigen!(a.m[idx]) + p.min_eig[idx] = 0.0 #minimum(fact.values) fill!(a.m[idx].data, 0.) for i in 1:length(fact.values) if fact.values[i] > 0. @@ -110,6 +115,16 @@ function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing return nothing end +function min_eig(a::AuxiliaryData, idx::Int, p::Params) + # if p.min_eig[idx] == -Inf + # @timeit "bad min eig" begin + # fact = eigen!(a.m[idx]) + # p.min_eig[idx] = minimum(fact.values) + # end + # end + return p.min_eig[idx] +end + function soc_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, opt::Options, p::Params)::Nothing for (idx, soc) in enumerate(cones.socone) soc_projection!(a.soc_v[idx], a.soc_s[idx]) diff --git a/src/residuals.jl b/src/residuals.jl index 70538fa..a3409a8 100644 --- a/src/residuals.jl +++ b/src/residuals.jl @@ -79,9 +79,14 @@ function soc_gap(v::ViewVector, s::ViewScalar) return norm(v, 2) - s[] end -function convergedrank(p::Params, cones::ConicSets, opt::Options)::Bool +function convergedrank(a, p::Params, cones::ConicSets, opt::Options)::Bool for (idx, sdp) in enumerate(cones.sdpcone) - if !(p.min_eig[idx] < opt.tol_psd || p.target_rank[idx] > opt.max_target_rank_krylov_eigs || sdp.sq_side < opt.min_size_krylov_eigs) + if !( + sdp.sq_side < opt.min_size_krylov_eigs || + p.target_rank[idx] > opt.max_target_rank_krylov_eigs || + min_eig(a, idx, p) < opt.tol_psd + ) + # @show min_eig(a, idx, p), -opt.tol_psd return false end end diff --git a/src/structs.jl b/src/structs.jl index f77d262..450b4b8 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -20,6 +20,7 @@ Base.@kwdef mutable struct Options log_freq::Int = 100 timer_verbose::Bool = false timer_file::Bool = false + disable_julia_logger = true # time options time_limit::Float64 = 3600_00. #100 hours @@ -73,6 +74,8 @@ Base.@kwdef mutable struct Options eigsolver_resid_seed::Int = 1234 # Arpack + # note that Arpack is Non-deterministic + # (https://github.com/mariohsouto/ProxSDP.jl/issues/69) arpack_tol::Float64 = 1e-10 #= 0: arpack random [usually faster - NON-DETERMINISTIC - slightly] diff --git a/test/jump_mimo.jl b/test/jump_mimo.jl index cc62d8f..97be138 100644 --- a/test/jump_mimo.jl +++ b/test/jump_mimo.jl @@ -30,8 +30,6 @@ function jump_mimo(solver, seed, n; verbose = false, test = false) objval = objective_value(model) stime = MOI.get(model, MOI.SolveTime()) - # @show tp = typeof(model.moi_backend.optimizer.model.optimizer) - # @show fieldnames(tp) rank = -1 try @show rank = model.moi_backend.optimizer.model.optimizer.sol.final_rank From 91fbbb3e5da86c1486ca3a2a5d8c4848f26d83df Mon Sep 17 00:00:00 2001 From: Joaquim Garcia Date: Tue, 22 Sep 2020 17:19:30 -0300 Subject: [PATCH 4/4] remove old code --- src/structs.jl | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/structs.jl b/src/structs.jl index 450b4b8..027eba2 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -116,29 +116,6 @@ Base.@kwdef mutable struct Options approx_norm::Bool = true end -function Options(args) - options = Options() - parse_args!(options, args) - return options -end - -function parse_args!(options, args) - for i in args - parse_arg!(options, i) - end - return nothing -end - -function parse_arg!(options::Options, arg) - fields = fieldnames(Options) - name = arg[1] - value = arg[2] - if name in fields - setfield!(options, name, value) - end - return nothing -end - mutable struct AffineSets n::Int # Size of primal variables p::Int # Number of linear equalities