9.1 Solving Linear Systems Involving the Basis Matrix

A linear optimization problem always has an optimal solution which is also a basic solution. In an optimal basic solution there are exactly \(m\) basic variables where \(m\) is the number of rows in the constraint matrix \(A\). Define

\[B \in \real^{m\times m}\]

as a matrix consisting of the columns of \(A\) corresponding to the basic variables. The basis matrix \(B\) is always non-singular, i.e.

\[\det (B) \neq 0\]

or, equivalently, \(B^{-1}\) exists. This implies that the linear systems

(1)\[B \bar{x} = w\]

and

(2)\[B^T \bar{x} = w\]

each have a unique solution for all \(w\).

MOSEK provides functions for solving the linear systems (1) and (2) for an arbitrary \(w\).

In the next sections we will show how to use MOSEK to

9.1.1 Basis identification

To use the solutions to (1) and (2) it is important to know how the basis matrix \(B\) is constructed.

Internally MOSEK employs the linear optimization problem

(3)\[\begin{split}\begin{array} {lccccl} \mbox{maximize} & & & c^T x & & \\ \mbox{subject to} & & & Ax - x^c & = & 0, \\ & l^x & \leq & x & \leq & u^x, \\ & l^c & \leq & x^c & \leq & u^c. \end{array}\end{split}\]

where

\[x^c \in \real^m \mbox{ and } x \in \real^n.\]

The basis matrix is constructed of \(m\) columns taken from

\[\left[ \begin{array} {cc} A & -I \end{array} \right].\]

If variable \(x_j\) is a basis variable, then the \(j\)-th column of \(A\), denoted \(a_{:,j}\), will appear in \(B\). Similarly, if \(x_i^c\) is a basis variable, then the \(i\)-th column of \(-I\) will appear in the basis. The ordering of the basis variables and therefore the ordering of the columns of \(B\) is arbitrary. The ordering of the basis variables may be retrieved by calling the function

    MSK_initbasissolve(task,basis);

This function initializes data structures for later use and returns the indexes of the basic variables in the array basis. The interpretation of the basis is as follows. If

\[\mathtt{basis}[i] < \mathtt{numcon},\]

then the \(i\)-th basis variable is \(x_i^c\). Moreover, the \(i\)-th column in \(B\) will be the \(i\)-th column of \(-I\). On the other hand if

\[\mathtt{basis}[i] \geq \mathtt{numcon},\]

then the \(i\)-th basis variable is the variable

\[x_{\mathtt{basis}[i]-\mathtt{numcon}}\]

and the \(i\)-th column of \(B\) is the column

\[A_{:,(\mathtt{basis}[i]-\mathtt{numcon})}.\]

For instance if \(\mathtt{basis}[0] = 4\) and \(\mathtt{numcon} = 5\), then since \(\mathtt{basis}[0]< \mathtt{numcon}\), the first basis variable is \(x_4^c\). Therefore, the first column of \(B\) is the fourth column of \(-I\). Similarly, if \(\mathtt{basis}[1] = 7\), then the second variable in the basis is \(x_{\mathtt{basis}[1]-\mathtt{numcon}} = x_2\). Hence, the second column of \(B\) is identical to \(a_{:,2}\).

An example

Consider the linear optimization problem:

(4)\[\begin{split}\begin{array} {lccl} \mbox{minimize} & x_0 + x_1 & & \\ \mbox{subject to} & x_0 + 2 x_1 & \leq & 2, \\ & x_0 + x_1 & \leq & 6, \\ & x_0, x_1 \geq 0. & & \end{array}\end{split}\]

Suppose a call to MSK_initbasissolve returns an array basis so that

basis[0] = 1,
basis[1] = 2.

Then the basis variables are \(x_1^c\) and \(x_0\) and the corresponding basis matrix \(B\) is

\[\begin{split}\left[ \begin{array} {cc} 0 & 1 \\ -1 & 1 \end{array} \right].\end{split}\]

Please note the ordering of the columns in \(B\) .

Listing 9.1 A program showing how to identify the basis.
#include "mosek.h"

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

