# SemiDefinite Programming¶

A semidefinite program (SDP) is an optimization problem (Wikipedia article Optimization_(mathematics)>) of the following form

$\begin{split}\min \sum_{i,j=1}^n C_{ij}X_{ij} & \qquad \text{(Dual problem)}\\ \text{Subject to:} & \sum_{i,j=1}^n A_{ijk}X_{ij} = b_k, \qquad k=1\dots m\\ &X \succeq 0\end{split}$

where the $$X_{ij}$$, $$1 \leq i,j \leq n$$ are $$n^2$$ variables satisfying the symmetry conditions $$x_{ij} = x_{ji}$$ for all $$i,j$$, the $$C_{ij}=C_{ji}$$, $$A_{ijk}=A_{kji}$$ and $$b_k$$ are real coefficients, and $$X$$ is positive semidefinite, i.e., all the eigenvalues of $$X$$ are nonnegative. The closely related dual problem of this one is the following, where we denote by $$A_k$$ the matrix $$(A_{kij})$$ and by $$C$$ the matrix $$(C_{ij})$$,

$\begin{split}\max \sum_k b_k x_k & \qquad \text{(Primal problem)}\\ \text{Subject to:} & \sum_k x_k A_k \preceq C.\end{split}$

Here $$(x_1,...,x_m)$$ is a vector of scalar variables. A wide variety of problems in optimization can be formulated in one of these two standard forms. Then, solvers are able to calculate an approximation to a solution. Here we refer to the latter problem as primal, and to the former problem as dual. The optimal value of the dual is always at least the optimal value of the primal, and usually (although not always) they are equal.

For instance, suppose you want to maximize $$x_1 - x_0$$ subject to

$\begin{split}\left( \begin{array}{cc} 1 & 2 \\ 2 & 3 \end{array} \right) x_0 + \left( \begin{array}{cc} 3 & 4 \\ 4 & 5 \end{array} \right) x_1 \preceq \left( \begin{array}{cc} 5 & 6 \\ 6 & 7 \end{array} \right),\quad \left( \begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array} \right) x_0 + \left( \begin{array}{cc} 2 & 2 \\ 2 & 2 \end{array} \right) x_1 \preceq \left( \begin{array}{cc} 3 & 3 \\ 3 & 3 \end{array} \right), \quad x_0\geq 0, x_1\geq 0.\end{split}$

An SDP can give you an answer to the problem above. Here is how it’s done:

1. You have to create an instance of SemidefiniteProgram.
2. Create a dictionary $$x$$ of integer variables via new_variable(), for example doing x = p.new_variable() if p is the name of the SDP instance.
3. Add those two matrix inequalities as inequality constraints via add_constraint().
4. Add another matrix inequality to specify nonnegativity of $$x$$.
5. Specify the objective function via set_objective(). In our case it is $$x_1 - x_0$$. If it is a pure constraint satisfaction problem, specify it as None.
6. To check if everything is set up correctly, you can print the problem via show.
7. Solve it and print the solution.

The following example shows all these steps:

sage: p = SemidefiniteProgram()
sage: x = p.new_variable()
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: c1 = matrix([[1.0, 0],[0,0]],sparse=True)
sage: c2 = matrix([[0.0, 0],[0,1]],sparse=True)
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: p.add_constraint(c1*x + c2*x >= matrix.zero(2,2,sparse=True))
sage: p.solver_parameter("show_progress", True)
sage: opt = p.solve()
pcost       dcost       gap    pres   dres   k/t
0: ...
...
Optimal solution found.
sage: print('Objective Value: {}'.format(N(opt,3)))
Objective Value: 1.0
sage: [N(x, 3) for x in sorted(p.get_values(x).values())]
[3.0e-8, 1.0]
sage: p.show()
Maximization:
x_0 - x_1
Constraints:
constraint_0: [3.0 4.0][4.0 5.0]x_0 + [1.0 2.0][2.0 3.0]x_1 <=  [5.0 6.0][6.0 7.0]
constraint_1: [2.0 2.0][2.0 2.0]x_0 + [1.0 1.0][1.0 1.0]x_1 <=  [3.0 3.0][3.0 3.0]
constraint_2: [ 0.0  0.0][ 0.0 -1.0]x_0 + [-1.0  0.0][ 0.0  0.0]x_1 <=  [0 0][0 0]
Variables:
x_0,  x_1


