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[1] - x[0])
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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)
sage: p.add_constraint(c1*x[0] + c2*x[1] >= matrix.zero(2,2,sparse=True))

sage: # needs cvxopt
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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(1)] - x[Integer(0)])
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> a3 = matrix([[Integer(5), RealNumber('6.')], [RealNumber('6.'), RealNumber('7.')]])
>>> b1 = matrix([[Integer(1), RealNumber('1.')], [RealNumber('1.'), RealNumber('1.')]])
>>> b2 = matrix([[Integer(2), RealNumber('2.')], [RealNumber('2.'), RealNumber('2.')]])
>>> b3 = matrix([[Integer(3), RealNumber('3.')], [RealNumber('3.'), RealNumber('3.')]])
>>> c1 = matrix([[RealNumber('1.0'), Integer(0)],[Integer(0),Integer(0)]],sparse=True)
>>> c2 = matrix([[RealNumber('0.0'), Integer(0)],[Integer(0),Integer(1)]],sparse=True)
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a3)
>>> p.add_constraint(b1*x[Integer(0)] + b2*x[Integer(1)] <= b3)
>>> p.add_constraint(c1*x[Integer(0)] + c2*x[Integer(1)] >= matrix.zero(Integer(2),Integer(2),sparse=True))

>>> # needs cvxopt
>>> p.solver_parameter("show_progress", True)
>>> opt = p.solve()
    pcost       dcost       gap    pres   dres   k/t
0: ...
...
Optimal solution found.
>>> print('Objective Value: {}'.format(N(opt,Integer(3))))
Objective Value: 1.0
>>> [N(x, Integer(3)) for x in sorted(p.get_values(x).values())]
[3.0e-8, 1.0]
>>> 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                                           # needs cvxopt
[ 85555.0 -85555.0]
[-85555.0  85555.0]

>>> from sage.all import *
>>> p.dual_variable(Integer(1))  # rel tol 2e-03                                           # needs cvxopt
[ 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()  # tol 8e-08                           # needs cvxopt
....:        + (p.dual_variable(1)*b3).trace())
0.0

>>> from sage.all import *
>>> opt - ((p.dual_variable(Integer(0))*a3).trace()  # tol 8e-08                           # needs cvxopt
...        + (p.dual_variable(Integer(1))*b3).trace())
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                          # needs cvxopt
0.0

>>> from sage.all import *
>>> (p.slack(Integer(0))*p.dual_variable(Integer(0))).trace()  # tol 2e-07                          # needs cvxopt
0.0

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

sage: # needs sage.graphs
sage: c = graphs.CycleGraph(7)
sage: c2 = c.distance_graph(2).adjacency_matrix()
sage: c3 = c.distance_graph(3).adjacency_matrix()
sage: p.<y> = SemidefiniteProgram()
sage: p.add_constraint((1/7)*matrix.identity(7)>=-y[0]*c2-y[1]*c3)
sage: p.set_objective(y[0]*(c2**2).trace()+y[1]*(c3**2).trace())
sage: x = p.solve(); x + 1                                                          # needs cvxopt
3.31766...

>>> from sage.all import *
>>> # needs sage.graphs
>>> c = graphs.CycleGraph(Integer(7))
>>> c2 = c.distance_graph(Integer(2)).adjacency_matrix()
>>> c3 = c.distance_graph(Integer(3)).adjacency_matrix()
>>> p = SemidefiniteProgram(names=('y',)); (y,) = p._first_ngens(1)
>>> p.add_constraint((Integer(1)/Integer(7))*matrix.identity(Integer(7))>=-y[Integer(0)]*c2-y[Integer(1)]*c3)
>>> p.set_objective(y[Integer(0)]*(c2**Integer(2)).trace()+y[Integer(1)]*(c3**Integer(2)).trace())
>>> x = p.solve(); x + Integer(1)                                                          # needs cvxopt
3.31766...

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

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

>>> from sage.all import *
>>> p.slack(Integer(0)).trace()   # tol 1e-14                                              # needs cvxopt sage.graphs
1.0

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

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

