Optimisation rapide de la fonction convexe "pathologique"

Question

J'ai un simple problème convexe J'essaie d'accélérer la solution de. Je résous le argmin (thêta) de

où thêta et rt est Nx1.

Je peux résoudre ce facilement avec 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]]) 

Mais pour une plus grande R cela devient trop lent, donc je suis en train une méthode basée sur gradient avec l » fmin_cgscipy:

goalfun est un fonction conviviale qui renvoie la valeur de la fonction et le gradient.

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 

assuré que la fonction et les gradients sont corrects:

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

Mais résoudre cela donne juste des ordures, quel que soit method, itérations, etc. (Les seules choses que les rendements Optimization terminated successfully est si x0 est pratiquement égale à la optimale thêta)

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

Ie Ce problème apparemment inoffensif, que l'on peut facilement traiter, s'avère complètement pathologique pour un solveur non-convexe. Ce problème est vraiment méchant, ou est-ce que je manque quelque chose? Quelle serait une alternative pour accélérer cela?

Answer 1

Je crois que le problème est qu'il est possible que theta soit tel que l'argument log devienne négatif. Il semble que vous ayez identifié ce problème et que goalfun ait retourné le n-uplet (100,100*ones(N)) dans ce cas, apparemment, comme une tentative heuristique de suggérer au solveur que cette "solution" n'est pas préférable. Cependant, une condition plus forte doit être imposée, c'est-à-dire que cette "solution" n'est pas possible . Bien sûr, cela peut être fait en fournissant des contraintes appropriées. (Fait intéressant, cvxpy semble gérer ce problème automatiquement.)

Voici un exemple d'exécution, sans se soucier de fournir des dérivés. Notez l'utilisation d'une estimation initiale possible 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]) 

Notez que quand je cours, je reçois un avertissement invalid value encountered in log, ce qui indique que, à un moment donné dans la recherche d'une valeur de theta est vérifiée qui satisfait à peine les contraintes. Cependant, le résultat est raisonnablement proche de cvxpy. Il serait intéressant de vérifier si la solution cvxpy change lorsque les contraintes sont explicitement imposées dans la formulation cvxpy.Problem.