寻找用户定义函数的局部最大值和最小值

Question

寻找用户定义函数的局部最大值和最小值

5

我的需求

我想要找到一个静止点列表，包括它们的值和位置，以及它们是极小值还是极大值。

我的函数长这样：

import numpy as np

def func(x,y):
  return (np.cos(x*10))**2 + (np.sin(y*10))**2

方法

以下是我考虑使用的方法：

在Mathematica上，我已经在做类似这样的事情。我对函数进行一次和两次微分。我查看一阶导数为0的点，计算它们的值和位置。然后，在那些位置处取二阶导数，检查它们是否是极小值或极大值。
我也想知道是否只需创建一个x和y中函数值的二维数组，并找到该数组的最大和最小值。但这需要我知道如何定义x和y的网格以可靠地捕获函数的行为。

对于后一种情况，我已经找到了一些方法，例如这个。

我只是想知道，在Python中哪种方法更具效率、速度、准确性甚至优雅？

- SuperCiocia

1

您遗漏了数值优化方法，这是介于这两种方法之间的一种方法。 - Davis Herring

在第二种方法中，如果我理解正确的话，您会将单调递增或递减函数的最后一个值解释为最大/最小值，这似乎不是您想要的。如果用户输入的函数具有这种单调区域，它也会被错误解释。对我来说，方法1或优化方法听起来最好。Python已经拥有了丰富的工具来实现这两种方法。 - atru

已经有一段时间了，但我认为你可以查看“最陡下降法”（steepest decent method），它相当简单易实现。此外，有很多Python程序可以实现它。 - atru

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

寻找静止点的列表，它们的值和位置，以及它们是极小值还是极大值。这通常是一个无法解决的问题。方法1（符号）适用于此，但对于复杂的函数，静止点没有符号解（一般情况下没有方法可以符号地解决两个方程组的问题）。对于像您的示例这样的简单函数，SymPy将很好地工作。以下是一个完整的示例，介绍如何通过Hessian矩阵的特征值来找到静止点并将其分类。

import sympy as sym
x, y = sym.symbols("x y")
f = sym.cos(x*10)**2 + sym.sin(y*10)**2
gradient = sym.derive_by_array(f, (x, y))
hessian = sym.Matrix(2, 2, sym.derive_by_array(gradient, (x, y)))

到目前为止，Hessian是一个2x2的符号矩阵：[[200*sin(10*x)**2 - 200*cos(10*x)**2, 0], [0, -200*sin(10*y)**2 + 200*cos(10*y)**2]]。接下来，我们通过将gradient等于零来找到静止点，并逐个将它们代入到Hessian中。

stationary_points = sym.solve(gradient, (x, y))
for p in stationary_points:
    value = f.subs({x: p[0], y: p[1]})
    hess = hessian.subs({x: p[0], y: p[1]})
    eigenvals = hess.eigenvals()
    if all(ev > 0 for ev in eigenvals):
        print("Local minimum at {} with value {}".format(p, value))
    elif all(ev < 0 for ev in eigenvals):
        print("Local maximum at {} with value {}".format(p, value))
    elif any(ev > 0 for ev in eigenvals) and any(ev < 0 for ev in eigenvals):
        print("Saddle point at {} with value {}".format(p, value))
    else:
        print("Could not classify the stationary point at {} with value {}".format(p, value))

最后一个从句是必要的，因为当Hessian矩阵只是半正定时，我们无法确定这是什么类型的稳定点（在(0, 0)处，x**2 + y**4和x**2 - y**4具有相同的Hessian矩阵，但行为不同）。输出：

Saddle point at (0, 0) with value 1
Local maximum at (0, pi/20) with value 2
Saddle point at (0, pi/10) with value 1
Local maximum at (0, 3*pi/20) with value 2
Local minimum at (pi/20, 0) with value 0
Saddle point at (pi/20, pi/20) with value 1
Local minimum at (pi/20, pi/10) with value 0
Saddle point at (pi/20, 3*pi/20) with value 1
Saddle point at (pi/10, 0) with value 1
Local maximum at (pi/10, pi/20) with value 2
Saddle point at (pi/10, pi/10) with value 1
Local maximum at (pi/10, 3*pi/20) with value 2
Local minimum at (3*pi/20, 0) with value 0
Saddle point at (3*pi/20, pi/20) with value 1
Local minimum at (3*pi/20, pi/10) with value 0
Saddle point at (3*pi/20, 3*pi/20) with value 1

显然，solve 没有找到所有的解（因为这些解是无限多的）。考虑使用 solve vs solveset，但在任何情况下，处理无限多的解都很困难。

使用 SciPy 进行数值优化

SciPy 提供了许多数值最小化例程，包括暴力搜索（也就是方法 2；通常速度非常慢）。这些是强大的方法，但需要考虑以下几点。

每次运行只能找到一个最小值。
用 -f 替换 f，你也可以找到一个最大值。
改变搜索的起始点（minimize 的参数 x0）可能会得到另一个最大值或最小值。但你永远不会知道是否还有其他未被发现的极值。
这些方法都无法找到鞍点。

混合策略

使用 lambdify 可以将符号表达式转换为 Python 函数，该函数可以传递给 SciPy 数值求解器。

from scipy.optimize import fsolve
grad = sym.lambdify((x, y), gradient)
fsolve(lambda v: grad(v[0], v[1]), (1, 2))

这将返回一些静止点，例如在这个例子中是[0.9424778, 2.04203522]。它是哪一个取决于初始猜测，这里是(1, 2)。通常（但并非总是），你会得到一个接近初始猜测的解。相比直接最小化方法，这种方法的优势在于可以检测到鞍点。然而，要找到所有的解仍然很困难，因为每次运行fsolve只能得到一个解。