>>> from sage.all import *
>>> # needs cvxopt
>>> p = SemidefiniteProgram(solver='cvxopt')
>>> p.base_ring()
Real Double Field
>>> x = p.new_variable()
>>> RealNumber('0.5') + Integer(3)/Integer(2)*x[Integer(1)]
0.5 + 1.5*x_0

For representing an SDP with exact data, there is another backend:

sage: from sage.numerical.backends.matrix_sdp_backend import MatrixSDPBackend
sage: p = SemidefiniteProgram(solver=MatrixSDPBackend, base_ring=QQ)
sage: p.base_ring()
Rational Field
sage: x = p.new_variable()
sage: 1/2 + 3/2 * x[1]
1/2 + 3/2*x_0

>>> from sage.all import *
>>> from sage.numerical.backends.matrix_sdp_backend import MatrixSDPBackend
>>> p = SemidefiniteProgram(solver=MatrixSDPBackend, base_ring=QQ)
>>> p.base_ring()
Rational Field
>>> x = p.new_variable()
>>> Integer(1)/Integer(2) + Integer(3)/Integer(2) * x[Integer(1)]
1/2 + 3/2*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[1] + 2*b[3]
5 + x_0 + 2*x_1

>>> from sage.all import *
>>> sdp = SemidefiniteProgram(names=('a', 'b',)); (a, b,) = sdp._first_ngens(2)
>>> a
SDPVariable
>>> Integer(5) + a[Integer(1)] + Integer(2)*b[Integer(3)]
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[2]
x_2 - 7*x_3
sage: sdp.show()
Maximization:

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

>>> from sage.all import *
>>> a[Integer(4), 'string', QQ]
x_2
>>> a[Integer(4), 'string', QQ] - Integer(7)*b[Integer(2)]
x_2 - 7*x_3
>>> sdp.show()
Maximization:
<BLANKLINE>
Constraints:
Variables:
  a[1],  b[3],  a[(4, 'string', Rational Field)],  b[2]

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[source]#

Bases: RuntimeError

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'...)

>>> from sage.all import *
>>> from sage.numerical.sdp import SDPSolverException
>>> SDPSolverException("Error")
SDPSolverException('Error'...)
class sage.numerical.sdp.SDPVariable[source]#

Bases: Element

SDPVariable is a variable used by the class SemidefiniteProgram.

Warning

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

items()[source]#

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

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> v = p.new_variable()
>>> p.set_objective(v[Integer(0)] + v[Integer(1)])
>>> sorted(v.items())
[(0, x_0), (1, x_1)]
keys()[source]#

Return the keys already defined in the dictionary.

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> v = p.new_variable()
>>> p.set_objective(v[Integer(0)] + v[Integer(1)])
>>> sorted(v.keys())
[0, 1]
values()[source]#

Return the symbolic variables associated to the current dictionary.

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> v = p.new_variable()
>>> p.set_objective(v[Integer(0)] + v[Integer(1)])
>>> sorted(v.values(), key=str)
[x_0, x_1]
class sage.numerical.sdp.SDPVariableParent[source]#

Bases: Parent

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[source]#

alias of SDPVariable

class sage.numerical.sdp.SemidefiniteProgram[source]#

