6.6 Semidefinite Optimization

Semidefinite optimization is a generalization of conic optimization, allowing the use of matrix variables belonging to the convex cone of positive semidefinite matrices

\[\PSD^r = \left\lbrace X \in \Symm^r: z^T X z \geq 0, \quad \forall z \in \real^r \right\rbrace,\]

where \(\Symm^r\) is the set of \(r \times r\) real-valued symmetric matrices.

MOSEK can solve semidefinite optimization problems of the form

\[\begin{split}\begin{array}{lccccll} \mbox{minimize} & & & \sum_{j=0}^{n-1} c_j x_j + \sum_{j=0}^{p-1} \left\langle \barC_j, \barX_j \right\rangle + c^f & & &\\ \mbox{subject to} & l_i^c & \leq & \sum_{j=0}^{n-1} a_{ij} x_j + \sum_{j=0}^{p-1} \left\langle \barA_{ij}, \barX_j \right\rangle & \leq & u_i^c, & i = 0, \ldots, m-1,\\ & l_j^x & \leq & x_j & \leq & u_j^x, & j = 0, \ldots, n-1,\\ & & & x \in \K, \barX_j \in \PSD^{r_j}, & & & j = 0, \ldots, p-1 \end{array}\end{split}\]

where the problem has \(p\) symmetric positive semidefinite variables \(\barX_j\in \PSD^{r_j}\) of dimension \(r_j\) with symmetric coefficient matrices \(\barC_j\in \Symm^{r_j}\) and \(\barA_{i,j}\in \Symm^{r_j}\). We use standard notation for the matrix inner product, i.e., for \(A,B\in \real^{m\times n}\) we have

\[\left\langle A,B \right\rangle := \sum_{i=0}^{m-1} \sum_{j=0}^{n-1} A_{ij} B_{ij}.\]

We demonstrate the setup of semidefinite variables and the matrices \(\barC\), \(\barA\) on the following examples:

6.6.1 Example SDO1

We consider the simple optimization problem with semidefinite and conic quadratic constraints:

