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:
Without loss of generality it is assumed that \(Q^o\) and \(Q^k\) are all symmetric because
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
with a negative semi-definite \(Q^k\) or of the form
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
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:
The matrix formulation of (6.42) has:
with the bounds:
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 MSK_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)):
qsubi[0] = 0; qsubj[0] = 0; qval[0] = 2.0;
qsubi[1] = 1; qsubj[1] = 1; qval[1] = 0.2;
qsubi[2] = 2; qsubj[2] = 0; qval[2] = -1.0;
qsubi[3] = 2; qsubj[3] = 2; qval[3] = 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:
r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
Source code
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */
#define NUMCON 1 /* Number of constraints. */
#define NUMVAR 3 /* Number of variables. */
#define NUMANZ 3 /* Number of non-zeros in A. */
#define NUMQNZ 4 /* Number of non-zeros in Q. */
static void MSKAPI printstr(void *handle,
const char str[])
{
printf("%s", str);
} /* printstr */
int main(int argc, const char *argv[])
{
double c[] = {0.0, -1.0, 0.0};
MSKboundkeye bkc[] = {MSK_BK_LO};
double blc[] = {1.0};
double buc[] = { +MSK_INFINITY };
MSKboundkeye bkx[] = {MSK_BK_LO,
MSK_BK_LO,
MSK_BK_LO
};
double blx[] = {0.0,
0.0,
0.0
};
double bux[] = { +MSK_INFINITY,
+MSK_INFINITY,
+MSK_INFINITY
};
MSKint32t aptrb[] = {0, 1, 2},
aptre[] = {1, 2, 3},
asub[] = {0, 0, 0};
double aval[] = {1.0, 1.0, 1.0};
MSKint32t qsubi[NUMQNZ];
MSKint32t qsubj[NUMQNZ];
double qval[NUMQNZ];
MSKint32t i, j;
double xx[NUMVAR];
MSKenv_t env = NULL;
MSKtask_t task = NULL;
MSKrescodee r;
/* Create the mosek environment. */
r = MSK_makeenv(&env, NULL);
if (r == MSK_RES_OK)
{
/* Create the optimization task. */
r = MSK_maketask(env, NUMCON, NUMVAR, &task);
if (r == MSK_RES_OK)
{
r = 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);
/* Optionally add a constant term to the objective. */
if (r == MSK_RES_OK)
r = MSK_putcfix(task, 0.0);
for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
{
/* Set the linear term c_j in the objective.*/
if (r == MSK_RES_OK)
r = MSK_putcj(task, j, c[j]);
/* Set the bounds on variable j.
blx[j] <= x_j <= bux[j] */
if (r == MSK_RES_OK)
r = MSK_putvarbound(task,
j, /* Index of variable.*/
bkx[j], /* Bound key.*/
blx[j], /* Numerical value of lower bound.*/
bux[j]); /* Numerical value of upper bound.*/
/* Input column j of A */
if (r == MSK_RES_OK)
r = MSK_putacol(task,
j, /* Variable (column) index.*/
aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
asub + aptrb[j], /* Pointer to row indexes of column j.*/
aval + aptrb[j]); /* Pointer to Values of column j.*/
}
/* 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.*/
if (r == MSK_RES_OK)
{
/*
* The lower triangular part of the Q
* matrix in the objective is specified.
*/
qsubi[0] = 0; qsubj[0] = 0; qval[0] = 2.0;
qsubi[1] = 1; qsubj[1] = 1; qval[1] = 0.2;
qsubi[2] = 2; qsubj[2] = 0; qval[2] = -1.0;
qsubi[3] = 2; qsubj[3] = 2; qval[3] = 2.0;
/* Input the Q for the objective. */
r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
}
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;
int j;
MSK_getsolsta(task, MSK_SOL_ITR, &solsta);
switch (solsta)
{
case MSK_SOL_STA_OPTIMAL:
MSK_getxx(task,
MSK_SOL_ITR, /* Request the interior solution. */
xx);
printf("Optimal primal solution\n");
for (j = 0; j < NUMVAR; ++j)
printf("x[%d]: %e\n", j, xx[j]);
break;
case MSK_SOL_STA_DUAL_INFEAS_CER:
case MSK_SOL_STA_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. Termination code: %d.\n", trmcode);
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);
}
}
MSK_deletetask(&task);
}
MSK_deleteenv(&env);
return (r);
} /* 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:
This is equivalent to
where
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 MSK_putqconk.
qsubi[0] = 0; qsubj[0] = 0; qval[0] = -2.0;
qsubi[1] = 1; qsubj[1] = 1; qval[1] = -2.0;
qsubi[2] = 2; qsubj[2] = 2; qval[2] = -0.2;
qsubi[3] = 2; qsubj[3] = 0; qval[3] = 0.2;
/* Put Q^0 in constraint with index 0. */
r = MSK_putqconk(task,
0,
4,
qsubi,
qsubj,
qval);
While MSK_putqconk adds quadratic terms to a specific constraint, it is also possible to input all quadratic terms in one chunk using the MSK_putqcon function.
Source code
here to download.¶#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */
#define NUMCON 1 /* Number of constraints. */
#define NUMVAR 3 /* Number of variables. */
#define NUMANZ 3 /* Number of non-zeros in A. */
#define NUMQNZ 4 /* Number of non-zeros in Q. */
static void MSKAPI printstr(void *handle,
const char str[])
{
printf("%s", str);
} /* printstr */
int main(int argc, const char *argv[])
{
MSKrescodee r;
double c[] = {0.0, -1.0, 0.0};
MSKboundkeye bkc[] = {MSK_BK_LO};
double blc[] = {1.0};
double buc[] = { +MSK_INFINITY};
MSKboundkeye bkx[] = {MSK_BK_LO,
MSK_BK_LO,
MSK_BK_LO
};
double blx[] = {0.0,
0.0,
0.0
};
double bux[] = { +MSK_INFINITY,
+MSK_INFINITY,
+MSK_INFINITY
};
MSKint32t aptrb[] = {0, 1, 2 },
aptre[] = {1, 2, 3},
asub[] = { 0, 0, 0};
double aval[] = { 1.0, 1.0, 1.0};
MSKint32t qsubi[NUMQNZ],
qsubj[NUMQNZ];
double qval[NUMQNZ];
MSKint32t j, i;
double xx[NUMVAR];
MSKenv_t env;
MSKtask_t task;
/* Create the mosek environment. */
r = MSK_makeenv(&env, NULL);
if (r == MSK_RES_OK)
{
/* Create the optimization task. */
r = MSK_maketask(env, NUMCON, NUMVAR, &task);
if (r == MSK_RES_OK)
{
r = 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);
/* Optionally add a constant term to the objective. */
if (r == MSK_RES_OK)
r = MSK_putcfix(task, 0.0);
for (j = 0; j < NUMVAR && r == MSK_RES_OK; ++j)
{
/* Set the linear term c_j in the objective.*/
if (r == MSK_RES_OK)
r = MSK_putcj(task, j, c[j]);
/* Set the bounds on variable j.
blx[j] <= x_j <= bux[j] */
if (r == MSK_RES_OK)
r = MSK_putvarbound(task,
j, /* Index of variable.*/
bkx[j], /* Bound key.*/
blx[j], /* Numerical value of lower bound.*/
bux[j]); /* Numerical value of upper bound.*/
/* Input column j of A */
if (r == MSK_RES_OK)
r = MSK_putacol(task,
j, /* Variable (column) index.*/
aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
asub + aptrb[j], /* Pointer to row indexes of column j.*/
aval + aptrb[j]); /* Pointer to Values of column j.*/
}
/* 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.*/
if (r == MSK_RES_OK)
{
/*
* The lower triangular part of the Q^o
* matrix in the objective is specified.
*/
qsubi[0] = 0; qsubj[0] = 0; qval[0] = 2.0;
qsubi[1] = 1; qsubj[1] = 1; qval[1] = 0.2;
qsubi[2] = 2; qsubj[2] = 0; qval[2] = -1.0;
qsubi[3] = 2; qsubj[3] = 2; qval[3] = 2.0;
/* Input the Q^o for the objective. */
r = MSK_putqobj(task, NUMQNZ, qsubi, qsubj, qval);
}
if (r == MSK_RES_OK)
{
/*
* The lower triangular part of the Q^0
* matrix in the first constraint is specified.
This corresponds to adding the term
- x_1^2 - x_2^2 - 0.1 x_3^2 + 0.2 x_1 x_3
*/
qsubi[0] = 0; qsubj[0] = 0; qval[0] = -2.0;
qsubi[1] = 1; qsubj[1] = 1; qval[1] = -2.0;
qsubi[2] = 2; qsubj[2] = 2; qval[2] = -0.2;
qsubi[3] = 2; qsubj[3] = 0; qval[3] = 0.2;
/* Put Q^0 in constraint with index 0. */
r = MSK_putqconk(task,
0,
4,
qsubi,
qsubj,
qval);
}
if (r == MSK_RES_OK)
r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);
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_LOG);
if (r == MSK_RES_OK)
{
MSKsolstae solsta;
int j;
MSK_getsolsta(task, MSK_SOL_ITR, &solsta);
switch (solsta)
{
case MSK_SOL_STA_OPTIMAL:
MSK_getxx(task,
MSK_SOL_ITR, /* Request the interior solution. */
xx);
printf("Optimal primal solution\n");
for (j = 0; j < NUMVAR; ++j)
printf("x[%d]: %e\n", j, xx[j]);
break;
case MSK_SOL_STA_DUAL_INFEAS_CER:
case MSK_SOL_STA_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. Termination code: %d.\n", trmcode);
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);
}
}
MSK_deletetask(&task);
}
MSK_deleteenv(&env);
return (r);
} /* main */