Bases: SageObject

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:

    • CVXOPT (solver="CVXOPT"). See the CVXOPT website.

    • If solver=None (default), the default solver is used (see default_sdp_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.

See also

EXAMPLES:

Computation of a basic Semidefinite Program:

sage: p = SemidefiniteProgram(maximization=False)
sage: x = p.new_variable()
sage: p.set_objective(x[0] - x[1])
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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)
sage: N(p.solve(), 2)                                                          # needs cvxopt
-3.0

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization=False)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(0)] - x[Integer(1)])
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> a3 = matrix([[Integer(5), RealNumber('6.')], [RealNumber('6.'), RealNumber('7.')]])
>>> b1 = matrix([[Integer(1), RealNumber('1.')], [RealNumber('1.'), RealNumber('1.')]])
>>> b2 = matrix([[Integer(2), RealNumber('2.')], [RealNumber('2.'), RealNumber('2.')]])
>>> b3 = matrix([[Integer(3), RealNumber('3.')], [RealNumber('3.'), RealNumber('3.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a3)
>>> p.add_constraint(b1*x[Integer(0)] + b2*x[Integer(1)] <= b3)
>>> N(p.solve(), Integer(2))                                                          # needs cvxopt
-3.0
add_constraint(linear_function, name=None)[source]#

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[1] + 5*x[2])
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(a1*x[1] + a2*x[2] <= a3)
sage: N(p.solve(), digits=3)                                                # needs cvxopt
16.2

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization=True)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(1)] + Integer(5)*x[Integer(2)])
>>> a1 = matrix([[Integer(1),Integer(2)],[Integer(2),Integer(3)]])
>>> a2 = matrix([[Integer(1),Integer(1)],[Integer(1),Integer(1)]])
>>> a3 = matrix([[Integer(1),-Integer(1)],[-Integer(1),Integer(1)]])
>>> p.add_constraint(a1*x[Integer(1)] + a2*x[Integer(2)] <= a3)
>>> N(p.solve(), digits=Integer(3))                                                # needs cvxopt
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[1] + 5*x[2])
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[1] + a2*x[2])
sage: N(p.solve(), digits=3)                                                # needs cvxopt
16.2

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization=True)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(1)] + Integer(5)*x[Integer(2)])
>>> a1 = matrix([[Integer(1),Integer(2)],[Integer(2),Integer(3)]])
>>> a2 = matrix([[Integer(1),Integer(1)],[Integer(1),Integer(1)]])
>>> a3 = matrix([[Integer(1),-Integer(1)],[-Integer(1),Integer(1)]])
>>> p.add_constraint(a3 >= a1*x[Integer(1)] + a2*x[Integer(2)])
>>> N(p.solve(), digits=Integer(3))                                                # needs cvxopt
16.2
base_ring()[source]#

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

>>> from sage.all import *
>>> p = SemidefiniteProgram(solver='cvxopt')
>>> p.base_ring()
Real Double Field
dual_variable(i, sparse=False)[source]#

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[0] - x[1])
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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)

sage: # needs cvxopt
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

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization=False)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(0)] - x[Integer(1)])
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> a3 = matrix([[Integer(5), RealNumber('6.')], [RealNumber('6.'), RealNumber('7.')]])
>>> b1 = matrix([[Integer(1), RealNumber('1.')], [RealNumber('1.'), RealNumber('1.')]])
>>> b2 = matrix([[Integer(2), RealNumber('2.')], [RealNumber('2.'), RealNumber('2.')]])
>>> b3 = matrix([[Integer(3), RealNumber('3.')], [RealNumber('3.'), RealNumber('3.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a3)
>>> p.add_constraint(b1*x[Integer(0)] + b2*x[Integer(1)] <= b3)

>>> # needs cvxopt
>>> p.solve()                                                              # tol 1e-08
-3.0
>>> x = p.get_values(x).values()
>>> -(a3*p.dual_variable(Integer(0))).trace() - (b3*p.dual_variable(Integer(1))).trace()     # tol 1e-07
-3.0

Dual variable is orthogonal to the slack

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

>>> from sage.all import *
>>> # needs cvxopt
>>> g = sum((p.slack(j)*p.dual_variable(j)).trace() for j in range(Integer(2))); g  # tol 1.2e-08
0.0
gen(i)[source]#

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]

>>> from sage.all import *
>>> sdp = SemidefiniteProgram()
>>> sdp.gen(Integer(0))
x_0
>>> [sdp.gen(i) for i in range(Integer(10))]
[x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9]
get_backend()[source]#

Return the backend instance used.

This might be useful when access 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[0] + a2*x[1] <= a1)
sage: b = p.get_backend()
sage: b.get_matrix()[0][0]
(
    [-1.0 -2.0]
-1, [-2.0 -3.0]
)

>>> from sage.all import *
>>> p = SemidefiniteProgram(solver="cvxopt")
>>> x = p.new_variable()
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a1)
>>> b = p.get_backend()
>>> b.get_matrix()[Integer(0)][Integer(0)]
(
    [-1.0 -2.0]
-1, [-2.0 -3.0]
)
get_values(*lists)[source]#

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[3] - x[5])
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[3] + a2*x[5] <= a3)
sage: p.add_constraint(b1*x[3] + b2*x[5] <= b3)
sage: N(p.solve(), 3)                                                       # needs cvxopt
-3.0

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