(6.11)\[\begin{split}\begin{array} {llcc} \mbox{minimize} & \left\langle \left[ \begin{array} {ccc} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{array} \right], \barX \right\rangle + x_0 & & \\ \mbox{subject to} & \left\langle \left[ \begin{array} {ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right], \barX \right\rangle + x_0 & = & 1, \\ & \left\langle \left[ \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{array} \right], \barX \right\rangle + x_1 + x_2 & = & 1/2, \\ & x_0 \geq \sqrt{{x_1}^2 + {x_2}^2}, & \barX \succeq 0, & \end{array}\end{split}\]

The problem description contains a 3-dimensional symmetric semidefinite variable which can be written explicitly as:

\[\begin{split}\barX = \left[ \begin{array} {ccc} \barX_{00} & \barX_{10} & \barX_{20} \\ \barX_{10} & \barX_{11} & \barX_{21} \\ \barX_{20} & \barX_{21} & \barX_{22} \end{array} \right] \in \PSD^3,\end{split}\]

and a conic quadratic variable \((x_0, x_1, x_2) \in \Q^3\). The objective is to minimize

\[2(\barX_{00} + \barX_{10} + \barX_{11} + \barX_{21} + \barX_{22}) + x_0,\]

subject to the two linear constraints

\[\begin{split}\begin{array}{ccc} \barX_{00} + \barX_{11} + \barX_{22} + x_0 & = & 1, \\ \barX_{00} + \barX_{11} + \barX_{22} + 2(\barX_{10} + \barX_{20} + \barX_{21}) + x_1 + x_2 & = & 1/2. \end{array}\end{split}\]

Setting up the linear and conic part

The linear and conic parts (constraints, variables, objective, cones) are set up using the methods described in the relevant tutorials; Sec. 6.1 (Linear Optimization), Sec. 6.3 (Conic Quadratic Optimization), Sec. 6.5 (Conic Exponential Optimization), Sec. 6.4 (Power Cone Optimization). Here we only discuss the aspects directly involving semidefinite variables.

Appending semidefinite variables

First, we need to declare the number of semidefinite variables in the problem, similarly to the number of linear variables and constraints. This is done with the function Task.appendbarvars.

      task.appendbarvars(dimbarvar);

Appending coefficient matrices

Coefficient matrices \(\barC_j\) and \(\barA_{ij}\) are constructed as weighted combinations of sparse symmetric matrices previously appended with the function Task.appendsparsesymmat.

        idx[0] = task.appendsparsesymmat(dimbarvar[0],
                                         barc_i,
                                         barc_j,
                                         barc_v);

The arguments specify the dimension of the symmetric matrix, followed by its description in the sparse triplet format. Only lower-triangular entries should be included. The function produces a unique index of the matrix just entered in the collection of all coefficient matrices defined by the user.

After one or more symmetric matrices have been created using Task.appendsparsesymmat, we can combine them to set up the objective matrix coefficient \(\barC_j\) using Task.putbarcj, which forms a linear combination of one or more symmetric matrices. In this example we form the objective matrix directly, i.e. as a weighted combination of a single symmetric matrix.

        task.putbarcj(0, idx, falpha);

Similarly, a constraint matrix coefficient \(\barA_{ij}\) is set up by the function Task.putbaraij.

        task.putbaraij(0, 0, idx, falpha);

Retrieving the solution

After the problem is solved, we read the solution using Task.getbarxj:

          task.getbarxj(mosek.soltype.itr,    /* Request the interior solution. */
                        0,
                        barx);

The function returns the half-vectorization of \(\barX_j\) (the lower triangular part stacked as a column vector), where the semidefinite variable index \(j\) is passed as an argument.

Source code

Listing 6.7 Source code solving problem (6.11). Click here to download.
package com.mosek.example;
import mosek.*;

public class sdo1 {
  public static void main(String[] argv) {
    int    numcon      = 2;  /* Number of constraints.              */
    int    numvar      = 3;  /* Number of conic quadratic variables */
    int    dimbarvar[] = {3};         /* Dimension of semidefinite cone */
    int    lenbarvar[] = {3 * (3 + 1) / 2}; /* Number of scalar SD variables  */

    mosek.boundkey bkc[] = { mosek.boundkey.fx,
                             mosek.boundkey.fx
                           };
    double[]     blc     = { 1.0, 0.5 };
    double[]     buc     = { 1.0, 0.5 };

    int[]        barc_i  = {0, 1, 1, 2, 2},
                 barc_j  = {0, 0, 1, 1, 2};
    double[]     barc_v  = {2.0, 1.0, 2.0, 1.0, 2.0};

    int[][]      asub    = {{0}, {1, 2}}; /* column subscripts of A */
    double[][]   aval    = {{1.0}, {1.0, 1.0}};

    int[][]      bara_i  = { {0,   1,   2},   {0,   1 ,  2,   1,   2,   2 } },
                 bara_j  = { {0,   1,   2},   {0,   0 ,  0,   1,   1,   2 } };
    
    double[][]   bara_v  = { {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}};
    int[]        conesub = { 0, 1, 2};


    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); }});

      /* 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);

      /* Append 'NUMBARVAR' semidefinite variables. */
      task.appendbarvars(dimbarvar);

      /* Optionally add a constant term to the objective. */
      task.putcfix(0.0);

      /* Set the linear term c_j in the objective.*/
      task.putcj(0, 1.0);

      for (int j = 0; j < numvar; ++j)
        task.putvarbound(j, mosek.boundkey.fr, -0.0, 0.0);

      /* Set the linear term barc_j in the objective.*/
      {
        long[] idx = new long[1];
        double[] falpha = { 1.0 };
        idx[0] = task.appendsparsesymmat(dimbarvar[0],
                                         barc_i,
                                         barc_j,
                                         barc_v);
        task.putbarcj(0, idx, falpha);
      }

      /* 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,           /* Index of constraint.*/
                         bkc[i],      /* Bound key.*/
                         blc[i],      /* Numerical value of lower bound.*/
                         buc[i]);     /* Numerical value of upper bound.*/

      /* Input A row by row */
      for (int i = 0; i < numcon; ++i)
        task.putarow(i,
                     asub[i],
                     aval[i]);

      /* Append the conic quadratic cone */
      task.appendcone(mosek.conetype.quad,
                      0.0,
                      conesub);

      /* Add the first row of barA */
      {
        long[] idx = new long[1];
        double[] falpha = {1.0};
        task.appendsparsesymmat(dimbarvar[0],
                                bara_i[0],
                                bara_j[0],
                                bara_v[0],
                                idx);

        task.putbaraij(0, 0, idx, falpha);
      }

      {
        long[] idx = new long[1];
        double[] falpha = {1.0};
        /* Add the second row of barA */
        task.appendsparsesymmat(dimbarvar[0],
                                bara_i[1],
                                bara_j[1],
                                bara_v[1],
                                idx);

        task.putbaraij(1, 0, idx, falpha);
      }

      /* Run optimizer */
      task.optimize();

      /* Print a summary containing information
         about the solution for debugging purposes*/
      task.solutionsummary (mosek.streamtype.msg);

      mosek.solsta[] solsta = new mosek.solsta[1];
      task.getsolsta (mosek.soltype.itr, solsta);

      switch (solsta[0]) {
        case optimal:
          double[] xx = new double[numvar];
          double[] barx = new double[lenbarvar[0]];

          task.getxx(mosek.soltype.itr, xx);
          task.getbarxj(mosek.soltype.itr,    /* Request the interior solution. */
                        0,
                        barx);
          System.out.println("Optimal primal solution");
          for (int i = 0; i < numvar; ++i)
            System.out.println("x[" + i + "]   : " + xx[i]);

          for (int i = 0; i < lenbarvar[0]; ++i)
            System.out.println("barx[" + i + "]: " + barx[i]);
          break;
        case dual_infeas_cer:
        case prim_infeas_cer:
          System.out.println("Primal or dual infeasibility certificate found.");
          break;
        case unknown:
          System.out.println("The status of the solution could not be determined.");
          break;
        default:
          System.out.println("Other solution status.");
          break;
      }
    }
  }
}

