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

Question

我想要的是

我想找到固定点的列表，它们的值和位置以及它们是最小值还是最大值。

我的功能看起来像：

import numpy as np

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

方法

这是我正在考虑使用的方法：

1. 我实际上已经在Mathematica上做了类似的事情。我一次又一次地区分函数。我看一下一阶导数为0的点，计算它们的值和位置。然后，我在那些位置取二阶导数，并检查它们是最小值还是最大值。

2. 我还想知道是否仅对x和y中的函数值进行2D数组化，并找到该数组的最大值和最小值。但这需要我知道如何精确地定义x和y网格，以可靠地捕获函数的行为

对于后一种情况下，我已经找到像某些方面这一个。

我只是想知道，哪种方法在效率，速度，准确性甚至Python的优雅方面更有意义？

Answer 1

找到驻点列表、它们的值和位置，以及它们是最小值还是最大值。

这通常是一个无法解决的问题。方法 1（符号）适合于此，但对于复杂函数，驻点没有符号解（没有符号求解一般两个方程组的方法）。

使用 SymPy 的符号解决方案

对于像您的示例这样的简单函数，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 是一个 2 × 2 的符号矩阵：[[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 矩阵只是半定时，我们无法判断那是哪种驻点（x**2 + y**4并且x**2 - y**4在 (0, 0) 处具有相同的 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没有找到所有的解决方案（有无限多个）。考虑求解与求解集，但无论如何，处理无限多个解都是困难的。

使用 SciPy 进行数值优化

SciPy 提供了很多数值最小化例程，包括蛮力（这是你的方法 2；通常它非常非常慢）。这些都是强大的方法，但请考虑以下几点。

1. 每次运行只会找到一个最小值。
2. 将 f 替换为 -f 您还可以找到最大值。
3. 更改搜索的起点（ 的参数 x0 minimize）可能会产生另一个最大值或最小值。尽管如此，您永远不会知道是否还有其他您尚未见过的极值。
4. 这些都找不到鞍点。

混合策略

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

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只能找到一个解决方案。