To return the optimal value of x[3]:

sage: N(p.get_values(x[3]),3)                                               # needs cvxopt
-1.0

>>> from sage.all import *
>>> N(p.get_values(x[Integer(3)]),Integer(3))                                               # needs cvxopt
-1.0

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

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

>>> from sage.all import *
>>> x_sol = p.get_values(x)                                               # needs cvxopt
>>> sorted(x_sol)                                                         # needs cvxopt
[3, 5]
linear_constraints_parent()[source]#

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> p.linear_constraints_parent()
Linear constraints over Real Double Field
linear_function(x)[source]#

Construct a new linear function.

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> p.linear_function({Integer(0):Integer(1)})
x_0
linear_functions_parent()[source]#

Return the parent for all linear functions.

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> p.linear_functions_parent()
Linear functions over Real Double Field
new_variable(name='')[source]#

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()

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> x = p.new_variable()

It behaves exactly as an usual dictionary would. It can use any key argument you may like, as x[5] 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.add_constraint(a1*x[0] + a1*x[3] <= 0)
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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> x = p.new_variable()
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a1*x[Integer(3)] <= Integer(0))
>>> p.show()
Maximization:
<BLANKLINE>
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()[source]#

Return 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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)
sage: p.add_constraint(b1*x[0] + a2*x[1] <= b3)
sage: p.number_of_constraints()
3

>>> from sage.all import *
>>> p = SemidefiniteProgram(solver = "cvxopt")
>>> x = p.new_variable()
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> a3 = matrix([[Integer(5), RealNumber('6.')], [RealNumber('6.'), RealNumber('7.')]])
>>> b1 = matrix([[Integer(1), RealNumber('1.')], [RealNumber('1.'), RealNumber('1.')]])
>>> b2 = matrix([[Integer(2), RealNumber('2.')], [RealNumber('2.'), RealNumber('2.')]])
>>> b3 = matrix([[Integer(3), RealNumber('3.')], [RealNumber('3.'), RealNumber('3.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a3)
>>> p.add_constraint(b1*x[Integer(0)] + b2*x[Integer(1)] <= b3)
>>> p.add_constraint(b1*x[Integer(0)] + a2*x[Integer(1)] <= b3)
>>> p.number_of_constraints()
3
number_of_variables()[source]#

Return the number of variables used so far.

EXAMPLES:

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

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> a = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> p.add_constraint(a*p[Integer(0)] - a*p[Integer(2)] <=  Integer(2)*a*p[Integer(4)]  )
>>> p.number_of_variables()
3
set_objective(obj)[source]#

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[1] + 5*x[2])
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(a1*x[1] + a2*x[2] <= a3)
sage: N(p.solve(), digits=3)                                                # needs cvxopt
16.2
sage: p.set_objective(None)
sage: _ = p.solve()                                                         # needs cvxopt

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization=True)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(1)] + Integer(5)*x[Integer(2)])
>>> a1 = matrix([[Integer(1),Integer(2)],[Integer(2),Integer(3)]])
>>> a2 = matrix([[Integer(1),Integer(1)],[Integer(1),Integer(1)]])
>>> a3 = matrix([[Integer(1),-Integer(1)],[-Integer(1),Integer(1)]])
>>> p.add_constraint(a1*x[Integer(1)] + a2*x[Integer(2)] <= a3)
>>> N(p.solve(), digits=Integer(3))                                                # needs cvxopt
16.2
>>> p.set_objective(None)
>>> _ = p.solve()                                                         # needs cvxopt
set_problem_name(name)[source]#

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 )

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> p.set_problem_name("Test program")
>>> p
Semidefinite Program "Test program" ( maximization, 0 variables, 0 constraints )
show()[source]#

