6.7 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
where \(\Symm^r\) is the set of \(r \times r\) real-valued symmetric matrices.
MOSEK can solve semidefinite optimization problems stated in the primal form,
where the problem has \(p\) symmetric positive semidefinite variables \(\barX_j\in \PSD^{r_j}\) of dimension \(r_j\). The symmetric coefficient matrices \(\barC_j\in \Symm^{r_j}\) and \(\barA_{i,j}\in \Symm^{r_j}\) are used to specify PSD terms in the linear objective and the linear constraints, respectively. The symmetric coefficient matrices \(\barF_{i,j}\in \Symm^{r_j}\) are used to specify PSD terms in the affine conic constraints. Note that \(q\) ((6.22)) is the total dimension of all the cones, i.e. \(q=\text{dim}(\K_1 \times \ldots \times \K_k)\), given there are \(k\) ACCs. We use standard notation for the matrix inner product, i.e., for \(A,B\in \real^{m\times n}\) we have
In addition to the primal form presented above, semidefinite problems can be expressed in their dual form. Constraints in this form are usually called linear matrix inequalities (LMIs). LMIs can be easily specified in MOSEK using the vectorized positive semidefinite cone which is defined as:
Vectorized semidefinite domain:
\[\PSD^{d,\mathrm{vec}} = \left\{(x_1,\ldots,x_{d(d+1)/2})\in \real^n~:~ \mathrm{sMat}(x)\in\PSD^d\right\},\]where \(n=d(d+1)/2\) and,
\[\begin{split}\mathrm{sMat}(x) = \left[\begin{array}{cccc}x_1 & x_2/\sqrt{2} & \cdots & x_{d}/\sqrt{2} \\ x_2/\sqrt{2} & x_{d+1} & \cdots & x_{2d-1}/\sqrt{2} \\ \cdots & \cdots & \cdots & \cdots \\ x_{d}/\sqrt{2} & x_{2d-1}/\sqrt{2} & \cdots & x_{d(d+1)/2}\end{array}\right],\end{split}\]or equivalently
\[\PSD^{d,\mathrm{vec}} = \left\{\mathrm{sVec}(X)~:~X\in\PSD^d\right\},\]where
\[\mathrm{sVec}(X) = (X_{11},\sqrt{2}X_{21},\ldots,\sqrt{2}X_{d1},X_{22},\sqrt{2}X_{32},\ldots,X_{dd}).\]
In other words, the domain consists of vectorizations of the lower-triangular part of a positive semidefinite matrix, with the non-diagonal elements additionally rescaled. LMIs can be expressed by restricting appropriate affine expressions to this cone type.
For other types of cones supported by MOSEK, see Sec. 15.9 (Supported domains) and the other tutorials in this chapter. Different cone types can appear together in one optimization problem.
We demonstrate the setup of semidefinite variables and their coefficient matrices in the following examples:
Sec. 6.7.1 (Example SDO1): A problem with one semidefinite variable and linear and conic constraints.
Sec. 6.7.2 (Example SDO2): A problem with two semidefinite variables with a linear constraint and bound.
Sec. 6.7.3 (Example SDO_LMI: Linear matrix inequalities and the vectorized semidefinite domain): A problem with linear matrix inequalities and the vectorized semidefinite domain.
6.7.1 Example SDO1¶
We consider the simple optimization problem with semidefinite and conic quadratic constraints:
The problem description contains a 3-dimensional symmetric semidefinite variable which can be written explicitly as:
and an affine conic constraint (ACC) \((x_0, x_1, x_2) \in \Q^3\). The objective is to minimize
subject to the two linear constraints
Setting up the linear and conic part
The linear and conic parts (constraints, variables, objective, ACC) are set up using the methods described in the relevant tutorials; Sec. 6.1 (Linear Optimization), Sec. 6.2 (From Linear to Conic 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.append_barvars.
task.append_barvars(&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.append_sparse_sym_mat.
let c_symmat_idx = task.append_sparse_sym_mat(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.append_sparse_sym_mat, we can combine them to set up the objective matrix coefficient \(\barC_j\) using Task.put_barc_j, 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.put_barc_j(0, &[c_symmat_idx], &[falpha])?;
Similarly, a constraint matrix coefficient \(\barA_{ij}\) is set up by the function Task.put_bara_ij.
task.put_bara_ij(0, 0, &[a_symmat_idx1][..], &[falpha][..])?;
Retrieving the solution
After the problem is solved, we read the solution using Task.get_barx_j:
let mut barx = vec![0.0,0.0,0.0,0.0,0.0,0.0];
task.get_barx_j(Soltype::ITR, /* Request the interior solution. */
0,
& mut 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
extern crate mosek;
use mosek::{Task,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
const NUMCON : usize = 2; /* Number of constraints. */
const NUMVAR : usize = 3; /* Number of conic quadratic variables */
fn main() -> Result<(),String>
{
let dimbarvar = vec![3]; /* Dimension of semidefinite cone */
let bkc = &[ mosek::Boundkey::FX, mosek::Boundkey::FX ];
let blc = &[ 1.0, 0.5 ];
let buc = &[ 1.0, 0.5 ];
let barc_i = &[0, 1, 1, 2, 2];
let barc_j = &[0, 0, 1, 1, 2];
let barc_v = &[2.0, 1.0, 2.0, 1.0, 2.0];
let aptrb = &[0, 1];
let aptre = &[1, 3];
let asub = &[0, 1, 2]; /* column subscripts of A */
let aval = &[1.0, 1.0, 1.0];
let bara_i = &[0, 1, 2, 0, 1, 2, 1, 2, 2];
let bara_j = &[0, 1, 2, 0, 0, 0, 1, 1, 2];
let bara_v = &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let falpha = 1.0;
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'NUMCON' empty constraints.
* The constraints will initially have no bounds. */
task.append_cons(NUMCON as i32)?;
/* Append 'NUMVAR' variables.
* The variables will initially be fixed at zero (x=0). */
task.append_vars(NUMVAR as i32)?;
/* Append 'NUMBARVAR' semidefinite variables. */
task.append_barvars(&dimbarvar[..])?;
/* Optionally add a constant term to the objective. */
task.put_cfix(0.0)?;
/* Set the linear term c_j in the objective.*/
task.put_c_j(0,1.0)?;
for j in 0..NUMVAR {
task.put_var_bound(j as i32,
mosek::Boundkey::FR,
-INF,
INF)?;
}
/* Set the linear term barc_j in the objective.*/
let c_symmat_idx = task.append_sparse_sym_mat(dimbarvar[0],
barc_i,
barc_j,
barc_v)?;
task.put_barc_j(0, &[c_symmat_idx], &[falpha])?;
for i in 0..NUMCON
{
/* Input A row by row */
task.put_a_row(i as i32,
& asub[aptrb[i]..aptre[i]],
& aval[aptrb[i]..aptre[i]])?;
/* Set the bounds on constraints.
* for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
task.put_con_bound(i as i32, /* Index of constraint.*/
bkc[i], /* Bound key.*/
blc[i], /* Numerical value of lower bound.*/
buc[i])?; /* Numerical value of upper bound.*/
}
{
/* Append the conic quadratic cone */
let afei = task.get_num_afe()?;
task.append_afes(3)?;
task.put_afe_f_entry_list(&[0,1,2],
&[0,1,2],
&[1.0,1.0,1.0])?;
let dom = task.append_quadratic_cone_domain(3)?;
task.append_acc_seq(dom,afei,&[0.0,0.0,0.0])?;
}
/* Add the first row of barA */
let a_symmat_idx1 =
task.append_sparse_sym_mat(dimbarvar[0],
& bara_i[..3],
& bara_j[..3],
& bara_v[..3])?;
task.put_bara_ij(0, 0, &[a_symmat_idx1][..], &[falpha][..])?;
/* Add the second row of barA */
let a_symmat_idx2 =
task.append_sparse_sym_mat(dimbarvar[0],
& bara_i[3..9],
& bara_j[3..9],
& bara_v[3..9])?;
task.put_bara_ij(1, 0, &[a_symmat_idx2][..], &[falpha][..])?;
let _trmcode = task.optimize()?;
task.write_data("sdo1.ptf")?;
/* Print a summary containing information
* about the solution for debugging purposes*/
task.solution_summary (Streamtype::MSG)?;
let solsta = task.get_sol_sta(Soltype::ITR)?;
match solsta
{
Solsta::OPTIMAL =>
{
let mut xx = vec![0.0,0.0,0.0];
task.get_xx(Soltype::ITR, /* Request the basic solution. */
& mut xx[..])?;
let mut barx = vec![0.0,0.0,0.0,0.0,0.0,0.0];
task.get_barx_j(Soltype::ITR, /* Request the interior solution. */
0,
& mut barx[..])?;
println!("Optimal primal solution");
for j in 0..NUMVAR as usize
{
println!("x[{}]: {}",j,xx[j]);
}
let n = dimbarvar[0] as usize;
for j in 0..n
{
for i in j..n
{
println!("barx[{},{}]: {}",i,j,barx[j*n+i-j*(j+1)/2]);
}
}
}
Solsta::DUAL_INFEAS_CER |
Solsta::PRIM_INFEAS_CER =>
{
println!("Primal or dual infeasibility certificate found.");
}
Solsta::UNKNOWN =>
{
/* If the solutions status is unknown, print the termination code
* indicating why the optimizer terminated prematurely. */
println!("The solution status is unknown.");
println!("The optimizer terminitated with code: {}",solsta);
}
_ =>
{
println!("Other solution status.");
}
}
Ok(())
} /* main */
6.7.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:
In our example \(\dim(\barX_1)=3\), \(\dim(\barX_2)=4\), \(b=23\), \(k=-3\) and
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.7.1 (Example SDO1).
Other than in Sec. 6.7.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.put_bara_block_triplet. 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.put_barc_block_triplet but only include \((j,k,l)\) and \(v\).
For explanations of other data structures used in the example see Sec. 6.7.1 (Example SDO1).
The code representing the above problem is shown below.
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
/* Input data */
let numcon : i32 = 2; /* Number of constraints. */
let dimbarvar : &[i32] = &[3, 4]; /* Dimension of semidefinite variables */
/* Objective coefficients concatenated */
let Cj : &[i32] = &[ 0, 0, 1, 1, 1, 1 ]; /* Which symmetric variable (j) */
let Ck : &[i32] = &[ 0, 2, 0, 1, 1, 2 ]; /* Which entry (k,l)->v */
let Cl : &[i32] = &[ 0, 2, 0, 0, 1, 2 ];
let Cv : &[f64] = &[ 1.0, 6.0, 1.0, -3.0, 2.0, 1.0 ];
/* Equality constraints coefficients concatenated */
let Ai : &[i32] = &[ 0, 0, 0, 0, 0, 0 ]; /* Which constraint (i = 0) */
let Aj : &[i32] = &[ 0, 0, 0, 1, 1, 1 ]; /* Which symmetric variable (j) */
let Ak : &[i32] = &[ 0, 2, 2, 1, 1, 3 ]; /* Which entry (k,l)->v */
let Al : &[i32] = &[ 0, 0, 2, 0, 1, 3 ];
let Av : &[f64] = &[ 1.0, 1.0, 2.0, 1.0, -1.0, -3.0 ];
/* The second constraint - one-term inequality */
let A2i : &[i32] = &[ 1 ]; /* Which constraint (i = 1) */
let A2j : &[i32] = &[ 1 ]; /* Which symmetric variable (j = 1) */
let A2k : &[i32] = &[ 1 ]; /* Which entry A(1,0) = A(0,1) = 0.5 */
let A2l : &[i32] = &[ 0 ];
let A2v : &[f64] = &[ 0.5 ];
let bkc = &[ mosek::Boundkey::FX,
mosek::Boundkey::UP ];
let blc = &[ 23.0, 0.0 ];
let buc = &[ 23.0, -3.0 ];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append numcon empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append numbarvar semidefinite variables. */
task.append_barvars(dimbarvar)?;
/* Set objective (6 nonzeros).*/
task.put_barc_block_triplet(Cj, Ck, Cl, Cv)?;
/* Set the equality constraint (6 nonzeros).*/
task.put_bara_block_triplet(Ai, Aj, Ak, Al, Av)?;
/* Set the inequality constraint (1 nonzero).*/
task.put_bara_block_triplet(A2i, A2j, A2k, A2l, A2v)?;
/* Set constraint bounds */
task.put_con_bound_slice(0, 2, bkc, blc, buc)?;
/* Run optimizer */
task.optimize()?;
task.solution_summary(Streamtype::MSG)?;
//mosek.solsta[] solsta = new mosek.solsta[1];
let solsta = task.get_sol_sta (Soltype::ITR)?;
match solsta {
Solsta::OPTIMAL => {
/* Retrieve the soution for all symmetric variables */
println!("Solution (lower triangular part vectorized):");
for (i,dimbarvari) in dimbarvar.iter().enumerate() {
//let dim = dimbarvar[i] * (dimbarvar[i] + 1) / 2;
let dim = dimbarvari * (dimbarvari+1)/2;
//double[] barx = new double[dim];
let mut barx : Vec<f64> = vec![0.0; dim as usize];
task.get_barx_j(Soltype::ITR, i as i32, barx.as_mut_slice())?;
println!("X{}: {:?}",i+1,barx);
// for (int j = 0; j < dim; ++j)
// System.out.print(barx[j] + " ");
// System.out.println();
}
},
Solsta::DUAL_INFEAS_CER|Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility certificate found."),
Solsta::UNKNOWN =>
println!("The status of the solution could not be determined."),
_ => println!("Other solution status.")
}
Ok(())
}
6.7.3 Example SDO_LMI: Linear matrix inequalities and the vectorized semidefinite domain¶
The standard form of a semidefinite problem is usually either based on semidefinite variables (primal form) or on linear matrix inequalities (dual form). However, MOSEK allows mixing of these two forms, as shown in (6.25)
The first affine expression is restricted to a linear domain and could also be modelled as a linear constraint (instead of an ACC). The lower triangular part of the linear matrix inequality (second constraint) can be vectorized and restricted to the Domaintype::SVEC_PSD_CONE. This allows us to express the constraints in (6.25) as the affine conic constraints shown in (6.26).
Vectorization of the LMI is performed as explained in Sec. 15.9 (Supported domains).
Setting up the linear part
The linear parts (objective, constraints, variables) and the semidefinite terms in the linear expressions are defined exactly as shown in the previous examples.
Setting up the affine conic constraints with semidefinite terms
To define the affine conic constraints, we first set up the affine expressions. The \(F\) matrix and the \(g\) vector are defined as usual. Additionally, we specify the coefficients for the semidefinite variables. The semidefinite coefficients shown in (6.26) are setup using the function Task.put_afe_barf_block_triplet.
task.put_afe_barf_block_triplet(barf_i, barf_j, barf_k, barf_l, barf_v)?;
These affine expressions are then included in their corresponding domains to construct the affine conic constraints. Lastly, the ACCs are appended to the task.
/* Append R+ domain and the corresponding ACC */
{
let dom = task.append_rplus_domain(1)?;
task.append_acc(dom, &[0], &[0.0])?;
}
/* Append SVEC_PSD domain and the corresponding ACC */
{
let dom = task.append_svec_psd_cone_domain(3)?;
task.append_acc(dom, &[1,2,3], &[0.0,0.0,0.0])?;
}
Source code
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numafe : i64 = 4; /* Number of affine expressions. */
let numvar : i32 = 2; /* Number of scalar variables */
let dimbarvar = &[2]; /* Dimension of semidefinite cone */
let lenbarvar = &[2 * (2 + 1) / 2]; /* Number of scalar SD variables */
let barc_j = &[0, 0];
let barc_k = &[0, 1];
let barc_l = &[0, 1];
let barc_v = &[1.0, 1.0];
let afeidx = &[0, 0, 1, 2, 2, 3];
let varidx = &[0, 1, 1, 0, 1, 0];
let f_val = &[-1.0, -1.0, 3.0, 2.0f64.sqrt(), 2.0f64.sqrt(), 3.0];
let g = &[0.0, -1.0, 0.0, -1.0];
let barf_i = &[0, 0];
let barf_j = &[0, 0];
let barf_k = &[0, 1];
let barf_l = &[0, 0];
let barf_v = &[0.0, 1.0];
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'NUMAFE' empty affine expressions. */
task.append_afes(numafe)?;
/* Append 'NUMVAR' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
/* Append 'NUMBARVAR' semidefinite variables. */
task.append_barvars(dimbarvar)?;
task.put_obj_sense(Objsense::MINIMIZE)?;
/* Optionally add a constant term to the objective. */
task.put_cfix(1.0)?;
/* Set the linear term c_j in the objective.*/
task.put_c_j(0, 1.0)?;
task.put_c_j(1, 1.0)?;
task.put_var_bound_slice_const(0,numvar, Boundkey::FR, -INF,INF)?;
/* Set the linear term barc_j in the objective.*/
task.put_barc_block_triplet(barc_j, barc_k, barc_l, barc_v)?;
/* Set up the affine conic constraints */
/* Construct the affine expressions */
/* F matrix */
task.put_afe_f_entry_list(afeidx, varidx, f_val)?;
/* g vector */
task.put_afe_g_slice(0, 4, g)?;
/* barF block triplets */
task.put_afe_barf_block_triplet(barf_i, barf_j, barf_k, barf_l, barf_v)?;
/* Append R+ domain and the corresponding ACC */
{
let dom = task.append_rplus_domain(1)?;
task.append_acc(dom, &[0], &[0.0])?;
}
/* Append SVEC_PSD domain and the corresponding ACC */
{
let dom = task.append_svec_psd_cone_domain(3)?;
task.append_acc(dom, &[1,2,3], &[0.0,0.0,0.0])?;
}
/* Run optimizer */
let _ = task.optimize()?;
/* Print a summary containing information
about the solution for debugging purposes*/
task.solution_summary (mosek::Streamtype::MSG)?;
let solsta = task.get_sol_sta(mosek::Soltype::ITR)?;
match solsta {
Solsta::OPTIMAL => {
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR,xx.as_mut_slice())?;
let mut barx = vec![0.0; lenbarvar[0]];
task.get_barx_j(Soltype::ITR, 0, barx.as_mut_slice())?; /* Request the interior solution. */
println!("Optimal primal solution");
println!(" x = {:?}",xx);
println!(" barx = {:?}",barx);
},
Solsta::DUAL_INFEAS_CER|Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility certificate found."),
Solsta::UNKNOWN =>
println!("The status of the solution could not be determined."),
_ =>
println!("Other solution status.")
}
Ok(())
}
