c, kktsolver[, Gl, hl[, Gs, hs[, A, b[, primalstart[, dualstart]]]]]) |
R[k]
of order m_k, returns a function for
solving the equation
'N'
]]]) evaluates the matrix-vector products
'N'
]]])
evaluates the linear mappings
'N'
]]])
evaluates the matrix vector products
The optimization problem
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]
The SDP
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]