Previous  Up  Next

Module Sdp

module Sdp: sig .. end

Common interface for SDP (Csdp or Moseksdp or Sdpa).

Primal-dual correspondence: the primal problem

max tr(C X)
    tr(A_1 X) = b_1
    .
    .
    .
    tr(A_m X) = b_m
    X psd

(X psd meaning X positive semi-definite) corresponds to the dual problem

min b^T y
    \sum_i y_i A_i - C = Z psd

where C, A_i and X are symmetric matrices, C, A_i and b_i are parameters whereas X (resp. y) is the primal (resp. dual) variable.

Block diagonal sparse and dense symmetric matrices.

type sparse_matrix = (int * int * float) list

Matrices. Sparse representation as triplet (i, j, x) meaning that the coefficient at line i >= 0 and column j >= 0 has value x. All forgotten coefficients are assumed to be 0.0. Since matrices are symmetric, only the lower triangular part (j <= i) must be given. No duplicates are allowed.

type matrix = float array array

Matrices. Dense representation, line by line. Must be symmetric. Type invariant: all lines have the same size.

type 'a block_diag = (int * 'a) list

Block diagonal matrices (sparse representation, forgetting null blocks). For instance, [(1, m1), (3, m2)] will be transformed into [m1; 0; m2]. No duplicates are allowed. There is no requirement for indices to be sorted.

val matrix_of_sparse : sparse_matrix -> matrix
val matrix_to_sparse : matrix -> sparse_matrix
val block_diag_of_sparse : sparse_matrix block_diag -> matrix block_diag
val block_diag_to_sparse : matrix block_diag -> sparse_matrix block_diag

SDP.

type solver = Sdp_default.solver =
| Csdp
| Mosek
| Sdpa
| SdpaGmp
| SdpaDd 
type options = {
   solver : solver; (*

default: see Sdp_default

*)
   verbose : int; (*

verbosity level, non negative integer, 0 (default) means no output (only for SDPA* and Mosek)

*)
   max_iteration : int; (*

maxIteration (default: 100)

*)
   stop_criterion : float; (*

epsilonStar and epsilonDash (default: 1.0E-7)

*)
   initial : float; (*

lambdaStar (default: 1.0E2)

*)
   precision : int; (*

precision (only for SDPA-GMP, default: 200)

*)
}

Options for calling SDP solvers (currently, options other than solver are only handled by Sdpa{,Gmp}).

val default : options

Default values above.

type 'a obj = 'a block_diag

Objective (matrix C).

type 'a constr =
| Eq of 'a block_diag * float
| Le of 'a block_diag * float
| Ge of 'a block_diag * float

Constraints (tr(A_i X) = b_i), Le (A_i, b_i) means tr(A_i X) <= b_i and will automatically be translated internally as tr(A_i X) + s_i = b_i by adding slack variables s_i >= 0. Similarly Ge (A_i, b_i) means tr(A_i X) >= b_i and gives tr(A_i X) - s_i = b_i with s_i >= 0.

val solve_sparse : ?options:options ->
       ?solver:solver ->
       ?init:matrix block_diag * float array * matrix block_diag ->
       sparse_matrix obj ->
       sparse_matrix constr list ->
       SdpRet.t * (float * float) *
       (matrix block_diag * float array * matrix block_diag)

solve_sparse obj constraints solves the SDP problem: max{ tr(obj X) | constraints, X psd }. If solver is provided, it will supersed the solver given in options. It returns both the primal and dual objective values and a witness for X (primal) and y and Z (dual). See above for details. There is no requirement for indices in obj and A_i to be present in every input. In case of success (or partial success), the block diagonal matrices returned for X and Z contain exactly the indices, sorted by increasing order, that appear in the objective or one of the constraints. Size of each diagonal block in X is the maximum size appearing for that block in the objective or one of the constraints. In case of success (or partial success), the array returned for y has the same size and same order than the input list of constraints.

val solve : ?options:options ->
       ?solver:solver ->
       ?init:matrix block_diag * float array * matrix block_diag ->
       matrix obj ->
       matrix constr list ->
       SdpRet.t * (float * float) *
       (matrix block_diag * float array * matrix block_diag)

Same as Sdp.solve_sparse with dense matrices as input. This can be more convenient for small problems. Otherwise, this is equivalent to Sdp.solve_sparse after conversion with Sdp.matrix_to_sparse. If an indice is present in one input, potential matrices for that indice in other inputs will be padded with 0 on right and bottom to have the same size.

Extended formulation.

Primal-dual correspondence: the primal problem

max c^T x + tr(C X)
    b_1^- <= a_1^T x + tr(A_1 X) <= b_1^+
    .
    .
    .
    b_m^- <= a_m^T x + tr(A_m X) <= b_m^+
    d_1^- <= x_1 <= d_1^+
    .
    .
    .
    d_n^- <= x_n <= d_n^+
    X psd

(X psd meaning X positive semi-definite) corresponds to the dual problem

min (b^+)^T s^+ - (b^-)^T s^- + (d^+)^T t^+ - (d^-)^T t^-
    y = s^+ - s^-
    \sum_i y_i A_i - C = Z psd
    [a_1,..., a_n] y - c + t^+ - t^- = 0
    s^+, s^-, t^+, t^- >= 0