Display 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[0] - x[1])
sage: p.add_constraint(a1*x[0] + a2*x[1]<= a3)
sage: p.show()
Maximization:
  hihi[0] - hihi[1]
Constraints:
  constraint_0: [1.0 2.0][2.0 3.0]hihi[0] + [2.0 3.0][3.0 4.0]hihi[1] <=  [3.0 4.0][4.0 5.0]
Variables:
   hihi[0],  hihi[1]

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> x = p.new_variable(name="hihi")
>>> a1 = matrix([[Integer(1),Integer(2)],[Integer(2),Integer(3)]])
>>> a2 = matrix([[Integer(2),Integer(3)],[Integer(3),Integer(4)]])
>>> a3 = matrix([[Integer(3),Integer(4)],[Integer(4),Integer(5)]])
>>> p.set_objective(x[Integer(0)] - x[Integer(1)])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)]<= a3)
>>> p.show()
Maximization:
  hihi[0] - hihi[1]
Constraints:
  constraint_0: [1.0 2.0][2.0 3.0]hihi[0] + [2.0 3.0][3.0 4.0]hihi[1] <=  [3.0 4.0][4.0 5.0]
Variables:
   hihi[0],  hihi[1]
slack(i, sparse=False)[source]#

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[0] - x[1])
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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)

sage: # needs cvxopt
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[0]*b1 + x[1]*b2 - b3 + B1       # tol 1e-09
[0.0 0.0]
[0.0 0.0]

>>> from sage.all import *
>>> p = SemidefiniteProgram(maximization = False)
>>> x = p.new_variable()
>>> p.set_objective(x[Integer(0)] - x[Integer(1)])
>>> a1 = matrix([[Integer(1), RealNumber('2.')], [RealNumber('2.'), RealNumber('3.')]])
>>> a2 = matrix([[Integer(3), RealNumber('4.')], [RealNumber('4.'), RealNumber('5.')]])
>>> a3 = matrix([[Integer(5), RealNumber('6.')], [RealNumber('6.'), RealNumber('7.')]])
>>> b1 = matrix([[Integer(1), RealNumber('1.')], [RealNumber('1.'), RealNumber('1.')]])
>>> b2 = matrix([[Integer(2), RealNumber('2.')], [RealNumber('2.'), RealNumber('2.')]])
>>> b3 = matrix([[Integer(3), RealNumber('3.')], [RealNumber('3.'), RealNumber('3.')]])
>>> p.add_constraint(a1*x[Integer(0)] + a2*x[Integer(1)] <= a3)
>>> p.add_constraint(b1*x[Integer(0)] + b2*x[Integer(1)] <= b3)

>>> # needs cvxopt
>>> p.solve()                         # tol 1e-08
-3.0
>>> B1 = p.slack(Integer(1)); B1               # tol 1e-08
[0.0 0.0]
[0.0 0.0]
>>> B1.is_positive_definite()
True
>>> x = sorted(p.get_values(x).values())
>>> x[Integer(0)]*b1 + x[Integer(1)]*b2 - b3 + B1       # tol 1e-09
[0.0 0.0]
[0.0 0.0]
solve(objective_only=False)[source]#

Solve 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)[source]#

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: # needs cvxopt
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[0] - x[1])
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[0] + a2*x[1] <= a3)
sage: p.add_constraint(b1*x[0] + b2*x[1] <= b3)
sage: N(p.solve(), 4)
     pcost       dcost       gap    pres   dres   k/t
 0:  1...
...
Optimal solution found.
-11.

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

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()

>>> from sage.all import *
>>> p = SemidefiniteProgram()
>>> v = p.new_variable()

The following command:

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

>>> from sage.all import *
>>> s = p.sum(v[i] for i in range(Integer(90)))

is much more efficient than:

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

>>> from sage.all import *
>>> s = sum(v[i] for i in range(Integer(90)))