Fast optimization of "pathological" convex function
<p>I have a simple convex problem I am trying to speed up the solution of. I am solving the argmin (<em>theta</em>) of</p>

<p><img src="http://latex.codecogs.com/gif.latex?-%5Csum_%7Bt%3D1%7D%5E%7BT%7D%5Clog%281%20&plus;%20%5Cboldsymbol%7B%5Ctheta%27%7D%20%5Cboldsymbol%7Br_t%7D%29" alt="eq"></p>

<p>where <em>theta</em> and <em>rt</em> is <em>Nx1</em>.</p>

<p>I can solve this easily with <code>cvxpy</code>
</p>

<pre><code>import numpy as np
from scipy.optimize import minimize
import cvxpy

np.random.seed(123)

T = 50
N = 5
R = np.random.uniform(-1, 1, size=(T, N))

cvtheta = cvxpy.Variable(N)
fn = -sum([cvxpy.log(1 + cvtheta.T * rt) for rt in R])

prob = cvxpy.Problem(cvxpy.Minimize(fn))
prob.solve()

prob.status
#'optimal'

prob.value
# -5.658335088091929

cvtheta.value
# matrix([[-0.82105079],
#         [-0.35475695],
#         [-0.41984643],
#         [ 0.66117397],
#         [ 0.46065358]])
</code></pre>

<p>But for a larger <code>R</code> this gets too slow, so I am trying a gradient based method with <code>scipy</code>'s <code>fmin_cg</code>:</p>

<p><code>goalfun</code> is a <code>scipy.minimize</code> friendly function that returns the function value and the gradient.</p>

<pre><code>def goalfun(theta, *args):
    R = args[0]
    N = R.shape[1]
    common = (1 + np.sum(theta * R, axis=1))**-1

    if np.any( common &lt; 0 ):
        return 1e2, 1e2 * np.ones(N)

    fun = np.sum(np.log(common))

    thetaprime = np.tile(theta, (N, 1)).T
    np.fill_diagonal(thetaprime, np.ones(N))
    grad = np.sum(np.dot(R, thetaprime) * common[:, None], axis=0)

    return fun, grad
</code></pre>

<p>Making sure the function and gradients are correct:</p>

<pre><code>goalfun(np.squeeze(np.asarray(cvtheta.value)), R)
# (-5.6583350819293603,
#  array([ -9.12423065e-09,  -3.36854633e-09,  -1.00983679e-08,
#          -1.49619901e-08,  -1.22987872e-08]))
</code></pre>

<p>But solving this just yields garbage, regardless of <code>method</code>, iterations, etc. (The only things that yields <code>Optimization terminated successfully</code>  is if <code>x0</code> is practically equal to the optimal <em>theta</em>)</p>

<pre><code>x0 = np.random.rand(R.shape[1])

minimize(fun=goalfun, x0=x0, args=R, jac=True, method='CG')
#   fun: 3.3690101669818775
#      jac: array([-11.07449021, -14.04017873, -13.38560561,  -5.60375334,  -2.89210078])
#  message: 'Desired error not necessarily achieved due to precision loss.'
#     nfev: 25
#      nit: 1
#     njev: 13
#   status: 2
#  success: False
#        x: array([ 0.00892177,  0.24404118,  0.51627475,  0.21119326, -0.00831957])
</code></pre>

<p>I.e. this seemingly innocuous problem that <code>cvxpy</code> handles with ease, turns out to be completely pathological for a non-convex solver. Is this problem really that nasty, or am I missing something? What would be an alternative to speed this up?</p>

<p>I believe the issue is that it is possible for <code>theta</code> to be such that the <code>log</code> argument becomes negative. It seems that you have identified this issue, and have <code>goalfun</code> return the tuple <code>(100,100*ones(N))</code> in this case, apparently, as a heuristic attempt to suggest the solver that this "solution" is not <em>preferable</em>. However, a stronger condition must be imposed, i.e., this "solution" is not <em>feasible</em>. Of course, this can be done by providing appropriate constraints. (Interestingly, <code>cvxpy</code> appears to handle this issue automatically.)  </p>

<p>Here is a sample run, without bothering with providing derivatives. Note the use of a feasible initial estimate <code>x0</code>.</p>

<pre><code>np.random.seed(123)

T = 50
N = 5
R = np.random.uniform(-1, 1, size=(T, N))

def goalfun(theta, *args):
    R = args[0]
    N = R.shape[1]
    common = (1 + np.sum(theta * R, axis=1))**-1

    return np.sum(np.log(common))

def con_fun(theta, *args):
    R = args[0]

    return 1+np.sum(theta * R, axis=1)


cons = ({'type': 'ineq', 'fun': lambda x: con_fun(x, R)})

x0 = np.zeros(R.shape[1])
minimize(fun=goalfun, x0=x0, args=R, constraints=cons)
</code></pre>

<blockquote>
<pre><code> fun: -5.658334806882614
 jac: array([ 0.0019, -0.0004, -0.0003,  0.0005, -0.0015,  0.    ])  message: 'Optimization terminated successfully.'
nfev: 92
 nit: 12
njev: 12   status: 0  success: True
   x: array([-0.8209, -0.3547, -0.4198,  0.6612,  0.4605])
</code></pre>
</blockquote>

<p>Note that when I run this, I get an <code>invalid value encountered in log</code> warning, indicating that at some point in the search a value of <code>theta</code> is checked which barely satisfies the constraints. However, the result is reasonably close to that of <code>cvxpy</code>. It would be interesting to check if the <code>cvxpy</code> solution changes when the constraints are explicitly imposed in the <code>cvxpy.Problem</code> formulation.</p>

python numpy scipy mathematical-optimization cvxopt