int main(int argc, const char *argv[])
{
  MSKenv_t     env;
  MSKtask_t    task;
  MSKint32t    numcon = 2, numvar = 2;
  double       c[]    = {1.0, 1.0};
  MSKint32t    ptrb[] = {0, 2},
               ptre[] = {2, 3};
  MSKint32t    asub[] = {0, 1, 0, 1};
  double       aval[] = {1.0, 1.0, 2.0, 1.0};
  MSKboundkeye bkc[]  = { MSK_BK_UP, MSK_BK_UP };
  double       blc[]  = { -MSK_INFINITY, -MSK_INFINITY };
  double       buc[]  = {2.0, 6.0};

  MSKboundkeye bkx[]  = { MSK_BK_LO, MSK_BK_LO };
  double       blx[]  = {0.0, 0.0};
  double       bux[]  = { +MSK_INFINITY, +MSK_INFINITY};
  MSKrescodee  r      = MSK_RES_OK;
  MSKint32t    i, nz;
  double       w[]    = {2.0, 6.0};
  MSKint32t    sub[]  = {0, 1};
  MSKint32t    *basis;

  if (r == MSK_RES_OK)
    r = MSK_makeenv(&env, NULL);

  if ( r == MSK_RES_OK )
    r = MSK_makeemptytask(env, &task);

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

  if ( r == MSK_RES_OK)
    r = MSK_inputdata(task, numcon, numvar, numcon, numvar,
                      c, 0.0,
                      ptrb, ptre, asub, aval, bkc, blc, buc, bkx, blx, bux);

  if (r == MSK_RES_OK)
    r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

  if (r == MSK_RES_OK)
    r = MSK_optimize(task);

  if (r == MSK_RES_OK)
    basis = MSK_calloctask(task, numcon, sizeof(MSKint32t));

  if (r == MSK_RES_OK)
    r = MSK_initbasissolve(task, basis);

  /* List basis variables corresponding to columns of B */
  for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
  {
    printf("basis[%d] = %d\n", i, basis[i]);
    if (basis[sub[i]] < numcon)
      printf ("Basis variable no %d is xc%d.\n", i, basis[i]);
    else
      printf ("Basis variable no %d is x%d.\n", i, basis[i] - numcon);
  }

  nz = 2;
  /* solve Bx = w */
  /* sub contains index of non-zeros in w.
     On return w contains the solution x and sub
     the index of the non-zeros in x.
   */
  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 0, &nz, sub, w);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to Bx = w:\n\n");

    /* Print solution and b. */

    for (i = 0; i < nz; ++i)
    {
      if (basis[sub[i]] < numcon)
        printf ("xc%d = %e\n", basis[sub[i]] , w[sub[i]] );
      else
        printf ("x%d = %e\n", basis[sub[i]] - numcon , w[sub[i]] );
    }
  }

  /* Solve B^T y = w */
  nz     = 1;     /* Only one element in sub is nonzero. */
  sub[0] = 1;     /* Only w[1] is nonzero. */
  w[0]   = 0.0;
  w[1]   = 1.0;

  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 1, &nz, sub, w);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to B^T y = w:\n\n");
    /* Print solution and y. */
    for (i = 0; i < nz; ++i)
      printf ("y%d = %e\n", sub[i] , w[sub[i]] );
  }

  return ( r );
}/* main */

In the example above the linear system is solved using the optimal basis for (4) and the original right-hand side of the problem. Thus the solution to the linear system is the optimal solution to the problem. When running the example program the following output is produced.

basis[0] = 1
Basis variable no 0 is xc1.
basis[1] = 2
Basis variable no 1 is x0.

Solution to Bx = b:

x0 = 2.000000e+00
xc1 = -4.000000e+00

Solution to B^Tx = c:

x1 = -1.000000e+00
x0 = 1.000000e+00

Please note that the ordering of the basis variables is

\[\begin{split}\left[ \begin{array} {c} x^c_1 \\ x_0 \end{array} \right]\end{split}\]

and thus the basis is given by:

\[\begin{split}B = \left[ \begin{array}{cc} 0 & 1 \\ -1 & 1 \end{array} \right]\end{split}\]

It can be verified that

\[\begin{split}\left[ \begin{array}{c} x^c_1 \\ x_0 \end{array} \right] = \left[ \begin{array} {c} -4 \\ 2 \end{array} \right]\end{split}\]

is a solution to

\[\begin{split}\left[ \begin{array}{cc} 0 & 1 \\ -1 & 1 \end{array} \right] \left[ \begin{array} {c} x^c_1 \\ x_0 \end{array} \right] = \left[ \begin{array} {c} 2 \\ 6 \end{array} \right].\end{split}\]

9.1.2 Solving arbitrary linear systems

MOSEK can be used to solve an arbitrary (rectangular) linear system

\[Ax = b\]

using the MSK_solvewithbasis function without optimizing the problem as in the previous example. This is done by setting up an \(A\) matrix in the task, setting all variables to basic and calling the MSK_solvewithbasis function with the \(b\) vector as input. The solution is returned by the function.

An example

Below we demonstrate how to solve the linear system

(5)\[\begin{split}\left[ \begin{array} {cc} 0 & 1 \\ -1 & 1 \end{array} \right] \left[ \begin{array} {c} x_0 \\ x_1 \end{array} \right] = \left[ \begin{array}{c} b_1 \\ b_2 \end{array} \right]\end{split}\]

with two inputs \(b=(1,-2)\) and \(b=(7,0)\) .

#include "mosek.h"

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


MSKrescodee put_a(MSKtask_t task,
                  double *aval,
                  MSKidxt *asub,
                  MSKidxt *ptrb,
                  MSKidxt *ptre,
                  int numvar,
                  MSKidxt *basis
                 )

