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

WIP: constrained optimization #50

Closed
wants to merge 29 commits into from
Closed

WIP: constrained optimization #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
680edc7
Add support for constraints and a interior-point methods
timholy Mar 8, 2014
809be96
Add constrained tests, and don't load problems unless running tests
timholy Mar 11, 2014
83a0bc3
Implement linear constraints
timholy Mar 14, 2014
f832562
Improve initial value of barrier coefficient
timholy Apr 2, 2014
5fabfe8
Update for deprecation of infs
timholy Apr 4, 2014
8f56008
Check type of function value, and fix type of initial t in interior-p…
timholy Apr 4, 2014
de95ece
Modify some asserts that were causing trouble for non-convex domains
timholy Apr 5, 2014
72f18ea
nelder-mead: check that at least one starting point has finite value
timholy May 23, 2014
048394a
Fix call to norm on high-dimensional objects
timholy Sep 11, 2014
4d536b7
fminbox: remove unused tol argument and add mu0
timholy Sep 11, 2014
7f48aeb
fixed some deprecation warnings thrown by julia 0.4. Used Compat pac…
Cody-G Feb 4, 2015
7ed8c2a
Merge pull request #99 from Cody-G/fix_deprecations
timholy Feb 14, 2015
0f35292
Fix ConstraintsNone parameter problem in step!
timholy Apr 9, 2015
03c5e04
Fixes for 0.4 deprecations
timholy Apr 9, 2015
4d1cdf4
Bump the default linesearchmax
timholy Apr 14, 2015
c0a3a3f
Julia 0.4 update
timholy Sep 23, 2015
5c1e65e
More 0.4 fixes
timholy Sep 23, 2015
c2b7cda
Guarantee that the initial direction is a descent direction for objec…
timholy Oct 21, 2015
92a1c08
Pass constraints through to minimizers
timholy Oct 21, 2015
ce8df0a
Automatically try a restart if we think we've converged
timholy Oct 22, 2015
28ce676
Ensure that the search direction is a descent direction
timholy Dec 1, 2015
dc4a4c9
interior_newton: increment t if the step converged
timholy Dec 18, 2015
67c9996
Fix returned function value when there are 0 parameters
timholy Dec 18, 2015
be3324a
Improvements to positive-definite newton step
timholy Dec 18, 2015
afccefc
Use PositiveFactorizations for the newton step
timholy Jan 26, 2016
dad89d2
Change convergence criteria in interior_newton to monitor solution, n…
timholy Jan 27, 2016
b66e6b6
dot->vecdot
timholy Jan 27, 2016
db76ca0
interior_newton: pass options, bail if encounter numeric instability
timholy Jan 27, 2016
2943c3f
interior_newton: also check for finiteness
timholy Jan 27, 2016
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
4 changes: 3 additions & 1 deletion REQUIRE
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
julia 0.2-
julia 0.4
Calculus
DualNumbers
Compat
PositiveFactorizations
28 changes: 21 additions & 7 deletions src/Optim.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,20 +1,28 @@
isdefined(Base, :__precompile__) && __precompile__()

module Optim
using Calculus
using Calculus, PositiveFactorizations
using Compat

import Base.dot,
Base.length,
import Base.length,
Base.push!,
Base.show,
Base.getindex,
Base.setindex!

export optimize,
interior,
linlsq,
DifferentiableFunction,
TwiceDifferentiableFunction
TwiceDifferentiableFunction,
ConstraintsBox

# Types
include("types.jl")

# Types for constrained optimization
include("constraints.jl")

# Automatic differentiation utilities
include("autodiff.jl")

Expand Down Expand Up @@ -65,11 +73,17 @@ module Optim
include("golden_section.jl")
include("brent.jl")

# Constrained optimization algorithms
include("interior.jl")

# End-User Facing Wrapper Functions
include("optimize.jl")

# Examples for testing
include(joinpath("problems", "unconstrained.jl"))

cgdescent(args...) = error("API has changed. Please use cg.")

# Tests
const basedir = dirname(Base.source_path())
const testpaths = [joinpath(basedir, "problems", "unconstrained.jl"),
joinpath(basedir, "problems", "constrained.jl")]

