Deep Dives

Machine Learning Optimizer in Constraint Parameter Space: The Interior Point Methods

A complete step-by-step case study with python from scratch

Table of Contents

· Motivation
· Problem Statement
· The Barrier Method
· Experiments
  ∘ Scenario 1
  ∘ Scenario 2
  ∘ Scenario 3
· Conclusion

Motivation

You are a decent trader in the stock market. You optimize your profit every day by some discreet strategy, and it has no limit: you can earn big in a day (positive profit) or lose so much that you tweet a joke out of it (negative profit). Since there are no constraints on the profit you optimize, there is a variety of strategies you can choose: short equity, pairs trading, heck even making a Geometric Brownian Motion model out of it.

Besides trading, you are also a great gardener in town. You have three types of plants and sell them at the closest traditional market every day: tomatoes, mint, and spinach. To obtain maximum profit, of course, you need to plant as many as possible. But there’s a catch: you have limited land and a certain trend of demands. So, you want effective land utilization and not sell too many at the same time to avoid unsold products. For example, you use 1/2 of your land for tomatoes, another 1/4 for spinach, 1/8 for mint, and the rest is left unoccupied. In other words, in this case, there are constraints to optimizing profit.

In real life, it’s easy to see that the trader scenario is less likely to occur than the gardener one. In mathematical formulation, the latter can be written as

where f and the functions cᵢ are all smooth, real-valued functions on a subset of ℝⁿ, and ℰ and ℐ are two finite sets of indices. We call f the objective function, while cᵢ, i ∈ ℰ are the equality constraints and cᵢ, i ∈ ℐ are the inequality constraints. A set of points or iterates x satisfying all cᵢ is called a feasible set and we denote it by Ω.

This formulation is called a constrained optimization problem. Now, how to solve it?

Problem Statement

As a concrete example, we will consider a simple constrained optimization problem as follows:

subject to constraints hᵢ(x₁, x₂):

In other words, we are minimizing a non-linear bivariate objective function f with constraints hᵢ. Neglecting hᵢ, it’s easy to verify that the minimizer is x* = (x₁*, x₂*) = (5, 6) with f(x₁, x₂) = 0. But with hᵢ, the minimizer is not as trivial as it seems. Here, we write x = (x₁, x₂) for convenience.

Import some libraries to get started.

Here’s f in code.

The Barrier Method

The Barrier Method is a part of Interior Point Methods, a class of algorithms that solve linear and nonlinear convex optimization problems, first introduced in 1948 by John von Neumann. However, the method was inefficient and slower in practice as compared to the Simplex method. Karmarkar improved the algorithm in 1984, and with further refinements, the Interior Point Methods have become extremely efficient and have found applications in virtually every field of continuous optimization.

Interior Point Methods typically solve the constrained convex optimization problem by applying Newton Method to a sequence of equality constrained problems. Barrier methods, as the name suggest, employ barrier functions to integrate inequality constraints into the objective function. Since we want to merge inequality constraints to the objective, the following conditions on barrier function ϕ seem natural:

• The barrier function should map the feasible set Ω to real space, i.e. ϕ : Ω → ℝ.
• ϕ(x) is integrated into the objective function and we intend to use Newton Method which is a gradient-based method, so the barrier function should be reasonably “smooth” (continuously differentiable).
• ϕ(x) should behave like a barrier and grow unbounded as x enters the boundary of the feasible set.
• The Hessian (second derivative for functions of scalars) of the barrier function needs to be positive definite, i.e. ∇²ϕ(x) ≻ 0, ∀ x ∈ Ω.

There are many examples of barrier functions. In particular, right now we will use a barrier function as follows.

where m is the number of constraints hᵢ(x) ≤ 0. It’s easy to see that this function satisfies all conditions explained above. Then, ϕ(x) can be integrated into the objective function f such that the optimization problem can be approached by a new one

In our case, this can be written as

where t > 0 is a parameter determined by the user. If t < ∞ then the solution of the new problem belongs to the interior of Ω, i.e. x ∈ Ω, and as t → ∞ the solution tends to the optimum, i.e. x(t) → x*.

Now, the gradient g and the Hessian H can be calculated as

with

After knowing all this, we are ready to create the Barrier Method algorithm:

1. Set an initial point x ∈ Ω, initial parameter t, and tolerance ε for stopping criterion. Here we will use ε = 1 × 10⁻⁵.
2. Do the following a-b-c loop until the stopping criterion is met. The stopping criterion used is the duality gap m/t ≤ ε, with m being the number of constraints hᵢ. In our case, m = 5.

a) Centering step: calculate

using Newton Method by doing this iteration

until the stopping criteria of the Newton Method are met. The criteria used are ǁxₖ₊₁ − xₖǁ₂ ≤ 10⁻⁵ or maximum iteration reaches 1000.

b) Update x := x*(t).

c) Increase t using the following equation

with ν is a parameter associated with the Barrier function ϕ. Here we use ν = 0.01.

