4.10 Example: Analytic Centering
The analytic centering problem is defined as
In the code below we solve the problem using Newton’s method. At each iteration
the Newton direction is computed by solving a positive definite set of linear
equations
(where A has rows aiT), and a suitable step size is determined by a backtracking line
search.
We use the level-3 BLAS function syrk() to form the Hessian matrix and the
LAPACK function posv() to solving the Newton system. The code can be further
optimized by replacing the matrix-vector products with the level-2 BLAS function
gemv().
from cvxopt import matrix, log, mul, div, blas, lapack, random
from math import sqrt
def acent(A,b):
"""
Returns the analytic center of A*x <= b.
We assume that b > 0 and the feasible set is bounded.
"""
MAXITERS = 100
ALPHA = 0.01
BETA = 0.5
TOL = 1e-8
m, n = A.size
x = matrix(0.0, (n,1))
H = matrix(0.0, (n,n))
for iter in xrange(MAXITERS):
# Gradient is g = A^T * (1./(b-A*x)).
d = (b - A*x)**-1
g = A.T * d
# Hessian is H = A^T * diag(d)^2 * A.
Asc = mul( d[:,n*[0]], A )
blas.syrk(Asc, H, trans=’T’)
# Newton step is v = -H^-1 * g.
v = -g
lapack.posv(H, v)
# Terminate if Newton decrement is less than TOL.
lam = blas.dot(g, v)
if sqrt(-lam) < TOL: return x
# Backtracking line search.
y = mul(A*v, d)
step = 1.0
while 1-step*max(y) < 0: step *= BETA
while True:
if -sum(log(1-step*y)) < ALPHA*step*lam: break
step *= BETA
x += step*v