end
2 changes: 1 addition & 1 deletion src/accelerated_gradient_descent.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ function accelerated_gradient_descent{T}(d::DifferentiableFunction,
return MultivariateOptimizationResults("Accelerated Gradient Descent",
initial_x,
x,
float64(f_x),
@compat(Float64(f_x)),
iteration,
iteration == iterations,
x_converged,
Expand Down
8 changes: 4 additions & 4 deletions src/bfgs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ macro bfgstrace()
end
end

function bfgs{T}(d::Union(DifferentiableFunction,
TwiceDifferentiableFunction),
function bfgs{T}(d::Union{DifferentiableFunction,
TwiceDifferentiableFunction},
initial_x::Vector{T};
initial_invH::Matrix = eye(length(initial_x)),
xtol::Real = 1e-32,
Expand Down Expand Up @@ -97,7 +97,7 @@ function bfgs{T}(d::Union(DifferentiableFunction,
# Increment the number of steps we've had to perform
iteration += 1

# Set the search direction
# Set the search direction
# Search direction is the negative gradient divided by the approximate Hessian
A_mul_B!(s, invH, gr)
for i in 1:n
Expand Down Expand Up @@ -178,7 +178,7 @@ function bfgs{T}(d::Union(DifferentiableFunction,
return MultivariateOptimizationResults("BFGS",
initial_x,
x,
float64(f_x),
@compat(Float64(f_x)),
iteration,
iteration == iterations,
x_converged,
Expand Down
14 changes: 7 additions & 7 deletions src/brent.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ macro brenttrace()
end
end

function brent{T <: FloatingPoint}(f::Function, x_lower::T, x_upper::T;
function brent{T <: AbstractFloat}(f::Function, x_lower::T, x_upper::T;
rel_tol::T = sqrt(eps(T)),
abs_tol::T = eps(T),
iterations::Integer = 1_000,
Expand All @@ -33,13 +33,13 @@ function brent{T <: FloatingPoint}(f::Function, x_lower::T, x_upper::T;
# Save for later
initial_lower = x_lower
initial_upper = x_upper

const golden_ratio::T = 0.5 * (3.0 - sqrt(5.0))

x_minimum = x_lower + golden_ratio*(x_upper-x_lower)
f_minimum = f(x_minimum)
f_calls = 1 # Number of calls to f

step = zero(T)
step_old = zero(T)

Expand All @@ -48,7 +48,7 @@ function brent{T <: FloatingPoint}(f::Function, x_lower::T, x_upper::T;

f_minimum_old = f_minimum
f_minimum_old_old = f_minimum

it = 0
converged = false

Expand All @@ -72,12 +72,12 @@ function brent{T <: FloatingPoint}(f::Function, x_lower::T, x_upper::T;
end

it += 1

if abs(step_old) > tolx
# Compute parabola interpolation
# x_minimum + p/q is the optimum of the parabola
# Also, q is guaranteed to be positive

r = (x_minimum - x_minimum_old) * (f_minimum - f_minimum_old_old)
q = (x_minimum - x_minimum_old_old) * (f_minimum - f_minimum_old)
p = (x_minimum - x_minimum_old_old) * q - (x_minimum - x_minimum_old) * r
Expand Down Expand Up @@ -150,7 +150,7 @@ function brent{T <: FloatingPoint}(f::Function, x_lower::T, x_upper::T;
initial_lower,
initial_upper,
x_minimum,
float64(f_minimum),
@compat(Float64(f_minimum)),
it,
converged,
rel_tol,
Expand Down
49 changes: 29 additions & 20 deletions src/cg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# Preconditioners
# * Empty preconditioner
cg_precondfwd(out::Array, P::Nothing, A::Array) = copy!(out, A)
cg_precondfwddot(A::Array, P::Nothing, B::Array) = dot(A, B)
cg_precondinvdot(A::Array, P::Nothing, B::Array) = dot(A, B)
cg_precondfwd(out::Array, P::Void, A::Array) = copy!(out, A)
cg_precondfwddot(A::Array, P::Void, B::Array) = vecdot(A, B)
cg_precondinvdot(A::Array, P::Void, B::Array) = vecdot(A, B)

# Diagonal preconditioner
function cg_precondfwd(out::Array, p::Vector, A::Array)
Expand Down Expand Up @@ -89,7 +89,7 @@ macro cgtrace()
dt["g(x)"] = copy(gr)
dt["Current step size"] = alpha
end
grnorm = norm(gr, Inf)
grnorm = norm(gr[:], Inf)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you use vecnorm(gr, Inf) here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah.

Thanks for the comment!

update!(tr,
iteration,
f_x,
Expand All @@ -101,9 +101,11 @@ macro cgtrace()
end
end

function cg{T}(df::Union(DifferentiableFunction,
TwiceDifferentiableFunction),
function cg{T}(df::Union{DifferentiableFunction,
TwiceDifferentiableFunction},
initial_x::Array{T};
constraints::AbstractConstraints = ConstraintsNone(),
interior::Bool = false,
xtol::Real = convert(T,1e-32),
ftol::Real = convert(T,1e-8),
grtol::Real = convert(T,1e-8),
Expand All @@ -117,7 +119,9 @@ function cg{T}(df::Union(DifferentiableFunction,
precondprep::Function = (P, x) -> nothing)

# Maintain current state in x and previous state in x_previous
x, x_previous = copy(initial_x), copy(initial_x)
x = copy(initial_x)
project!(x, constraints)
x_previous = copy(x)

# Count the total number of iterations
iteration = 0
Expand Down Expand Up @@ -146,7 +150,7 @@ function cg{T}(df::Union(DifferentiableFunction,
# Store f(x) in f_x
f_x = df.fg!(x, gr)
@assert typeof(f_x) == T
f_x_previous = nan(T)
f_x_previous = convert(T,NaN)
f_calls, g_calls = f_calls + 1, g_calls + 1
copy!(gr_previous, gr)

Expand Down Expand Up @@ -191,12 +195,12 @@ function cg{T}(df::Union(DifferentiableFunction,
iteration += 1

# Reset the search direction if it becomes corrupted
dphi0 = dot(gr, s)
dphi0 = vecdot(gr, s)
if dphi0 >= 0
for i in 1:n
@inbounds s[i] = -gr[i]
end
dphi0 = dot(gr, s)
dphi0 = vecdot(gr, s)
if dphi0 < 0
break
end
Expand All @@ -208,23 +212,28 @@ function cg{T}(df::Union(DifferentiableFunction,
@assert typeof(dphi0) == T
push!(lsr, zero(T), f_x, dphi0)

alphamax = interior ? toedge(x, s, constraints) : convert(T,Inf)

# Pick the initial step size (HZ #I1-I2)
alpha, mayterminate, f_update, g_update =
alphatry(alpha, df, x, s, x_ls, gr_ls, lsr)
alphatry(alpha, df, x, s, x_ls, gr_ls, lsr, constraints, alphamax)
f_calls, g_calls = f_calls + f_update, g_calls + g_update

if alpha == zero(T)
x_converged = true
break
end

# Determine the distance of movement along the search line
alpha, f_update, g_update =
linesearch!(df, x, s, x_ls, gr_ls, lsr, alpha, mayterminate)
linesearch!(df, x, s, x_ls, gr_ls, lsr, alpha, mayterminate, constraints, alphamax)
f_calls, g_calls = f_calls + f_update, g_calls + g_update

# Maintain a record of previous position
copy!(x_previous, x)

# Update current position
for i in 1:n
@inbounds x[i] = x[i] + alpha * s[i]
end
step!(x, x, s, alpha, constraints)

# Maintain a record of the previous gradient
copy!(gr_previous, gr)
Expand Down Expand Up @@ -254,14 +263,14 @@ function cg{T}(df::Union(DifferentiableFunction,
# Calculate the beta factor (HZ2012)
precondprep(P, x)
dPd = cg_precondinvdot(s, P, s)
etak::T = eta * dot(s, gr_previous) / dPd
etak::T = eta * vecdot(s, gr_previous) / dPd
for i in 1:n
@inbounds y[i] = gr[i] - gr_previous[i]
end
ydots = dot(y, s)
ydots = vecdot(y, s)
cg_precondfwd(pgr, P, gr)
betak = (dot(y, pgr) - cg_precondfwddot(y, P, y) *
dot(gr, s) / ydots) / ydots
betak = (vecdot(y, pgr) - cg_precondfwddot(y, P, y) *
vecdot(gr, s) / ydots) / ydots
beta = max(betak, etak)
for i in 1:n
@inbounds s[i] = beta * s[i] - pgr[i]
Expand All @@ -273,7 +282,7 @@ function cg{T}(df::Union(DifferentiableFunction,
return MultivariateOptimizationResults("Conjugate Gradient",
initial_x,
x,
float64(f_x),
@compat(Float64(f_x)),
iteration,
iteration == iterations,
x_converged,
Expand Down
Loading