Most solvers, e.g. the default Sage SDP solver CVXOPT, solve simultaneously the pair of primal and dual problems. Thus we can get the optimizer $$X$$ of the dual problem as follows, as diagonal blocks, one per each constraint, via dual_variable(). E.g.:

sage: p.dual_variable(1)  # rel tol 2e-03
[ 85555.0 -85555.0]
[-85555.0  85555.0]


We can see that the optimal value of the dual is equal (up to numerical noise) to $$opt$$.:

sage: opt-((p.dual_variable(0)*a3).trace()+(p.dual_variable(1)*b3).trace())  # tol 8e-08
0.0


Dual variable blocks at optimality are orthogonal to “slack variables”, that is, matrices $$C-\sum_k x_k A_k$$, cf. (Primal problem) above, available via slack(). E.g.:

sage: (p.slack(0)*p.dual_variable(0)).trace()       # tol 2e-07
0.0


More interesting example, the Lovasz theta of the 7-gon:

sage: c=graphs.CycleGraph(7)
sage: p.<y>=SemidefiniteProgram()
sage: p.set_objective(y*(c2**2).trace()+y*(c3**2).trace())
sage: x=p.solve(); x+1
3.31766...


Unlike in the previous example, the slack variable is very far from 0:

sage: p.slack(0).trace()        # tol 1e-14
1.0


The default CVXOPT backend computes with the Real Double Field, for example:

sage: p = SemidefiniteProgram(solver='cvxopt')
sage: p.base_ring()
Real Double Field
sage: x = p.new_variable()
sage: 0.5 + 3/2*x
0.5 + 1.5*x_0


## Linear Variables and Expressions¶

To make your code more readable, you can construct SDPVariable objects that can be arbitrarily named and indexed. Internally, this is then translated back to the $$x_i$$ variables. For example:

sage: sdp.<a,b> = SemidefiniteProgram()
sage: a
SDPVariable
sage: 5 + a + 2*b
5 + x_0 + 2*x_1


Indices can be any object, not necessarily integers. Multi-indices are also allowed:

sage: a[4, 'string', QQ]
x_2
sage: a[4, 'string', QQ] - 7*b
x_2 - 7*x_3
sage: sdp.show()
Maximization:

Constraints:
Variables:
a,  b,  a[(4, 'string', Rational Field)],  b


## Index of functions and methods¶

Below are listed the methods of SemidefiniteProgram. This module also implements the SDPSolverException exception, as well as the SDPVariable class.

 add_constraint() Adds a constraint to the SemidefiniteProgram base_ring() Return the base ring dual_variable() Return optimal dual variable block get_backend() Return the backend instance used get_values() Return values found by the previous call to solve() linear_constraints_parent() Return the parent for all linear constraints linear_function() Construct a new linear function linear_functions_parent() Return the parent for all linear functions new_variable() Return an instance of SDPVariable associated to the SemidefiniteProgram number_of_constraints() Return the number of constraints assigned so far number_of_variables() Return the number of variables used so far set_objective() Set the objective of the SemidefiniteProgram set_problem_name() Set the name of the SemidefiniteProgram slack() Return the slack variable block at the optimum show() Display the SemidefiniteProgram in a human-readable way solve() Solve the SemidefiniteProgram solver_parameter() Return or define a solver parameter sum() Efficiently compute the sum of a sequence of LinearFunction elements

AUTHORS:

• Ingolfur Edvardsson (2014/08): added extension for exact computation
• Dima Pasechnik (2014-) : supervision, minor fixes, duality
exception sage.numerical.sdp.SDPSolverException

Exception raised when the solver fails.

SDPSolverException is the exception raised when the solver fails.

EXAMPLES:

sage: from sage.numerical.sdp import SDPSolverException
sage: SDPSolverException("Error")
SDPSolverException('Error'...)

class sage.numerical.sdp.SDPVariable

SDPVariable is a variable used by the class SemidefiniteProgram.

Warning

