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 ), 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.base import matrix, log, mul, div
from cvxopt import 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