numpy.linalg.求解右侧超过三维

Question

我正在尝试求解一个具有 3x3 矩阵a和任意形状(3, ...)的右侧b的方程组。如果b有一个或两个维度，numpy.linalg.solve就可以了。不过，它细分为更多维度：

import numpy
a = numpy.random.rand(3, 3)
b = numpy.random.rand(3)
numpy.linalg.solve(a, b)  # okay
b = numpy.random.rand(3, 4)
numpy.linalg.solve(a, b)  # okay
b = numpy.random.rand(3, 4, 5)
numpy.linalg.solve(a, b)  # ERR

ValueError: solve: Input operand 1 has a mismatch in its core 
dimension 0, with gufunc signature (m,m),(m,n)->(m,n) (size 5 is 
different from 3)

我本来期望一个形状(3, 4, 5)的输出数组sol其解决方案对应于右侧b[:, i, j] sol[:, i, j]。

关于如何最好地解决此问题的任何提示？

Answer 1

暂时将b整形为(3, 20)，求解线性系统，然后将结果数组整形为b(3， 4， 5(的原始形状：

In [34]: a = numpy.random.rand(3, 3)
In [35]: b = numpy.random.rand(3, 4, 5)
In [36]: x = numpy.linalg.solve(a, b.reshape(b.shape[0], -1)).reshape(b.shape)

或

使用 np.swapaxes 将b的第一个轴与第二个轴交换，求解线性系统，然后恢复轴：

In [58]: x = np.swapaxes(np.linalg.solve(a, np.swapaxes(b, 0, 1)), 0, 1)

---

健全性检查：

In [38]: np.einsum('ij,jkl', a, x)
Out[38]: 
array([[[ 0.44859955,  0.22967928,  0.74336067,  0.47440575,  0.53798895],
        [ 0.80045696,  0.54138958,  0.89870834,  0.56862419,  0.28217437],
        [ 0.02093982,  0.78534718,  0.77208236,  0.41568151,  0.95100661],
        [ 0.03820421,  0.47067312,  0.71928294,  0.30852615,  0.64454321]],
       [[ 0.31757072,  0.30527186,  0.36768759,  0.95869289,  0.86601996],
        [ 0.60616508,  0.69927063,  0.53470332,  0.88906606,  0.76066344],
        [ 0.95411847,  0.51116677,  0.29338398,  0.04418815,  0.96210206],
        [ 0.23449429,  0.64159963,  0.7732404 ,  0.4314741 ,  0.81279619]],
       [[ 0.6399571 ,  0.57640652,  0.0186913 ,  0.66304489,  0.83372239],
        [ 0.28426522,  0.62367363,  0.37163699,  0.78217433,  0.90573787],
        [ 0.91066088,  0.06699638,  0.43079394,  0.00263537,  0.399102  ],
        [ 0.17711441,  0.48724858,  0.05526752,  0.34251648,  0.94059739]]])
In [39]: b
Out[39]: 
array([[[ 0.44859955,  0.22967928,  0.74336067,  0.47440575,  0.53798895],
        [ 0.80045696,  0.54138958,  0.89870834,  0.56862419,  0.28217437],
        [ 0.02093982,  0.78534718,  0.77208236,  0.41568151,  0.95100661],
        [ 0.03820421,  0.47067312,  0.71928294,  0.30852615,  0.64454321]],
       [[ 0.31757072,  0.30527186,  0.36768759,  0.95869289,  0.86601996],
        [ 0.60616508,  0.69927063,  0.53470332,  0.88906606,  0.76066344],
        [ 0.95411847,  0.51116677,  0.29338398,  0.04418815,  0.96210206],
        [ 0.23449429,  0.64159963,  0.7732404 ,  0.4314741 ,  0.81279619]],
       [[ 0.6399571 ,  0.57640652,  0.0186913 ,  0.66304489,  0.83372239],
        [ 0.28426522,  0.62367363,  0.37163699,  0.78217433,  0.90573787],
        [ 0.91066088,  0.06699638,  0.43079394,  0.00263537,  0.399102  ],
        [ 0.17711441,  0.48724858,  0.05526752,  0.34251648,  0.94059739]]])

使用 np.allclose()，这样您就不必手动浏览数字和检查，特别是对于大型数组：

In [32]: b_ = np.einsum('ij,jkl', a, x)
In [33]: np.allclose(b, b_)
Out[33]: True

Answer 2

我想

补充一点，手册明确指出：

a ： (...， M， M( array_like

系数矩阵。

b ： {(...， M，

(， (...， M， K(}，

array_like

纵坐标或"因变量"值。

因此，最后一个维度之前的值必须与 a (M( 的最后两个维度相同。除此之外，它的行为符合您的预期 - 可以生成更多维度，返回与B相同维度的结果。通过这种方式，自然计算Ax=B的解，并自动转换维度 - 只需求解许多具有维度 (M，K( 的解的方程组，并将它们嵌入外部维度。在您的情况下，3在开头而不是中间会混淆算法。具有 3 个维度的示例;

>>> a=np.random.rand(9).reshape(3,3)
>>> b=np.random.rand(12).reshape(2,3,2)
>>> np.linalg.solve(a,b)
array([[[-0.63673083,  0.57508091],
        [ 0.87653408,  0.46092677],
        [ 0.61128222, -0.19641607]],
       [[-0.91645601,  1.30939652],
        [ 0.83591936, -0.17006344],
        [ 0.19086912,  0.29082206]]])