¿Cuál es la forma más rápida de minimizar una función en python?

Así que tengo el siguiente problema para minimizar. Tengo un vector w que necesito encontrar para minimizar la siguiente función:

 import numpy as np from scipy.optimize import minimize matrix = np.array([[1.0, 1.5, -2.], [0.5, 3.0, 2.5], [1.0, 0.25, 0.75]]) def fct(x): return x.dot(matrix).dot(x) x0 = np.ones(3) / 3 cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1.0}) bnds = [(0, 1)] * 3 w = minimize(fct, x0, method='SLSQP', bounds=bnds, constraints=cons)['x'] 

Elegí el method='SLSQP' porque parece que es el único que permite bounds y constraints . Mi problema es que tendré que hacer un bucle de mi solución en selecciones múltiples, así que estoy tratando de ganar algo de velocidad aquí ¿Mi solución es la más rápida usando un optimizador o habría otras soluciones más rápidas? Gracias.

Introducción

En general, el enfoque más rápido siempre será el más adaptado al problema.

Como todos los algoritmos de optimización dentro de scipy.minimize son bastante generales, siempre habrá métodos más rápidos, obteniendo el rendimiento de las características especiales de su problema. Será un compromiso, cuánto análisis y trabajo se realiza para obtener un rendimiento.

Es importante tener en cuenta que, por ejemplo, el SLSQP es un algoritmo que puede abordar problemas no convexos, en cuyo caso la convergencia a algún óptimo local está garantizada (ignorando problemas numéricos dentro de la implementación; lo que siempre es un problema posible) .

Este poder tiene un precio: SLSQP será menos rápido y menos robusto en comparación con los algoritmos que están diseñados específicamente para problemas convexos (e incluso dentro de los problemas convexos, aunque todos son polinómicamente solubles, hay más fáciles como Progtwigción Lineal y más difíciles como Progtwigción Semidefinita ).

Análisis del problema

Como se indicó en los comentarios anteriores, para una matriz general indefinida M , este problema no es convexo (con una alta probabilidad; no estoy dando una prueba formal), lo que significa que no hay un enfoque general factible sin más suposiciones (ignorar análisis especial ya que algunos problemas no convexos pueden resolverse globalmente en tiempo polinomial).

Esto significa:

  • Cada algoritmo de optimización dentro de Scipy garantizará, como máximo, un óptimo local, que podría ser arbitrariamente malo en comparación con el óptimo global.

Supuesto: M es positivo-definido / negativo-definido

Si asumimos que la matriz M es positiva-definida o negativa-definida, pero no indefinida, este es un problema de optimización convexo. Como parece estar interesado en este caso, aquí hay algunas observaciones y enfoques.

Esto significa:

  • SLSQP está garantizado para converger al óptimo global
  • Puede usar solucionadores diseñados específicamente para problemas de optimización convexos
    • Solucionadores comerciales: Gurobi, CPLEX, Mosek
    • Solucionadores de código abierto: ECOS, SCS

Código de ejemplo usando Python + cvxpy + ecos / scs

No hay un solucionador de optimización convexo especial, excepto linprog, que es para la progtwigción lineal y, por lo tanto, no puede abordar este problema.

Sin embargo, hay otras alternativas, como se mencionó anteriormente, y hay muchas rutas posibles para usarlas.

Aquí presentaré uno de los más sencillos:

  • cvxpy se utiliza para la formulación modelo
    • ¡Se comprobará automáticamente que este problema es convexo!
      • (La construcción de modelos y la inferencia de convexidad pueden ser costosas)
  • ecos
    • Solucionador de propósito general para muchos problemas de optimización convexos
      • Basado en el método de punto interior
  • scs
    • Solucionador de propósito general para muchos problemas de optimización convexos
      • Basado en el método de dirección alterna de multiplicadores (ADMM)
    • Soporta dos enfoques diferentes para resolver ecuaciones lineales:
      • directo (factorización basada)
      • indirecto (basado en gradiente conjugado)
        • Soporte de GPU para este ya que se trata de productos vectoriales matriciales
    • Soluciones menos precisas en general en comparación con ECOS, pero a menudo mucho más rápidas

Código de ejemplo:

  • Ejemplo usando:
    • Matriz 1000×1000
    • Solver: SCS
      • modo indirecto
      • UPC
      • Multiproceso (automáticamente si BLAS está disponible)

Código:

 import time import numpy as np from cvxpy import * # Convex-Opt """ Create some random pos-def matrix """ N = 1000 matrix_ = np.random.normal(size=(N,N)) matrix = np.dot(matrix_, matrix_.T) """ CVXPY-based Convex-Opt """ print('\ncvxpy\n') x = Variable(N) constraints = [x >= 0, x <= 1, sum(x) == 1] objective = Minimize(quad_form(x, matrix)) problem = Problem(objective, constraints) time_start = time.perf_counter() problem.solve(solver=SCS, use_indirect=True, verbose=True) # or: solver=ECOS time_end = time.perf_counter() print(problem.value) print('cvxpy (modelling) + ecos/scs (solving) used (secs): ', time_end - time_start) 

