SciPy：使用多个b向量进行逐元素非负最小二乘

Question

我需要解决一个线性问题 Ax = b，通过最小二乘法来得到 x。而且 x 的所有元素都必须是非负的，所以我使用了 scipy.optimize.nnls（文档可以在这里找到）。

问题是，我需要多次解决这个问题，使用同一个 A 矩阵和多个 b 向量。我有一个三维的 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]

这段代码是完全可以运行的，但感觉不太优雅。更重要的是，它可能会非常慢（我实际的代码使用的是非常大的数据集，并且会循环几千次，随机改变参数，所以效率很重要）。

不久前，我问过一个非常相似的问题，关于逐元素的矩阵乘法。我了解到 np.einsum，在很多情况下都非常有用。我希望能找到一个类似的函数来解决最小二乘法的问题，但一直没能找到。如果有人知道可能有效的函数，或者有其他高效的解决方案，真是太感谢了！

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)

不过，如果你使用例如岭回归或普通最小二乘法（可能在你的情况下不太有用），那么就可以通过矩阵乘法得到一个简单的解决方案，所有的操作都可以通过一次重塑和矩阵乘法完成，把多目标的问题处理放到C语言层面。