Skip to content

Commit

Permalink
fix newton step and add consts
Browse files Browse the repository at this point in the history
  • Loading branch information
@joaquimg
joaquimg committed Jan 12, 2025
1 parent be87589 commit 1cd433e
Showing 1 changed file with 43 additions and 20 deletions.
63 changes: 43 additions & 20 deletions src/projections.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,11 @@ function projection_on_set(d::DefaultDistance, v::AbstractVector{T}, ::MOI.DualE
return SVector{3,T}(v[1] + p[1], v[2] + p[2], v[3] + p[3])
end

const DEFAULT_POWER_CONE_MAX_ITERS_NEWTON = 100
const DEFAULT_POWER_CONE_MAX_ITERS_BISSECTION = 1_000
const DEFAULT_POWER_CONE_TOL_CONV = 1e-7
const DEFAULT_POWER_CONE_TOL_IN_CONE = 1e-10

"""
projection_on_set(::DefaultDistance, v::AbstractVector{T}, ::MOI.PowerCone; [max_iters_newton = 3, max_iters_bisection = 1000]) where {T}
Expand All @@ -401,8 +406,8 @@ function projection_on_set(
::DefaultDistance,
v::AbstractVector{T},
s::MOI.PowerCone;
max_iters_newton = 3,
max_iters_bisection = 1000,
max_iters_newton = DEFAULT_POWER_CONE_MAX_ITERS_NEWTON,
max_iters_bisection = DEFAULT_POWER_CONE_MAX_ITERS_BISSECTION,
) where {T}
_check_dimension(v, s)

Expand All @@ -413,25 +418,25 @@ function projection_on_set(
# if in polar cone Ko = -K*
return zeros(T, 3)
end
if abs(v[3]) <= 1e-10
if abs(v[3]) <= DEFAULT_POWER_CONE_TOL_IN_CONE
return [max(v[1],0), max(v[2],0), 0]
end

_, proj4 = _solve_system_power_cone(
v,
s,
max_iters_newton = max_iters_newton,
max_iters_bisection=max_iters_bisection,
max_iters_bisection = max_iters_bisection,
)
return proj4
end

function _in_pow_cone(v::AbstractVector{T}, cone::MOI.PowerCone; tol=1e-10) where {T}
function _in_pow_cone(v::AbstractVector{T}, cone::MOI.PowerCone; tol=DEFAULT_POWER_CONE_TOL_IN_CONE) where {T}
α = cone.exponent
return v[1] >= 0 && v[2] >= 0 && tol + v[1]^α * v[2]^(1-α) >= abs(v[3])
end

function _in_pow_cone(v::AbstractVector{T}, cone::MOI.DualPowerCone; tol=1e-10) where {T}
function _in_pow_cone(v::AbstractVector{T}, cone::MOI.DualPowerCone; tol=DEFAULT_POWER_CONE_TOL_IN_CONE) where {T}
α = cone.exponent
return (
v[1] >= 0 && v[2] >= 0 && tol + (v[1])^α * (v[2])^(1-α) >= α^α * (1-α)^(1-α) * abs(v[3])
Expand Down Expand Up @@ -468,7 +473,12 @@ function _power_cone_system(r, x, y, z, α)
)
end

function _solve_system_power_cone_bisection(v::AbstractVector{T}, s::MOI.PowerCone; max_iters=10_000, tol=1e-5) where {T}
function _solve_system_power_cone_bisection(
v::AbstractVector{T},
s::MOI.PowerCone;
max_iters=DEFAULT_POWER_CONE_MAX_ITERS_BISSECTION,
tol=DEFAULT_POWER_CONE_TOL_CONV
) where {T}
x, y, z = v
α = s.exponent
# Φ is positive for r = 0 and negative for r = |z| so we can just bisect
Expand All @@ -492,7 +502,12 @@ function _solve_system_power_cone_bisection(v::AbstractVector{T}, s::MOI.PowerCo
return r
end

function _solve_system_power_cone_newton(v::AbstractVector{T}, s::MOI.PowerCone; max_iters=10_000, tol=1e-5) where {T}
function _solve_system_power_cone_newton(
v::AbstractVector{T},
s::MOI.PowerCone;
max_iters=DEFAULT_POWER_CONE_MAX_ITERS_NEWTON,
tol=DEFAULT_POWER_CONE_TOL_CONV
) where {T}
x, y, z = v
α = s.exponent
# Solve with Newton method
Expand All @@ -516,7 +531,7 @@ function _solve_system_power_cone_newton(v::AbstractVector{T}, s::MOI.PowerCone;
dΦ = dΦ_dr(r, Φ, px, py, dpx, dpy)

# Newton step, bounded to interval
r = min(max(r + Φ/dΦ, 0), abs(z))
r = min(max(r - Φ/dΦ, 0), abs(z))
end
return r
end
Expand All @@ -536,22 +551,25 @@ References:
function _solve_system_power_cone(
v::AbstractVector{T},
s::MOI.PowerCone;
max_iters_newton = 3,
max_iters_bisection = 1000,
tol=1e-5,
max_iters_newton = DEFAULT_POWER_CONE_MAX_ITERS_NEWTON,
max_iters_bisection = DEFAULT_POWER_CONE_MAX_ITERS_BISSECTION,
tol=DEFAULT_POWER_CONE_TOL_CONV,
) where {T}
x, y, z = v
α = s.exponent
# We first try to quickly get an accurate solution with Newton:
r = _solve_system_power_cone_newton(v, s; max_iters = max_iters_newton, tol)
r_newton = _solve_system_power_cone_newton(v, s; max_iters = max_iters_newton, tol)
Φ_newton = !isnan(r_newton) ? _power_cone_system(r_newton, x, y, z, α) : NaN
# When the optimal solution `r` is close to the boundaries `0` and `|z|`,
# Newton has a tendency to overshoot and hit the boundary. Once the boundary
# is hit, it will compute `NaN` so we'll need bisection to get an answer.
if isnan(r) || abs(_power_cone_system(r, x, y, z, α)) > tol
r = _solve_system_power_cone_bisection(v, s; max_iters = max_iters_bisection, tol)
println(r)
Φ = _power_cone_system(r, x, y, z, α)
if abs(Φ) > tol

Φ_bissection = Inf
if isnan(r_newton) || abs(Φ_newton) > tol
r_bissection = _solve_system_power_cone_bisection(v, s; max_iters = max_iters_bisection, tol)
println(r_bissection)
Φ_bissection = _power_cone_system(r_bissection, x, y, z, α)
if abs(Φ_bissection) > tol
# This happens for instance for
# `_solve_system_power_cone([-10, 10, 1e-3], MOI.PowerCone(0.15))`
# The value of `r` is found by the bisection to be between
Expand All @@ -560,9 +578,14 @@ function _solve_system_power_cone(
# 0.001
# It's not possible to do better for `Float64` but the tolerance requested
# by the user is not satisfied so we still warn
@warn("Error `$(abs(Φ)) > $tol` after maximum iterations hit for projection of $v onto $s")
@warn("Error `$(abs(Φ_bissection)) > $tol` after maximum iterations hit for projection of $v onto $s")
end
end
r = if isnan(r_newton) || abs(Φ_newton) > abs(Φ_bissection)
r_bissection
else
r_newton
end
px = _power_cone_system_factor(x, α, z, r)
py = _power_cone_system_factor(y, 1 - α, z, r)
return r, [px / T(2), py / T(2), sign(z) * r]
Expand Down Expand Up @@ -842,7 +865,7 @@ function projection_gradient_on_set(::DefaultDistance, v::AbstractVector{T}, s::
# if in polar cone Ko = -K*
return zeros(T, 3, 3)
end
if abs(v[3]) <= 1e-10
if abs(v[3]) <= DEFAULT_POWER_CONE_TOL_IN_CONE
return _pow_cone_∇proj_case_3(v, s)
end

Expand Down

0 comments on commit 1cd433e

Please sign in to comment.