where C, A_i and X are symmetric matrices, C, A_i, a_i, b_i and d are parameters whereas x, X, y, s^+, s^-, t^+ and t^- are variables.

Note that the simple formulation on top is a particular case of this for n = 0 and b^- = b^+.

The bounds b_i^- and d_j^- (resp. b_i^+ and d_j^+) can be neg_infinity (resp. infinity), in which case, the corresponding s or t in the dual is 0 (and the corresponding term disappears from the dual objective).

type vector = (int * float) list

No duplicates are allowed.

type 'a obj_ext = vector * 'a block_diag

Objective (vector c and matrix C).

type 'a constr_ext = vector * 'a block_diag * float * float

Constraints. (a_i, A_i, b_i_m, b_i_p) encodes the constraint b_i_m <= a_i^T x + tr(A_i X) <= b_i_p. b_i_m can be neg_infinity and b_i_p can be infinity.

type bounds = (int * float * float) list

Bounds on scalar variables. Each (i, lb, ub) in the list encodes the constraint lb <= x_i <= ub. lb can be neg_infinity and ub can be infinity. No duplicated indices i are allowed.

val solve_ext_sparse : ?options:options ->
       ?solver:solver ->
       sparse_matrix obj_ext ->
       sparse_matrix constr_ext list ->
       bounds ->
       SdpRet.t * (float * float) *
       (vector * matrix block_diag * float array *
        matrix block_diag)

solve_ext_sparse obj constraints bounds solves the SDP problem: max{ obj | constraints, bounds, X psd }. If solver is provided, it will supersed the solver given in options. It returns both the primal and dual objective values and a witness for (x, X) (primal) and (y, Z) (dual). See above for details. Variables x_i not bounded in bounds are assumed to be unbounded (i.e., bounded by neg_infinity and infinity). Bounds in bounds about variables x_i not appearing in the objective or constraints may be ignored. There is no requirement for indices in obj and a_i, A_i to be present in every input. In case of success (or partial success), the vector (resp. block diagonal matrix) returned for x (resp. X, Z) contains exactly the indices, sorted by increasing order, that appear in the linear (resp. matrix) part of the objective or one of the constraints. Size of each diagonal block in X or Z is the maximum size appearing for that block in the objective or one of the constraints. In case of success (or partial success), the array returned for y has the same size and same order than the input list of constraints.

val solve_ext : ?options:options ->
       ?solver:solver ->
       matrix obj_ext ->
       matrix constr_ext list ->
       bounds ->
       SdpRet.t * (float * float) *
       (vector * matrix block_diag * float array *
        matrix block_diag)

Same as Sdp.solve_ext_sparse with dense matrices as input. This can be more convenient for small problems. Otherwise, this is equivalent to Sdp.solve_ext_sparse after conversion with Sdp.matrix_to_sparse. If an indice is present in one input, potential matrices for that indice in other inputs will be padded with 0 on right and bottom to have the same size.

Printing functions.

val pp_sparse_matrix : Stdlib.Format.formatter -> sparse_matrix -> unit
val pp_matrix : Stdlib.Format.formatter -> matrix -> unit
val pp_block_diag : (Stdlib.Format.formatter -> 'a -> unit) ->
       Stdlib.Format.formatter -> 'a block_diag -> unit
val pp_obj : (Stdlib.Format.formatter -> 'a -> unit) ->
       Stdlib.Format.formatter -> 'a obj -> unit
val pp_constr : (Stdlib.Format.formatter -> 'a -> unit) ->
       Stdlib.Format.formatter -> 'a constr -> unit
val pp_sparse : Stdlib.Format.formatter ->
       sparse_matrix obj * sparse_matrix constr list -> unit
val pp : Stdlib.Format.formatter ->
       matrix obj * matrix constr list -> unit
val pp_vector : Stdlib.Format.formatter -> vector -> unit
val pp_obj_ext : (Stdlib.Format.formatter -> 'a -> unit) ->
       Stdlib.Format.formatter -> 'a obj_ext -> unit
val pp_constr_ext : (Stdlib.Format.formatter -> 'a -> unit) ->
       Stdlib.Format.formatter -> 'a constr_ext -> unit
val pp_bounds : Stdlib.Format.formatter -> bounds -> unit
val pp_ext_sparse : Stdlib.Format.formatter ->
       sparse_matrix obj_ext * sparse_matrix constr_ext list *
       bounds -> unit
val pp_ext : Stdlib.Format.formatter ->
       matrix obj_ext * matrix constr_ext list * bounds -> unit
val pp_ext_sparse_sedumi : Stdlib.Format.formatter ->
       sparse_matrix obj_ext * sparse_matrix constr_ext list *
       bounds -> unit
val pp_ext_sedumi : Stdlib.Format.formatter ->
       matrix obj_ext * matrix constr_ext list * bounds -> unit

Miscellaneous functions.

val pfeas_stop_crit : ?options:options -> ?solver:solver -> float list -> float

pfeas_stop_crit (b_1,..., b_m) returns eps used as primal feasibility error stopping criteria by solver solver on a problem with given b_i (i.e., if the solver successfully terminates, |tr(A_i X) - b_i| <= eps should hold for each constraint i).