矩阵求逆 (3,3) python - 硬编码与 numpy.linalg.inv

对于大量矩阵，我需要计算定义为的距离度量：

尽管我确实知道强烈建议不要使用矩阵求逆，但我没有找到解决方法。因此，我尝试通过对矩阵求逆进行硬编码来提高性能，因为所有矩阵的大小均为 (3,3)。

我预计这至少会是一个微小的改进，但事实并非如此。

为什么 numpy.linalg.inv 比这种硬编码矩阵求逆更快/性能更高？

此外，我有什么替代方案来改善这个瓶颈？

def inversion(m):    
    m1, m2, m3, m4, m5, m6, m7, m8, m9 = m.flatten()
    determinant = m1*m5*m9 + m4*m8*m3 + m7*m2*m6 - m1*m6*m8 - m3*m5*m7 - m2*m4*m9  
    return np.array([[m5*m9-m6*m8, m3*m8-m2*m9, m2*m6-m3*m5],
                     [m6*m7-m4*m9, m1*m9-m3*m7, m3*m4-m1*m6],
                     [m4*m8-m5*m7, m2*m7-m1*m8, m1*m5-m2*m4]])/determinant

与随机 3*3 矩阵进行时序比较：

%timeit np.linalg.inv(a)

100000 个循环，3 次最佳：每个循环 12.5 µs

%timeit inversion(a)

100000 个循环，3 次最佳：每个循环 13.9 µs

密切相关但完全不重复的是这个post https://codereview.stackexchange.com/questions/154926/efficient-extraction-of-patch-features-over-an-image关于代码审查，解释了背景和整个功能。

EDIT:正如 @Divakar 在评论中建议的那样，m.ravel() 而不是 m.flatten() 稍微改进了反转，因此时序比较现在产生：

numpy - 100000 个循环，3 个循环中最好的：每个循环 12.6 µs

硬编码 - 100000 个循环，3 次最佳：每个循环 12.8 µs

尽管差距正在缩小，但硬编码的速度仍然较慢。为何如此？

这是一个简单的优化，节省了 9 次乘法和 3 次减法

def inversion(m):    
    m1, m2, m3, m4, m5, m6, m7, m8, m9 = m.ravel()
    inv = np.array([[m5*m9-m6*m8, m3*m8-m2*m9, m2*m6-m3*m5],
                    [m6*m7-m4*m9, m1*m9-m3*m7, m3*m4-m1*m6],
                    [m4*m8-m5*m7, m2*m7-m1*m8, m1*m5-m2*m4]])
    return inv / np.dot(inv[0], m[:, 0])

你可以通过一次性完成整个跟踪来挤出更多的操作（如果我数正确的话，还有 24 次乘法）：

def det(m):
   m1, m2, m3, m4, m5, m6, m7, m8, m9 = m.ravel()
   return np.dot(m[:, 0], [m5*m9-m6*m8, m3*m8-m2*m9, m2*m6-m3*m5])
   # or try m1*(m5*m9-m6*m8) + m4*(m3*m8-m2*m9) + m7*(m2*m6-m3*m5)
   # probably the fastest would be to inline the two calls to det
   # I'm not doing it here because of readability but you should try it

def dist(m, n):
   m1, m2, m3, m4, m5, m6, m7, m8, m9 = m.ravel()
   n1, n2, n3, n4, n5, n6, n7, n8, n9 = n.ravel()
   return 0.5 * np.dot(
       m.ravel()/det(m) + n.ravel()/det(n),
       [m5*n9-m6*n8, m6*n7-m4*n9, m4*n8-m5*n7, n3*m8-n2*m9, n1*m9-n3*m7,
        n2*m7-n1*m8, m2*n6-m3*n5, m3*n4-m1*n6, m1*n5-m2*n4])

好的，这是内联版本：

import numpy as np
from timeit import timeit

def dist(m, n):
   m1, m2, m3, m4, m5, m6, m7, m8, m9 = m.ravel()
   n1, n2, n3, n4, n5, n6, n7, n8, n9 = n.ravel()
   return 0.5 * np.dot(
       m.ravel()/(m1*(m5*m9-m6*m8) + m4*(m3*m8-m2*m9) + m7*(m2*m6-m3*m5))
       + n.ravel()/(n1*(n5*n9-n6*n8) + n4*(n3*n8-n2*n9) + n7*(n2*n6-n3*n5)),
       [m5*n9-m6*n8, m6*n7-m4*n9, m4*n8-m5*n7, n3*m8-n2*m9, n1*m9-n3*m7,
        n2*m7-n1*m8, m2*n6-m3*n5, m3*n4-m1*n6, m1*n5-m2*n4])

def dist_np(m, n):
    return 0.5 * np.diag(np.linalg.inv(m)@n + np.linalg.inv(n)@m).sum()

for i in range(3):
    A, B = np.random.random((2,3,3))
    print(dist(A, B), dist_np(A, B))
    print('pp     ', timeit('f(A,B)', number=10000, globals={'f':dist, 'A':A, 'B':B}))
    print('numpy  ', timeit('f(A,B)', number=10000, globals={'f':dist_np, 'A':A, 'B':B}))

prints:

2.20109953156 2.20109953156
pp      0.13215381593909115
numpy   0.4334693900309503
7.50799877993 7.50799877993
pp      0.13934064202476293
numpy   0.32861811900511384
-0.780284449609 -0.780284449609
pp      0.1258618349675089
numpy   0.3110764700686559

请注意，您可以通过使用该函数的矢量化版本进行批处理来节省大量资金。该测试计算两批 100 个矩阵之间的所有 10,000 个成对距离：

def dist(m, n):
    m = np.moveaxis(np.reshape(m, m.shape[:-2] + (-1,)), -1, 0)
    n = np.moveaxis(np.reshape(n, n.shape[:-2] + (-1,)), -1, 0)
    m1, m2, m3, m4, m5, m6, m7, m8, m9 = m
    n1, n2, n3, n4, n5, n6, n7, n8, n9 = n
    return 0.5 * np.einsum("i...,i...->...",
        m/(m1*(m5*m9-m6*m8) + m4*(m3*m8-m2*m9) + m7*(m2*m6-m3*m5))
        + n/(n1*(n5*n9-n6*n8) + n4*(n3*n8-n2*n9) + n7*(n2*n6-n3*n5)),
        [m5*n9-m6*n8, m6*n7-m4*n9, m4*n8-m5*n7, n3*m8-n2*m9, n1*m9-n3*m7,
         n2*m7-n1*m8, m2*n6-m3*n5, m3*n4-m1*n6, m1*n5-m2*n4])

def dist_np(m, n):
    return 0.5 * (np.linalg.inv(m)@n + np.linalg.inv(n)@m)[..., np.arange(3), np.arange(3)].sum(axis=-1)

for i in range(3):
    A = np.random.random((100,1,3,3))
    B = np.random.random((1,100,3,3))
    print(np.allclose(dist(A, B), dist_np(A, B)))
    print('pp     ', timeit('f(A,B)', number=100, globals={'f':dist, 'A':A, 'B':B}))
    print('numpy  ', timeit('f(A,B)', number=100, globals={'f':dist_np, 'A':A, 'B':B}))

prints:

True
pp      0.14652886800467968
numpy   1.5294789629988372
True
pp      0.1482033939100802
numpy   1.6455406049499288
True
pp      0.1279512889450416
numpy   1.370200254023075