在获得帮助之后，我一直尝试在脚本中实现here，但是我无法聪明地运行它。
我需要为4072x3080图像的每个像素使用此算法，整个过程大约需要1h30，因此我尝试以某种方式强制执行此操作，但出现此错误：
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.
最佳答案
func中的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]使用broadcasting生成形状为(4072*3080, 47)的2D数组。 2D数组f * np.exp(-g*x[:, None])，然后也将形状也是
(4072*3080, 47)的形状在最后一个轴（即axis=-1）上相加。留下形状为
(4072*3080,)的一维数组。 fsolve求解x并返回形状为(4072*3080,)的一维数组。res = res.reshape(Kshape)将解决方案重塑为形状为(4072, 3080)。