You should not instantiate this class directly. Instead, use SemidefiniteProgram.new_variable().

items()

Return the pairs (keys,value) contained in the dictionary.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: v = p.new_variable()
sage: p.set_objective(v + v)
sage: sorted(v.items())
[(0, x_0), (1, x_1)]

keys()

Return the keys already defined in the dictionary.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: v = p.new_variable()
sage: p.set_objective(v + v)
sage: sorted(v.keys())
[0, 1]

values()

Return the symbolic variables associated to the current dictionary.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: v = p.new_variable()
sage: p.set_objective(v + v)
sage: sorted(v.values(), key=str)
[x_0, x_1]

class sage.numerical.sdp.SDPVariableParent

Parent for SDPVariable.

Warning

This class is for internal use. You should not instantiate it yourself. Use SemidefiniteProgram.new_variable() to generate sdp variables.

Element

alias of SDPVariable

class sage.numerical.sdp.SemidefiniteProgram

The SemidefiniteProgram class is the link between Sage, semidefinite programming (SDP) and semidefinite programming solvers.

A Semidefinite Programming (SDP) consists of variables, linear constraints on these variables, and an objective function which is to be maximised or minimised under these constraints.

See the Wikipedia article Semidefinite_programming for further information on semidefinite programming, and the SDP module for its use in Sage.

INPUT:

• solver – selects a solver:
• maximization
• When set to True (default), the SemidefiniteProgram is defined as a maximization.
• When set to False, the SemidefiniteProgram is defined as a minimization.

EXAMPLES:

Computation of a basic Semidefinite Program:

sage: p = SemidefiniteProgram(solver = "cvxopt", maximization=False)
sage: x = p.new_variable()
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: N(p.solve(), 2)
-3.0

add_constraint(linear_function, name=None)

Adds a constraint to the SemidefiniteProgram.

INPUT:

• linear_function – Two different types of arguments are possible:
• A linear function. In this case, arguments min or max have to be specified.
• A linear constraint of the form A <= B, A >= B, A <= B <= C, A >= B >= C or A == B. In this case, arguments min and max will be ignored.
• name – A name for the constraint.

EXAMPLES:

Let’s solve the following semidefinite program:

\begin{split}\begin{aligned} \text{maximize} &\qquad x + 5y \qquad \\ \text{subject to} &\qquad \left( \begin{array}{cc} 1 & 2 \\ 2 & 3 \end{array} \right) x + \left( \begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array} \right) y \preceq \left( \begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right) \end{aligned}\end{split}

This SDP can be solved as follows:

sage: p = SemidefiniteProgram(maximization=True)
sage: x = p.new_variable()
sage: p.set_objective(x + 5*x)
sage: a1 = matrix([[1,2],[2,3]])
sage: a2 = matrix([[1,1],[1,1]])
sage: a3 = matrix([[1,-1],[-1,1]])
sage: N(p.solve(),digits=3)
16.2


One can also define double-bounds or equality using the symbol >= or ==:

sage: p = SemidefiniteProgram(maximization=True)
sage: x = p.new_variable()
sage: p.set_objective(x + 5*x)
sage: a1 = matrix([[1,2],[2,3]])
sage: a2 = matrix([[1,1],[1,1]])
sage: a3 = matrix([[1,-1],[-1,1]])
sage: p.add_constraint(a3 >= a1*x + a2*x)
sage: N(p.solve(),digits=3)
16.2

base_ring()

Return the base ring.

OUTPUT:

A ring. The coefficients that the chosen solver supports.

EXAMPLES:

sage: p = SemidefiniteProgram(solver='cvxopt')
sage: p.base_ring()
Real Double Field

dual_variable(i, sparse=False)

The $$i$$-th dual variable.

Available after self.solve() is called, otherwise the result is undefined.

INPUT:

• index (integer) – the constraint’s id

OUTPUT:

The matrix of the $$i$$-th dual variable.

EXAMPLES:

Dual objective value is the same as the primal one:

sage: p = SemidefiniteProgram(maximization = False)
sage: x = p.new_variable()
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: p.solve()                                                         # tol 1e-08
-3.0
sage: x = p.get_values(x).values()
sage: -(a3*p.dual_variable(0)).trace()-(b3*p.dual_variable(1)).trace()  # tol 1e-07
-3.0


