6.10 Quadratic Optimization

MOSEK can solve quadratic and quadratically constrained problems, as long as they are convex. This class of problems can be formulated as follows:

(6.40)\[\begin{split}\begin{array}{lrcccll} \mbox{minimize} & & & \half x^T Q^o x + c^T x + c^f & & & \\ \mbox{subject to} & l_k^c & \leq & \half x^T Q^k x + \sum_{j=0}^{n-1} a_{k,j} x_j & \leq & u_k^c, & k =0,\ldots ,m-1, \\ & l_j^x & \leq & x_j & \leq & u_j^x, & j=0,\ldots ,n-1. \end{array}\end{split}\]

Without loss of generality it is assumed that \(Q^o\) and \(Q^k\) are all symmetric because

\[x^T Q x = \half x^T(Q+Q^T)x.\]

This implies that a non-symmetric \(Q\) can be replaced by the symmetric matrix \(\half(Q+Q^T)\).

The problem is required to be convex. More precisely, the matrix \(Q^o\) must be positive semi-definite and the \(k\)th constraint must be of the form

(6.41)\[ l_k^c \leq \half x^T Q^k x + \sum_{j=0}^{n-1} a_{k,j} x_j\]

with a negative semi-definite \(Q^k\) or of the form

\[\half x^T Q^k x + \sum_{j=0}^{n-1} a_{k,j} x_j \leq u_k^c.\]

with a positive semi-definite \(Q^k\). This implies that quadratic equalities are not allowed. Specifying a non-convex problem will result in an error when the optimizer is called.

A matrix is positive semidefinite if all the eigenvalues of \(Q\) are nonnegative. An alternative statement of the positive semidefinite requirement is

\[x^T Q x \geq 0, \quad \forall x.\]

If the convexity (i.e. semidefiniteness) conditions are not met MOSEK will not produce reliable results or work at all.

6.10.1 Example: Quadratic Objective

We look at a small problem with linear constraints and quadratic objective:

(6.42)\[\begin{split}\begin{array}{lll} \mbox{minimize} & & x_1^2 + 0.1 x_2^2 + x_3^2 - x_1 x_3 - x_2 \\ \mbox{subject to} & 1 \leq & x_1 + x_2 + x_3 \\ & 0 \leq & x. \end{array}\end{split}\]

The matrix formulation of (6.42) has:

\[\begin{split}Q^o = \left[ \begin{array}{ccc} 2 & 0 & -1\\ 0 & 0.2 & 0\\ -1 & 0 & 2 \end{array} \right], c = \left[ \begin{array}c 0\\ -1\\ 0 \end{array} \right], A = \left[ \begin{array} {ccc} 1 & 1 & 1 \end{array} \right],\end{split}\]

with the bounds:

\[\begin{split}l^c = 1, u^c = \infty , l^x = \left[ \begin{array}c 0 \\ 0 \\ 0 \end{array} \right] \mbox{ and } u^x = \left[ \begin{array} c \infty \\ \infty \\ \infty \end{array} \right]\end{split}\]

Please note the explicit \(\half\) in the objective function of (6.40) which implies that diagonal elements must be doubled in \(Q\), i.e. \(Q_{11}=2\) even though \(1\) is the coefficient in front of \(x_1^2\) in (6.42).

Setting up the linear part

The linear parts (constraints, variables, objective) are set up using exactly the same methods as for linear problems, and we refer to Sec. 6.1 (Linear Optimization) for all the details. The same applies to technical aspects such as defining an optimization task, retrieving the solution and so on.

Setting up the quadratic objective

The quadratic objective is specified using the function Task.putqobj. Since \(Q^o\) is symmetric only the lower triangular part of \(Q^o\) is inputted. In fact entries from above the diagonal may not appear in the input.

The lower triangular part of the matrix \(Q^o\) is specified using an unordered sparse triplet format (for details, see Sec. 15.1.4 (Matrix Formats)):

      int[]    qsubi = {0,   1,    2,   2  };
      int[]    qsubj = {0,   1,    0,   2  };
      double[] qval =  {2.0, 0.2, -1.0, 2.0};

Please note that

  • only non-zero elements are specified (any element not specified is 0 by definition),

  • the order of the non-zero elements is insignificant, and

  • only the lower triangular part should be specified.

Finally, this definition of \(Q^o\) is loaded into the task:

      task.putqobj(qsubi, qsubj, qval);

Source code

Listing 6.18 Source code implementing problem (6.42). Click here to download.
package com.mosek.example;
import mosek.*;

public class qo1 {
  static final int numcon = 1;   /* Number of constraints.             */
  static final int numvar = 3;   /* Number of variables.               */
  static final int NUMANZ = 3;   /* Number of numzeros in A.           */
  static final int NUMQNZ = 4;   /* Number of nonzeros in Q.           */

