在矩阵上使用 scipy.fsolve
Using scipy.fsolve on a matrix
在获得帮助后 我一直在尝试在我的脚本中实现它，但我 运行 巧妙地失败了。
我需要对 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.
optimize.fsolve中的func可以接受一维向量，但不能接受二维数组。
因此，即使 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) 的二维数组。二维数组f * np.exp(-g*x[:, None]),
这也是形状 (4072*3080, 47)，然后在最后一个轴上求和（即 axis=-1）。
这将留下形状为 (4072*3080,) 的一维数组。 fsolve 求解 x 和 returns 形状为 (4072*3080,).
的一维数组
res = res.reshape(Kshape) 将解重塑为 (4072, 3080).
将 n 标量求根问题的解决方案转换为矢量化 n 维求根问题通常不是一个好主意。标量求解器做出的奇异决策通常在所有输入变体上并不统一，并且系统的任何通用求解器都会尝试计算或近似雅可比行列式。虽然我们知道此设置是对角线，但将其传达给求解器可能会很复杂。此外，如前所述，当对矢量化问题统一采取步长、线搜索等决策时，几乎可以保证对每个标量问题都是次优的。
最好的策略是减少要解决的标量问题的数量，或者至少使用一维性质来获得反函数的良好近似，然后只需要很少的迭代就可以完善解决方案。
K
中少量整数值的情况
您想在密切相关的输入上多次调用一个复杂的函数（迭代反演），甚至可能只取有限数量的值。
更快计算的第一个想法是首先确定 K 的值范围，为每个实际存在的值计算 x 的值，然后将这些值分配回去图像数组。例如，如果 K 仅采用从 0 到 255 的整数值，则生成一个数组 x_vals，其中 func(x_vals[k],f,g,0)=k 允许将结果数组作为 x_vals[K]一拉
x_vals = np.array([0.1,0.2,0.3,0.4])
K = np.array([[1,2],[0,2],[3,3]]);
x=x_vals[K]; print x
array([[ 0.2, 0.3],
[ 0.1, 0.3],
[ 0.4, 0.4]])
K中非整数值范围相对较小的情况
如果K在一些相对（相对于f中的值）较小的范围内包含大量（非整数）值，它仍然可以为计算解决方案提供很大的改进对于一些采样 K_samples = np.linspace(np.amin(K), np.amax(K), 256) 以便 func(x_samples[k],f,g,0)=K_samples[k])`
然后解决方案或至少非常好的近似使用插值
获得进一步迭代细化
x = np.interp(K, K_samples, x_samples)
使用 func 的前向求值
如果向量f和g都是正值，那么指数之和是一个单调下降函数，其反函数可以通过简单的查找得到table 和它的线性插值。对于指数和函数的值，一个得到
sum(f)*exp(-max(g)*x) <= K <= sum(f)*exp(-min(g)*x)
这样就可以计算 x 的范围为
- log(max(K)/sum(f)) / max(g) <= x <= - log(min(K)/sum(f)) / min(g)
使用这些上限和下限生成一个数组x_samples 通过对func 的前向评估得到相应的K_samples。如上所述进行线性插值以获得反函数的近似值。
您可以通过创建一个接受 N 个输入和 N 个输出的函数来向量化问题，其中 N 是像素数。这涉及展平输入图像并将其视为一维数组。在此设置中，输入独立于输出，因此雅可比矩阵是对角线的。因为 fsolve 计算雅可比矩阵的完全近似值，所以您最终会 运行 内存不足 (MemoryError)。相反，您可以将 scipy.optimize.root 与 method='diagbroyden' 一起使用，它通过仅跟踪对角线雅可比来使用近似值：
import numpy as np
import scipy.optimize as optimize
def func(x, f, g, K):
return np.sum(f * np.exp(-g*x[:, None]), axis=1) - K
np.random.seed(123)
f = np.random.uniform(size=(47,))
g = np.random.uniform(size=f.shape)
img = np.random.uniform(size=(4072, 3328)).ravel()
K = func(img, f, g, 0)
res = optimize.root(func, method='diagbroyden', x0=0.5*np.ones(img.size), args=(f, g, K))
print('Success:', res.success)
print('Message:', res.message)
assert np.allclose(img, res.x)
但是，使用此方法您无法利用可为您的特定函数计算的解析导数。
在获得帮助后
我需要对 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.
optimize.fsolve中的func可以接受一维向量，但不能接受二维数组。
因此，即使 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) 的二维数组。二维数组f * np.exp(-g*x[:, None]),
这也是形状 (4072*3080, 47)，然后在最后一个轴上求和（即 axis=-1）。
这将留下形状为 (4072*3080,) 的一维数组。 fsolve 求解 x 和 returns 形状为 (4072*3080,).
res = res.reshape(Kshape) 将解重塑为 (4072, 3080).
将 n 标量求根问题的解决方案转换为矢量化 n 维求根问题通常不是一个好主意。标量求解器做出的奇异决策通常在所有输入变体上并不统一，并且系统的任何通用求解器都会尝试计算或近似雅可比行列式。虽然我们知道此设置是对角线，但将其传达给求解器可能会很复杂。此外，如前所述，当对矢量化问题统一采取步长、线搜索等决策时，几乎可以保证对每个标量问题都是次优的。
最好的策略是减少要解决的标量问题的数量，或者至少使用一维性质来获得反函数的良好近似，然后只需要很少的迭代就可以完善解决方案。
K
中少量整数值的情况您想在密切相关的输入上多次调用一个复杂的函数（迭代反演），甚至可能只取有限数量的值。
更快计算的第一个想法是首先确定 K 的值范围，为每个实际存在的值计算 x 的值，然后将这些值分配回去图像数组。例如，如果 K 仅采用从 0 到 255 的整数值，则生成一个数组 x_vals，其中 func(x_vals[k],f,g,0)=k 允许将结果数组作为 x_vals[K]一拉
x_vals = np.array([0.1,0.2,0.3,0.4])
K = np.array([[1,2],[0,2],[3,3]]);
x=x_vals[K]; print x
array([[ 0.2, 0.3],
[ 0.1, 0.3],
[ 0.4, 0.4]])
K中非整数值范围相对较小的情况
如果K在一些相对（相对于f中的值）较小的范围内包含大量（非整数）值，它仍然可以为计算解决方案提供很大的改进对于一些采样 K_samples = np.linspace(np.amin(K), np.amax(K), 256) 以便 func(x_samples[k],f,g,0)=K_samples[k])`
然后解决方案或至少非常好的近似使用插值
获得进一步迭代细化x = np.interp(K, K_samples, x_samples)
使用 func 的前向求值
如果向量f和g都是正值，那么指数之和是一个单调下降函数，其反函数可以通过简单的查找得到table 和它的线性插值。对于指数和函数的值，一个得到
sum(f)*exp(-max(g)*x) <= K <= sum(f)*exp(-min(g)*x)
这样就可以计算 x 的范围为
- log(max(K)/sum(f)) / max(g) <= x <= - log(min(K)/sum(f)) / min(g)
使用这些上限和下限生成一个数组x_samples 通过对func 的前向评估得到相应的K_samples。如上所述进行线性插值以获得反函数的近似值。
您可以通过创建一个接受 N 个输入和 N 个输出的函数来向量化问题，其中 N 是像素数。这涉及展平输入图像并将其视为一维数组。在此设置中，输入独立于输出，因此雅可比矩阵是对角线的。因为 fsolve 计算雅可比矩阵的完全近似值，所以您最终会 运行 内存不足 (MemoryError)。相反，您可以将 scipy.optimize.root 与 method='diagbroyden' 一起使用，它通过仅跟踪对角线雅可比来使用近似值：
import numpy as np
import scipy.optimize as optimize
def func(x, f, g, K):
return np.sum(f * np.exp(-g*x[:, None]), axis=1) - K
np.random.seed(123)
f = np.random.uniform(size=(47,))
g = np.random.uniform(size=f.shape)
img = np.random.uniform(size=(4072, 3328)).ravel()
K = func(img, f, g, 0)
res = optimize.root(func, method='diagbroyden', x0=0.5*np.ones(img.size), args=(f, g, K))
print('Success:', res.success)
print('Message:', res.message)
assert np.allclose(img, res.x)
但是，使用此方法您无法利用可为您的特定函数计算的解析导数。