6.6.2 Example SDO2

We now demonstrate how to define more than one semidefinite variable using the following problem with two matrix variables and two types of constraints:

(6.12)\[\begin{split}\begin{array}{lrll} \mbox{minimize} & \langle C_1,\barX_1\rangle + \langle C_2,\barX_2\rangle & & \\ \mbox{subject to} & \langle A_1,\barX_1\rangle + \langle A_2,\barX_2\rangle & = & b, \\ & (\barX_2)_{01} & \leq & k, \\ & \barX_1, \barX_2 & \succeq & 0. \end{array}\end{split}\]

In our example \(\dim(\barX_1)=3\), \(\dim(\barX_2)=4\), \(b=23\), \(k=-3\) and

\[\begin{split}C_1= \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 6 \end{array}\right], A_1= \left[\begin{array}{ccc} 1 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 2 \end{array}\right],\end{split}\]
\[\begin{split}C_2= \left[\begin{array}{cccc} 1 & -3 & 0 & 0\\ -3 & 2 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array}\right], A_2= \left[\begin{array}{cccc} 0 & 1 & 0 & 0\\ 1 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -3 \\ \end{array}\right],\end{split}\]

are constant symmetric matrices.

Note that this problem does not contain any scalar variables, but they could be added in the same fashion as in Sec. 6.6.1 (Example SDO1).

Other than in Sec. 6.6.1 (Example SDO1) we don’t append coefficient matrices separately but we directly input all nonzeros in each constraint and all nonzeros in the objective at once. Every term of the form \((\barA_{i,j})_{k,l}(\barX_j)_{k,l}\) is determined by four indices \((i,j,k,l)\) and a coefficient value \(v=(\barA_{i,j})_{k,l}\). Here \(i\) is the number of the constraint in which the term appears, \(j\) is the index of the semidefinite variable it involves and \((k,l)\) is the position in that variable. This data is passed in the call to Task.putbarablocktriplet. Note that only the lower triangular part should be specified explicitly, that is one always has \(k\geq l\). Semidefinite terms \((\barC_j)_{k,l}(\barX_j)_{k,l}\) of the objective are specified in the same way in Task.putbarcblocktriplet but only include \((j,k,l)\) and \(v\).

For explanations of other data structures used in the example see Sec. 6.6.1 (Example SDO1).

The code representing the above problem is shown below.

Listing 6.8 Implementation of model (6.12). Click here to download.
    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); }});

      /* Append numcon empty constraints.
         The constraints will initially have no bounds. */
      task.appendcons(numcon);

      /* Append numbarvar semidefinite variables. */
      task.appendbarvars(dimbarvar);

      /* Set objective (6 nonzeros).*/
      task.putbarcblocktriplet(6, Cj, Ck, Cl, Cv);

      /* Set the equality constraint (6 nonzeros).*/
      task.putbarablocktriplet(6, Ai, Aj, Ak, Al, Av);

      /* Set the inequality constraint (1 nonzero).*/
      task.putbarablocktriplet(1, A2i, A2j, A2k, A2l, A2v);

      /* Set constraint bounds */
      task.putconboundslice(0, 2, bkc, blc, buc);

      /* Run optimizer */
      task.optimize();
      task.solutionsummary(mosek.streamtype.msg);

      mosek.solsta[] solsta = new mosek.solsta[1];
      task.getsolsta (mosek.soltype.itr, solsta);

      switch (solsta[0]) {
        case optimal:

          /* Retrieve the soution for all symmetric variables */
          System.out.println("Solution (lower triangular part vectorized):");
          for(int i = 0; i < numbarvar; i++) {
            int dim = dimbarvar[i] * (dimbarvar[i] + 1) / 2;
            double[] barx = new double[dim];

            task.getbarxj(mosek.soltype.itr, i, barx);

            System.out.print("X" + (i+1) + ": ");
            for (int j = 0; j < dim; ++j)
              System.out.print(barx[j] + " ");
            System.out.println();
          }

          break;
        case dual_infeas_cer:
        case prim_infeas_cer:
          System.out.println("Primal or dual infeasibility certificate found.");
          break;
        case unknown:
          System.out.println("The status of the solution could not be determined.");
          break;
        default:
          System.out.println("Other solution status.");
          break;
      }
    }
  }
}