Linearly Constrained Optimization
Mathwrist Code Example Series
(January 20, 2023)
Abstract
This document gives some code examples on how to use the linearly constrained optimization features in Mathwrist’s C++ Numerical Programming Library (NPL).
1 Box Constraints
This example illustrates minimizing a quadratic objective subject to simple bound constraints , .
Here, we use two simple bound optimization solvers (modified Newton and Quasi-Newton methods) for general smooth functions.
1
constexpr size_t N = 100;
3
// Objective function:
4
//
5
// number of variables 400
6
// f(x) = \sum_{i=0}^{n-1} x_i^2 -
7
// \sum_{i=0}^{n-2}x_i x_{i+1} - x_0 - x_{n-1}
9
// Set up the Hessian matrix
10
Matrix H(N, N);
11
{
12
for (size_t i = 0; i < N; ++i)
13
{
14
H(i, i) = 2;
15
if (i > 0)
16
H(i, i - 1) = -1;
17
if (i + 1 < N)
18
H(i, i + 1) = -1;
19
}
20
};
22
struct Fx : public FunctionND
23
{
24
const Matrix& m_H;
25
explicit Fx(const Matrix& H) : m_H(H) {}
27
// Clone is not used in this example.
28
FunctionND* clone() const { return 0; }
30
double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const
31
{
32
Matrix Hx = m_H * x;
33
double retval = 0.5 * x.inner_product(Hx);
34
retval -= x(0, 0) + x(N - 1, 0);
36
if (_g != 0)
37
{
38
*_g = Hx;
40
// Gradient wrt x_0 and x_{n-1} has extra term -1
41
(*_g)(0, 0) -= 1;
42
(*_g)(N-1, 0) -= 1;
43
}
45
if (_H != 0)
46
*_H = m_H;
48
return retval;
49
}
50
size_t dimension_size() const { return N; }
51
};
53
// simple bounds: 0 <= x_i <= 0.9, i < n-1.
54
ftl::Buffer<Bound> x_bounds(N);
55
x_bounds.fill_n(0, N-1, Bound(0, 0.9));
57
Fx fx(H);
59
GS_SimpleBoundNewton newton(fx, x_bounds, false);
60
newton.set_max_iter(1000);
62
// Initial guess x = 0
63
Matrix x(N, 1);
65
// solve by modified Newton.
66
double y = newton(x);
67
assert(std::abs(y + 0.9925) < 1.e-8);
69
for (size_t i = 0; i < N-1; ++i)
70
{
71
assert(x_bounds[i].violation(x(i, 0), 1.e-10) == 0);
72
}
74
// Quasi Newton solver
75
GS_SimpleBoundQuasiNewton quasi(fx, x_bounds);
77
// Initial guess x = 0
78
x.zeros();
80
// solve by discrete Newton.
81
y = newton(x);
82
assert(std::abs(y + 0.9925) < 1.e-8);
84
for (size_t i = 0; i < N-1; ++i)
85
{
86
assert(x_bounds[i].violation(x(i, 0), 1.e-10) == 0);
87
}
2 General Linear Constraints
This example illustrates minimizing an indefinite quadratic objective function subject to linear and simple bound constraints.
, subject to
Here, we use active-set based modified Newton and Quasi-Newton methods for general smooth objective function. Alternatively, we can use QP active set method as well.
1
// 2 general inear equality constraint wrt 3 variables.
2
//
3
// 0.05 x_1 + 0.1x_2 + 0.15x_3 = 0.1
4
// x_1 + x_2 + x_3 = 1
5
LinearProg::Constraints constr(2, 3);
6
constr.A(0, 0) = 0.05;
7
constr.A(0, 1) = 0.1;
8
constr.A(0, 2) = 0.15;
9
constr.A(1, 0) = 1;
10
constr.A(1, 1) = 1;
11
constr.A(1, 2) = 1;
12
constr.general_bounds[0] = Bound(0.1, 0.1);
13
constr.general_bounds[1] = Bound(1, 1);
15
// simple bounds x1 >= 0, x2 >= 0, x3 >= 0
16
constr.simple_bounds[0] = Bound::lower(0);
17
constr.simple_bounds[1] = Bound::lower(0);
18
constr.simple_bounds[2] = Bound::lower(0);
20
// Objective function
21
// f(x) = 0.2 x_1^2 + 0.21 x_2^2 + 0.28 x_3^2 +
22
// 0.3x_1 x_2 - 0.5 x_1 x_3 + 0.18 x_2 x_3
23
struct Fx : public FunctionND
24
{
25
// Clone is not used in this example.
26
FunctionND* clone() const { return 0; }
28
double eval(const Matrix& x, Matrix* _g=0, Matrix* _H=0) const
29
{
30
double x1 = x(0, 0);
31
double x2 = x(1, 0);
32
double x3 = x(2, 0);
34
double retval = 0.2 * x1 * x1 + 0.21 * x2 * x2 +
35
0.28 * x3 * x3 + 0.3 * x1 * x2 -
36
0.5 * x1 * x3 + 0.18 * x2 * x3;
38
if (_g != 0)
39
{
40
(*_g)(0, 0) = 0.4 * x1 + 0.3 * x2 - 0.5 * x3;
41
(*_g)(1, 0) = 0.42 * x2 + 0.3 * x1 + 0.18 * x3;
42
(*_g)(2, 0) = 0.56 * x3 - 0.5 * x1 + 0.18 * x2;
43
}
45
if (_H != 0)
46
{
47
(*_H)(0, 0) = 0.4;
48
(*_H)(0, 1) = (*_H)(1,0) = 0.3;
49
(*_H)(0, 2) = (*_H)(2, 0) = -0.5;
51
(*_H)(1, 1) = 0.42;
52
(*_H)(1, 2) = (*_H)(2, 1) = 0.18;
54
(*_H)(2, 2) = 0.56;
55
}
57
return retval;
58
}
59
size_t dimension_size() const { return 3; }
60
};
62
Fx fx;
64
Matrix x(3, 1);
66
GS_ActiveSetNewton newton(fx, constr, true);
67
double y = newton(x);
69
assert(std::abs(y + 0.005) < 1.e-10);
70
assert(std::abs(x(0, 0) - 0.5) < 1.e-10);
71
assert(std::abs(x(1, 0)) < 1.e-10);
72
assert(std::abs(x(2, 0) - 0.5) < 1.e-10);
74
GS_ActiveSetQuasiNewton quasi(fx, constr);
75
x.zeros();
76
y = quasi(x);
78
assert(std::abs(y + 0.005) < 1.e-10);
79
assert(std::abs(x(0, 0) - 0.5) < 1.e-10);
80
assert(std::abs(x(1, 0)) < 1.e-10);
81
assert(std::abs(x(2, 0) - 0.5) < 1.e-10);