Dual variable is orthogonal to the slack

sage: g = sum((p.slack(j)*p.dual_variable(j)).trace() for j in range(2)); g # tol 1.2e-08
0.0

gen(i)

Return the linear variable $$x_i$$.

EXAMPLES:

sage: sdp = SemidefiniteProgram()
sage: sdp.gen(0)
x_0
sage: [sdp.gen(i) for i in range(10)]
[x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9]

get_backend()

Returns the backend instance used.

This might be useful when acces to additional functions provided by the backend is needed.

EXAMPLES:

This example prints a matrix coefficient:

sage: p = SemidefiniteProgram(solver="cvxopt")
sage: x = p.new_variable()
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: p.add_constraint(a1*x + a2*x <= a1)
sage: b = p.get_backend()
sage: b.get_matrix()
(
[-1.0 -2.0]
-1, [-2.0 -3.0]
)

get_values(*lists)

Return values found by the previous call to solve().

INPUT:

• Any instance of SDPVariable (or one of its elements), or lists of them.

OUTPUT:

• Each instance of SDPVariable is replaced by a dictionary containing the numerical values found for each corresponding variable in the instance.
• Each element of an instance of a SDPVariable is replaced by its corresponding numerical value.

EXAMPLES:

sage: p = SemidefiniteProgram(solver = "cvxopt", maximization = False)
sage: x = p.new_variable()
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: N(p.solve(),3)
-3.0


To return the optimal value of x:

sage: N(p.get_values(x),3)
-1.0


To get a dictionary identical to x containing optimal values for the corresponding variables

sage: x_sol = p.get_values(x)
sage: sorted(x_sol)
[3, 5]


Obviously, it also works with variables of higher dimension:

sage: x_sol = p.get_values(x)

linear_constraints_parent()

Return the parent for all linear constraints.

See linear_functions for more details.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: p.linear_constraints_parent()
Linear constraints over Real Double Field

linear_function(x)

Construct a new linear function.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: p.linear_function({0:1})
x_0

linear_functions_parent()

Return the parent for all linear functions.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: p.linear_functions_parent()
Linear functions over Real Double Field

new_variable(name='')

Returns an instance of SDPVariable associated to the current instance of SemidefiniteProgram.

A new variable x is defined by:

sage: p = SemidefiniteProgram()
sage: x = p.new_variable()


It behaves exactly as an usual dictionary would. It can use any key argument you may like, as x or x["b"], and has methods items() and keys().

INPUT:

• dim – integer. Defines the dimension of the dictionary. If x has dimension $$2$$, its fields will be of the form x[key1][key2]. Deprecated.
• name – string. Associates a name to the variable.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: x = p.new_variable()
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: p.show()
Maximization:

Constraints:
constraint_0: [1.0 2.0][2.0 3.0]x_0 + [1.0 2.0][2.0 3.0]x_1 <=  [0 0][0 0]
Variables:
x_0,  x_1

number_of_constraints()

Returns the number of constraints assigned so far.

EXAMPLES:

sage: p = SemidefiniteProgram(solver = "cvxopt")
sage: x = p.new_variable()
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: p.add_constraint(b1*x + a2*x <= b3)
sage: p.number_of_constraints()
3

number_of_variables()

Returns the number of variables used so far.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: a = matrix([[1, 2.], [2., 3.]])
sage: p.add_constraint(a*p - a*p <=  2*a*p  )
sage: p.number_of_variables()
3

set_objective(obj)

Sets the objective of the SemidefiniteProgram.

INPUT:

• obj – A semidefinite function to be optimized. ( can also be set to None or 0 when just looking for a feasible solution )

EXAMPLES:

Let’s solve the following semidefinite program:

\begin{split}\begin{aligned} \text{maximize} &\qquad x + 5y \qquad \\ \text{subject to} &\qquad \left( \begin{array}{cc} 1 & 2 \\ 2 & 3 \end{array} \right) x + \left( \begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array} \right) y \preceq \left( \begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array} \right) \end{aligned}\end{split}

This SDP can be solved as follows:

