python - 在矩阵上使用 scipy.fsolve

Question

根据我得到的帮助here我一直在尝试在我的脚本中实现它，但我未能巧妙地运行它。

我需要对 4072x3080 图像的每个像素使用此算法，整个过程大约需要 1 小时 30 秒，所以我尝试以某种方式强制它，但出现此错误:

ValueError                                Traceback (most recent call last)
<ipython-input-12-99c1f41dbba7> in <module>()
----> 1 res = scipy.optimize.fsolve(func, x0=np.ones((K.shape[0], K.shape[1])), args=(f[:,None], g[:,None], K))

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    146                'diag': diag}
    147 
--> 148     res = _root_hybr(func, x0, args, jac=fprime, **options)
    149     if full_output:
    150         x = res['x']

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    212     if not isinstance(args, tuple):
    213         args = (args,)
--> 214     shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
    215     if epsfcn is None:
    216         epsfcn = finfo(dtype).eps

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     25 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     26                 output_shape=None):
---> 27     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     28     if (output_shape is not None) and (shape(res) != output_shape):
     29         if (output_shape[0] != 1):

<ipython-input-7-911c817cb57d> in func(x, f, g, K)
      1 def func(x, f, g, K):
----> 2     return np.sum(f * np.exp(-g*x), axis=0) - K
      3 
      4 
      5 def derivative(x, f, g, K):

ValueError: operands could not be broadcast together with shapes (13551616,) (4072,3328) 

这是我一直在尝试的代码:

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x), axis=0) - K


def derivative(x, f, g, K):
    return np.sum(-g*f * np.exp(-g*x), axis=0)

+

res = scipy.optimize.fsolve(func, x0=np.ones((K.shape[0], K.shape[1])), args=(f[:,None], g[:,None], K))

f 和 g 都是 (47,) 数组，其中 K 是 ( 4072, 3328) 图片

否则，缓慢的过程会这样进行:(但无论如何，这个过程都会因导数而失败。

res = np.ones((mbn.shape[0],mbn.shape[1]))
for i_x in range(0,mbn.shape[0]):
    if i_x%10 == 0:
        print i_x/4070 
    for i_y in range(0,mbn.shape[1]):
        res[i_x,i_y] = scipy.optimize.fsolve(func, x0=1, args=(f[:], g[:], K[i_x,i_y]) )

如果我尝试将慢速方法与导数一起使用，这将是错误

    ---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: object of too small depth for desired array

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-8-3587dcccfd93> in <module>()
      4         print i_x/4070
      5     for i_y in range(0,mbn.shape[1]):
----> 6         res[i_x,i_y] = scipy.optimize.fsolve(func, fprime=derivative, x0=1, args=(f[:], g[:], K[i_x,i_y]) )

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    146                'diag': diag}
    147 
--> 148     res = _root_hybr(func, x0, args, jac=fprime, **options)
    149     if full_output:
    150         x = res['x']

/*/python2.7/site-packages/scipy/optimize/minpack.pyc in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    232         with _MINPACK_LOCK:
    233             retval = _minpack._hybrj(func, Dfun, x0, args, 1,
--> 234                                      col_deriv, xtol, maxfev, factor, diag)
    235 
    236     x, status = retval[0], retval[-1]

error: Result from function call is not a proper array of floats.

Answer 1

optimize.fsolve 中的 func可以接受一维向量，但不能接受二维数组。因此，即使 K 和 x 是二维的，对于此计算，我们有必要将它们 reshape 为一维数组。

Kshape = K.shape
K = K.ravel()

然后在调用optimize.fsolve之后，您可以将结果重新整形为二维:

res = optimize.fsolve(func, x0=np.ones(K.shape).ravel(), args=(f, g, K))
res = res.reshape(Kshape)

<小时/>

然后，您可以通过这样编写计算来避免双重 for 循环:

import numpy as np
import scipy.optimize as optimize

np.random.seed(123)

def func(x, f, g, K):
    return np.sum(f * np.exp(-g*x[:, None]), axis=-1) - K


def derivative(x, f, g, K):
    return np.sum(-g*f * np.exp(-g*x[:, None]), axis=-1)


f = np.random.uniform(size=(47,))
g = np.random.uniform(size=f.shape)
K = np.random.uniform(size=(4072,3080))
Kshape = K.shape
K = K.ravel()

res = optimize.fsolve(func, x0=np.ones(K.shape).ravel(), args=(f, g, K))
res = res.reshape(Kshape)
print(res)

请注意，g*x[:, None] 使用 broadcasting生成形状为 (4072*3080, 47) 的二维数组。二维数组 f * np.exp(-g*x[:, None])，其形状也是 (4072*3080, 47)，然后在最后一个轴(即 axis=-1)上求和。这会留下一个形状为 (4072*3080,) 的一维数组。 fsolve 求解 x 并返回形状为 (4072*3080,) 的一维数组。

res = res.reshape(Kshape) 将解 reshape 为形状 (4072, 3080)。