我需要解决线性问题Ax = b
,使用最小二乘法获得x
。 x
的所有元素都必须是非负数,因此我使用的是 scipy.optimize.nnls
(文档 here)。
麻烦的是,我需要用一个 A
矩阵和许多 b
向量多次解决这个问题。我有一个 3d numpy ndarray
,其中沿轴 0
的向量是 b
向量,而其他两个轴对应于空间中的点。我希望将所有 x
向量输出到相应的数组,以便保留每个答案的空间信息。
问题的第一步看起来像这样:
A = np.random.rand(5,3)
b_array = B = np.random.rand(5,100,100)
x_array = np.zeros((3,100,100))
for i in range(100):
for j in range(100):
x_array[:,i,j] = sp.optimize.nnls(A, b_array[:,i,j])[0]
此代码功能完美,但感觉完全不优雅。更重要的是,它可能会非常慢(我的实际代码使用非常大的数据集并且循环了数千次随机参数更改,因此效率很重要)。
不久前我问过this very similar question关于逐元素矩阵乘法。我被介绍到 np.einsum
,它在许多情况下都被证明非常有用。我曾希望最小二乘解有类似的功能,但一直找不到任何东西。如果有人知道可能有效的函数,或者有效/Python 方式解决此问题的替代方法,将不胜感激!
最佳答案
NNLS 没有封闭形式的解决方案,除了为设计矩阵共享内存外,同时处理这些问题不会获得算法加速。尽管将多目标功能降低到 C 级别可能会导致一些加速,但看起来 scipy 实现一次只支持一个目标,因此循环看起来是这里唯一的选择。这个问题是并行的,所以你可以使用例如joblib
并行化循环如下
from joblib import Parallel, delayed
from itertools import product
from scipy.optimize import nnls
results = Parallel(n_jobs=10)(delayed(nnls)(A, b_array[:,i,j])[0]
for i, j in product(range(100), range(100)))
x_array = np.array(results).reshape(100, 100, -1).transpose(2, 0, 1)
但是,如果您使用例如岭回归或 OLS(可能对您的情况没有用),然后解决方案是可通过矩阵乘法获得的封闭形式,并且所有事情都可以在一次 reshape 和矩阵乘法中完成,将问题的多目标方面推向 C 级处理.
关于python - SciPy:使用 b 向量数组的逐元素非负最小二乘法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26718222/