时间戳

首页 > 解决方案 > 在矩阵上使用 scipy.fsolve

python - 在矩阵上使用 scipy.fsolve

问题描述

在我在这里得到的帮助之后，我一直在尝试在我的脚本中实现它，但我未能巧妙地运行它。

我需要对 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.

标签： pythonpython-2.7scipynumericnumerical-methods

解决方案

funcin可以接受一维optimize.fsolve向量，但不能接受二维数组。所以即使K和x是二维的，对于这个计算，我们有必要将它们重塑为一维数组。

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

然后在调用之后optimize.fsolve，您可以将结果重新整形为再次 2D：

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]使用广播来生成 shape 的 2D 数组(4072*3080, 47)。然后在最后一个轴（即）上对f * np.exp(-g*x[:, None])同样具有 shape的 2D 数组求和。这留下了 shape 的一维数组。求解并返回 shape 的一维数组。(4072*3080, 47)axis=-1(4072*3080,)fsolvex(4072*3080,)

res = res.reshape(Kshape)将解决方案重塑为具有形状(4072, 3080)。

推荐阅读