向量化 LU 分解求解多个 b答案

【问题标题】：Vectorized LU decomposition solve for multiple b向量化 LU 分解求解多个 b
【发布时间】：2020-09-07 20:57:19
【问题描述】：

我正在使用 scipy 的 LU 分解对形状为 (n,n) 的二次矩阵 A 进行预处理，然后一遍又一遍地求解形状为 (...,n) 的多个右手边 B。但是scipy.linalg.lu_solve 只接受 b 的向量，而不接受像 (m,n) 或 (k,m,n) 这样的矩阵。如何包装lu_solve 以处理(...,n) 形状的参数？ Numpy 的linalg.solve 将接受多个b，但不允许分离的LU 因子和求解操作。

【问题讨论】：

标签： numpy scipy numpy-ndarray

【解决方案1】：

lu_solve的文档中没有提到，但实际上b可以包含多个向量。如果A 具有形状(n, n)，则b 可以具有形状(n, m)。例如，

In [44]: A                                                                                                           
Out[44]: 
array([[ 1.01,  0.02, -0.01],
       [ 0.02,  1.04, -0.02],
       [-0.01, -0.02,  1.01]])

In [45]: b                                                                                                           
Out[45]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

In [46]: lu = lu_factor(A)                                                                                           

In [47]: x = lu_solve(lu, b)                                                                                         

In [48]: x                                                                                                           
Out[48]: 
array([[ 0.        ,  0.98113208,  1.96226415,  2.94339623],
       [ 4.        ,  4.96226415,  5.9245283 ,  6.88679245],
       [ 8.        ,  9.01886792, 10.03773585, 11.05660377]])

In [49]: A.dot(x)                                                                                                    
Out[49]: 
array([[ 0.,  1.,  2.,  3.],
       [ 4.,  5.,  6.,  7.],
       [ 8.,  9., 10., 11.]])

高维b 必须具有形状(n, ...)。请注意，对于多于二维的形状，用A.dot(x) 测试结果将不起作用，因为x 的形状将与 NumPy 的矩阵乘法不兼容。例如，这里B 的形状为(3, 2, 5)：

In [40]: A                                                                      
Out[40]: 
array([[ 1.01,  0.02, -0.01],
       [ 0.02,  1.04, -0.02],
       [-0.01, -0.02,  1.01]])

In [41]: B = np.random.rand(3, 2, 5)                                            

In [42]: lu = lu_factor(A)                                                      

In [43]: x = lu_solve(lu, B)                                                    

In [44]: x.shape                                                                
Out[44]: (3, 2, 5)

In [45]: xx = np.moveaxis(x, 0, 1)                                              

In [46]: np.allclose(A.dot(xx), B)                                              
Out[46]: True

【讨论】：

``` 导入 scipy.linalg，numpy 作为 np n = 10 A = np.random.rand(n,n) LU_and_piv = scipy.linalg.lu_factor(A) B = np.random。 rand(2,5,n) X = scipy.linalg.lu_solve(LU_and_piv, B) ``` 抛出 ValueError：尺寸不兼容。我正在使用 scipy 1.5.0
我明白了。我将如何包装 scipy 的 lu_solve 以使其更像 numpy 的求解，其中 B 的形状必须为 (...,n) 并且输出类似于 (...,n)
在包装器中，您可以使用b = np.moveaxis(B, -1, 0) 之类的东西，调用lu_solve 和b 以获取x，然后返回np.moveaxis(x, 0, -1)。