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]].接下来,我们通过将梯度等于零来找到静止点,并将它们逐个插入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具有相同的Hessian at(0, 0)但不同的行为).输出:
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 vs solveset,但无论如何,处理无限多的解决方案很难.
使用SciPy进行数值优化
SciPy提供了很多numerical minimization routines,包括brute force(这是你的方法2;通常它非常慢).这些是强大的方法,但请考虑这些要点.
>每次运行只会找到一个最小值.
>用-f替换f也可以找到最大值.
>更改搜索的起点(最小化的参数x0)可能会产生另一个最大值或最小值.尽管如此,你永远不会知道还有其他你还没见过的极值.
>这些都不会找到马鞍点.
混合策略
使用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只会产生一个.