8.6 Exploiting Structure in LPs and SDPs

The solvers lp() and sdp() are interfaces to a common function conelp(), which can also be called directly. When calling conelp(), the user must provide functions for evaluating the constraint functions and for solving the linear equations (KKT equations) that are solved in each iteration of the algorithm. This is useful for LPs and SDPs that possess some interesting structure that makes it possible to solve the KKT equations fast.

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.

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]