  public static void main (String[] args) {
    // Since the value infinity is never used, we define
    // 'infinity' symbolic purposes only
    double infinity = 0;
    double[] c = {0.0, -1.0, 0.0};

    mosek.boundkey[] bkc  = { mosek.boundkey.lo };
    double[] blc = {1.0};
    double[] buc = {infinity};

    mosek.boundkey[] bkx  = { mosek.boundkey.lo,
                              mosek.boundkey.lo,
                              mosek.boundkey.lo
                            };
    double[] blx = {0.0,
                    0.0,
                    0.0
                   };
    double[] bux = {infinity,
                    infinity,
                    infinity
                   };

    int[][]    asub  = { {0},   {0}, {0} };
    double[][] aval  = { {1.0}, {1.0}, {1.0} };
    double[] xx   = new double[numvar];

    try (Env  env  = new Env();
         Task task = new Task(env, 0, 0)) {
      // Directs the log task stream to the user specified
      // method task_msg_obj.stream
      task.set_Stream(
        mosek.streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});
      /* Give MOSEK an estimate of the size of the input data.
      This is done to increase the speed of inputting data.
      However, it is optional. */
      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      task.appendcons(numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      task.appendvars(numvar);

      for (int j = 0; j < numvar; ++j) {
        /* Set the linear term c_j in the objective.*/
        task.putcj(j, c[j]);
        /* Set the bounds on variable j.
           blx[j] <= x_j <= bux[j] */
        task.putvarbound(j, bkx[j], blx[j], bux[j]);
        /* Input column j of A */
        task.putacol(j,                     /* Variable (column) index.*/
                     asub[j],               /* Row index of non-zeros in column j.*/
                     aval[j]);              /* Non-zero Values of column j. */
      }
      /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
      for (int i = 0; i < numcon; ++i)
        task.putconbound(i, bkc[i], blc[i], buc[i]);

      /*
       The lower triangular part of the Q
       matrix in the objective is specified.
      */

      int[]    qsubi = {0,   1,    2,   2  };
      int[]    qsubj = {0,   1,    0,   2  };
      double[] qval =  {2.0, 0.2, -1.0, 2.0};

      /* Input the Q for the objective. */

      task.putqobj(qsubi, qsubj, qval);

      /* Solve the problem */
      mosek.rescode r = task.optimize();
      System.out.println (" Mosek warning:" + r.toString());
      // Print a summary containing information
      //   about the solution for debugging purposes
      task.solutionsummary(mosek.streamtype.msg);

      mosek.solsta solsta[] = new mosek.solsta[1];
      /* Get status information about the solution */
      task.getsolsta(mosek.soltype.itr, solsta);

      /* Get the solution */
      task.getxx(mosek.soltype.itr, // Interior solution.
                 xx);

      switch (solsta[0]) {
        case optimal:
          System.out.println("Optimal primal solution\n");
          for (int j = 0; j < numvar; ++j)
            System.out.println ("x[" + j + "]:" + xx[j]);
          break;
        case dual_infeas_cer:
        case prim_infeas_cer:
          System.out.println("Primal or dual infeasibility\n");
          break;
        case unknown:
          System.out.println("Unknown solution status.\n");
          break;
        default:
          System.out.println("Other solution status");
          break;
      }
    }
    catch (mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.toString());
      throw e;
    }
  } /* Main */
}

6.10.2 Example: Quadratic constraints

In this section we show how to solve a problem with quadratic constraints. Please note that quadratic constraints are subject to the convexity requirement (6.41).

Consider the problem:

\[\begin{split}\begin{array}{lcccl} \mbox{minimize} & & & x_1^2 + 0.1 x_2^2 + x_3^2 - x_1 x_3 - x_2 & \\ \mbox{subject to} & 1 & \leq & x_1 + x_2 + x_3 - x_1^2 - x_2^2 - 0.1 x_3^2 + 0.2 x_1 x_3, & \\ & & & x \geq 0. & \end{array}\end{split}\]

This is equivalent to

(6.43)\[\begin{split}\begin{array}{lccl} \mbox{minimize} & \half x^T Q^o x + c^T x & & \\ \mbox{subject to} & \half x^T Q^0 x + A x & \geq & b, \\ & x\geq 0, \end{array}\end{split}\]

where

\[\begin{split}Q^o = \left[ \begin{array}{ccc} 2 & 0 & -1 \\ 0 & 0.2 & 0 \\ -1 & 0 & 2 \end{array} \right], c = \left[ \begin{array}{ccc} 0 &-1 & 0 \end{array} \right]^T, A = \left[ \begin{array}{ccc} 1 & 1 & 1 \end{array} \right], b = 1.\end{split}\]
\[\begin{split}Q^0 = \left[ \begin{array}{ccc} -2 & 0 & 0.2 \\ 0 & -2 & 0 \\ 0.2 & 0 & -0.2 \end{array} \right].\end{split}\]

The linear parts and quadratic objective are set up the way described in the previous tutorial.

Setting up quadratic constraints

