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

Question

我目前有一个 Fortran 函数，我想使用 SciPy 使用 Ctypes 包装它来优化它。这可能吗？也许我在实施过程中做错了什么。例如，假设我有：

cost.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

然后尝试找到一个（局部）最小值：

cost.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 return 值调用 f2 时，"minimimum" 只是在 (1, 1) 的初始猜测中给出。如果我使用 "Python" return 值调用 f2，优化会找到正确的最小值。

我已经尝试重新定义 sin_2_cos 以获取维度 (2) 数组输入，但看到了类似的行为。也许我需要用最小化直接调用 sin_2_cos （但我将如何为参数指定 c_float ）？任何想法表示赞赏！

编辑：对于下面的评论，请注意两条评论的 print(f2(...)) 行产生了预期值。因此，我相信通过 Python f2 函数可以正确调用 Fortran 函数。

Answer 1

您的 Fortran 代码使用单精度浮点值（即 32 位浮点数）。（在 C 中，float 是单精度，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 代码转换为使用双精度：

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。

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

Wrong result from scipy.optimize.minimize used on Fortran function using ctypes

python

fortran

ctypes

scipy

cost.f90

cost.py