反转或点kxnxn矩阵的快速方法 [英] fast way to invert or dot kxnxn matrix

问题描述

是否有使用numpy计算kxnxn矩阵的逆的快速方法(在每个k切片上计算逆)?换句话说，有一种方法可以向量化以下代码:

>>>from numpy.linalg import inv
>>>a-random(4*2*2).reshape(4,2,2)
>>>b=a.copy()
>>>for k in range(len(a)):
>>>    b[k,:,:] = inv(a[k,:,:])

解决方案

首先要求逆.我已经研究了np.linalg.tensorinv和np.linalg.tensorsolve.

不幸的是，我认为tensorinv不会给您想要的东西.它需要数组为正方形".这不包括您要执行的操作，因为它们对正方形的定义是np.prod(a[:i]) == np.prod(a[i:])，其中i是0、1或2(通常是数组的轴之一)；这可以作为tensorinv的第三个参数ind给出.这意味着，如果您有长度为M的NxN个矩阵的通用数组，则需要具有(对于i = 1)NxN == NxM，通常这是不正确的(在您的示例中实际上是正确的，但无论如何它都无法给出正确的答案).

现在，tensorsolve也许可能会有所帮助.但是，在将a矩阵数组作为第一个参数传递给tensorsolve之前，这将涉及一些繁重的构造工作.因为我们希望b是矩阵数组方程" a*b = 1(其中1是单位矩阵的数组)的解，而1的形状与a和，我们不能简单地将您在上面定义的a作为tensorsolve的第一个参数提供.相反，它必须是形状为(M，N，N，M，N，N)或(M，N，N，N，M，N)或(M，N，N，N，N，N，M的形状的数组).这是必须的，因为tensorsolve将在这最后三个轴上与b相乘，并对它们求和，以使结果(函数的第二个自变量)再次具有形状(M，N，N).

其次，关于点积(您的标题表明，这也是您的问题的一部分).这是非常可行的.两种选择.

首先:詹姆斯·汉斯曼(James Hensman)的这篇博文提供了一些好的建议.

第二:为了清晰起见，我个人更喜欢使用np.einsum.例如:

a=np.random.random((7,2,2))
b=np.random.random((7,2,2))
np.einsum('ijk,ikl->ijl', a,b)

这将对数组a和b中的所有7个矩阵"进行矩阵相乘.它似乎比上面博客文章中的数组方法慢大约2倍，但是仍然比示例中使用for循环快大约70倍.实际上，对于较大的数组(例如10000个5x5矩阵)，einsum方法似乎要快一些(不确定为什么).

希望有帮助.

Is there a fast way to calculate the inverse of a kxnxn matrix using numpy (the inverse being calculated at each k-slice)? In other words, is there a way to vectorize the following code:

>>>from numpy.linalg import inv
>>>a-random(4*2*2).reshape(4,2,2)
>>>b=a.copy()
>>>for k in range(len(a)):
>>>    b[k,:,:] = inv(a[k,:,:])

解决方案

First about getting the inverse. I have looked into both np.linalg.tensorinv and np.linalg.tensorsolve.

I think unfortunately tensorinv will not give you what you want. It needs the array to be "square". This excludes what you want to do, because their definition of square is that np.prod(a[:i]) == np.prod(a[i:]) where i is 0, 1 or 2 (one of the axes of the array in general); this can be given as the third argument ind of tensorinv. This means that if you have a general array of NxN matrices of length M, you need to have e.g. (for i = 1) NxN == NxM, which is not true in general (it actually is true in your example, but it does not give the correct answer anyway).

Now, maybe something is possible with tensorsolve. This would however involve some heavy construction work on the a matrix-array before it is passed as the first argument to tensorsolve. Because we would want b to be the solution of the "matrix-array equation" a*b = 1 (where 1 is an array of identity matrices) and 1 would have the same shape as a and b, we cannot simply supply the a you defined above as the first argument to tensorsolve. Rather, it needs to be an array with shape (M,N,N,M,N,N) or (M,N,N,N,M,N) or (M,N,N,N,N,M). This is necessary, because tensorsolve would multiply with b over these last three axes and also sum over them so that the result (the second argument to the function) is again of shape (M,N,N).

Then secondly, about dot products (your title suggests that's also part of your question). This is very doable. Two options.

First: this blog post by James Hensman gives some good suggestions.

Second: I personally like using np.einsum better for clarity. E.g.:

a=np.random.random((7,2,2))
b=np.random.random((7,2,2))
np.einsum('ijk,ikl->ijl', a,b)

This will matrix-multiply all 7 "matrices" in arrays a and b. It seems to be about 2 times slower than the array-method from the blog post above, but it's still about 70 times faster than using a for loop as in your example. In fact, with larger arrays (e.g. 10000 5x5 matrices) the einsum method seems to be slightly faster (not sure why).

Hope that helps.

这篇关于反转或点kxnxn矩阵的快速方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！