#!/usr/bin/python # The SDP example of section 8.7 (Exploiting structure). from cvxopt import blas, lapack, solvers, matrix, normal 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 is a matrix. blas.scal(beta, y) blas.axpy(x, y, alpha = -alpha, incy = n+1) else: # x is a matrix; y is a vector. blas.scal(beta, y) blas.axpy(x, y, alpha = -alpha, incx = n+1) def cngrnc(r, x, alpha = 1.0): """ Congruence transformation x := alpha * r'*x*r. r and x are square matrices. """ # Scale diagonal of x by 1/2. x[::n+1] *= 0.5 # a := tril(x)*r a = +r tx = matrix(x, (n,n)) blas.trmm(tx, a, side='L') # x := alpha*(a*r' + r*a') blas.syr2k(r, a, tx, trans = 'T', alpha = alpha) x[:] = tx[:] def Fkkt(W): rti = W['rti'][0] # t = rti*rti' as a nonsymmetric matrix. t = matrix(0.0, (n,n)) blas.gemm(rti, rti, t, transB = 'T') # Cholesky factorization of tsq = t.*t. tsq = t**2 lapack.potrf(tsq) def f(x, y, z): """ Solve -diag(z) = bx -diag(x) - inv(rti*rti') * z * inv(rti*rti') = bs On entry, x and z contain bx and bs. On exit, they contain the solution, with z scaled (inv(rti)'*z*inv(rti) is returned instead of z). We first solve ((rti*rti') .* (rti*rti')) * x = bx - diag(t*bs*t) and take z = -rti' * (diag(x) + bs) * rti. """ # tbst := t * zs * t = t * bs * t tbst = matrix(z, (n,n)) cngrnc(t, tbst) # x := x - diag(tbst) = bx - diag(rti*rti' * bs * rti*rti') x -= tbst[::n+1] # x := (t.*t)^{-1} * x = (t.*t)^{-1} * (bx - diag(t*bs*t)) lapack.potrs(tsq, x) # z := z + diag(x) = bs + diag(x) z[::n+1] += x # z := -rti' * z * rti = -rti' * (diag(x) + bs) * rti cngrnc(rti, z, alpha = -1.0) return f c = matrix(1.0, (n,1)) # Initial feasible x: x = 1.0 - min(lambda(w)). lmbda = matrix(0.0, (n,1)) lapack.syevx(+w, lmbda, range='I', il=1, iu=1) x0 = matrix(-lmbda[0]+1.0, (n,1)) s0 = +w s0[::n+1] += x0 # Initial feasible z is identity. z0 = matrix(0.0, (n,n)) z0[::n+1] = 1.0 dims = {'l': 0, 'q': [], 's': [n]} sol = solvers.conelp(c, Fs, w[:], dims, kktsolver = Fkkt, primalstart = {'x': x0, 's': s0[:]}, dualstart = {'z': z0[:]}) return sol['x'], matrix(sol['z'], (n,n)) n = 100 w = normal(n,n) x, z = mcsdp(w)