How to do constrained optimization

Ronny Bergmann

This tutorial is a short introduction to using solvers for constrained optimization in Manopt.jl.

Introduction

A constrained optimization problem is given by

\[\tag{P} \begin{align*} \operatorname*{arg\,min}_{p∈\mathcal M} & f(p)\\ \text{such that} &\quad g(p) \leq 0\\ &\quad h(p) = 0,\\ \end{align*}\]

where $f: \mathcal M → ℝ$ is a cost function, and $g: \mathcal M → ℝ^m$ and $h: \mathcal M → ℝ^n$ are the inequality and equality constraints, respectively. The $\leq$ and $=$ in (P) are meant element-wise.

This can be seen as a balance between moving constraints into the geometry of a manifold $\mathcal M$ and keeping some, since they can be handled well in algorithms, see [BH19], [LB19] for details.

using Distributions, LinearAlgebra, Manifolds, Manopt, RandomRandom.seed!(42);

In this tutorial we want to look at different ways to specify the problem and its implications. We start with specifying an example problem to illustrate the different available forms.

We consider the problem of a Nonnegative PCA, cf. Section 5.1.2 in [LB19]

Let $v_0 ∈ ℝ^d$, $\lVert v_0 \rVert=1$, be a given spike signal, that is a signal that is sparse with only $s=\lfloor δd \rfloor$ nonzero entries.

\[Z = \sqrt{σ} v_0v_0^{\mathrm{T}}+N,\]

where $\sigma$ is a signal-to-noise ratio and $N$ is a matrix with random entries, distributed with zero mean and standard deviation $1/d$ on the off-diagonal and $2/d$ on the diagonal

d = 150; # dimension of v0σ = 0.1^2; # SNRδ = 0.1; sp = Int(floor(δ * d)); # SparsityS = sample(1:d, sp; replace = false);v0 = [i ∈ S ? 1 / sqrt(sp) : 0.0 for i in 1:d];N = rand(Normal(0, 1 / d), (d, d)); N[diagind(N, 0)] .= rand(Normal(0, 2 / d), d);Z = sqrt(σ) * v0 * transpose(v0) + N;

In order to recover $v_0$ we consider the constrained optimization problem on the sphere $\mathbb S^{d-1}$ given by

\[\begin{align*} \operatorname*{arg\,min}_{p∈\mathbb S^{d-1}} & -p^{\mathrm{T}}Zp\\ \text{such that} &\quad p \geq 0\\ \end{align*}\]

or in the previous notation $f(p) = -p^{\mathrm{T}}Zp$ and $g(p) = -p$. We first initialize the manifold under consideration

M = Sphere(d - 1)

Sphere(149)

A first augmented Lagrangian run

We first define $f$ and $g$ as usual functions

f(M, p) = -transpose(p) * Z * p;g(M, p) = -p;

since $f$ is a function defined in the embedding $ℝ^d$ as well, we obtain its gradient by projection.

grad_f(M, p) = project(M, p, -transpose(Z) * p - Z * p);

For the constraints this is a little more involved, since each function $g_i=g(p)_i=-p_i$ has to return its own gradient. These are again in the embedding just $\operatorname{grad} g_i(p) = -e_i$ the $i$ th unit vector. We can project these again onto the tangent space at $p$:

grad_g(M, p) = project.(    Ref(M), Ref(p), [[i == j ? -1.0 : 0.0 for j in 1:d] for i in 1:d]);

We further start in a random point:

p0 = rand(M);

Let’s verify a few things for the initial point

f(M, p0)

0.005667399180991248

Let’s also check how much the function $g$ is positive:

maximum(g(M, p0))

0.17885478285466855

Now as a first method we can just call the Augmented Lagrangian Method with a simple call:

@time v1 = augmented_Lagrangian_method(    M, f, grad_f, p0; g = g, grad_g = grad_g,    debug = [:Iteration, :Cost, :Stop, " | ", (:Change, "Δp : %1.5e"), 20, "\n"],    stopping_criterion = StopAfterIteration(300) | (        StopWhenSmallerOrEqual(:ϵ, 1.0e-5) & StopWhenChangeLess(M, 1.0e-8)    ));

Initial f(x): 0.005667 | 
# 20    f(x): -0.123557 | Δp : 1.04367e-07
# 40    f(x): -0.123557 | Δp : 2.19364e-10
# 60    f(x): -0.123557 | Δp : 7.59446e-11
The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-5).
At iteration 67 the algorithm performed a step with a change (7.594532846062582e-11) less than 1.0471285480508967e-5.
  9.913634 seconds (22.70 M allocations: 1.684 GiB, 5.56% gc time, 94.44% compilation time)

Now the function value is lower and the point is nearly within the constraints, namely up to numerical inaccuracies

f(M, v1)

-0.12355672067040904

maximum(g(M, v1))

1.0787741898413087e-11

A faster augmented Lagrangian run

Now this is a little slow, so we can modify two things:

1. Gradients should be evaluated in place, so for example

grad_f!(M, X, p) = project!(M, X, p, -transpose(Z) * p - Z * p);

1. The constraints are currently always evaluated all together, since the function grad_g always returns a vector of gradients. We first change the constraints function into a vector of functions. We further change the gradient both into a vector of gradient functions $\operatorname{grad} g_i,i=1,\ldots,d$, as well as gradients that are computed in place.