Besides executing the algorithm, note that the BarrierMethod function above also prints out the current iterate x, the value of the objective function f(x), and the duality gap m/t.

Experiments

This is where the fun begins. To visualize the process of the optimization, we will first make a plot_all function. This function plots the result from BarrierMethod function with the following details:

1. The yellow shaded area indicates the feasible set
2. The black area indicates the non-feasible set
3. The magenta curve is the iteration path
4. The cyan point is the solution found by the algorithm

The convergence of iteration should depend on the initial choice of x and t. For that, we pick three choices as follows:

1. Initial values x = (0.50, 0.75) and t = 0.1
2. Initial values x = (0.50, 0.75) and t = 1
3. Initial values x = (1.50, 2.00) and t = 0.1

Scenario 1: Initial values x = (0.50, 0.75) and t = 0.1

Initial condition: x = [0.5  0.75], f(x) = 47.8125

Iteration: 1 	 x = [1.0202 1.1642], f(x) = 39.2236, gap = 28.2609
Iteration: 2 	 x = [1.1088 1.1778], f(x) = 38.3945, gap = 15.9735
Iteration: 3 	 x = [1.2357 1.1815], f(x) = 37.3873, gap = 9.0285
Iteration: 4 	 x = [1.3893 1.1683], f(x) = 36.3827, gap = 5.1031
Iteration: 5 	 x = [1.5441 1.1408], f(x) = 35.5548, gap = 2.8843
Iteration: 6 	 x = [1.6786 1.1079], f(x) = 34.9641, gap = 1.6303
Iteration: 7 	 x = [1.7838 1.0771], f(x) = 34.5792, gap = 0.9215
Iteration: 8 	 x = [1.8603 1.0519], f(x) = 34.3412, gap = 0.5208
Iteration: 9 	 x = [1.9129 1.0333], f(x) = 34.1983, gap = 0.2944
Iteration: 10 	 x = [1.9473 1.0205], f(x) = 34.1142, gap = 0.1664
Iteration: 11 	 x = [1.9689 1.0122], f(x) = 34.0654, gap = 0.0940
Iteration: 12 	 x = [1.9819 1.0072], f(x) = 34.0372, gap = 0.0532
Iteration: 13 	 x = [1.9896 1.0041], f(x) = 34.0211, gap = 0.0300
Iteration: 14 	 x = [1.9941 1.0024], f(x) = 34.0120, gap = 0.0170
Iteration: 15 	 x = [1.9966 1.0013], f(x) = 34.0068, gap = 0.0096
Iteration: 16 	 x = [1.9981 1.0008], f(x) = 34.0038, gap = 0.0054
Iteration: 17 	 x = [1.9989 1.0004], f(x) = 34.0022, gap = 0.0031
Iteration: 18 	 x = [1.9994 1.0002], f(x) = 34.0012, gap = 0.0017
Iteration: 19 	 x = [1.9997 1.0001], f(x) = 34.0007, gap = 0.0010
Iteration: 20 	 x = [1.9998 1.0001], f(x) = 34.0004, gap = 0.0006
Iteration: 21 	 x = [1.9999 1.0000], f(x) = 34.0002, gap = 0.0003
Iteration: 22 	 x = [1.9999 1.0000], f(x) = 34.0001, gap = 0.0002
Iteration: 23 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0001
Iteration: 24 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0001
Iteration: 25 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0000
Iteration: 26 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0000
Iteration: 27 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0000
Iteration: 28 	 x = [2.0000 1.0000], f(x) = 34.0000, gap = 0.0000

The algorithm found a solution x* = (2, 1) after 28 iterations by minimizing the duality gap such that m/t ≤ ε. The value of the function is minimized at f(x*) = 34. It can be confirmed from the contour plot that x* = (2, 1) is indeed the only solution. Note that the iteration path is not going smooth to x* because of the centering step.

Scenario 2: Initial values x = (0.50, 0.75) and t = 1

Initial condition: x = [0.5  0.75], f(x) = 47.8125