Ejemplo de salida:

 cvxpy ---------------------------------------------------------------------------- SCS v1.2.6 - Splitting Conic Solver (c) Brendan O'Donoghue, Stanford University, 2012-2016 ---------------------------------------------------------------------------- Lin-sys: sparse-indirect, nnz in A = 1003002, CG tol ~ 1/iter^(2.00) eps = 1.00e-03, alpha = 1.50, max_iters = 2500, normalize = 1, scale = 1.00 Variables n = 1001, constraints m = 3003 Cones: primal zero / dual free vars: 1 linear vars: 2000 soc vars: 1002, soc blks: 1 Setup time: 6.76e-02s ---------------------------------------------------------------------------- Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s) ---------------------------------------------------------------------------- 0| inf inf -nan -inf -inf inf 1.32e-01 100| 1.54e-02 1.48e-04 7.63e-01 -5.31e+00 -4.28e+01 1.10e-11 1.15e+00 200| 1.53e-02 1.10e-04 7.61e-01 -3.87e+00 -3.17e+01 1.08e-11 1.95e+00 300| 1.53e-02 7.25e-05 7.55e-01 -2.47e+00 -2.08e+01 1.07e-11 2.79e+00 400| 1.53e-02 3.61e-05 7.39e-01 -1.11e+00 -1.03e+01 1.06e-11 3.61e+00 500| 7.64e-03 2.55e-04 1.09e-01 -2.01e-01 -6.32e-02 1.05e-11 4.64e+00 560| 7.71e-06 4.24e-06 8.61e-04 2.17e-01 2.16e-01 1.05e-11 5.70e+00 ---------------------------------------------------------------------------- Status: Solved Timing: Solve time: 5.70e+00s Lin-sys: avg # CG iterations: 1.71, avg solve time: 9.98e-03s Cones: avg projection time: 3.97e-06s ---------------------------------------------------------------------------- Error metrics: dist(s, K) = 5.1560e-16, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 2.4992e-17 |Ax + s - b|_2 / (1 + |b|_2) = 7.7108e-06 |A'y + c|_2 / (1 + |c|_2) = 4.2390e-06 |c'x + b'y| / (1 + |c'x| + |b'y|) = 8.6091e-04 ---------------------------------------------------------------------------- c'x = 0.2169, -b'y = 0.2157 ============================================================================ 0.21689554805292935 cvxpy (modelling) + ecos/scs (solving) used (secs): 7.105745473999832 

Ejemplo extra: 5000x5000 usa ~ 9 minutos (sin ajustes de parámetros).

Algunas pequeñas observaciones adicionales:

  • SCS es más rápido que ECOS (esperado)
  • SCS / ECOS, ambos más rápidos que los ingenuos (no dan jacobi-matrix) SLSQP (para todo al menos cada N> = 100); más rápido = órdenes de magnitud para N grande
  • Revisé la equivalencia de este método en comparación con SLSQP para pequeños ejemplos

Basándome en los comentarios de Pylang, calculé el jacobiano de mi función que lleva a la siguiente función:

 def fct_deriv(x): return 2 * matrix.dot(x) 

El problema de optimización se convierte en el siguiente

 minimize(fct, x0, method='SLSQP', jac=fct_deriv, bounds=bnds, constraints=cons)['x'] 

Sin embargo, esa solución no permite agregar el Hessian ya que el método SLSQP no lo permite. Existen otros métodos de optimización, pero SLSQP es el único que acepta límites y restricciones al mismo tiempo (lo que es fundamental para mi problema de optimización).

Vea a continuación el código completo:

 import numpy as np from scipy.optimize import minimize matrix = np.array([[1.0, 1.5, -2.], [0.5, 3.0, 2.5], [1.0, 0.25, 0.75]]) def fct(x): return x.dot(matrix).dot(x) def fct_deriv(x): return 2 * matrix.dot(x) x0 = np.ones(3) / 3 cons = ({'type': 'eq', 'fun': lambda x: x.sum() - 1.0}) bnds = [(0, 1)] * 3 w = minimize(fct, x0, method='SLSQP', jac=fct_deriv, bounds=bnds, constraints=cons)['x'] 

Editado (agregado el jacobiano de la restricción):

 cons2 = ({'type': 'eq', 'fun': lambda x: x.sum() - 1.0, 'jac': lambda x: np.ones_like(x)}) w = minimize(fct, x0, method='SLSQP', jac=fct_deriv, bounds=bnds, constraints=cons2)['x']