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

11 浏览
0 Comments

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

我的需求

\n我想要找到静止点的列表，包括它们的值和位置，以及它们是极小值还是极大值。\n我的函数如下：\n

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

\n

方法

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

    \n

  1. 实际上，我已经在Mathematica上做了类似的工作。我对函数进行一次和两次求导。我查看一阶导数为0的点，计算它们的值和位置。然后我在这些位置处对二阶导数进行计算，以确定它们是极小值还是极大值。

    2. \n

  2. 我还在思考是否只需创建一个二维数组来存储函数在x和y上的值，并找到该数组的最大值和最小值。但这需要我知道如何精确定义x和y的网格以可靠地捕捉函数的行为。

    3. \n

\n对于后一种情况，我已经找到了一些方法，比如这个。\n我只是想知道，哪种方法在Python中更加高效、快速、准确，甚至更加优雅？

0

1 答案

0 Comments

用户定义函数的局部极大值和极小值的查找是一个常见的问题。然而，对于复杂函数，很难找到解析解。解决这个问题的一种方法是使用符号计算方法，但是对于一般的二元方程组，目前还没有符号解法。因此，需要采用数值优化方法来解决这个问题。

对于简单的函数，可以使用SymPy库来进行符号计算。下面是一个使用SymPy找到函数的驻点并根据海森矩阵的特征值对其进行分类的完整示例：

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)))
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))

这段代码先求解梯度为零的驻点，然后将驻点逐个代入海森矩阵中进行分类。需要注意的是，当海森矩阵只是半正定的时候，无法确定驻点的性质，需要进行额外处理。

然而，使用符号计算方法很难找到所有的解（因为解有无穷多个）。可以尝试使用SymPy的`solveset`函数来处理无穷多解的情况，但是处理起来会比较困难。

另一种解决方法是使用数值优化方法，例如SciPy库提供的各种数值最小化算法，包括`brute force`方法。这些方法可以找到一个最小值，通过将函数取负可以找到最大值。但是需要注意的是，每次运行只能找到一个最值，而且初始搜索点的选择可能会导致找到不同的最值。此外，这些方法不能找到鞍点。

为了综合利用符号计算和数值优化方法，可以使用SymPy的`lambdify`函数将符号表达式转化为可以传递给SciPy数值优化方法的Python函数。下面是一个示例代码：

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

这段代码将梯度函数转化为可以传递给`fsolve`方法的函数，并通过初始猜测得到一个驻点。需要注意的是，每次运行`fsolve`只能得到一个解，因此无法找到所有解。

用户定义函数的局部极值点的查找是一个复杂的问题。可以使用符号计算方法找到一些解析解，也可以使用数值优化方法找到数值解。通过综合利用这两种方法，可以得到更全面的结果。

0