Quadratic Programming
In this chapter we turn the LP problem from Chapter 10 into a Quadratic Programming (QP) problem, and the first MIP model from Chapter 11 into a Mixed Integer Quadratic Programming (MIQP) problem. The chapter shows how to
- define quadratic objective functions,
- incrementally define and solve problems.
Chapter 7 shows how to formulate and solve the same examples with Mosel and in Chapter 17 the QP problem is input and solved directly with Xpress Optimizer.
Problem description
The investor may also look at his portfolio selection problem from a different angle: instead of maximizing the estimated return and limiting the portion of high-risk investments he now wishes to minimize the risk whilst obtaining a certain target yield. He adopts the Markowitz idea of getting estimates of the variance/covariance matrix of estimated returns on the securities. (For example, hardware and software company worths tend to move together, but are oppositely correlated with the success of theatrical production, as people go to the theater more when they have become bored with playing with their new computers and computer games.) The return on theatrical productions are highly variable, whereas the treasury bill yield is certain.
Question 1: Which investment strategy should the investor adopt to minimize the variance subject to getting some specified minimum target yield?
Question 2: Which is the least variance investment strategy if the investor wants to choose at most four different securities (again subject to getting some specified minimum target yield)?
The first question leads us to a Quadratic Programming problem, that is, a Mathematical Programming problem with a quadratic objective function and linear constraints. The second question necessitates the introduction of discrete variables to count the number of securities, and so we obtain a Mixed Integer Quadratic Programming problem. The two cases will be discussed separately in the following two sections.
QP
To adapt the model developed in Chapter 2 to the new way of looking at the problem, we need to make the following changes:
- New objective function: mean variance instead of total return.
- The risk-related constraint disappears.
- Addition of a new constraint: target yield.
The new objective function is the mean variance of the portfolio, namely:
where VARst is the variance/covariance matrix of all shares. This is a quadratic objective function (an objective function becomes quadratic either when a variable is squared, e.g., frac12, or when two variables are multiplied together, e.g., frac1 · frac2).
The target yield constraint can be written as follows:
The limit on the North-American shares as well as the requirement to spend all the money, and the upper bounds on the fraction invested into every share are retained. We therefore obtain the following complete mathematical model formulation:
∑s ∈ NA fracs ≥ 0.5
∑s ∈ SHARES fracs = 1
∑s ∈ SHARES RETs·fracs ≥ TARGET
∀ s ∈ SHARES: 0 ≤ fracs ≤ 0.3
Implementation with BCL
The estimated returns and the variance/covariance matrix are given in the data file foliocppqp.dat:
! trs haw thr tel brw hgw car bnk sof elc 0.1 0 0 0 0 0 0 0 0 0 ! treasury 0 19 -2 4 1 1 1 0.5 10 5 ! hardware 0 -2 28 1 2 1 1 0 -2 -1 ! theater 0 4 1 22 0 1 2 0 3 4 ! telecom 0 1 2 0 4 -1.5 -2 -1 1 1 ! brewery 0 1 1 1 -1.5 3.5 2 0.5 1 1.5 ! highways 0 1 1 2 -2 2 5 0.5 1 2.5 ! cars 0 0.5 0 0 -1 0.5 0.5 1 0.5 0.5 ! bank 0 10 -2 3 1 1 1 0.5 25 8 ! software 0 5 -1 4 1 1.5 2.5 0.5 8 16 ! electronics
We may read this datafile with the function XPRBreadarrlinecb: all comments preceded by ! and also empty lines are skipped. We read an entire line at once indicating the format of an entry (`g') and the separator (any number of spaces or tabulations).
For the definition of the objective function we now use a quadratic expression (equally represented by the class XPRBexpr). Since we now wish to minimize the problem, we use the default optimization sense setting and optimization as a continuous problem is again started with the method lpOptimize (with empty string argument indicating the default algorithm).
#include <iostream>
#include "xprb_cpp.h"
using namespace std;
using namespace ::dashoptimization;
#define DATAFILE "foliocppqp.dat"
#define TARGET 9 // Target yield
#define MAXNUM 4 // Max. number of different assets
#define NSHARES 10 // Number of shares
#define NNA 4 // Number of North-American shares
double RET[] = {5,17,26,12,8,9,7,6,31,21}; // Estimated return in investment
int NA[] = {0,1,2,3}; // Shares issued in N.-America
double VAR[NSHARES][NSHARES]; // Variance/covariance matrix of
// estimated returns
int main(int argc, char **argv)
{
int s,t;
XPRBprob p("FolioQP"); // Initialize a new problem in BCL
XPRBexpr Na,Return,Cap,Num,Variance;
XPRBvar frac[NSHARES]; // Fraction of capital used per share
FILE *datafile;
// Read `VAR' data from file
datafile=fopen(DATAFILE,"r");
for(s=0;s<NSHARES;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 200, "g ", VAR[s], NSHARES);
fclose(datafile);
// Create the decision variables
for(s=0;s<NSHARES;s++)
frac[s] = p.newVar(XPRBnewname("frac(%d)",s+1), XPRB_PL, 0, 0.3);
// Objective: mean variance
for(s=0;s<NSHARES;s++)
for(t=0;t<NSHARES;t++) Variance += VAR[s][t]*frac[s]*frac[t];
p.setObj(Variance); // Set the objective function
// Minimum amount of North-American values
for(s=0;s<NNA;s++) Na += frac[NA[s]];
p.newCtr(Na >= 0.5);
// Spend all the capital
for(s=0;s<NSHARES;s++) Cap += frac[s];
p.newCtr(Cap == 1);
// Target yield
for(s=0;s<NSHARES;s++) Return += RET[s]*frac[s];
p.newCtr(Return >= TARGET);
// Solve the problem
p.lpOptimize("");
// Solution printing
cout << "With a target of " << TARGET << " minimum variance is " <<
p.getObjVal() << endl;
for(s=0;s<NSHARES;s++)
cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
return 0;
}This program produces the following solution output with a dual-core processor (notice that the default algorithm for solving QP problems is Newton-Barrier, not the Simplex as in all previous examples):
Reading Problem FolioQP
Problem Statistics
3 ( 0 spare) rows
10 ( 0 spare) structural columns
24 ( 0 spare) non-zero elements
76 quadratic elements
Global Statistics
0 entities 0 sets 0 set members
Minimizing QP FolioQP
Original problem has:
3 rows 10 cols 24 elements
76 qobjelem
Presolved problem has:
3 rows 10 cols 24 elements
76 qobjelem
Barrier cache sizes : L1=32K L2=8192K
Using AVX support
Cores per CPU (CORESPERCPU): 8
Barrier starts, using up to 8 threads, 4 cores
Matrix ordering - Dense cols.: 9 NZ(L): 92 Flops: 584
Its P.inf D.inf U.inf Primal obj. Dual obj. Compl.
0 1.90e+001 1.85e+002 3.70e+000 8.7840000e+002 -1.1784000e+003 4.5e+003
1 1.69e-001 1.58e+000 3.29e-002 7.1810240e+000 -2.7042733e+002 3.1e+002
2 3.31e-003 1.48e-002 6.45e-004 5.1672666e+000 -1.2127681e+001 1.7e+001
3 6.12e-007 2.66e-015 2.78e-017 1.5558934e+000 -4.8803143e+000 6.4e+000
4 9.71e-017 1.39e-015 2.78e-017 7.2498306e-001 1.4062618e-001 5.8e-001
5 3.64e-017 1.01e-015 5.55e-017 5.6634270e-001 5.2690100e-001 3.9e-002
6 3.97e-017 7.15e-016 5.55e-017 5.5894833e-001 5.5591229e-001 3.0e-003
7 1.22e-016 1.33e-015 5.55e-017 5.5760205e-001 5.5726378e-001 3.4e-004
8 9.18e-017 1.68e-015 5.55e-017 5.5741308e-001 5.5738503e-001 2.8e-005
9 2.36e-016 3.29e-016 5.55e-017 5.5739403e-001 5.5739328e-001 7.5e-007
Barrier method finished in 0 seconds
Uncrunching matrix
Optimal solution found
Barrier solved problem
9 barrier iterations in 0s
Final objective : 5.573940299651456e-01
Max primal violation (abs / rel) : 7.347e-17 / 7.347e-17
Max dual violation (abs / rel) : 0.0 / 0.0
Max complementarity viol. (abs / rel) : 6.075e-07 / 1.012e-07
All values within tolerances
With a target of 9 minimum variance is 0.557394
0: 30%
1: 7.15401%
2: 7.38237%
3: 5.46362%
4: 12.6561%
5: 5.91283%
6: 0.333491%
7: 29.9979%
8: 1.0997%
9: 6.97039e-06%MIQP
We now wish to express the fact that at most a given number MAXNUM of different assets may be selected into the portfolio, subject to all other constraints of the previous QP model. In Chapter 11 we have already seen how this can be done, namely by introducing an additional set of binary decision variables buys that are linked logically to the continuous variables:
Through this relation, a variable buys will be at 1 if a fraction fracs greater than 0 is selected into the portfolio. If, however, buys equals 0, then fracs must also be 0.
To limit the number of different shares in the portfolio, we then define the following constraint:
Implementation with BCL
We may modify the previous QP model or simply append the following lines to the program of the previous section, just after the solution printing: the problem is then solved once as a QP and once as a MIQP in a single program run.
XPRBvar buy[NSHARES]; // 1 if asset is in portfolio, 0 otherwise
// Create the decision variables
for(s=0;s<NSHARES;s++)
buy[s] = p.newVar(XPRBnewname("buy(%d)",s+1), XPRB_BV);
// Limit the total number of assets
for(s=0;s<NSHARES;s++) Num += buy[s];
p.newCtr(Num <= MAXNUM);
// Linking the variables
for(s=0;s<NSHARES;s++) p.newCtr(frac[s] <= buy[s]);
// Solve the problem
p.mipOptimize("");
// Solution printing
cout << "With a target of " << TARGET << " and at most " << MAXNUM <<
" assets, minimum variance is " << p.getObjVal() << endl;
for(s=0;s<NSHARES;s++)
cout << s << ": " << frac[s].getSol()*100 << "% (" << buy[s].getSol()
<< ")" << endl;When executing the MIQP model, we obtain the following solution output:
Reading Problem FolioQP
Problem Statistics
14 ( 514 spare) rows
20 ( 0 spare) structural columns
54 ( 5056 spare) non-zero elements
76 quadratic elements
Global Statistics
10 entities 0 sets 0 set members
Minimizing MIQP FolioQP
Original problem has:
14 rows 20 cols 54 elements 10 globals
76 qobjelem
Presolved problem has:
14 rows 20 cols 54 elements 10 globals
76 qobjelem
LP relaxation tightened
Will try to keep branch and bound tree memory usage below 14.8Gb
Crash basis containing 0 structural columns created
Its Obj Value S Ninf Nneg Sum Inf Time
0 .000000 p 1 4 .100000 0
8 .000000 p 0 0 .000000 0
8 4.609000 p 0 0 .000000 0
Its Obj Value S Nsft Nneg Dual Inf Time
27 .557393 QP 0 0 .000000 0
QP solution found
Optimal solution found
Primal solved problem
27 simplex iterations in 0s
Final objective : 5.573934108103899e-01
Max primal violation (abs / rel) : 1.804e-16 / 1.804e-16
Max dual violation (abs / rel) : 1.776e-15 / 1.776e-15
Max complementarity viol. (abs / rel) : 2.670e-16 / 1.007e-16
All values within tolerances
Starting root cutting & heuristics
Its Type BestSoln BestBound Sols Add Del Gap GInf Time
a 4.094715 .557393 1 86.39
b 1.839000 .557393 2 69.69
q 1.825619 .557393 3 69.47
k 1.419003 .557393 4 60.72
1 K 1.419003 .557393 4 3 0 60.72 0
2 K 1.419003 .557393 4 9 2 60.72 0
3 K 1.419003 .557393 4 7 6 60.72 0
4 K 1.419003 .557393 4 5 6 60.72 0
5 K 1.419003 .557393 4 11 5 60.72 0
6 K 1.419003 .558122 4 8 10 60.67 0
7 K 1.419003 .570670 4 11 9 59.78 0
8 K 1.419003 .570670 4 5 12 59.78 0
9 K 1.419003 .583638 4 5 3 58.87 0
10 K 1.419003 .612496 4 4 0 56.84 0
11 K 1.419003 .618043 4 6 5 56.45 0
12 K 1.419003 .620360 4 8 4 56.28 0
13 K 1.419003 .620360 4 0 6 56.28 0
Heuristic search started
Heuristic search stopped
Cuts in the matrix : 14
Cut elements in the matrix : 138
Starting tree search.
Deterministic mode with up to 8 running threads and up to 16 tasks.
Node BestSoln BestBound Sols Active Depth Gap GInf Time
a 9 1.248762 .919281 5 2 5 26.38
*** Search completed *** Time: 0 Nodes: 15
Number of integer feasible solutions found is 5
Best integer solution found is 1.248762
Best bound is 1.248752
Uncrunching matrix
With a target of 9 and at most 4 assets, minimum variance is 1.24876
0: 30% (1)
1: 20% (1)
2: 0% (0)
3: 0% (0)
4: 23.8095% (1)
5: 26.1905% (1)
6: 0% (0)
7: 0% (0)
8: 0% (0)
9: 0% (0)The log of the Branch-and-Bound search tells us this time that 5 integer feasible solutions have been found (all by the MIP heuristics) and a total of 15 nodes have been enumerated to complete the search.With the additional constraint on the number of different assets the minimum variance is more than twice as large as in the QP problem.