SciPy:使用b向量数组的元素方式非负最小二乘

Question

我需要解决线性问题Ax = b,x使用最小二乘法获得.所有元素x必须是非负的,所以我使用scipy.optimize.nnls(文档在这里).

麻烦的是,我需要用单个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]

这段代码非常实用,但感觉完全不优雅.更重要的是,它可能会非常慢(我的实际代码使用非常大的数据集,并且通过随机参数更改循环数千次,因此效率很重要).

不久之后,我问了一个关于元素矩阵乘法的非常相似的问题.我被介绍过np.einsum,在许多情况下证明它非常有用.我原本希望最小二乘解决方案有类似的功能,但一直找不到任何东西.如果有人知道可能有效的功能,或者有效/热情地解决这个问题的替代方法,那将非常感激!

Answer 1

NNLS没有封闭形式的解决方案,除了为设计矩阵共享内存之外,没有通过一起处理问题可以获得算法加速.虽然将多目标功能降低到C级可能会导致一些加速,但看起来scipy实现一次只支持一个目标,因此循环看起来像这里唯一的选择.这个问题非常平行,因此您可以使用eg 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级治疗.