Tips for custom kkt solvers

[JTK]

I've been using a custom KKT solver to solve a conelp problem with
16000 variables, 32000 constraints. I usually get solutions in several
minutes on an aging desktop. I'm posting a sketch of some
observations which made this possible. All that follows also applies
without modification to coneqp problems. We assume the problem (as
passed to cvxopt) has no equality constraints.

Storing matrices in memory is not feasible, but for my problem,
evaluation of G and G.T requires only O(n log n) flops, so iterative
solvers are an option.

The problem we need to solve in the custom KKT solver has the form
(A B)(x)=
(a)
(C D)(y)=
(b)

Solve this using the Schur complement. First solve for x in
(A-BD^{-1}C)x = a-BD^{-1}b
(*)

With x in hand, solve for y:
y= D^{-1}(b-Cx) (**)

The matrix A-BD^{-1}C = P+G.T*W^{-T}*W^{-1}*G, is guaranteed to be
nonnegative definite. Moreover, D^{-1} is directly available from W.
In short, solving (*) is, in principle, efficiently doable using
conjugate gradient, minres or whatever.

An undocumented function (and therefore used with some risk) is
misc.scale(x,W,trans=t, inverse=i)
It seems to be an efficient method of evaluating W*x, W^{-1}*x, W^{T}
*x, W^{-T}*x. I found misc.scale was faster than any custom
evaluations I cobbled together.

If there is any interest, I'll post a toy example or two.

[Joachim Dahl]

Dear Jeffery,

please share your findings... it is great to see people making advanced
use of
CVXOPT like this.

Best regards
Joachim

JTK skrev:

[JTK]

Okay -- I'll make two separate posts -- the code below can serve as a
template, the next post contains code for the same example, but solved
using conjugate-gradient iteration on the positive definite matrix P
+G.T*W'*W*G.

--Jeff

from cvxopt import solvers, matrix, normal, setseed
from cvxopt.lapack import *
from cvxopt.blas import *
import cvxopt.misc as misc
import numpy as np

"""
Create random matrices c, B, Gq for the problem

minimize c.T*x + x.T*B*x/2

subject to nrm2(Gq*x) <= 1.0
x >= 0.0
"""

# problem size
n =25
# size of Gq
gq=15
# size of Gl
gl=n

setseed(1)
c=normal(n,1)
B=normal(n,n)
B=B.T*B

"""
linear portion of
G, b stored in Gl, bl
quad portion are in Gq, bq
"""
Gl=-matrix(np.eye(gl), tc='d')
hl= matrix(0., (gl,1))

Gq=normal(gq,n)
Gq[0,:]=0.
hq=matrix(0., (gq,1))
hq[0]=1

h=matrix([hl,hq])
G=matrix([Gl,Gq])

dims = {'l': gl, 'q': [gq], 's': []}




"""
y:= alpha*G*x + beta*y, trans=='N'
y:= alpha*G.T*x + beta*y, otherwise
"""
def Gfun(x, y, alpha = 1.0, beta = 0.0, trans = 'N'):
gemv(G,x,y,alpha=alpha, beta=beta,trans=trans)

def Fkkt(W):
"""
to solve block matrix
(A B)(x)=(a)
(C D)(y)=(b)

First solve for x:

(A-BD^{-1}C)x = a-BD^{-1}b (*)

(A-BD^{-1}C) is positive definite

Next solve for y:
Cx+Dy=b so y= D^{-1}(b-Cx) (**)
"""
def f(bx, by, bz):
n=len(bx)
m=len(bz)
a = matrix(0. ,(n,1))
b = matrix(0. ,(m,1))

"""
Get the right hand side of (*)
"""
copy(bz, b, m)
misc.scale(b, W, inverse='I', trans='T')
misc.scale(b, W, inverse='I', trans='N')
Gfun(b, a, trans='T')
axpy(bx, a)

"""
solve (*) using favorite solver
"""
"""
begin CUSTOM SOLVER
"""
H=+G
misc.scale(H, W, inverse='I', trans='T')
misc.scale(H, W, inverse='I', trans='N')
H=B+G.T*H
sysv(H,a)
copy(a,bx)
"""
end CUSTOM SOLVER
"""

"""
now solve (**)
"""
Gfun(bx,b)
scal(-1, b)
axpy(bz, b)
misc.scale(b, W, inverse='I', trans='T')
scal(-1, b)
copy(b, bz)
return f

