python - 使用 ctypes 在 Fortran 函数上使用 scipy.optimize.minimize 的错误结果

Question

我目前有一个 Fortran 函数，我想使用 SciPy 对其进行优化，并使用 Ctypes 对其进行包装。这可能吗？也许我在实现过程中做错了什么。例如，假设我有:

成本.f90

module cost_fn
  use iso_c_binding, only: c_float
  implicit none

contains

  function sin_2_cos(x,y) bind(c)
      real(c_float) :: x, y, sin_2_cos
      sin_2_cos = sin(x)**2 * cos(y)
  end function sin_2_cos

end module cost_fn

我编译的:

gfortran -fPIC -shared -g -o cost.so cost.f90

然后尝试通过以下方式找到(局部)最小值:

成本.py

#!/usr/bin/env python

from ctypes import *
import numpy as np
import scipy.optimize as sopt

cost = cdll.LoadLibrary('./cost.so')

cost.sin_2_cos.argtypes = [POINTER(c_float), POINTER(c_float)]
cost.sin_2_cos.restype = c_float

def f2(x):
    return cost.sin_2_cos(c_float(x[0]), c_float(x[1]))
    # return np.sin(x[0])**2 * np.cos(x[1])

# print(f2( [1, 1] ))
# print(f2( [0.5 * np.pi, np.pi] ))

print( sopt.minimize( f2, (1.0, 1.0), options={'disp': True}, tol=1e-8) )

我期望局部最小值 f2(pi/2, pi) = -1。当我使用 cost.sin_2_cos 返回值调用 f2 时，“最小值”仅在初始猜测 (1, 1) 处给出。如果我使用“Python”返回值调用 f2，optimize 会找到正确的最小值。

我尝试重新定义 sin_2_cos 以采用维度(2)数组输入，但看到了类似的行为。也许我需要直接使用 minimize 调用 sin_2_cos (但是我将如何为参数指定 c_float )？任何想法表示赞赏!

编辑:对于下面的注释，请注意两个注释的 print(f2(...)) 行生成预期值。因此，我相信 Fortran 函数正在通过 Python f2 函数正确调用。

Answer 1

您的 Fortran 代码使用单精度浮点值(即 32 位 float )。 (在 C 语言中，float 是单精度，double 是 double 。)scipy.optimize.minimize() 使用的默认方法使用有限差异来近似函数的导数。也就是说，为了估计导数 f'(x0)，它会计算 (f(x0+eps) - f(x0))/eps，其中 eps 是步长。用于计算有限差分的默认步长约为 1.49e-08。不幸的是，这个值小于值 1 周围的单精度值的间距。因此，当最小化器将 eps 添加到 1 时，结果仍然是 1。这意味着函数在同一点计算，有限差分结果为 0。这是最小值的条件，因此求解器决定完成。

求解器选项eps设置有限差分步长。将其设置为大于 1.19e-7 的值。例如，如果我使用 options={'disp': True, 'eps': 2e-6}，我会得到一个解决方案。

顺便说一句，您可以使用 numpy.finfo() 找到该值 1.19e-7:

In [4]: np.finfo(np.float32(1.0)).eps
Out[4]: 1.1920929e-07

如果您在 minimize() 函数中使用选项 method='nelder-mead'，您也会得到一个解决方案。该方法不依赖于有限差分。

最后，您可以将 Fortran 代码转换为使用 double :

module cost_fn
  use iso_c_binding, only: c_double
  implicit none

contains

  function sin_2_cos(x,y) bind(c)
      real(c_double) :: x, y, sin_2_cos
      sin_2_cos = sin(x)**2 * cos(y)
  end function sin_2_cos

end module cost_fn

然后更改 Python 代码以使用 ctypes.c_double 而不是 ctypes.c_float。