在python中,最小化函数的最快方法是什么?

2024-05-14 00:14:23 发布

您现在位置:Python中文网/ 问答频道 /正文

所以我有以下问题要最小化。我需要找到一个向量w,以便最小化以下函数:

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

我选择method='SLSQP'是因为它似乎是唯一允许boundsconstraints的。我的问题是我将不得不循环我的多个选择的解决方案,所以我试图在这里获得一些速度。我的解决方案是使用优化器的最快的,还是有其他更快的解决方案?谢谢您。


Tags: importnp解决方案向量matrixmethoddotconstraints
2条回答

基于pylang注释,我计算了函数的jacobian,得到了以下函数:

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

优化问题如下

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

但是,该解决方案不允许添加Hessian,因为SLSQP方法不允许添加Hessian。其他优化方法也存在,但SLSQP是唯一同时接受边界和约束的方法(这是我优化问题的核心)。

完整代码见下文:

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

已编辑(添加了约束的雅可比):

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

导言

一般来说,最快的方法总是最适合问题的。

由于scipy.minimize中的所有优化算法都非常通用,因此 始终采用更快的方法,从问题的特殊特性中获得性能。 这将是一个权衡,需要做多少分析和工作来获得绩效。

需要注意的是,例如SLSQP是一个算法,它能够 处理非凸问题,在这种情况下,保证收敛到某些局部最优 (忽略实现中的数字问题;这总是一个可能的问题)。

这种能力是有代价的:与算法相比,SLSQP的速度和健壮性都会降低 专门为凸问题(甚至在凸问题中,尽管 它们都是多项式可解的,有简单的线性规划,也有困难的 作为半定义编程)。

问题分析

如上所述,对于一些一般的不定矩阵M,这个问题 是非凸的(概率很高;我没有给出形式证明),这意味着, 如果没有进一步的假设(忽略 特殊分析是一些非凸问题可以在多项式时间内全局求解)。

这意味着:

  • scipy中的每一个优化算法,最多只能保证一个局部最优,这可能 与全局最优值相比任意糟糕

假设:M为正定/负定

如果我们假设矩阵M是正定或负定的,但不是不定的,这是一个凸优化问题。由于你似乎对这个案例感兴趣,这里有一些评论和方法。

这意味着:

  • SLSQP保证收敛到全局最优
  • 可以使用专门为凸优化问题设计的解算器
    • 商业偿付机构:Gurobi、CPLEX、Mosek
    • 开源解决方案:ECOS、SCS

使用Python+cvxpy+ecos/scs的示例代码

除了用于线性规划的linprog,没有特殊的凸优化求解器 因此无法解决这个问题。

不过,如上所述,还有其他选择,有许多可能的路线 使用它们。

在这里,我将介绍一个最简单的:

  • cvxpy用于模型公式化
    • 它会自动证明这个问题是凸的!
      • (模型建立和凸性推理可能代价高昂)
  • ecos
    • 多凸优化问题的通用解
      • 基于内点法
  • scs
    • 多凸优化问题的通用解
      • 基于乘法器交替方向法(ADMM)
    • 支持求解线性方程的两种不同方法:
      • 直接(基于因子分解)
      • 间接(基于共轭梯度)
        • GPU支持这个,因为它是关于矩阵向量积的
    • 一般来说,与ECOS相比,它的精度较低,但通常要快得多

示例代码:

  • 示例使用:
    • 1000x1000矩阵
    • 解算器:SCS
      • 间接模式
      • 中央处理器
      • 多线程(如果BLAS可用,则自动)

代码:

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)

示例输出:

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

额外示例:5000x5000使用~9分钟(不调整参数)。

一些小小的额外评论:

  • SCS比ECOS快(预期)
  • SCS/ECOS都比naive快(不给出jacobi矩阵)SLSQP(至少每N>;=100);更快=大N的数量级
  • 我检查了这个方法与SLSQP在小例子中的等价性

相关问题 更多 >