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 |