Быстрая оптимизация "патологической" выпуклой функции

Question

У меня есть простая выпуклая проблема, которую я пытаюсь ускорить решение. Я решаю аргмин (тета)

где тета и rt - это Nx1.

Я могу решить это легко с cvxpy

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]])

Но для большего R это становится слишком медленным, поэтому я пытаюсь метод на основе градиента с scipy"s fmin_cg:

goalfun это scipy.minimize дружественная функция, которая возвращает значение функции и градиент.

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

    if np.any( common < 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

Убедиться, что функция и градиенты верны:

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]))

Но решение этой проблемы просто приводит к мусору, независимо от method, итерации и т. д. (единственное, что дает Optimization terminated successfully если x0 практически равен оптимальной тэте)

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])

Т.е. эта, казалось бы, безобидная проблема cvxpy обрабатывает с легкостью, оказывается совершенно патологическим для невыпуклого решателя. Эта проблема действительно такая неприятная, или я что-то упустил? Что может быть альтернативой, чтобы ускорить это?

Answer 1

Я считаю, что проблема в том, что это возможно для theta быть таким, чтобы log аргумент становится отрицательным. Похоже, что вы определили эту проблему, и есть goalfun вернуть кортеж (100,100*ones(N)) в данном случае, по-видимому, в качестве эвристической попытки предложить решателю, что это "решение" не является предпочтительным. Однако должно быть наложено более сильное условие, т. Е. Это "решение" неосуществимо. Конечно, это может быть сделано путем предоставления соответствующих ограничений. (Что интересно, cvxpy кажется, чтобы решить эту проблему автоматически.)

Вот примерный прогон, не заботясь о предоставлении производных. Обратите внимание на использование возможной начальной оценки x0,

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)

 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])

Обратите внимание, что когда я запускаю это, я получаю invalid value encountered in log предупреждение, указывающее, что в какой-то момент в поиске значение theta проверяется, что едва удовлетворяет ограничениям. Тем не менее, результат достаточно близок к cvxpy, Было бы интересно проверить, если cvxpy решение меняется, когда ограничения явно наложены в cvxpy.Problem рецептура.