All of the examples can be found in Jupyter notebook form here.

Convex Optimization in Julia

Madeleine Udell | ISMP 2015

Convex.jl team

• Convex.jl: Madeleine Udell, Karanveer Mohan, David Zeng, Jenny Hong

Collaborators/Inspiration:

• CVX: Michael Grant, Stephen Boyd
• CVXPY: Steven Diamond, Eric Chu, Stephen Boyd
• JuliaOpt: Miles Lubin, Iain Dunning, Joey Huchette

# initial package installation

# Make the Convex.jl module available
using Convex, SparseArrays, LinearAlgebra
using SCS # first order splitting conic solver [O'Donoghue et al., 2014]

# Generate random problem data
m = 50;  n = 100
A = randn(m, n)
x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector
noise = .1*randn(m)    # gaussian noise
b = A*x♮ + noise      # noisy linear observations

# Create a (column vector) variable of size n.
x = Variable(n)

# nonnegative elastic net with regularization
λ = 1
μ = 1
problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)

# Solve the problem by calling solve!
solve!(problem, () -> SCS.Optimizer(verbose=0))

println("problem status is ", problem.status) # :Optimal, :Infeasible, :Unbounded etc.
println("optimal value is ", problem.optval)

problem status is OPTIMAL
optimal value is 32.81676511626289

using Interact, Plots
# Interact.WebIO.install_jupyter_nbextension() # might be helpful if you see `WebIO` warnings in Jupyter
@manipulate throttle=.1 for λ=0:.1:5, μ=0:.1:5
    global A
    problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)
    solve!(problem, () -> SCS.Optimizer(verbose=0))
    histogram(evaluate(x), xlims=(0,3.5), label="x")
end

Quick convex prototyping

Variables

# Scalar variable
x = Variable()

Variable
size: (1, 1)
sign: real
vexity: affine
id: 806…655

# (Column) vector variable
y = Variable(4)

Variable
size: (4, 1)
sign: real
vexity: affine
id: 661…933

# Matrix variable
Z = Variable(4, 4)

Variable
size: (4, 4)
sign: real
vexity: affine
id: 715…088

Expressions

Convex.jl allows you to use a wide variety of functions on variables and on expressions to form new expressions.

x + 2x

+ (affine; real)
├─ real variable (id: 806…655)
└─ * (affine; real)
   ├─ 2
   └─ real variable (id: 806…655)

e = y[1] + logdet(Z) + sqrt(x) + minimum(y)

+ (concave; real)
├─ index (affine; real)
│  └─ 4-element real variable (id: 661…933)
├─ logdet (concave; real)
│  └─ 4×4 real variable (id: 715…088)
├─ geomean (concave; positive)
│  ├─ real variable (id: 806…655)
│  └─ [1.0]
└─ minimum (concave; real)
   └─ 4-element real variable (id: 661…933)

Examine the expression tree

e.children[2]

logdet (concave; real)
└─ 4×4 real variable (id: 715…088)

Constraints

A constraint is convex if convex combinations of feasible points are also feasible. Equivalently, feasible sets are convex sets.

In other words, convex constraints are of the form

• convexExpr <= 0
• concaveExpr >= 0
• affineExpr == 0

x <= 0

<= constraint (affine)
├─ real variable (id: 806…655)
└─ 0

square(x) <= sum(y)

<= constraint (convex)
├─ qol_elem (convex; positive)
│  ├─ real variable (id: 806…655)
│  └─ [1.0]
└─ sum (affine; real)
   └─ 4-element real variable (id: 661…933)

M = Z
for i = 1:length(y)
    global M += rand(size(Z)...)*y[i]
end
M ⪰ 0

sdp constraint (affine)
└─ + (affine; real)
   ├─ 4×4 real variable (id: 715…088)
   ├─ * (affine; real)
   │  ├─ 4×4 Array{Float64,2}
   │  └─ index (affine; real)
   │     └─ …
   ├─ * (affine; real)
   │  ├─ 4×4 Array{Float64,2}
   │  └─ index (affine; real)
   │     └─ …
   ├─ * (affine; real)
   │  ├─ 4×4 Array{Float64,2}
   │  └─ index (affine; real)
   │     └─ …
   └─ * (affine; real)
      ├─ 4×4 Array{Float64,2}
      └─ index (affine; real)
         └─ …

Problems

x = Variable()
y = Variable(4)
objective = 2*x + 1 - sqrt(sum(y))
constraint = x >= maximum(y)
p = minimize(objective, constraint)

minimize
└─ + (convex; real)
   ├─ * (affine; real)
   │  ├─ 2
   │  └─ real variable (id: 100…590)
   ├─ 1
   └─ - (convex; negative)
      └─ geomean (concave; positive)
         ├─ …
         └─ …
subject to
└─ >= constraint (convex)
   ├─ real variable (id: 100…590)
   └─ maximum (convex; real)
      └─ 4-element real variable (id: 153…790)

status: `solve!` not called yet

# solve the problem
solve!(p, () -> SCS.Optimizer(verbose=0))
p.status

OPTIMAL::TerminationStatusCode = 1

x.value

0.2500000000385289

# can evaluate expressions directly
evaluate(objective)

0.5000000000686586

Pass to solver

call a MathProgBase solver suited for your problem class

• see the list of Convex.jl operations to find which cones you're using
• see the list of solvers for an up-to-date list of solvers and which cones they support

to solve problem using a different solver, just import the solver package and pass the solver to the solve! method: eg

using Mosek
solve!(p, Mosek.Optimizer)

Warmstart

# Generate random problem data
m = 50;  n = 100
A = randn(m, n)
x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector
noise = .1*randn(m)    # gaussian noise
b = A*x♮ + noise      # noisy linear observations

# Create a (column vector) variable of size n.
x = Variable(n)

# nonnegative elastic net with regularization
λ = 1
μ = 1
problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)
@time solve!(problem, () -> SCS.Optimizer(verbose=0))
λ = 1.5
@time solve!(problem, () -> SCS.Optimizer(verbose=0), warmstart = true)

  0.213972 seconds (13.18 k allocations: 3.374 MiB)
  0.986449 seconds (1.60 M allocations: 83.633 MiB, 9.11% gc time)

DCP examples

# affine
x = Variable(4)
y = Variable(2)
sum(x) + y[2]

+ (affine; real)
├─ sum (affine; real)
│  └─ 4-element real variable (id: 204…209)
└─ index (affine; real)
   └─ 2-element real variable (id: 867…939)

2*maximum(x) + 4*sum(y) - sqrt(y[1] + x[1]) - 7 * minimum(x[2:4])

+ (convex; real)
├─ * (convex; real)
│  ├─ 2
│  └─ maximum (convex; real)
│     └─ 4-element real variable (id: 204…209)
├─ * (affine; real)
│  ├─ 4
│  └─ sum (affine; real)
│     └─ 2-element real variable (id: 867…939)
├─ - (convex; negative)
│  └─ geomean (concave; positive)
│     ├─ + (affine; real)
│     │  ├─ …
│     │  └─ …
│     └─ [1.0]
└─ - (convex; real)
   └─ * (concave; real)
      ├─ 7
      └─ minimum (concave; real)
         └─ …

# not dcp compliant
log(x) + square(x)

+ (Convex.NotDcp; real)
├─ log (concave; real)
│  └─ 4-element real variable (id: 204…209)
└─ qol_elem (convex; positive)
   ├─ 4-element real variable (id: 204…209)
   └─ 4×1 Array{Float64,2}

# $f$ is convex increasing and $g$ is convex
square(pos(x))

qol_elem (convex; positive)
├─ max (convex; positive)
│  ├─ 4-element real variable (id: 204…209)
│  └─ 0
└─ 4×1 Array{Float64,2}

# $f$ is convex decreasing and $g$ is concave
invpos(sqrt(x))

qol_elem (convex; positive)
├─ 4×1 Array{Float64,2}
└─ geomean (concave; positive)
   ├─ 4-element real variable (id: 204…209)
   └─ 4×1 Array{Float64,2}

# $f$ is concave increasing and $g$ is concave
sqrt(sqrt(x))

geomean (concave; positive)
├─ geomean (concave; positive)
│  ├─ 4-element real variable (id: 204…209)
│  └─ 4×1 Array{Float64,2}
└─ 4×1 Array{Float64,2}

This page was generated using Literate.jl.