Coverage for local/lib/python2.7/site-packages/sage/numerical/sdp.pyx : 80%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
r""" SemiDefinite Programming
A semidefinite program (:wikipedia:`SDP <Semidefinite_programming>`) is an optimization problem (:wikipedia:`Optimization_(mathematics)>`) of the following form
.. MATH::
\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
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})`,
.. MATH::
\max \sum_k b_k x_k & \qquad \text{(Primal problem)}\\ \text{Subject to:} & \sum_k x_k A_k \preceq C.
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
.. MATH::
\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.
An SDP can give you an answer to the problem above. Here is how it's done:
#. You have to create an instance of :class:`SemidefiniteProgram`. #. Create a dictionary `x` of integer variables via :meth:`~SemidefiniteProgram.new_variable`, for example doing ``x = p.new_variable()`` if ``p`` is the name of the SDP instance. #. Add those two matrix inequalities as inequality constraints via :meth:`~SemidefiniteProgram.add_constraint`. #. Add another matrix inequality to specify nonnegativity of `x`. #. Specify the objective function via :meth:`~SemidefiniteProgram.set_objective`. In our case it is `x_1 - x_0`. If it is a pure constraint satisfaction problem, specify it as ``None``. #. To check if everything is set up correctly, you can print the problem via :meth:`show <sage.numerical.sdp.SemidefiniteProgram.show>`. #. :meth:`Solve <sage.numerical.sdp.SemidefiniteProgram.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: 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(round(opt,3))) Objective Value: 1.0 sage: [round(x,3) for x in p.get_values(x).values()] [0.0, 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 :meth:`~SemidefiniteProgram.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 :meth:`~SemidefiniteProgram.slack`. E.g.::
sage: (p.slack(0)*p.dual_variable(0)).trace() # tol 2e-07 0.0
More interesting example, the :func:`Lovasz theta <sage.graphs.lovasz_theta.lovasz_theta>` of the 7-gon::
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 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[1] 0.5 + 1.5*x_0
Linear Variables and Expressions --------------------------------
To make your code more readable, you can construct :class:`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
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: <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 :class:`SemidefiniteProgram`. This module also implements the :class:`SDPSolverException` exception, as well as the :class:`SDPVariable` class.
.. csv-table:: :class: contentstable :widths: 30, 70 :delim: |
:meth:`~SemidefiniteProgram.add_constraint` | Adds a constraint to the ``SemidefiniteProgram`` :meth:`~SemidefiniteProgram.base_ring` | Return the base ring :meth:`~SemidefiniteProgram.dual_variable` | Return optimal dual variable block :meth:`~SemidefiniteProgram.get_backend` | Return the backend instance used :meth:`~SemidefiniteProgram.get_values` | Return values found by the previous call to ``solve()`` :meth:`~SemidefiniteProgram.linear_constraints_parent` | Return the parent for all linear constraints :meth:`~SemidefiniteProgram.linear_function` | Construct a new linear function :meth:`~SemidefiniteProgram.linear_functions_parent` | Return the parent for all linear functions :meth:`~SemidefiniteProgram.new_variable` | Return an instance of ``SDPVariable`` associated to the ``SemidefiniteProgram`` :meth:`~SemidefiniteProgram.number_of_constraints` | Return the number of constraints assigned so far :meth:`~SemidefiniteProgram.number_of_variables` | Return the number of variables used so far :meth:`~SemidefiniteProgram.set_objective` | Set the objective of the ``SemidefiniteProgram`` :meth:`~SemidefiniteProgram.set_problem_name` | Set the name of the ``SemidefiniteProgram`` :meth:`~SemidefiniteProgram.slack` | Return the slack variable block at the optimum :meth:`~SemidefiniteProgram.show` | Display the ``SemidefiniteProgram`` in a human-readable way :meth:`~SemidefiniteProgram.solve` | Solve the ``SemidefiniteProgram`` :meth:`~SemidefiniteProgram.solver_parameter` | Return or define a solver parameter :meth:`~SemidefiniteProgram.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
""" #***************************************************************************** # Copyright (C) 2014 Ingolfur Edvardsson <ingolfured@gmail.com> # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #***************************************************************************** from __future__ import print_function
from sage.structure.parent cimport Parent from sage.structure.element cimport Element
cdef class SemidefiniteProgram(SageObject): r""" 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:`Semidefinite_programming` for further information on semidefinite programming, and the :mod:`SDP module <sage.numerical.sdp>` for its use in Sage.
INPUT:
- ``solver`` -- selects a solver:
- CVXOPT (``solver="CVXOPT"``). See the `CVXOPT <http://www.cvxopt.org/>`_ website.
- If ``solver=None`` (default), the default solver is used (see :func:`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.
.. SEEALSO::
- :func:`default_sdp_solver` -- Returns/Sets the default SDP solver.
EXAMPLES:
Computation of a basic Semidefinite Program::
sage: p = SemidefiniteProgram(solver = "cvxopt", 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: round(p.solve(), 2) -3.0 """
def __init__(self, solver=None, maximization=True, r""" Constructor for the ``SemidefiniteProgram`` class.
INPUT:
- ``solver`` -- the following solvers should be available through this class:
- CVXOPT (``solver="CVXOPT"``). See the `CVXOPT <http://www.cvxopt.org/>`_ web site.
-If ``solver=None`` (default), the default solver is used (see ``default_sdp_solver`` method.
- ``maximization``
- When set to ``True`` (default), the ``SemidefiniteProgram`` is defined as a maximization. - When set to ``False``, the ``SemidefiniteProgram`` is defined as a minimization.
- ``names`` -- list/tuple/iterable of string. Default names of the SDP variables. Used to enable the ``sdp.<x> = SemidefiniteProgram()`` syntax.
.. SEEALSO::
- :meth:`default_sdp_solver` -- Returns/Sets the default SDP solver.
EXAMPLES::
sage: p = SemidefiniteProgram(maximization=True)
"""
# Associates an index to the variables
def linear_functions_parent(self): """ Return the parent for all linear functions.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: p.linear_functions_parent() Linear functions over Real Double Field """
def linear_constraints_parent(self): """ Return the parent for all linear constraints.
See :mod:`~sage.numerical.linear_functions` for more details.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: p.linear_constraints_parent() Linear constraints over Real Double Field """
def __call__(self, x): """ Construct a new linear function.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: p.linear_function({0:1}) x_0 """
def _repr_(self): r""" Returns a short description of the ``SemidefiniteProgram``.
EXAMPLES::
sage: p = SemidefiniteProgram() 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] + b2*x[1] <= a1) sage: print(p) Semidefinite Program ( maximization, 2 variables, 3 constraints ) """
def __getitem__(self, v): r""" Returns the symbolic variable corresponding to the key from a default dictionary.
It returns the element asked, and otherwise creates it. If necessary, it also creates the default dictionary.
This method lets the user define LinearProgram without having to define independent dictionaries when it is not necessary for him.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: p.set_objective(p['x'] + p['z']) sage: p['x'] x_0 """
def base_ring(self): """ 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 """
def set_problem_name(self,name): r""" 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 ) """
def new_variable(self, name=""): r""" 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[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: <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 """
def _first_ngens(self, n): """ Construct the first `n` SDPVariables.
This method is used for the generater syntax (see below). You probably shouldn't use it for anything else.
INPUT:
- ``n`` -- integer. The number of variables to construct.
OUTPUT:
A tuple of not necessarily positive :class:`SDPVariable` instances.
EXAMPLES::
sage: sdp.<a,b> = SemidefiniteProgram() sage: a[0] + b[2] x_0 + x_1 sage: sdp.show() Maximization: <BLANKLINE> Constraints: Variables: a[0], b[2] """
def gen(self, 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] """
cpdef int number_of_constraints(self): r""" 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[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 """
cpdef int number_of_variables(self): r""" 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[0] - a*p[2] <= 2*a*p[4] ) sage: p.number_of_variables() 3 """
def show(self): r""" 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[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] """ cdef int i, j
# inv_variables associates a SDPVariable object to an id
# varid_name associates variables id to names
##### Sense and objective function print("+ {}".format(b.obj_constant_term)) print("- {}".format(-b.obj_constant_term))
##### Constraints # Constraint's name else: continue
##### Variables
def get_values(self, *lists): r""" 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: round(p.solve(),3) -3.0
To return the optimal value of ``x[3]``::
sage: round(p.get_values(x[3]),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: x_sol.keys() [3, 5]
Obviously, it also works with variables of higher dimension::
sage: x_sol = p.get_values(x)
""" if len(l) == 1: val.append([self.get_values(l[0])]) else: c = [] [c.append(self.get_values(ll)) for ll in l] val.append(c) #val.append(self._values[l])
else: return val
def set_objective(self,obj): r""" 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:
.. MATH::
\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}
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: round(p.solve(),5) 16.2 sage: p.set_objective(None) sage: _ = p.solve() """
# If the objective is None, or a constant, we want to remember # that the objective function has been defined ( the user did not # forget it ). In some SDO problems, you just want a feasible solution # and do not care about any function being optimal. cdef int i
else:
def add_constraint(self, linear_function, name=None): r""" 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:
.. MATH::
\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}
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: round(p.solve(),5) 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: round(p.solve(),5) 16.2
TESTS:
Complex constraints::
sage: p = SemidefiniteProgram(solver = "cvxopt") sage: b = p.new_variable() sage: a1 = matrix([[1,2],[2,3]]) sage: a2 = matrix([[1,-2],[-2,4]]) sage: p.add_constraint(a1*b[8] - a1*b[15] <= a2*b[8]) sage: p.show() Maximization: <BLANKLINE> Constraints: constraint_0: [ 0.0 4.0][ 4.0 -1.0]x_0 + [-1.0 -2.0][-2.0 -3.0]x_1 <= [0 0][0 0] Variables: x_0, x_1
Empty constraint::
sage: p=SemidefiniteProgram() sage: p.add_constraint(sum([]))
"""
else:
else: raise ValueError('argument must be a linear function or constraint, got '+str(linear_function))
def solve(self, objective_only=False): r""" 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.
TESTS:
The SDP from the header of this module::
sage: p = SemidefiniteProgram(solver = "cvxopt", 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., 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: round(p.solve(),4) -11.0 sage: x = p.get_values(x) sage: round(x[0],4) -8.0 sage: round(x[1],4) 3.0 """
cpdef dual_variable(self, int 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[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: 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
TESTS::
sage: p.dual_variable(7) ... Traceback (most recent call last): ... IndexError: list index out of range """
cpdef slack(self, int 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[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: 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 = 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]
TESTS::
sage: p.slack(7) ... Traceback (most recent call last): ... IndexError: list index out of range """
def solver_parameter(self, 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[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: round(p.solve(),4) pcost dcost gap pres dres k/t 0: 1... ... Optimal solution found. -11.0 """ else:
cpdef sum(self, L): r""" Efficiently computes the sum of a sequence of :class:`~sage.numerical.linear_functions.LinearFunction` elements.
INPUT:
- ``L`` -- list of :class:`~sage.numerical.linear_functions.LinearFunction` instances.
.. 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)) """
def get_backend(self): r""" 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[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] ) """
r""" 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',)
TESTS:
No solution::
sage: p=SemidefiniteProgram(solver="cvxopt") sage: x=p.new_variable() sage: p.set_objective(x[0]) sage: a = matrix([[1,2],[2,4]]) sage: b = matrix([[1,9],[9,4]]) sage: p.add_constraint( a*x[0] == b ) sage: p.solve() ... Traceback (most recent call last): ... SDPSolverException: ...
The value of the exception::
sage: from sage.numerical.sdp import SDPSolverException sage: e = SDPSolverException("Error") sage: print(e) Error """ pass
cdef class SDPVariable(Element): r""" ``SDPVariable`` is a variable used by the class ``SemidefiniteProgram``.
.. warning::
You should not instantiate this class directly. Instead, use :meth:`SemidefiniteProgram.new_variable`. """
def __init__(self, parent, sdp, name): r""" Constructor for ``SDPVariable``.
INPUT:
- ``parent`` -- :class:`SDPVariableParent`. The parent of the SDP variable.
- ``sdp`` -- :class:`SemidefiniteProgram`. The underlying linear program.
- ``name`` -- A name for the ``SDPVariable``.
- ``lower_bound``, ``upper_bound`` -- lower bound and upper bound on the variable. Set to ``None`` to indicate that the variable is unbounded.
For more informations, see the method ``SemidefiniteProgram.new_variable``.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: p.new_variable() SDPVariable """
def __getitem__(self, i): r""" Returns the symbolic variable corresponding to the key.
Returns the element asked, otherwise creates it.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: v = p.new_variable() sage: p.set_objective(v[0] + v[1]) sage: v[0] x_0
""" cdef int j
def _repr_(self): r""" Returns a representation of self.
EXAMPLES::
sage: p=SemidefiniteProgram() sage: v=p.new_variable() sage: v SDPVariable """
def keys(self): r""" 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: list(v.keys()) [0, 1] """
def items(self): r""" 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: list(v.items()) [(0, x_0), (1, x_1)] """
def values(self): r""" 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: list(v.values()) [x_0, x_1] """
cdef _matrix_rmul_impl(self, m): """ Implement the action of a matrix multiplying from the right. """
cdef _matrix_lmul_impl(self, m): """ Implement the action of a matrix multiplying from the left. """
cpdef _acted_upon_(self, mat, bint self_on_left): """ Act with matrices on SDPVariables.
EXAMPLES::
sage: p = SemidefiniteProgram() sage: v = p.new_variable() sage: m = matrix([[1,2], [3,4]]) sage: v * m (1.0, 2.0)*x_0 + (3.0, 4.0)*x_1 sage: m * v (1.0, 3.0)*x_0 + (2.0, 4.0)*x_1 """
cdef class SDPVariableParent(Parent): """ Parent for :class:`SDPVariable`.
.. warning::
This class is for internal use. You should not instantiate it yourself. Use :meth:`SemidefiniteProgram.new_variable` to generate sdp variables. """
def _repr_(self): r""" Return representation of self.
OUTPUT:
String.
EXAMPLES::
sage: sdp.<v> = SemidefiniteProgram() sage: v.parent() Parent of SDPVariables """
def _an_element_(self): """ Construct a SDP variable.
OUTPUT:
This is required for the coercion framework. We raise a ``TypeError`` to abort search for any coercion to another parent for binary operations. The only interesting operations involving :class:`SDPVariable` elements are actions by matrices.
EXAMPLES::
sage: sdp.<x> = SemidefiniteProgram() sage: parent = x.parent() sage: parent.an_element() # indirect doctest Traceback (most recent call last): ... TypeError: disallow coercion """
def _element_constructor_(self, sdp, name=""): """ The Element constructor
INPUT/OUTPUT:
See :meth:`SDPVariable.__init__`.
EXAMPLES::
sage: sdp = SemidefiniteProgram() sage: sdp.new_variable() # indirect doctest SDPVariable """
|