sol1=solvers.coneqp(B, c, Gfun, h, dims, kktsolver=Fkkt)
x1=sol1['x']
z1=sol1['z']
sol2=solvers.coneqp(B, c, G, h, dims)
x2=sol2['x']
z2=sol2['z']
print ""
print "Error of custom FKKT solution: ", nrm2(x1-x2)
print "Error of custom FKKT dual: ", nrm2(z1-z2)
print "Is nrm2(Gq*x1) <= 1?"
print " nrm2(Gq*x1)=", nrm2(Gq*x1)
print "Is min(x1)>=0?"
print " min(x1)=", min(x1)

[JTK]

The code below solves the KKT matrix using only matrix-vector
products. A conjugate gradient algorithm is included. The problem
considered is the same as in the previous post.

--Jeff

from cvxopt import solvers, matrix, normal, setseed
from cvxopt.lapack import *
from cvxopt.blas import *
import cvxopt.misc as misc
import numpy as np

from math import sqrt, floor

# use for diagnostics
checkProgress=True

congGrad(A,x,b,tol)

attempts to solve A*x==b using iteration

uses x as initial starting value
places the solution of A*x==b into x

we require the syntax
A(x,y)
to perform
y:=A*x

see appendix B of
http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
"""
def conjGrad(A, x, b, tol):
n=len(b)
m=floor(sqrt(n))

#r = b - A(x)
r=matrix(0., (n,1))
d=matrix(0., (n,1))
q=matrix(0., (n, 1))
A(x,r)
scal(-1,r)
axpy(b, r)
copy(r,d)
delnew = dotu(r,r)
del0 = delnew
threshold=del0*tol**2

for i in range(n):
if(delnew<threshold):
break
A(d,q)
alpha = delnew/dotu(d, q)
#x[:]=x+alpha*d
axpy(d,x,alpha)

if( i % m or i==0):
#r = r-alpha*q
axpy(q,r, -alpha)

else:
#r = b-A(x)
A(x,r)
scal(-1, r)
axpy(b, r)

delold = delnew
delnew = dotu(r,r)
if(delnew<threshold):
#r = b-A(x)
A(x,r)
scal(-1,r)
axpy(b, r)
delnew = dotu(r,r)

beta = delnew/delold
#d = s+beta*d
scal(beta,d)
axpy(r,d)

"""
y:= alpha*G*x + beta*y, trans=='N'
y:= alpha*G.T*x + beta*y, otherwise
"""
def Gfun(x, y, alpha = 1.0, beta = 0.0, trans = 'N'):
gemv(G,x,y,alpha=alpha, beta=beta,trans=trans)

def Fkkt(W):
"""

y:=(P+G'*W'*W*G)*x
"""
_y=matrix(0., (gl+gq,1))
def D(x,y):
Gfun(x, _y)
misc.scale(_y, W, inverse='I', trans='T')
misc.scale(_y, W, inverse='I', trans='N')
Gfun(_y, y, trans='T')
y[:]=B*x+y

"""
to solve block matrix
(A B)(x)=(a)
(C D)(y)=(b)

First solve for x:
(A-BD^{-1}C)x = a-BD^{-1}b (*)

(A-BD^{-1}C) is positive definite

Next solve for y:
Cx+Dy=b so y= D^{-1}(b-Cx) (**)
"""
def f(bx, by, bz):
n=len(bx)
m=len(bz)
a = matrix(0. ,(n,1))
b = matrix(0. ,(m,1))

x0= matrix(0. ,(n,1))

"""
Get the right hand side of (*)
"""
copy(bz, b, m)
misc.scale(b, W, inverse='I', trans='T')
misc.scale(b, W, inverse='I', trans='N')
Gfun(b, a, trans='T')
axpy(bx, a)
"""
solve (*) using favorite solver
"""

"""
begin CUSTOM SOLVER

use conjugate gradient method
"""
# initialize parameters
itnmax=10
algtol=1e-6
errtol=algtol*10
err=+x0
solnerr=1.
itn=0
while solnerr>errtol and itn<itnmax:
conjGrad(D, x0, a, tol=algtol)
D(x0, err)
solnerr=nrm2(err-a)
itn+=1
"""
check the quality of solution
"""
if checkProgress:
print "Solution error: ",nrm2(err-a), "iteration: ",
itn

copy(x0, bx)