SMCP — Python extension for Sparse Matrix Cone Programs¶

Installation with pip¶

SMCP can be installed via pip using the following command:

pip install smcp

Installation from source¶

The package requires Python version 2.7 or newer. To build the package from source, the Python header files and libraries must be installed, as well as the core binaries. SMCP also requires the Python extension modules CVXOPT 1.1.9 or later and CHOMPACK 2.3.2 or later.

The source package is available here. The extension can be built and installed as follows:

python setup.py install
python example.py

Overview¶

Let \(\mathbf{S}^n\) denote the set of symmetric matrices of order \(n\), and let and let \(C \bullet X\) denote the standard inner product on \(\mathbf{S}^n\).

\(\mathbf{S}_V^n\) is the set of symmetric matrices of order \(n\) and with sparsity pattern \(V\), i.e., \(X \in \mathbf{S}_V^n\) if and only if \(X_{ij} = 0\) for all \((i,j) \neq V\). We will assume that all diagonal elements are included in \(V\). The projection \(Y\) of \(X\) on the subspace \(\mathbf{S}_V^n\) is denoted \(Y=P_V(X)\), i.e., \(Y_{ij} = X_{ij}\) if \((i,j) \in V\) and \(Y_{ij} = 0\) otherwise. The inequality signs \(\succeq\) and \(\succ\) denote matrix inequality. We define \(|V|\) as the number of nonzeros in the lower triangular part of \(V\), and \(\mathrm{nnz}(A_i)\) denotes the number of nonzeros in the matrix A_i.

\(\mathbf{S}_{V,+}^n\) and \(\mathbf{S}_{V,++}^n\) are the sets of positive semidefinite and positive definite matrices in \(\mathbf{S}_V^n\), and similarly, \(\mathbf{S}_{V,c+}^n = \{ P_V(X)\ |\ X \succeq 0 \}\) and \(\mathbf{S}_{V,c++}^n = \{ P_V(X)\ |\ X \succ 0 \}\) are the sets of matrices in \(\mathbf{S}_V^n\) that have a positive semidefinite completion and a positive definite completion, respectively. We denote with \(\succeq_{\rm c}\) and \(\succ_{\rm c}\) matrix inequality with respect to the cone \(\mathbf{S}_{V,c+}^n\).

SMCP solves a pair of primal and dual linear cone programs:

The variables are \(X \in \mathbf{S}_V^n\), \(S \in \mathbf{S}_V^n\), and \(y \in \mathbf{R}^m\), and the problem data are the matrix \(C \in \mathbf{S}_V^n\), the vector \(b \in \mathbf{R}^m\), and \(A_i \in \mathbf{S}_V^n, i=1,\ldots,m\).

Compositions of cones are handled implicitly by defining a block diagonal sparsity pattern \(V\). Dense blocks and general sparse blocks correspond to standard positive semidefinite matrix constraints, diagonal blocks corresponds to linear inequality constraints, and second-order cone constraints can be embedded in an LMI with an “arrow pattern”, i.e,

The chordal SDP solvers¶

smcp.solvers.chordalsolver_feas(A, b[, primalstart[, dualstart[, scaling='primal'[, kktsolver='chol']]]])¶

Solves the pair of cone programs (P) and (D) using a feasible start interior-point method. If no primal and/or dual feasible starting point is specified, the algorithm tries to find a feasible starting point based on some simple heuristics. An exception is raised if no starting point can be found. In this case a Phase I problem must be solved, or the (experimental) infeasible start interior-point method chordalsolver_esd can be used.

The columns of the sparse matrix A are vectors of length \(n^2\) and the \(m+1\) columns of A are:

\[\left[ \mathbf{vec}(C) \ \mathbf{vec}(A_1) \ \cdots \ \mathbf{vec}(A_m) \right].\]

Only the rows of A corresponding to the lower triangular elements of the aggregate sparsity pattern \(V\) are accessed.

The optional argument primalstart is a dictionary with the key x which can be used to specify an initial value for the primal variable \(X\). Similarly, the optional argument dualstart must be a dictionary with keys s and y.

The optional argument scaling takes one of the values 'primal' (default) or 'dual'.

The optional argument kktsolver is used to specify the KKT solver. Possible values include:

'chol' (default)

solves the KKT system via a Cholesky factorization of the Schur complement

'qr'

solves the KKT system via a QR factorization

The solver returns a dictionary with the following keys:

'primal objective', 'dual objective'

primal objective value and dual objective value.

'primal infeasibility', 'dual infeasibility'

residual norms of primal and dual infeasibility.

'x', 's', and 'y'

primal and dual variables.

'iterations'

number of iterations.

'cputime', 'time'

total cputime and real time.

'gap'

duality gap.

'relative gap'

relative duality gap.

'status'

• has the value 'optimal' if

  \[\frac{\| b-\mathcal{A}(X) \|_2}{\max\{1,\|b\|_2 \}} \leq \epsilon_{\rm feas}, \qquad \frac{\|\mathcal{A}^{\rm adj}(y) + S - C\|_F}{\max\{1,\|C\|_F \}} \leq \epsilon_{\rm feas}, \qquad X \succ_{\rm c} 0,\qquad S \succ 0,\]

  and

  \[X\bullet S \leq \epsilon_{\rm abs} \quad \text{or} \quad \left( \min\{C\bullet X, -b^Ty\} \leq 0 , \frac{X\bullet S}{-\min\{C\bullet X,-b^Ty \}} \leq \epsilon_{\rm rel} \right).\]

• has the value 'unknown' otherwise.

The following options can be set using the dictionary smcp.solvers.options:

'delta' (default: 0.9)

a positive constant between 0 and 1; an approximate tangent direction is computed when the Newton decrement is less than delta.

'eta' (default: None)

None or a positive float. If 'eta' is a positive number, a step in the approximate tangent direction is taken such that

\[\Omega(X+\alpha \Delta X, S + \alpha \Delta S) \approx \eta\]

where \(\Omega(X,S)\) is the proximity function

\[\Omega(X,S) = \phi_{\rm c}(X) + \phi(S) + n\cdot \log\frac{X\bullet S}{n} + n.\]

If 'eta' is None, the step length \(\alpha\) in the approximate tangent direction is computed as

\[\begin{split}\alpha_p &= \arg \max \{ \alpha \in (0,1] \,|\, X + \alpha \Delta X \succeq_{\rm c} 0 \} \\ \alpha_d &= \arg \max \{ \alpha \in (0,1] \,|\, S + \alpha \Delta S \succeq 0 \} \\ \alpha &= \texttt{step}\cdot\min(\alpha_p,\alpha_d)\end{split}\]

where \(\texttt{step}\) is the value of the option 'step' (default: 0.98).

'prediction' (default: True)

True or False. This option is effective only when 'eta' is None. If 'prediction' is True, a step in the approximate tangent direction is never taken but only used to predict the duality gap. If 'prediction' is False, a step in the approximate tangent direction is taken.

'step' (default: 0.98)

positive float between 0 and 1.

'lifting' (default: True)

True or False; determines whether or not to apply lifting before taking a step in the approximate tangent direction.

'show_progress'

True or False; turns the output to the screen on or off (default: True).

'maxiters'

maximum number of iterations (default: 100).

'abstol'

absolute accuracy (default: 1e-6).

'reltol'

relative accuracy (default: 1e-6).

'feastol'

tolerance for feasibility conditions (default: 1e-8).

'refinement'

number of iterative refinement steps when solving KKT equations (default: 1).

'cholmod'

use Cholmod’s AMD embedding (defaults: False).

'dimacs'

report DIMACS error measures (default: True).

smcp.solvers.chordalsolver_esd(A, b[, primalstart[, dualstart[, scaling='primal'[, kktsolver='chol']]]])¶

Solves the pair of cone programs (P) and (D) using an extended self-dual embedding. This solver is currently experimental.

The columns of the sparse matrix A are vectors of length \(n^2\) and the \(m+1\) columns of A are:

\[\left[ \mathbf{vec}(C) \ \mathbf{vec}(A_1) \ \cdots \ \mathbf{vec}(A_m) \right].\]

Only the rows of A corresponding to the lower triangular elements of the aggregate sparsity pattern \(V\) are accessed.

The optional argument primalstart is a dictionary with the key x which can be used to specify an initial value for the primal variable \(X\). Similarly, the optional argument dualstart must be a dictionary with keys s and y.

The optional argument scaling takes one of the values 'primal' (default) or 'dual'.

The optional argument kktsolver is used to specify the KKT solver. Possible values include:

'chol' (default)

solves the KKT system via a Cholesky factorization of the Schur complement

'qr'

solves the KKT system via a QR factorization

The solver returns a dictionary with the following keys:

'primal objective', 'dual objective'

primal objective value and dual objective value.

'primal infeasibility', 'dual infeasibility'

residual norms of primal and dual infeasibility.

'x', 's', and 'y'

primal and dual variables.

'iterations'

number of iterations.

'cputime', 'time'

total cputime and real time.

'gap'

duality gap.

'relative gap'

relative duality gap.

'status'

• has the value 'optimal' if

  \[\frac{(1/\tau) \| \tau b-\mathcal{A}(X) \|_2}{\max\{1,\|b\|_2 \}} \leq \epsilon_{\rm feas}, \qquad \frac{(1/\tau) \|\mathcal{A}^{\rm adj}(y) + S - \tau C\|_F}{\max\{1,\|C\|_F \}} \leq \epsilon_{\rm feas}, \qquad X \succ_{\rm c} 0,\qquad S \succ 0,\]

  and

  \[\frac{X\bullet S}{\tau^2} \leq \epsilon_{\rm abs} \quad \text{or} \quad \left( \min\{C\bullet X, -b^Ty\} \leq 0 , \frac{(1/\tau) X\bullet S}{-\min\{C\bullet X,-b^Ty \}} \leq \epsilon_{\rm rel} \right).\]

• has the value 'primal infeasible' if

  \[b^Ty = 1, \qquad \frac{\| \mathcal{A}^{\rm adj}(y) + S \|_F}{\max\{1,\|C\|_F\}} \leq \epsilon_{\rm feas}, \qquad S \succ 0.\]

• has the value 'dual infeasible' if

  \[C\bullet X = -1, \qquad \frac{\| \mathcal{A}(X) \|_2}{\max\{1,\|b\|_2\}} \leq \epsilon_{\rm feas}, \qquad X \succ_{\rm c} 0.\]

• has the value 'unknown' if maximum number iterations is reached or if a numerical error is encountered.

The following options can be set using the dictionary smcp.solvers.options:

'show_progress'

True or False; turns the output to the screen on or off (default: True).

'maxiters'

maximum number of iterations (default: 100).

'abstol'

absolute accuracy (default: 1e-6).

'reltol'

relative accuracy (default: 1e-6).

'feastol'

tolerance for feasibility conditions (default: 1e-8).

'refinement'

number of iterative refinement steps when solving KKT equations (default: 1).

'cholmod'

use Cholmod’s AMD embedding (defaults: False).

'dimacs'

report DIMACS error measures (default: True).

Solver interfaces¶

The following functions implement CVXOPT-like interfaces to the experimental solver chordalsolver_esd.

smcp.solvers.conelp(c, G, h[, dims[, kktsolver='chol']])¶

Interface to chordalsolver_esd.

smcp.solvers.lp(c, G, h[, kktsolver='chol'])¶

Interface to conelp; see CVXOPT documentation for more information.

smcp.solvers.socp(c[, Gl, hl[, Gq, hq[, kktsolver='chol']]])¶

Interface to conelp; see CVXOPT documentation for more information.

smcp.solvers.sdp(c[, Gl, hl[, Gs, hs[, kktsolver='chol']]])¶

Interface to conelp; see CVXOPT documentation for more information.

The SDP object¶

class SDP(filename)¶

Class for SDP problems. Simplifies reading and writing SDP data files and includes a wrapper for chordalsolver_esd.

The constructor accepts sparse SDPA data files (extension ‘dat-s’) and data files created with the save method (extension ‘pkl’). Data files compressed with Bzip2 can also be read (extensions ‘dat-s.bz2’ and ‘pkl.bz2’).

m¶

Number of constraints.

n¶

Order of semidefinite variable.

A¶

Problem data: sparse matrix of size \(n^2 \times (m+1)\) with columns \(\mathbf{vec}(C),\mathbf{vec}(A_1),\ldots,\mathbf{vec}(A_m)\). Only the lower triangular elements of \(C,A_1,\ldots,A_m\) are stored.

b¶

Problems data: vector of length \(m\).

V¶

Sparse matrix with aggregate sparsity pattern (lower triangle).

nnz¶

Number of nonzero elements in lower triangle of aggregate sparsity pattern.

nnzs¶

Vector with number of nonzero elements in lower triangle of \(A_0,\ldots,A_m\).

nzcols¶

Vector with number of nonzero columns in \(A_1,\ldots,A_m\).

issparse¶

True if the number of nonzeros is less than \(0.5 \cdot n(n+1)/2\), otherwise false.

ischordal¶

True if aggregate sparsity pattern is chordal, otherwise false.

get_A(i)¶

Returns the \(i\)‘th coeffiecient matrix \(A_i\) (\(0\leq i \leq m\)) as a sparse matrix. Only lower triangular elements are stored.

write_sdpa([fname[, compress=False]])¶

Writes SDP data to SDPA sparse data file. The extension ‘dat-s’ is automatically added to the filename. The method is an interface to sdpa_write.

If compress is true, the data file is compressed with Bzip2 and ‘bz2’ is appended to the filename.

save([fname[, compress=False]])¶

Writes SDP data to file using cPickle. The extension ‘pkl’ is automatically added to the filename.

If compress is true, the data file is compressed with Bzip2 and ‘bz2’ is appended to the filename.

solve_feas([scaling='primal'[, kktsolver='chol'[, primalstart[, dualstart]]]])¶

Interface to the feasible start solver chordalsolver_feas. Returns dictionary with solution.

solve_phase1([kktsolver='chol'[, M=1e5]])¶

Solves a Phase I problem to find a feasible (primal) starting point:

\[\begin{split}\begin{array}{ll} \mbox{minimize} & s \\ \mbox{subject to} & A_i \bullet X = b_i, \quad i=i,\ldots,m \\ & \mathbf{tr}(X) \leq M \\ & X + (s-\epsilon)I \succeq_{\rm c} 0,\, s \geq 0 \end{array}\end{split}\]

The variables are \(X \in \mathbf{S}_V^n\) and \(s \in \mathbf{R}\), and \(\epsilon \in \mathbf{R}_{++}\) is a small constant. If \(s^\star < \epsilon\), the method returns \(X^\star\) (which is a strictly feasible starting point in the original problem) and a dictionary (with information about the Phase I problem). If \(s >= \epsilon\) the method returns (None, None).

solve_esd([scaling='primal'[, kktsolver='chol'[, primalstart[, dualstart]]]])¶

Interface to chordalsolver_esd. Returns dictionary with solution.

solve_cvxopt([primalstart[, dualstart]])¶

Interface to cvxopt.solvers.sdp(). Returns dictionary with solution. (Note that this simple interface does not yet specify block structure properly.)

The following example demostrates how to load and solve a problem from an SDPA sparse data file:

>>> from smcp import SDP
>>> P = SDP('qpG11.dat-s')
>>> print P
<SDP: n=1600, m=800, nnz=3200> qpG11
>>> sol = P.solve_feas(kktsolver='chol')
>>> print sol['primal objective']
-2448.6588977
>>> print sol['dual objective']
-2448.65913565
>>> print sol['gap']
0.00023794772363
>>> print sol['relative gap']
9.71747121876e-08

Auxiliary routines¶

smcp.completion(X)¶

Computes the maximum determinant positive definite completion of a sparse matrix X.

Example:

>>> from smcp import mtxnorm_SDP, completion
>>> P = mtxnorm_SDP(p=10,q=2,r=10)
>>> sol = P.solve_feas(kktsolver='chol')
>>> X = completion(sol['x'])

smcp.misc.ind2sub(siz, ind)¶

Converts indices to subscripts.

Parameters:	siz (integer) – matrix order ind (matrix) – vector with indices
Returns:	matrix I with row subscripts and matrix J with column subscripts

smcp.misc.sub2ind(siz, I, J)¶

Converts subscripts to indices.

Parameters:	siz (integer tuple) – matrix size I (matrix) – row subscripts J (matrix) – column subscripts
Returns:	matrix with indices

smcp.misc.sdpa_read(file_obj)¶

Reads data from sparse SDPA data file (file extension: ‘dat-s’). A description of the sparse SDPA file format can be found in the document SDPLIB/FORMAT and in the SDPA User’s Manual.

Example:

>>> f = open('qpG11.dat-s')
>>> A, b, blockstruct = smcp.misc.sdpa_read(f)
>>> f.close()

smcp.misc.sdpa_readhead(file_obj)¶

Reads header from sparse SDPA data file and returns the order \(n\), the number of constraints \(m\), and a vector with block sizes.

Example:

>>> f = open('qpG11.dat-s')
>>> n, m, blockstruct = smcp.misc.sdpa_readhead(f)
>>> f.close()

smcp.misc.sdpa_write(file_obj, A, b, blockstruct)¶

Writes SDP data to sparse SDPA file.

Example:

>>> f = open('my_data_file.dat-s','w')
>>> smcp.misc.sdpa_write(f,A,b,blockstruct)
>>> f.close()

Analysis routines¶

smcp.analysis.embed_SDP(P[, order[, cholmod]])¶

Computes chordal embedding and returns SDP object with chordal sparsity pattern.

Parameters:	P (SDP) – SDP object with problem data order (string) – ‘AMD’ (default) or ‘METIS’ cholmod (boolean) – use Cholmod to compute embedding (default is false)
Returns:	SDP object with chordal sparsity

Note that CVXOPT must be compiled and linked to METIS in order to use the METIS ordering.

The following routines require Matplotlib:

smcp.analysis.spy(P[, i[, file[, scale]]])¶

Plots aggregate sparsity pattern of SDP object P or sparsity pattern of \(A_i\).

Parameters:	P (SDP) – SDP object with problem data i (integer) – index between 0 and m file (string) – saves plot to file scale (float) – downsamples plot

smcp.analysis.clique_hist(P)¶

Plots clique histogram if P.ischordal is true, and otherwise an exception is raised.

Parameters:	P (SDP) – SDP object with problem data

smcp.analysis.nnz_hist(P)¶

Plots histogram of number of nonzeros in lower triangle of \(A_1,\ldots,A_m\).

Parameters:	P (SDP) – SDP object with problem data

Random problem generators¶

class mtxnorm_SDP(p, q, r[, density[, seed]])¶

Inherits from SDP class.

Generates random data \(F_i,G \in \mathbf{R}^{p\times q}\) for the matrix norm minimization problem

\[\begin{array}{ll} \mbox{minimize} & \lVert F(z) + G \rVert_2 \end{array}\]

with the variable \(z\in \mathbf{R}^r\) and where \(F(z) = z_1F_1 + \cdots z_r F_r\). The problem is cast as an equivalent SDP:

\[\begin{split}\begin{array}{ll} \mbox{minimize} & t \\ \mbox{subject to} & \begin{bmatrix} tI & (F(z)+G)^T \\ F(z)+G & tI \end{bmatrix} \succeq 0. \end{array}\end{split}\]

The sparsity of \(F_i\) can optionally be chosen by specifying the parameter density which must be a float between 0 and 1 (default is 1 which corresponds to dense matrices).

Example:

>>> from smcp import mtxnorm_SDP
>>> P = mtxnorm_SDP(p=200,q=10,r=200)
>>> print P
<SDP: n=210, m=201, nnz=2210> mtxnorm_p200_q10_r200
>>> sol = P.solve_feas(kktsolver='qr')

class base.band_SDP(n, m, bw[, seed])¶

Generates random SDP with band sparsity and m constraints, of order n, and with bandwidth bw (bw=0 corresponds to a diagonal, bw=1 is tridiagonal etc.). Returns SDP object. The optional parameter seed sets the random number generator seed.

Example:

>>> from smcp import band_SDP
>>> P = band_SDP(n=100,m=100,bw=2,seed=10)
>>> print P
<SDP: n=100, m=100, nnz=297> band_n100_m100_bw2
>>> X,p1sol = P.solve_phase1(kktsolver='qr')
>>> P.solve_feas(kktsolver='qr',primalstart={'x':X})
>>> print sol['primal objective'],sol['dual objective']
31.2212701455 31.2212398351

class base.rand_SDP(V, m[, density[, seed]])¶

Generates random SDP with sparsity pattern V and m constraints. Returns SDP object.

The sparsity of \(A_i\) can optionally be chosen by specifying the parameter density which must be a float between 0 and 1 (default is 1 which corresponds to dense matrices).

SDP solvers¶

The following interior-point solvers were used in our experiments:

• Method M1 (SMCP 0.3a, feasible start solver with kktsolver=’chol’)
• Method M1c (SMCP 0.3a, feasible start solver with kktsolver=’chol’ and solvers.options[‘cholmod’]=True)
• Method M2 (SMCP 0.3a, feasible start solver with kktsolver=’qr’)
• CSDP 6.0.1
• DSDP 5.8
• SDPA 7.3.1
• SDPA-C 6.2.1 (binary dist.)
• SDPT3 4.0b (64-bit Matlab)
• SeDuMi 1.2 (64-bit Matlab)

Error measures¶

We report DIMACS error measures when available. The six error measures are defined as:

Here \(\|C\|_{\max} = \max_{i,j} |C_{ij}|\), and \(C\bullet X = \mathbf{tr}(C^TX)\).

Note that \(\epsilon_2(X,y,S)=0\) and \(\epsilon_4(X,y,S)=0\) since all iterates \((X,y,S)\) satisfy \(X \in \mathbf{S}_{V,c++}^n\) and \(S \in \mathbf{S}_{V,++}^n\).

Experimental setup¶

The following experiments were conducted on a desktop computer with an Intel Core 2 Quad Q6600 CPU (2.4 GHz), 4 GB of RAM, and running Ubuntu 9.10 (64 bit).

The problem instances used in the experiments are available for download here and the SDPLIB problems are available here.

We use the least-norm solution to the set of equations \(A_i \bullet X, i=1,\ldots,m,\) as starting point when it is strictly feasible, and otherwise we solve the phase I problem

Here \(\epsilon\) is a small constant, and the constraint \(\mathbf{tr}(X) \leq M\) is added to bound the feasible set.

SDPs with band structure¶

We consider a family of SDP instances where the data matrices \(C,A_1,\ldots,A_m\) are of order \(n\) and banded with a common bandwidth \(2w+1\).

Matrix norm minimization¶

We consider the matrix norm minimization problem

where \(F(x) = x_1 F_1 + \cdots + x_r F_r\) and \(G,F_i \in \mathbf{R}^{p \times q}\) are the problem data. The problem can be formulated as an SDP:

This SDP has dimensions \(m=r+1\) and \(n = p + q\). We generate \(G\) as a dense \(p\times q\) matrix, and the matrices \(F_i\) are generated such that the number of nonzero entries in each matrix is given by \(\max (1,dpq)\) where the parameter \(d \in [0,1]\) determines sparsity. The locations of nonzero entries in \(F_i\) are random. Thus, the problem family is parameterized by the tuple \((p,q,r,d)\).

Overlapping cliqes¶

We consider a family of SDPs which have an aggregate sparsity patterns \(V\) with \(l\) cliques of order \(N\). The cliques are given by

where \(u\) (\(0 \leq u \leq N-1\)) is the overlap between neighboring cliques. The sparsity pattern is illustrated below:

Note that \(u=0\) corresponds to a block diagonal sparsity pattern and \(u=N-1\) corresponds to a band pattern.

SDPLIB problems¶

The following experiment is based on problem instances from SDPLIB.

Nonchordal sparsity patterns¶

The following problems are based on sparsity patterns from the University of Florida Sparse Matrix Collection (UFSMC). We use as problem identifier the name rsX where X is the ID number of a sparsity pattern from UFSMC. Each problem instance has \(m\) constraints and the number of nonzeros in the lower triangle of \(A_i\) is \(\max\{1,\mathrm{round}(0.005 |V|)\}\) where \(|V|\) is the number of nonzeros in the lower triangle of the aggregate sparsity pattern \(V\), and \(C\) has \(|V|\) nonzeros.