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

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 中函数值的二维数组，然后找到该数组的最大值和最小值。但这需要我知道如何精细地定义 x 和 y 网格以可靠地捕获函数的行为

对于后一种情况，我已经找到了一些方法，比如 this one .

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

Answer 1

find a list of the stationary points, of their values and locations, and of whether they are minima or maxima.

这通常是一个无法解决的问题。方法 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 没有找到所有 解(其中有无穷多个)。考虑 solve vs solveset但无论如何，处理无限多的解是很困难的。

使用 SciPy 进行数值优化

SciPy 提供了很多 numerical minimization routines , 包括 brute force (这是您的方法 2；通常它非常非常慢)。这些都是强大的方法，但请考虑以下几点。

1. 每次运行只会找到一个最小值。
2. 将 f 替换为 -f 也可以找到最大值。
3. 更改搜索的起点(minimize 的参数 x0)可能会产生另一个最大值或最小值。不过，您永远不会知道是否还有其他您还没有看到的极值。
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 只会出现一个。