g2 = [(M, p) -> -p[i] for i in 1:d];grad_g2! = [    (M, X, p) -> project!(M, X, p, [i == j ? -1.0 : 0.0 for j in 1:d]) for i in 1:d];

We obtain

@time v2 = augmented_Lagrangian_method(    M, f, grad_f!, p0; g = g2, grad_g = grad_g2!, evaluation = InplaceEvaluation(),    debug = [:Iteration, :Cost, :Stop, " | ", (:Change, "Δp : %1.5e"), 20, "\n"],    stopping_criterion = StopAfterIteration(300) | (        StopWhenSmallerOrEqual(:ϵ, 1.0e-5) & StopWhenChangeLess(M, 1.0e-8)    ));

Initial f(x): 0.005667 | 
# 20    f(x): -0.123557 | Δp : 1.04367e-07
# 40    f(x): -0.123557 | Δp : 2.19364e-10
# 60    f(x): -0.123557 | Δp : 7.59446e-11
The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-5).
At iteration 67 the algorithm performed a step with a change (7.594532846062582e-11) less than 1.0471285480508967e-5.
  2.852497 seconds (7.26 M allocations: 640.897 MiB, 2.44% gc time, 96.36% compilation time)

As a technical remark: note that (by default) the change to InplaceEvaluations affects both the constrained solver as well as the inner solver of the subproblem in each iteration.

f(M, v2)

-0.12355672067040904

maximum(g(M, v2))

1.0787741898413087e-11

These are very similar to the previous values, but the solver took much less time and fewer memory allocations.

Exact penalty method

As a second solver, we have the Exact Penalty Method, which currently is available with two smoothing variants. These make the inner problem smooth, so that it can be solved with a solver for smooth optimization, by default again quasi-Newton: LogarithmicSumOfExponentials and LinearQuadraticHuber. We compare both here as well. The first smoothing technique is the default, so we can just call

@time v3 = exact_penalty_method(    M, f, grad_f!, p0; g = g2, grad_g = grad_g2!, evaluation = InplaceEvaluation(),    debug = [:Iteration, :Cost, :Stop, " | ", :Change, 50, "\n"],);

Initial f(x): 0.005667 | 
# 50    f(x): -0.122792 | Last Change: 0.001668
# 100   f(x): -0.123555 | Last Change: 0.000006
The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
At iteration 102 the algorithm performed a step with a change (3.024488503740816e-7) less than 1.0e-6.
  2.912256 seconds (13.14 M allocations: 3.684 GiB, 6.91% gc time, 78.79% compilation time)

We obtain a similar cost value as for the Augmented Lagrangian Solver from before, but here the constraint is actually fulfilled and not just numerically “on the boundary”.

f(M, v3)

-0.12355544268449428

maximum(g(M, v3))

-3.5897980609997717e-6

The second smoothing technique is often beneficial, when we have a lot of constraints (in the previously mentioned vectorial manner), since we can avoid several gradient evaluations for the constraint functions here. This leads to a faster iteration time.

@time v4 = exact_penalty_method(    M, f, grad_f!, p0; g = g2, grad_g = grad_g2!,    evaluation = InplaceEvaluation(),    smoothing = LinearQuadraticHuber(),    debug = [:Iteration, :Cost, :Stop, " | ", :Change, 50, "\n"],);

Initial f(x): 0.005667 | 
# 50    f(x): -0.123559 | Last Change: 0.000004
# 100   f(x): -0.123557 | Last Change: 0.000000
The value of the variable (ϵ) is smaller than or equal to its threshold (1.0e-6).
At iteration 100 the algorithm performed a step with a change (1.1298698721628557e-8) less than 1.0715193052376023e-6.
  2.447675 seconds (9.58 M allocations: 1.735 GiB, 5.87% gc time, 86.68% compilation time)

For the result we see the same behavior as for the other smoothing.

f(M, v4)

-0.1235566792814103

maximum(g(M, v4))

3.0266225123116394e-8

Comparing to the unconstrained solver

We can compare this to the global optimum on the sphere, which is the unconstrained optimization problem, where we can just use Quasi Newton.

Note that this is much faster, since every iteration of the constrained solvers above performs a full quasi-Newton run as its sub solver.

@time w1 = quasi_Newton(    M, f, grad_f!, p0; evaluation = InplaceEvaluation());

  1.259024 seconds (3.23 M allocations: 172.658 MiB, 1.76% gc time, 99.57% compilation time)

f(M, w1)

-0.13990874034056583

But for sure the constraints here are not fulfilled and we have quite positive entries in $g(w_1)$

maximum(g(M, w1))

0.1180320073974675

Literature

Status `~/work/Manopt.jl/Manopt.jl/tutorials/Project.toml`
  [47edcb42] ADTypes v1.24.0
  [6e4b80f9] BenchmarkTools v1.8.0
  [5ae59095] Colors v0.13.1
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [31c24e10] Distributions v0.25.131
  [26cc04aa] FiniteDifferences v0.12.34
  [f6369f11] ForwardDiff v1.4.5
  [8ac3fa9e] LRUCache v1.6.2
  [af67fdf4] ManifoldDiff v0.4.5
  [1cead3c2] Manifolds v0.11.29
  [3362f125] ManifoldsBase v2.5.1
  [0fc0a36d] Manopt v0.6.7 `.`
  [91a5bcdd] Plots v1.41.7
  [731186ca] RecursiveArrayTools v4.5.1
  [37e2e46d] LinearAlgebra v1.12.0
  [9a3f8284] Random v1.11.0