Iteration: 1 	 x = [5.3475 6.2329], f(x) = 0.1750, gap = 2.8261
Iteration: 2 	 x = [5.2030 6.1341], f(x) = 0.0592, gap = 1.5974
Iteration: 3 	 x = [5.1171 6.0767], f(x) = 0.0196, gap = 0.9029
Iteration: 4 	 x = [5.0670 6.0436], f(x) = 0.0064, gap = 0.5103
Iteration: 5 	 x = [5.0381 6.0247], f(x) = 0.0021, gap = 0.2884
Iteration: 6 	 x = [5.0216 6.0140], f(x) = 0.0007, gap = 0.1630
Iteration: 7 	 x = [5.0123 6.0079], f(x) = 0.0002, gap = 0.0921
Iteration: 8 	 x = [5.0069 6.0045], f(x) = 0.0001, gap = 0.0521
Iteration: 9 	 x = [5.0039 6.0025], f(x) = 0.0000, gap = 0.0294
Iteration: 10 	 x = [5.0022 6.0014], f(x) = 0.0000, gap = 0.0166
Iteration: 11 	 x = [5.0013 6.0008], f(x) = 0.0000, gap = 0.0094
Iteration: 12 	 x = [5.0007 6.0005], f(x) = 0.0000, gap = 0.0053
Iteration: 13 	 x = [5.0004 6.0003], f(x) = 0.0000, gap = 0.0030
Iteration: 14 	 x = [5.0002 6.0001], f(x) = 0.0000, gap = 0.0017
Iteration: 15 	 x = [5.0001 6.0001], f(x) = 0.0000, gap = 0.0010
Iteration: 16 	 x = [5.0001 6.0000], f(x) = 0.0000, gap = 0.0005
Iteration: 17 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0003
Iteration: 18 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0002
Iteration: 19 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0001
Iteration: 20 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0001
Iteration: 21 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0000
Iteration: 22 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0000
Iteration: 23 	 x = [5.0000 6.0000], f(x) = 0.0000, gap = 0.0000

A solution is found to be x* = (5, 6) with f(x) = 0 after 23 iterations. However, since the initial t is too big, at some point in the iteration process, the temporary iterate x broke through the feasible set, and the iterations that follow try to find a solution outside the feasible set, which we know is (5, 6).

Scenario 3: Initial values x = (1.50, 2.00) and t = 0.1

Initial condition: x = [1.5 2. ], f(x) = 28.2500

Iteration: 1 	 x = [1.4440 8.0101], f(x) = 16.6853, gap = 28.2609
Iteration: 2 	 x = [1.5565 7.2585], f(x) = 13.4414, gap = 15.9735
Iteration: 3 	 x = [1.6742 6.7642], f(x) = 11.6450, gap = 9.0285
Iteration: 4 	 x = [1.7785 6.4526], f(x) = 10.5827, gap = 5.1031
Iteration: 5 	 x = [1.8587 6.2633], f(x) = 9.9374, gap = 2.8843
Iteration: 6 	 x = [1.9138 6.1514], f(x) = 9.5476, gap = 1.6303
Iteration: 7 	 x = [1.9490 6.0864], f(x) = 9.3161, gap = 0.9215
Iteration: 8 	 x = [1.9704 6.0491], f(x) = 9.1810, gap = 0.5208
Iteration: 9 	 x = [1.9830 6.0279], f(x) = 9.1031, gap = 0.2944
Iteration: 10 	 x = [1.9903 6.0158], f(x) = 9.0585, gap = 0.1664
Iteration: 11 	 x = [1.9945 6.0089], f(x) = 9.0332, gap = 0.0940
Iteration: 12 	 x = [1.9969 6.0050], f(x) = 9.0188, gap = 0.0532
Iteration: 13 	 x = [1.9982 6.0029], f(x) = 9.0106, gap = 0.0300
Iteration: 14 	 x = [1.9990 6.0016], f(x) = 9.0060, gap = 0.0170
Iteration: 15 	 x = [1.9994 6.0009], f(x) = 9.0034, gap = 0.0096
Iteration: 16 	 x = [1.9997 6.0005], f(x) = 9.0019, gap = 0.0054
Iteration: 17 	 x = [1.9998 6.0003], f(x) = 9.0011, gap = 0.0031
Iteration: 18 	 x = [1.9999 6.0002], f(x) = 9.0006, gap = 0.0017
Iteration: 19 	 x = [1.9999 6.0001], f(x) = 9.0003, gap = 0.0010
Iteration: 20 	 x = [2.0000 6.0001], f(x) = 9.0001, gap = 0.0006
Iteration: 21 	 x = [2.0000 6.0000], f(x) = 9.0001, gap = 0.0003
Iteration: 22 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0002
Iteration: 23 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0001
Iteration: 24 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0001
Iteration: 25 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0000
Iteration: 26 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0000
Iteration: 27 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0000
Iteration: 28 	 x = [2.0000 6.0000], f(x) = 9.0000, gap = 0.0000

This is an example of wrong practice in implementing the Barrier Method. Note that the initial x is not in the feasible set, which causes the iterations to produce iterates that violate the constraints. Hence, we end up with the solution x* which we don’t want in the first place.

Conclusion

The Barrier Method is an effective algorithm for optimization problems with constraints. It integrates inequality constraints to the objective function using a barrier function. The method uses the gradient and hessian of the objective function, and the Newton method as the centering step. We’ve seen the Barrier Method in action and found out that the initial iterate x and parameter t are crucial for it to work the way we desire.

🔥 Hi there! If you enjoy this story and want to support me as a writer, consider becoming a member. For only $5 a month, you’ll get unlimited access to all stories on Medium. If you sign up using my link, I’ll earn a small commission.

🔖 Want to know more about how classical machine learning models work and how they optimize their parameters? Or an example of MLOps megaprojects? What about cherry-picked top-notch articles of all time? Continue reading: