SciPy：使用 b 向量数组的逐元素非负最小二乘

Question

我需要解决线性问题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 解决这个问题的替代方法，那将不胜感激！

Answer 1

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（在您的情况下可能没有用），然后解决方案是可通过矩阵乘法获得的封闭形式，并且一切都可以通过一次整形和矩阵乘法完成，将问题的多目标方面推到C级处理.