8.6 Exploiting Structure in LPs and SDPs

lp()

sdp()

conelp()

conelp( c, kktsolver[, Gl, hl[, Gs, hs[, A, b[, primalstart[, dualstart]]]]])

Solves the pair of primal and dual SDPs (8.2). The arguments c, hl, hs, b, primalstart, dualstart have the same meaning as in sdp(). The arguments kktsolver, Gl, Gs, A are functions that must handle the following calling sequences.

kktsolver(d, R) with d a positive dense real matrix of size (ml,1), and R a list of N square dense real matrices R[k] of order m_k, returns a function for solving the equation

$\begin{eqnarray*} A^T u_y + G_\mathrm{l}^T u_{z_\mathrm{l}} + G_\mathrm{s}^T(... ...} R^{-1} u_{z_\mathrm{s}} R^{-T} R^{-1} & = & b_{z_\mathrm{s}}. \end{eqnarray*}$

The function created by "f = kktsolver(d, R)" will be called as "f(bx, by, bzl, bzs)". On entry, bx, by, bzl and bzs contain the righthand side. On exit, they should contain the solution of the KKT system, with uzl and uzs scaled:

$\begin{displaymath} b_x := u_x, \qquad b_y := u_y, \qquad b_{z_\mathrm{l}} :=... ...}, \qquad b_{z_\mathrm{s}} := R^{-1} u_{z_\mathrm{s}} R^{-T}. \end{displaymath}$
Gl(x, y[, alpha=1.0[, beta=0.0[, trans='N']]]) evaluates the matrix-vector products

$\begin{displaymath} y := \alpha G_\mathrm{l}x + \beta y \quad (\mathrm{trans} ... ...thrm{l}^T x + \beta y \quad (\mathrm{trans} = \mathrm{'T'}). \end{displaymath}$
Gs(x, y[, alpha=1.0[, beta=0.0[, trans='N']]]) evaluates the linear mappings

$\begin{displaymath} y := \alpha G_\mathrm{s}(x) + \beta y \quad (\mathrm{trans... ...hrm{s}^T(x) + \beta y \quad (\mathrm{trans} = \mathrm{'T'}). \end{displaymath}$
A(x, y[, alpha=1.0[, beta=0.0[, trans='N']]]) evaluates the matrix vector products

$\begin{displaymath} y := \alpha Ax + \beta y \quad (\mathrm{trans} = \mathrm{'... ...\alpha A^Tx + \beta y \quad (\mathrm{trans} = \mathrm{'T'}). \end{displaymath}$

Example: 1-norm approximation

The optimization problem

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & \Vert Pu-q\Vert _1 \end{array}\end{displaymath}$

can be formulated as an LP

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & {\bf 1}^T v \\ \mbox{subject to} & -v \preceq Pu - q \preceq v. \end{array}\end{displaymath}$

By exploiting the structure in the inequalities, the cost of an iteration of an interior-point method can be reduced to the cost of least-squares problem of the same dimensions. (See section 11.8.2 in the book Convex Optimization.) The code belows taks advantage of this fact.

from cvxopt import base, blas, lapack, solvers
from cvxopt.base import matrix, spmatrix, mul, div

def l1(P, q):
    """

    Returns the solution u, w of the l1 approximation problem

        (primal) minimize    ||P*u - q||_1       
    
        (dual)   maximize    q'*w
                 subject to  P'*w = 0
                             ||w||_infty <= 1.
    """

    m, n = P.size

    # Solve equivalent LP 
    #
    #     minimize    [0; 1]' * [u; v]
    #     subject to  [P, -I; -P, -I] * [u; v] <= [q; -q]
    #
    #     maximize    -[q; -q]' * z 
    #     subject to  [P', -P']*z  = 0
    #                 [-I, -I]*z + 1 = 0 
    #                 z >= 0 
    
    c = matrix(n*[0.0] + m*[1.0])
    h = matrix([q, -q])

    def Fi(x, y, alpha=1.0, beta=0.0, trans='N'):    
        if trans=='N':
            # y := alpha * [P, -I; -P, -I] * x + beta*y
            u = P*x[:n]
            y[:m] = alpha * ( u - x[n:]) + beta*y[:m]
            y[m:] = alpha * (-u - x[n:]) + beta*y[m:]

        else:
            # y := alpha * [P', -P'; -I, -I] * x + beta*y
            y[:n] =  alpha * P.T * (x[:m] - x[m:]) + beta*y[:n]
            y[n:] = -alpha * (x[:m] + x[m:]) + beta*y[n:]


    def kktsolver(d, R): 

        # Returns a function f(x,y,zl,zs) that solves
        #
        # [ 0  0  P'      -P'      ] [ x[:n] ]   [ bx[:n]  ]
        # [ 0  0 -I       -I       ] [ x[n:] ]   [ bx[n:]  ]
        # [ P -I -D1^{-1}  0       ] [ zl[:m]] = [ bzl[:m] ]
        # [-P -I  0       -D2^{-1} ] [ zl[m:]]   [ bzl[m:] ]
        #
        # where D1 = diag(d[:m])^2, D2 = diag(d[m:])^2.
        #
        # On entry bx, bzl are stored in x, zl.
        # On exit x, zl contain the solution, with zl scaled: zl./d is
        # returned instead of zl. 

        # Factor A = 4*P'*D*P where D = d1.*d2 ./(d1+d2) and
        # d1 = d[:m].^2, d2 = d[m:].^2.

        d1, d2 = d[:m]**2, d[m:]**2
        D = div( mul(d1,d2), d1+d2 )  
        A = P.T * spmatrix(4*D, range(m), range(m)) * P
        lapack.potrf(A)

        def f(x, y, zl, zs):

            # Solve for x[:n]:
            #
            #    A*x[:n] = bx[:n] + P' * ( ((D1-D2)*(D1+D2)^{-1})*bx[n:]
            #              + (2*D1*D2*(D1+D2)^{-1}) * (bzl[:m] - bzl[m:]) ).
            x[:n] += P.T * ( mul(div(d1-d2, d1+d2), x[n:]) + mul(2*D, zl[:m]-zl[m:]) )
            lapack.potrs(A, x)

            # x[n:] := (D1+D2)^{-1} * (bx[n:] - D1*bzl[:m] - D2*bzl[m:] + (D1-D2)*P*x[:n])
            u = P*x[:n]
            x[n:] =  div(x[n:] - mul(d1, zl[:m]) - mul(d2, zl[m:]) + mul(d1-d2, u), d1+d2)

            # z[:m] := d1[:m] .* ( P*x[:n] - x[n:] - bzl[:m])
            # z[m:] := d2[m:] .* (-P*x[:n] - x[n:] - bzl[m:]) 
            zl[:m] = mul(d[:m],  u-x[n:]-zl[:m])
            zl[m:] = mul(d[m:], -u-x[n:]-zl[m:])

        return f

    sol = solvers.conelp(c, kktsolver, Gl=Fi, hl=h) 
    return sol['x'][:n],  sol['zl'][m:] - sol['zl'][:m]

Example: SDP with diagonal linear term

The SDP

$\begin{displaymath} \begin{array}{ll} \mbox{minimize} & {\bf 1}^T x \\ \mbox{subject to} & W + \mbox{\bf diag}\,(x) \succeq 0 \end{array} \end{displaymath}$

can be solved efficiently by exploiting properties of the diag operator.

from cvxopt import base, blas, lapack, solvers
from cvxopt.base import matrix

def mcsdp(w):
    """
    Returns solution x, z to 

        (primal)  minimize    sum(x)
                  subject to  w + diag(x) >= 0

        (dual)    maximize    -tr(w*z)
                  subject to  diag(z) = 1
                              z >= 0.
    """

    n = w.size[0]

    def Fs(x, y, alpha=1.0, beta=0.0, trans='N'):
        """
            y := alpha*(-diag(x)) + beta*y.   
        """
        if trans=='N':
            # x is a vector; y[0] is a matrix.
            y[0] *= beta
            y[0][::n+1] -= alpha * x
        else:   
            # x[0] is a matrix; y is a vector.
            y *= beta
            y -= alpha * x[::n+1] 
	 

    def cngrnc(r, x, alpha=1.0):
        """
        Congruence transformation

	    x := alpha * r'*x*r.

        r and x are square 'd' matrices.  
        """

        # Scale diagonal of x by 1/2.  
        x[::n+1] *= 0.5
    
        # a := tril(x)*r 
        a = +r
        blas.trmm(x, a, side='L')

        # x := alpha*(a*r' + r*a') 
        blas.syr2k(r, a, x, trans='T', alpha=alpha)


    def kktsolver(d, r):

        # t = r*r' as a nonsymmetric matrix.
        t = matrix(0.0, (n,n))
        blas.gemm(r[0], r[0], t, transB='T') 

        # Cholesky factorization of tsq = t.*t.
        tsq = t**2
	lapack.potrf(tsq)

	def f(x, y, zl, zs):
            """
            Solve
                          -diag(zs)               = bx
                -diag(x) - inv(r*r')*zs*inv(r*r') = bs.

            On entry, x and zs contain bx and bs.  
            On exit, they contain the solution, with zs scaled 
            (inv(r)'*zs*inv(r) is returned instead of zs).

            We solve 

                ((r*r') .* (r*r')) * x = bx - diag(t*bs*t)

            and take zs  = -r' * (diag(x) + bs) * r.
            """

            # tbst := t * zs * t = t * bs * t
            tbst = +zs[0]
            cngrnc(t, tbst) 

            # x := x - diag(tbst) = bx - diag(r*r' * bs * r*r')
            x -= tbst[::n+1]

            # x := (t.*t)^{-1} * x = (t.*t)^{-1} * (bx - diag(t*bs*t))
            lapack.potrs(tsq, x)

            # zs := zs + diag(x) = bs + diag(x)
            zs[0][::n+1] += x 

            # zs := -r' * zs * r = -r' * (diag(x) + bs) * r 
            cngrnc(r[0], zs[0], alpha=-1.0)

	return f

    c = matrix(1.0, (n,1))
    sol = solvers.conelp(c, kktsolver, Gs=Fs, hs=[w]) 
    return sol['x'], sol['zs'][0]