To add quadratic terms to the constraints we use the function Task.putqconk.

      int[]    qsubi = {0,   1,    2,   2  };
      int[]    qsubj = {0,   1,    2,   0  };
      double[] qval =  { -2.0, -2.0, -0.2, 0.2};

      /* put Q^0 in constraint with index 0. */

      task.putqconk (0,
                     qsubi,
                     qsubj,
                     qval);

While Task.putqconk adds quadratic terms to a specific constraint, it is also possible to input all quadratic terms in one chunk using the Task.putqcon function.

Source code

Listing 6.19 Implementation of the quadratically constrained problem (6.43). Click here to download.
package com.mosek.example;
import mosek.*;

public class qcqo1 {
  static final int numcon = 1;   /* Number of constraints.             */
  static final int numvar = 3;   /* Number of variables.               */
  static final int NUMANZ = 3;   /* Number of numzeros in A.           */
  static final int NUMQNZ = 4;   /* Number of nonzeros in Q.           */

  public static void main (String[] args) {
    // Since the value infinity is never used, we define
    // 'infinity' symbolic purposes only
    double infinity = 0;
    double[] c = {0.0, -1.0, 0.0};

    mosek.boundkey[]    bkc  = {mosek.boundkey.lo};
    double[] blc = {1.0};
    double[] buc = {infinity};

    mosek.boundkey[]  bkx
    = {mosek.boundkey.lo,
       mosek.boundkey.lo,
       mosek.boundkey.lo
      };
    double[] blx = {0.0,
                    0.0,
                    0.0
                   };
    double[] bux = {infinity,
                    infinity,
                    infinity
                   };

    int[][]    asub  = { {0},   {0},   {0} };
    double[][] aval  = { {1.0}, {1.0}, {1.0} };

    double[] xx   = new double[numvar];


    try (mosek.Env  env  = new mosek.Env();
         mosek.Task task = new mosek.Task(env, 0, 0)) {
      // Directs the log task stream to the user specified
      // method task_msg_obj.stream
      task.set_Stream(
        mosek.streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});

      /* Give MOSEK an estimate of the size of the input data.
      This is done to increase the speed of inputting data.
      However, it is optional. */
      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      task.appendcons(numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      task.appendvars(numvar);

      for (int j = 0; j < numvar; ++j) {
        /* Set the linear term c_j in the objective.*/
        task.putcj(j, c[j]);
        /* Set the bounds on variable j.
         blx[j] <= x_j <= bux[j] */
        task.putvarbound(j, bkx[j], blx[j], bux[j]);
        /* Input column j of A */
        task.putacol(j,                     /* Variable (column) index.*/
                     asub[j],               /* Row index of non-zeros in column j.*/
                     aval[j]);              /* Non-zero Values of column j. */
      }
      /* Set the bounds on constraints.
      for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
      for (int i = 0; i < numcon; ++i)
        task.putconbound(i, bkc[i], blc[i], buc[i]);
      /*
       * The lower triangular part of the Q
       * matrix in the objective is specified.
       */

      int[]   qosubi = { 0,   1,   2,    2 };
      int[]   qosubj = { 0,   1,   0,    2 };
      double[] qoval = { 2.0, 0.2, -1.0, 2.0 };

      /* Input the Q for the objective. */

      task.putqobj(qosubi, qosubj, qoval);

      /*
       * The lower triangular part of the Q^0
       * matrix in the first constraint is specified.
       * This corresponds to adding the term
       * x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2
       */

      int[]    qsubi = {0,   1,    2,   2  };
      int[]    qsubj = {0,   1,    2,   0  };
      double[] qval =  { -2.0, -2.0, -0.2, 0.2};

      /* put Q^0 in constraint with index 0. */

      task.putqconk (0,
                     qsubi,
                     qsubj,
                     qval);

      task.putobjsense(mosek.objsense.minimize);

      /* Solve the problem */

      try {
        mosek.rescode termcode = task.optimize();
      } catch (mosek.Warning e) {
        System.out.println (" Mosek warning:");
        System.out.println (e.toString ());
      }
      // Print a summary containing information
      //   about the solution for debugging purposes
      task.solutionsummary(mosek.streamtype.msg);

      mosek.solsta solsta[] = new mosek.solsta[1];
      /* Get status information about the solution */
      task.getsolsta(mosek.soltype.itr, solsta);

      task.getxx(mosek.soltype.itr, // Interior solution.
                 xx);
      switch (solsta[0]) {
        case optimal:
          System.out.println("Optimal primal solution\n");
          for (int j = 0; j < numvar; ++j)
            System.out.println ("x[" + j + "]:" + xx[j]);
          break;
        case dual_infeas_cer:
        case prim_infeas_cer:
          System.out.println("Primal or dual infeasibility.\n");
          break;
        case unknown:
          System.out.println("Unknown solution status.\n");
          break;
        default:
          System.out.println("Other solution status");
          break;
      }


    }
    catch (mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.msg);
      throw e;
    }
  } /* Main */
}