6.4 Semidefinite Optimization

Semidefinite optimization is a generalization of conic quadratic 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}.\]

6.4.1 Example SDO1

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

(1)\[\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 quadratic part

The linear and quadratic parts (constraints, variables, objective, cones) are set up using the methods described in the relevant tutorials; Sec. 6.1 (Linear Optimization) and Sec. 6.3 (Conic Quadratic 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 MSK_appendbarvars.

        r = MSK_appendbarvars(task, NUMBARVAR, 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 MSK_appendsparsesymmat.

        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   5,
                                   barc_i,
                                   barc_j,
                                   barc_v,
                                   &idx);

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 MSK_appendsparsesymmat, we can combine them to set up the objective matrix coefficient \(\barC_j\) using MSK_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.

        r = MSK_putbarcj(task, 0, 1, &idx, &falpha);

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

        r = MSK_putbaraij(task, 0, 0, 1, &idx, &falpha);

Retrieving the solution

After the problem is solved, we read the solution using MSK_getbarxj:

              MSK_getbarxj(task,
                           MSK_SOL_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 8 Source code solving problem (1). Click here to download.
#include <stdio.h>

#include "mosek.h"    /* Include the MOSEK definition file.  */

#define NUMCON    2   /* Number of constraints.              */
#define NUMVAR    3   /* Number of conic quadratic variables */
#define NUMANZ    3   /* Number of non-zeros in A            */
#define NUMBARVAR 1   /* Number of semidefinite variables    */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  MSKint32t    DIMBARVAR[] = {3};         /* Dimension of semidefinite cone */
  MSKint64t    LENBARVAR[] = {3 * (3 + 1) / 2}; /* Number of scalar SD variables  */

  MSKboundkeye bkc[] = { MSK_BK_FX, MSK_BK_FX };
  double       blc[] = { 1.0, 0.5 };
  double       buc[] = { 1.0, 0.5 };

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

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

  MSKint32t    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};
  MSKint32t    conesub[] = {0, 1, 2};

  MSKint32t    i, j;
  MSKint64t    idx;
  double       falpha = 1.0;

  MSKrealt     *xx;
  MSKrealt     *barx;
  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if ( r == MSK_RES_OK )
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, NUMCON, 0, &task);

    if ( r == MSK_RES_OK )
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'NUMCON' empty constraints.
       The constraints will initially have no bounds. */
      if ( r == MSK_RES_OK )
        r = MSK_appendcons(task, NUMCON);

      /* Append 'NUMVAR' variables.
      The variables will initially be fixed at zero (x=0). */
      if ( r == MSK_RES_OK )
        r = MSK_appendvars(task, NUMVAR);

      /* Append 'NUMBARVAR' semidefinite variables. */
      if ( r == MSK_RES_OK ) {
        r = MSK_appendbarvars(task, NUMBARVAR, DIMBARVAR);
      }

      /* Optionally add a constant term to the objective. */
      if ( r == MSK_RES_OK )
        r = MSK_putcfix(task, 0.0);

      /* Set the linear term c_j in the objective.*/
      if ( r == MSK_RES_OK )
        r = MSK_putcj(task, 0, 1.0);

      for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
        r = MSK_putvarbound(task,
                            j,
                            MSK_BK_FR,
                            -MSK_INFINITY,
                            MSK_INFINITY);

      /* Set the linear term barc_j in the objective.*/
      if ( r == MSK_RES_OK )
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   5,
                                   barc_i,
                                   barc_j,
                                   barc_v,
                                   &idx);

      if ( r == MSK_RES_OK )
        r = MSK_putbarcj(task, 0, 1, &idx, &falpha);

      /* Set the bounds on constraints.
        for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            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 (i = 0; i < NUMCON && r == MSK_RES_OK; ++i)
        r = MSK_putarow(task,
                        i,
                        aptre[i] - aptrb[i],
                        asub     + aptrb[i],
                        aval     + aptrb[i]);

      /* Append the conic quadratic cone */
      if ( r == MSK_RES_OK )
        r = MSK_appendcone(task,
                           MSK_CT_QUAD,
                           0.0,
                           3,
                           conesub);

      /* Add the first row of barA */
      if ( r == MSK_RES_OK )
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   3,
                                   bara_i,
                                   bara_j,
                                   bara_v,
                                   &idx);

      if ( r == MSK_RES_OK )
        r = MSK_putbaraij(task, 0, 0, 1, &idx, &falpha);

      /* Add the second row of barA */
      if ( r == MSK_RES_OK )
        r = MSK_appendsparsesymmat(task,
                                   DIMBARVAR[0],
                                   6,
                                   bara_i + 3,
                                   bara_j + 3,
                                   bara_v + 3,
                                   &idx);

      if ( r == MSK_RES_OK )
        r = MSK_putbaraij(task, 1, 0, 1, &idx, &falpha);

      if ( r == MSK_RES_OK )
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);

        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary (task, MSK_STREAM_MSG);

        if ( r == MSK_RES_OK )
        {
          MSKsolstae solsta;

          MSK_getsolsta (task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
            case MSK_SOL_STA_NEAR_OPTIMAL:
              xx   = (MSKrealt *) MSK_calloctask(task, NUMVAR, sizeof(MSKrealt));
              barx = (MSKrealt *) MSK_calloctask(task, LENBARVAR[0], sizeof(MSKrealt));

              MSK_getxx(task,
                        MSK_SOL_ITR,
                        xx);
              MSK_getbarxj(task,
                           MSK_SOL_ITR,    /* Request the interior solution. */
                           0,
                           barx);

              printf("Optimal primal solution\n");
              for (i = 0; i < NUMVAR; ++i)
                printf("x[%d]   : % e\n", i, xx[i]);

              for (i = 0; i < LENBARVAR[0]; ++i)
                printf("barx[%d]: % e\n", i, barx[i]);

              MSK_freetask(task, xx);
              MSK_freetask(task, barx);

              break;

            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
            case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
            case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;

            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined.\n");
              break;

            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc (r,
                         symname,
                         desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return ( r );
} /* main */