{
  MSKrescodee r = MSK_RES_OK;
  int i;
  MSKstakeye *skx = NULL , *skc = NULL;

  skx = (MSKstakeye *) calloc(numvar, sizeof(MSKstakeye));
  if (skx == NULL && numvar)
    r = MSK_RES_ERR_SPACE;

  skc = (MSKstakeye *) calloc(numvar, sizeof(MSKstakeye));
  if (skc == NULL && numvar)
    r = MSK_RES_ERR_SPACE;

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
  {
    skx[i] = MSK_SK_BAS;
    skc[i] = MSK_SK_FIX;
  }

  /* Create a coefficient matrix and right hand
     side with the data from the linear system */
  if (r == MSK_RES_OK)
    r = MSK_appendvars(task, numvar);

  if (r == MSK_RES_OK)
    r = MSK_appendcons(task, numvar);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putacol(task, i, ptre[i] - ptrb[i], asub + ptrb[i], aval + ptrb[i]);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putconbound(task, i, MSK_BK_FX, 0, 0);

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

  /* Allocate space for the solution and set status to unknown */

  if (r == MSK_RES_OK)
    r = MSK_deletesolution(task, MSK_SOL_BAS);

  /* Define a basic solution by specifying
     status keys for variables & constraints. */
  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putsolutioni (
          task,
          MSK_ACC_VAR,
          i,
          MSK_SOL_BAS,
          skx[i],
          0.0,
          0.0,
          0.0,
          0.0);

  for (i = 0; i < numvar && r == MSK_RES_OK; ++i)
    r = MSK_putsolutioni (
          task,
          MSK_ACC_CON,
          i,
          MSK_SOL_BAS,
          skc[i],
          0.0,
          0.0,
          0.0,
          0.0);

  if (r == MSK_RES_OK)
    r = MSK_initbasissolve(task, basis);

  free (skx);
  free (skc);

  return ( r );

}

#define NUMCON 2
#define NUMVAR 2

int main(int argc, const char *argv[])
{
  MSKenv_t  env;
  MSKtask_t task;
  MSKrescodee r = MSK_RES_OK;
  MSKintt   numvar = NUMCON;
  MSKintt   numcon = NUMVAR;   /* we must have numvar == numcon */
  int       i, nz;
  double    aval[] = { -1.0, 1.0, 1.0};
  MSKidxt   asub[] = {1, 0, 1};
  MSKidxt   ptrb[] = {0, 1};
  MSKidxt   ptre[] = {1, 3};

  MSKidxt   bsub[NUMCON];
  double    b[NUMCON];

  MSKidxt   *basis = NULL;

  if (r == MSK_RES_OK)
    r = MSK_makeenv(&env, NULL);

  if ( r == MSK_RES_OK )
    r = MSK_makeemptytask(env, &task);

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

  basis = (MSKidxt *) calloc(numcon, sizeof(MSKidxt));
  if ( basis == NULL && numvar)
    r = MSK_RES_ERR_SPACE;

  /* Put A matrix and factor A.
     Call this function only once for a given task. */
  if (r == MSK_RES_OK)
    r = put_a( task,
               aval,
               asub,
               ptrb,
               ptre,
               numvar,
               basis
             );

  /* now solve rhs */
  b[0] = 1;
  b[1] = -2;
  bsub[0] = 0;
  bsub[1] = 1;
  nz = 2;

  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 0, &nz, bsub, b);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to Bx = b:\n\n");
    /* Print solution and show correspondents
       to original variables in the problem */
    for (i = 0; i < nz; ++i)
    {
      if (basis[bsub[i]] < numcon)
        printf("This should never happen\n");
      else
        printf ("x%d = %e\n", basis[bsub[i]] - numcon , b[bsub[i]] );
    }
  }

  b[0] = 7;
  bsub[0] = 0;
  nz = 1;

  if (r == MSK_RES_OK)
    r = MSK_solvewithbasis(task, 0, &nz, bsub, b);

  if (r == MSK_RES_OK)
  {
    printf("\nSolution to Bx = b:\n\n");
    /* Print solution and show correspondents
       to original variables in the problem */
    for (i = 0; i < nz; ++i)
    {
      if (basis[bsub[i]] < numcon)
        printf("This should never happen\n");
      else
        printf ("x%d = %e\n", basis[bsub[i]] - numcon , b[bsub[i]] );
    }
  }

  free (basis);
  return r;
}

The most important step in the above example is the definition of the basic solution, where we define the status key for each variable. The actual values of the variables are not important and can be selected arbitrarily, so we set them to zero. All variables corresponding to columns in the linear system we want to solve are set to basic and the slack variables for the constraints, which are all non-basic, are set to their bound.

The program produces the output:

Solution to Bx = b:

x1 = 1
x0 = 3

Solution to Bx = b:

x1 = 7
x0 = 7