Thursday, April 13, 2017

QP as LP: cutting planes

In (1) I described a simple linearization scheme for a QP model we we can solve it as a straight LP. For a simple (convex) quadratic function \(f(x)\), instead of solving:

\[ \min\> f(x)=x^2\]

we solve

\[\begin{align} \min\>&z\\
                              &z \ge a_i x + b_i\end{align}\]

In this post I do things slightly different: instead of adding the linear inequalities ahead of time, we add them one at the time based on the previously found optimal point. This approach is called a cutting plane technique (2).

Example: portfolio model

We consider again the simple portfolio model:

\[\boxed{\begin{align} \min \>& \sum_t w_t^2\\
                              & w_t = \sum_i r’_{i,t} x_i\\
                              & \sum_i \mu_i x_i \ge R\\
                              & \sum_i x_i = 1\\
                              &x_i\ge 0\end{align}}
\]

The model is initially linearized as an LP model:

\[\boxed{\begin{align}\min \>& \sum_t z_t\\
                              & w_t = \sum_i r’_{i,t} x_i\\
                              & \sum_i \mu_i x_i \ge R\\
                              & \sum_i x_i = 1\\
                              &x_i, z_t\ge 0\end{align}}
\]

Later on, during the CP algorithm, we will add linear cuts of the form

\[z_t \ge a_t w_t + b_t\]

The algorithm would look like:

  1. \(k:=1\)
  2. Solve model LP, let \(w^*_t\) be the optimal values for \(w_t\).
  3. if \(k=\)MAXIT: STOP
  4. Add the constraint \(z_t \ge a_t w_t + b_t\) where \(a_t=2 w^*_t\) and \(b_t=-(w^*_t)^2\). Note that we add one cut for each \(w_t\) here (our dataset has \(T=717\) time periods). 
  5. \(k:=k+1\)
  6. go to step 2

Here we stop simply when a certain number of iterations MAXIT has been exceeded. That can be refined by stopping when the objective does not seem to change much anymore. Another optimization could be to only add cuts that are different from the ones added before (for some \(t\) we may converge quicker than for others).

The algorithm converges very fast:

----    118 PARAMETER report2 

           obj(LP)         w'w

iter1               1.02907546
iter2   0.21577087  0.46879016
iter3   0.33210537  0.48432099
iter4   0.37143835  0.41397603
iter5   0.40990988  0.41362293
iter6   0.41109694  0.41222051
iter7   0.41183602  0.41212331
iter8   0.41192720  0.41204267
iter9   0.41196991  0.41204112
iter10  0.41197620  0.41203782

Note that the optimal QP solution has an objective: 0.41202816.

This is pretty good performance especially because the dataset is not very small (it has 717 time periods and 83 stocks).

Here is a picture of the cuts introduced for the first element \(z_1=w_1^2\):

animation

A combined approach

We can even combine the two methods:

  1. Start with a coarse-grained (i.e. cheap) initial set of cuts. In (1) we used \(n=10\) inequalities per \(w_t\). For this experiment I reduced this to \(n=5\).
  2. Then apply our cutting plane algorithm. Instead of MAXIT=10 iterations we now do 5 iterations.

This yields even faster convergence:

----    142 PARAMETER report3 

          obj(LP)         w'w

iter1  0.32921509  0.41401219
iter2  0.40812966  0.41471624
iter3  0.41017228  0.41211351
iter4  0.41188869  0.41208158
iter5  0.41195624  0.41203588

References
  1. QP as LP: piecewise linear functions, http://yetanothermathprogrammingconsultant.blogspot.com/2017/04/qp-as-lp-piecewise-linear-functions.html
  2. J. E. Kelley, Jr. "The Cutting-Plane Method for Solving Convex Programs." J. Soc. Indust. Appl. Math. Vol. 8, No. 4, pp. 703-712, December, 1960.
  3. Cardinality Constrained Portfolio Optimization: MIQP without MIQP Solver, http://yetanothermathprogrammingconsultant.blogspot.com/2016/02/cardinality-constrained-portfolio.html. Here this cutting plane method is applied to a MIQP model (not strictly “allowed” as this is no longer convex, but useful as a heuristic).