I have manually translated and adapted the QuadProg solver for quadratic programming problems made by Berwin A. Turlach. His code was originally published under the GNU Library License, which has now been superseded by the GNU Lesser License. This adapted version honors the original work and is thus distributed under the same license.
Introduction
Despite the name, the terms linear or quadratic programming have little resemblance to the set of activities most people now know as programming. Those terms usually usually refers to a specific set of function optimization methods, i.e. methods which can be used to determine the maximum or minimum points of special kinds of functions under a given number of solution constraints. For example, suppose we would like to determine the minimum value of the function:
f(x, y) = 2x + y + 4
Under the constraints that x and y must be non-negative (i.e. either positive or zero). This may seem fairly simple and trivial, but remember that practical linear programming problems may have hundreds or even thousands of variables and possibly million constraints.
When the problem to be solved involves a quadratic function instead of a linear function, but still presents linear constraints, this problem can be cast as a quadratic programming problem. Quadratic functions are polynomial functions in each each term may have at most a total degree of 2. For example, consider the function
f(x, y, z) = 2x² + 5xy + y² - z² + x – 5.
Now let's check the sum of the degrees for each variable on the polynomial terms. We start by writing the missing terms of the polynomial
f(x, y, z) = 2 x2y0z0 + 5 x1y1z0 + 2 x0y2z0 - x0y0z2 + x1y0z0 - 5 x0y0z0
and then proceed to check the sum of the degrees at each term. In the first term, 2+0+0 = 2. For the second, 1+1+0 = 2, and so on. Those functions have a nice property that they can be expressed in a matrix form
f(x) = 1/2 xT Q x + cTx.
Here, x and c are vectors. The matrix Q is a symmetric matrix specifying how the variables combine in the quadratic terms of the function. If this matrix is positive definite, then the function is convex, and the optimization has a single, unique optimum (we say it has a global optimum). The coefficients c specify the linear terms of our function.
Source code
The available source code is based on a translation of the Fortran code written by Berwin A. Turlach. However, some modifications have been made. Fortran uses column-major ordering for matrices, meaning that matrices are stored in memory in sequential order of column elements. Almost all other languages use row-major ordering, including C, C++, Java and C#. In order to improve data locality, I have modified the code to use the transpose of the original matrices D and A. I have also modified the QP formulation adopted in the Goldfarb and Idnani paper to reflect the form presented in the introduction.
This code is part of the Accord.NET Framework. However, the version available in this blog post will most likely not be the most recently, updated, fixed and enhanced version of the code. For the latest version, be sure to download the latest version of the framework on the project site or through a NuGet package.
Using the code
The first step in solving a quadratic programming problem is, well, specifying the problem. To specify a quadratic programming problem, one would need two components: a matrix D describing the relationship between the quadratic terms, and a vector d describing the linear terms. Perhaps this would work better with an example.
Suppose we are trying to solve a minimization problem. Given a function, the goal in such problems is to find the correct set of function arguments which would result in the minimum possible value for the function. An example of a quadratic minimization problem is given below:
However, note that this problem involves a set of constraints. The required solution for this minimization problem is required to lie in the interval specified by the constraints. More specifically, any x and y pair candidate for being a minimal of the function must respect the relations x - y = 5 and x >= 10. Thus, instead of lying in the unconstrained minimum of the function surface shown above, the solution lies slightly off the center of the surface. This is an obvious easy problem to solve manually, but it will fit for this demonstration.
As it can be seen (and also live demonstrated by asking Wolfram Alpha) the solution lies on the point (10,5), and the constrained minimum of the function is given by 170. So, now that we know what a quadratic programming problem looks like, how can we actually solve it?
Specifying the objective function
The first step in solving a quadratic programming problem is to specify the objective function. Using this code, there are three ways to specify it. Each of them has their own advantages and disadvantages.
1. Manually specifying the QP matrix.
This is the most common approach for numerical software, and probably the most cumbersome for the user. The problem matrix has to be specified manually. This matrix is sometimes denoted Q, D or H as it actually denotes the Hessian matrix for the problem.
The matrix Q is used to describe the quadratic terms of our problem. It is a n x n matrix, in which n corresponds to the number of variables in our problem, covering all possible combinations of variables. Recall our example given on the start of this section. We have 2 variables, x and y. Thus, our matrix Q is 2 x 2. The possible combinations for x and y are expressed in the table below.
To form our matrix Q, we can take all coefficients associated with each pair mentioned on the table above. The diagonal elements should also be multiplied by two (this is actually because the matrix is the Hessian matrix of the problem: it is the matrix of all second-order derivatives for the function. Since we have only at most quadratic terms, the elementary power rule of derivation “drops” the ² from the x² and y² terms – I think a mathematician would hit me with a stick for explaining it like this, but it serves well for a quick, non-technical explanation).
Remember our quadratic terms were 2x² - 1xy + 4y². Writing the terms on their proper position and differentiating, we have:
As it can be seen, the matrix is also symmetric (and often, but not always, positive definite). The next step, more trivial, is to write a vector d containing the linear terms. The linear terms are –5x –6y, and thus our vector d can be given by:
Therefore our C# code can be created like this:
double[,] Q =
{
{ +4, -1 },
{ -1, +8 },
};
double[] d =
{
-5,
-6
};
2. Using lambda expressions
This approach is a bit more intuitive and less error prone. However, it involves lambdas functions and some people find it hard to follow them. Another disadvantage is that we will lose the edit & continue debugging ability of visual studio. The advantage is that the compiler may catch some obvious syntax errors automatically.
double x = 0, y = 0;
QuadraticObjectiveFunction f = // 2x² - xy + 4y² - 5x - 6y
new QuadraticObjectiveFunction(() => 2*(x*x) - (x*y) + 4*(y*y) - 5*x - 6*y);
Note that the x and y variables could have been initialized to any value. They are only used as symbols, and not used in any computations.
3. Using text strings
This approach is more intuitive but a bit more error prone. The function can be specified using strings, as in a standard mathematical formula. However, since all we have are strings, there is no way to enforce static, compile time checking.
QuadraticObjectiveFunction f =
new QuadraticObjectiveFunction("2x² - xy + 4y² - 5x - 6y");
Couldn’t be easier.
Specifying the constraints
The next step in specifying a quadratic programming problem is to specify the constraints. The constraints can be specified in almost the same way as the objective function.
1. Manually specifying the constraints matrix
The first option is to manually specify the constraints matrix A and vector b. The constraint matrix expresses the way the variables should be combined when compared to corresponding value on vector b. It is possible to specify either equality constraints or inequality constraints. The formulation used in this code is slightly different from the one used in Turlach’s original implementation. The constraints are supposed to be in the form:
A1 x = b1
A2 x = b2
This means that each line of matrix A expresses a possible combination of variables x which should be compared to the corresponding line of b. An integer variable m can be specified to denote how many of the first rows of matrix A should be treated as equalities rather than inequalities. Recall that in our example the constraints are given by 1x -1y = 5 and 1x = 10. Lets write this down in a tabular form:
#
|
x
|
y
|
?
|
b
|
q1
|
1
|
-1
|
=
|
5
|
q2
|
1
|
0
|
=
|
10
|
Thus our matrix A and vector b can be specified as:
And not forgetting that m = 1, because the first constraint is actually an equality.
2. Using classes and objects
A more natural way to specify constraints is using the classes and objects of the Accord.NET Framework. The LinearConstraint class allows one to specify a single constraint using an object-oriented approach. It doesn’t have the most intuitive usage on earth, but has much more expressiveness. It can also be read aloud, it that adds anything! :-)
List<LinearConstraint> list = new List<LinearConstraint>();
list.Add(new LinearConstraint(numberOfVariables: 1)
{
VariablesAtIndices = new[] { 0 }, // index 0 (x)
ShouldBe = ConstraintType.GreaterThanOrEqualTo,
Value = 10
});
list.Add(new LinearConstraint(numberOfVariables: 2)
{
VariablesAtIndices = new int[] { 0, 1 }, // index 0 (x) and index 1 (y)
CombinedAs = new double[] { 1, -1 }, // when combined as 1x -1y
ShouldBe = ConstraintType.EqualTo,
Value = 5
});
The specification is centered around the notion that variables are numbered and have an associated index. For example, x is the zero-th variable of the problem. Thus x has an index of 0 and y has an index of 1. So for example, reading aloud the last constraint, it is possible to express how the variables at indices 0 and 1, when combined as 1x and –1y, should be equal to value 5.
2. Using lambda expressions
A more intuitive way to express constraints is again using lambda expressions. And again the problems are the same: some people find it hard to follow and we lose edit & continue.
var constraints = new List<LinearConstraint>();
constraints.Add(new LinearConstraint(f, () => x - y == 5));
constraints.Add(new LinearConstraint(f, () => x >= 10));
3. Using text strings
Same as above, but with strings.
var constraints = new List<LinearConstraint>();
constraints.Add(new LinearConstraint(f, "x - y = 5"));
constraints.Add(new LinearConstraint(f, "x >= 10"));
Finally, creating and solving the problem
Once we have specified what do we want, we can now ask the code for a solution. In case we have opted for manually specifying the matrix A, vector b and integer m, we can use:
// Create the optimization problem
var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables:2, A, b, m);
In case we have opted for creating a list of constraints instead, we can use:
// Create our optimization problem
var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables: 2, constraints: list);
After the solver object has been created, we can call Minimize() to solve the problem. In case we have opted for manually specifying Q and d, we can use:
// Attempt to solve the problem
double minimumValue = solver.Minimize(Q, d);
And in case we have opted for creating a QuadraticObjectiveFunction object, we can use:
// Attempt to solve the problem
double minimumValue = target.Minimize(f);
In either case, the solution will be available in the Solution property of the solver object, and will be given by:
double value = solver.Value; // f(x,y) = 170
double x = solver.Solution[0]; // x = 10
double y = solver.Solution[1]; // y = 5
Sample application
The Accord.NET Framework now includes a sample application demonstrating the use of the Goldfarb-Idnani Quadratic Programming Solver. It can be downloaded at the Accord.NET Framework site, and also comes together with recent versions of the framework (> 2.6).
Solver sample application included in the Accord.NET Framework.
Remarks
Because the code has been translated by hand (in contrast of using automatic translators such as f2c) there could be potential bugs in the code. I have tested the code behavior against R’s quadprog package and still didn’t find errors. But this does not mean the code is bug-free. As always, as is the case of everything else in this blog, this code is published in good faith, but I can not guarantee the correctness of everything. Please read the disclaimer for more information.
References
B. A. Turlach. QuadProg: Quadratic Programming Routines (1998). Website, available at http://school.maths.uwa.edu.au/~berwin/software/quadprog.html.
D. Goldfarb and A. Idnani. Dual and Primal-Dual Methods for Solving Strictly Convex Quadratic Programs (1982). In J. P. Hennart (ed.), Numerical Analysis, Springer-Verlag, Berlin, pages 226–239.
D. Goldfarb and A. Idnani. A numerically stable dual method for solving strictly convex quadratic programs (1983). Mathematical Programming, 27, 1–33.
Wikipedia contributors, Quadratic programming. In Wikipedia, The Free Encyclopedia. Available on: http://en.wikipedia.org/wiki/Quadratic_programming.
Wikipedia contributors, Quadratic form. In Wikipedia, The Free Encyclopedia. Available on: http://en.wikipedia.org/wiki/Quadratic_form.
Wikipedia contributors, Linear programming. In Wikipedia, The Free Encyclopedia. Available on: http://en.wikipedia.org/wiki/Linear_programming.
- Wikipedia contributors, Mathematical optimization. In Wikipedia, The Free Encyclopedia. Available on: http://en.wikipedia.org/wiki/Mathematical_optimization.
Wolfram Alpha LLC. Wolfram|Alpha, 2012