sage: p = SemidefiniteProgram(maximization=True)
sage: x = p.new_variable()
sage: p.set_objective(x + 5*x)
sage: a1 = matrix([[1,2],[2,3]])
sage: a2 = matrix([[1,1],[1,1]])
sage: a3 = matrix([[1,-1],[-1,1]])
sage: N(p.solve(),digits=3)
16.2
sage: p.set_objective(None)
sage: _ = p.solve()

set_problem_name(name)

Sets the name of the SemidefiniteProgram.

INPUT:

• name – A string representing the name of the SemidefiniteProgram.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: p.set_problem_name("Test program")
sage: p
Semidefinite Program "Test program" ( maximization, 0 variables, 0 constraints )

show()

Displays the SemidefiniteProgram in a human-readable way.

EXAMPLES:

When constraints and variables have names

sage: p = SemidefiniteProgram()
sage: x = p.new_variable(name="hihi")
sage: a1 = matrix([[1,2],[2,3]])
sage: a2 = matrix([[2,3],[3,4]])
sage: a3 = matrix([[3,4],[4,5]])
sage: p.set_objective(x - x)
sage: p.show()
Maximization:
hihi - hihi
Constraints:
constraint_0: [1.0 2.0][2.0 3.0]hihi + [2.0 3.0][3.0 4.0]hihi <=  [3.0 4.0][4.0 5.0]
Variables:
hihi,  hihi

slack(i, sparse=False)

Slack of the $$i$$-th constraint

Available after self.solve() is called, otherwise the result is undefined

INPUT:

• index (integer) – the constraint’s id.

OUTPUT:

The matrix of the slack of the $$i$$-th constraint

EXAMPLES:

sage: p = SemidefiniteProgram(maximization = False)
sage: x = p.new_variable()
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 5.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 2.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: p.solve()                         # tol 1e-08
-3.0
sage: B1 = p.slack(1); B1               # tol 1e-08
[0.0 0.0]
[0.0 0.0]
sage: B1.is_positive_definite()
True
sage: x = sorted(p.get_values(x).values())
sage: x*b1 + x*b2 - b3 + B1       # tol 1e-09
[0.0 0.0]
[0.0 0.0]

solve(objective_only=False)

Solves the SemidefiniteProgram.

INPUT:

• objective_only – Boolean variable.
• When set to True, only the objective function is returned.
• When set to False (default), the optimal numerical values are stored (takes computational time).

OUTPUT:

The optimal value taken by the objective function.

solver_parameter(name, value=None)

Return or define a solver parameter.

The solver parameters are by essence solver-specific, which means their meaning heavily depends on the solver used.

(If you do not know which solver you are using, then you are using CVXOPT).

INPUT:

• name (string) – the parameter
• value – the parameter’s value if it is to be defined, or None (default) to obtain its current value.

EXAMPLES:

sage: p.<x> = SemidefiniteProgram(solver = "cvxopt", maximization = False)
sage: p.solver_parameter("show_progress", True)
sage: p.solver_parameter("show_progress")
True
sage: p.set_objective(x - x)
sage: a1 = matrix([[1, 2.], [2., 3.]])
sage: a2 = matrix([[3, 4.], [4., 2.]])
sage: a3 = matrix([[5, 6.], [6., 7.]])
sage: b1 = matrix([[1, 1.], [1., 1.]])
sage: b2 = matrix([[2, 2.], [2., 1.]])
sage: b3 = matrix([[3, 3.], [3., 3.]])
sage: p.add_constraint(a1*x + a2*x <= a3)
sage: p.add_constraint(b1*x + b2*x <= b3)
sage: N(p.solve(),4)
pcost       dcost       gap    pres   dres   k/t
0:  1...
...
Optimal solution found.
-11.

sum(L)

Efficiently computes the sum of a sequence of LinearFunction elements.

INPUT:

Note

The use of the regular sum function is not recommended as it is much less efficient than this one.

EXAMPLES:

sage: p = SemidefiniteProgram()
sage: v = p.new_variable()


The following command:

sage: s = p.sum(v[i] for i in range(90))


is much more efficient than:

sage: s = sum(v[i] for i in range(90))