Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Multiple centrality corrections #75

Merged
merged 3 commits into from
Dec 16, 2020
Merged

Multiple centrality corrections #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
77a9899
Re-organize MPC code
mtanneau Dec 15, 2020
66e428b
Add Gondzio corrections to MPC; set default max corrections to 3.
mtanneau Dec 15, 2020
5c3903d
Cosmetic changes
mtanneau Dec 15, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions src/IPM/HSD/step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
α = max_step_length(pt, Δ)
γ = (one(T) - α)^2 * min(one(T) - α, params.GammaMin)
η = one(T) - γ

# Mehrotra corrector
@timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0,
# Right-hand side of Newton system
Expand All @@ -104,8 +104,8 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
ncor += 1

# Compute extra-corrector
αc = compute_higher_corrector!(Δc,
hsd, γ,
αc = compute_higher_corrector!(Δc,
hsd, γ,
hx, hy, h0,
Δ, α_, params.CentralityOutlierThreshold
)
Expand All @@ -121,7 +121,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
Δ.κ = Δc.κ
α = αc
end

if αc < T(11 // 10) * α_
break
end
Expand All @@ -130,7 +130,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
# # not enough improvement, step correcting
# break
# end

end

# Update current iterate
Expand All @@ -144,7 +144,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
pt.τ += α * Δ.τ
pt.κ += α * Δ.κ
update_mu!(pt)

return nothing
end

Expand Down Expand Up @@ -213,11 +213,11 @@ function solve_newton_system!(Δ::Point{T, Tv},

# II. Recover Δτ, Δx, Δy
# Compute Δτ
@timeit hsd.timer "ξg_" ξg_ = (ξg + ξtk / pt.τ
@timeit hsd.timer "ξg_" ξg_ = (ξg + ξtk / pt.τ
- dot((ξxzl ./ pt.xl) .* dat.lflag, dat.l .* dat.lflag) # l'(Xl)^-1 * ξxzl
+ dot((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag)
- dot(((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag)
- dot(((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag) #
+ dot((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag)
- dot(((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag)
- dot(((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag) #
)

@timeit hsd.timer "Δτ" Δ.τ = (
Expand Down Expand Up @@ -276,7 +276,7 @@ function max_step_length(x::Vector{T}, dx::Vector{T}) where{T}
n == size(dx, 1) || throw(DimensionMismatch())
a = T(Inf)

#=@inbounds=# for i in Base.OneTo(n)
@inbounds for i in 1:n
if dx[i] < zero(T)
if (-x[i] / dx[i]) < a
a = (-x[i] / dx[i])
Expand All @@ -299,7 +299,7 @@ function max_step_length(pt::Point{T, Tv}, δ::Point{T, Tv}) where{T, Tv<:Abstra

at = δ.τ < zero(T) ? (-pt.τ / δ.τ) : oneunit(T)
ak = δ.κ < zero(T) ? (-pt.κ / δ.κ) : oneunit(T)

α = min(one(T), axl, axu, azl, azu, at, ak)

return α
Expand All @@ -316,7 +316,7 @@ Requires the solution of one Newton system.
# Arguments
- `Δc`: Corrected search direction, modified in-place
- `hsd`: The HSD optimizer
- `γ`:
- `γ`:
- `hx, hy, h0`: Terms obtained from the preliminary augmented system solve
- `Δ`: Current predictor direction
- `α`: Maximum step length in predictor direction
Expand Down Expand Up @@ -398,4 +398,4 @@ function compute_higher_corrector!(Δc::Point{T, Tv},
# Compute corrected step-length
αc = max_step_length(pt, Δc)
return αc
end
end
52 changes: 32 additions & 20 deletions src/IPM/MPC/MPC.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
primal_status::SolutionStatus # Primal solution status
dual_status::SolutionStatus # Dual solution status


primal_objective::T # Primal bound: c'x
dual_objective::T # Dual bound: b'y + l' zl - u'zu

Expand All @@ -28,27 +27,34 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
pt::Point{T, Tv} # Current primal-dual iterate
res::Residuals{T, Tv} # Residuals at current iterate

Δ::Point{T, Tv} # Predictor direction
Δ::Point{T, Tv} # Predictor
Δc::Point{T, Tv} # Corrector

# Step sizes
αp::T
αd::T

# Newton system RHS

ξp::Tv
ξl::Tv
ξu::Tv
ξd::Tv
ξxzl::Tv
ξxzu::Tv

# KKT solver
kkt::Tk
regP::Tv # Primal regularization
regD::Tv # Dual regularization

function MPC(
dat::IPMData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T}
) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}}

m, n = dat.nrow, dat.ncol
p = sum(dat.lflag) + sum(dat.uflag)

# Allocate some memory
# Working memory
pt = Point{T, Tv}(m, n, p, hflag=false)
res = Residuals(
tzeros(Tv, m), tzeros(Tv, n), tzeros(Tv, n),
Expand All @@ -58,6 +64,14 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
Δ = Point{T, Tv}(m, n, p, hflag=false)
Δc = Point{T, Tv}(m, n, p, hflag=false)

# Newton RHS
ξp = tzeros(Tv, m)
ξl = tzeros(Tv, n)
ξu = tzeros(Tv, n)
ξd = tzeros(Tv, n)
ξxzl = tzeros(Tv, n)
ξxzu = tzeros(Tv, n)

# Initial regularizations
regP = tones(Tv, n)
regD = tones(Tv, m)
Expand All @@ -70,6 +84,7 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T}
T(Inf), T(-Inf),
TimerOutput(),
pt, res, Δ, Δc, zero(T), zero(T),
ξp, ξl, ξu, ξd, ξxzl, ξxzu,
kkt, regP, regD
)
end
Expand Down Expand Up @@ -160,15 +175,15 @@ function update_solver_status!(mpc::MPC{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
else
mpc.dual_status = Sln_Unknown
end

# Check for optimal solution
if ρp <= ϵp && ρd <= ϵd && ρg <= ϵg
mpc.primal_status = Sln_Optimal
mpc.dual_status = Sln_Optimal
mpc.solver_status = Trm_Optimal
return nothing
end

# TODO: Primal/Dual infeasibility detection
# Check for infeasibility certificates
if max(
Expand Down Expand Up @@ -257,12 +272,12 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
if params.OutputLevel > 0
# Display log
@printf "%4d" mpc.niter

# Objectives
ϵ = dat.objsense ? one(T) : -one(T)
@printf " %+14.7e" ϵ * mpc.primal_objective
@printf " %+14.7e" ϵ * mpc.dual_objective

# Residuals
@printf " %8.2e" max(mpc.res.rp_nrm, mpc.res.rl_nrm, mpc.res.ru_nrm)
@printf " %8.2e" mpc.res.rd_nrm
Expand Down Expand Up @@ -295,14 +310,14 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
|| mpc.solver_status == Trm_DualInfeasible
)
break
elseif mpc.niter >= params.IterationsLimit
elseif mpc.niter >= params.IterationsLimit
mpc.solver_status = Trm_IterationLimit
break
elseif ttot >= params.TimeLimit
mpc.solver_status = Trm_TimeLimit
break
end


# TODO: step
# For now, include the factorization in the step function
Expand All @@ -314,7 +329,7 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}
if isa(err, PosDefException) || isa(err, SingularException)
# Numerical trouble while computing the factorization
mpc.solver_status = Trm_NumericalProblem

elseif isa(err, OutOfMemoryError)
# Out of memory
mpc.solver_status = Trm_MemoryLimit
Expand All @@ -328,16 +343,13 @@ function ipm_optimize!(mpc::MPC{T}, params::IPMOptions{T}) where{T}

break
end

mpc.niter += 1

end

# TODO: print message based on termination status
params.OutputLevel > 0 && println("Solver exited with status $((mpc.solver_status))")

return nothing

end

function compute_starting_point(mpc::MPC{T}) where{T}
Expand All @@ -351,7 +363,7 @@ function compute_starting_point(mpc::MPC{T}) where{T}
# Get initial iterate
KKT.solve!(zeros(T, n), pt.y, mpc.kkt, false .* mpc.dat.b, mpc.dat.c) # For y
KKT.solve!(pt.x, zeros(T, m), mpc.kkt, mpc.dat.b, false .* mpc.dat.c) # For x

# I. Recover positive primal-dual coordinates
δx = one(T) + max(
zero(T),
Expand All @@ -375,11 +387,11 @@ function compute_starting_point(mpc::MPC{T}) where{T}
=#
pt.zl .= ( z ./ (dat.lflag + dat.uflag)) .* dat.lflag
pt.zu .= (-z ./ (dat.lflag + dat.uflag)) .* dat.uflag

δz = one(T) + max(zero(T), (-3 // 2) * minimum(pt.zl), (-3 // 2) * minimum(pt.zu))
pt.zl[dat.lflag] .+= δz
pt.zu[dat.uflag] .+= δz

mpc.pt.τ = one(T)
mpc.pt.κ = zero(T)

Expand All @@ -397,4 +409,4 @@ function compute_starting_point(mpc::MPC{T}) where{T}
update_mu!(mpc.pt)

return nothing
end
end
Loading