4.10 Example: Analytic Centering

The analytic centering problem is defined as

           ∑
minimize  -   mi=1log(bi - aTi x).

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

  T           -2                  -1
A  diag(b- Ax)  Av = - diag(b- Ax)  1

(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