python - 向量化python中矩阵中所有组合的角度计算

Question

查看下面的编辑

我有一个相当大的矩阵 x (3xn, n >> 1000)，其中关于每一列的关系的信息有限。 从这个矩阵 x 中，我需要找到最大的角度和两列的对应索引。 目前我正在使用两个 for 循环，这需要很长时间，请参见下面的代码。

我知道我可以向量化向量的范数，例如通过 np.linalg.norm(x, axis=0)。 另外，我找到了一些关于如何向量化点积的信息 [here] . 但是，这只会计算每列之间的点积，而不是第 1 列和第 3 列之间的点积。

有没有办法优化甚至向量化问题？

import numpy as np

# Calculate the max angle between two vectors
len_x = 1000
x = np.random.rand(3,len_x) # 3xn matrix

ang_max = 0 # maximum angle
n_max   = 0 # n-index of maximum angle
m_max   = 0 # m-index of maximum angle

for n in range(len_x - 1):
     vn = np.array([x[0][n], \
                    x[1][n], \
                    x[2][n]])
     for m in range(n + 1, len_x):
         vm = np.array([x[0][m], \
                        x[1][m], \
                        x[2][m]])

         ang = np.arccos(np.dot(vn,vm)/(np.linalg.norm(vn)*np.linalg.norm(vm)))

         if ang > ang_max:
             ang_max = ang
             n_max = n
             m_max = m

编辑:

使用 np.einsum 我能够将内部循环矢量化。 这将执行时间缩短了 44 倍，至少在我的测试用例中是这样(从 31.98 秒减少到 0.725 秒)。 这是一个显着的改进，但我认为对外部循环进行矢量化可能会以类似的因素提供另一个加速。 因此，我不会将问题标记为已回答，因为外循环也未矢量化。

import numpy as np

# Calculate the max angle between two vectors
len_x = 1000
x = np.random.rand(3,len_x) # 3xn matrix

ang_max = 0 # maximum angle
n_max   = 0 # n-index of maximum angle
m_max   = 0 # m-index of maximum angle

vm      = np.array(x)
vm_norm = np.linalg.norm(vm, axis=0)

for n in range(len_x - 1):
    vn = np.array([float(x[n]) - xs, \
                   float(y[n]) - ys, \
                   float(z[n]) - zs])
    ang = np.arccos(np.divide(np.einsum('ij,ij->j', vn[:,None], vm[:,:,0]),(np.linalg.norm(vn)*vm_norm)[:,0]))
    maxang = max(ang)

    if maxang > ang_max:
        ang_max = maxang
        n_max = n
        m_max = np.where(ang == maxang)[0][0]

Answer 1

原来，我要找的其实叫[Gramian Matrix] ，一个矩阵，其元素是内积的所有可能组合。在我的例子中，转换为以下代码:

import numpy as np

# Calculate the max angle between two vectors
len_x = 1000
x = np.random.rand(3,len_x) # 3xn matrix

ang_max = 0 # maximum angle
n_max   = 0 # n-index of maximum angle
m_max   = 0 # m-index of maximum angle

vm      = np.array(x)
vm_norm = np.linalg.norm(vm, axis=0)

ang = np.arccos(np.divide(vm.T.dot(vm), (vm_norm.dot(vm_norm))))

ang_max = np.amax(ang)
argw    = np.argwhere(ang == np.amax(ang))
n_max   = argw[0, 0]
m_max   = argw[0, 1]

该代码让我的性能又提高了大约 4 倍，将执行时间缩短到约 0.2 秒。

现在，由于内积是对称的，vm.dot(vm) 的一半计算是多余的。我不确定如何在不增加太多复杂性的情况下进一步优化这种边缘情况。但是，我有一种强烈的感觉，即 numpy 在某种程度上足够聪明，可以检测冗余并自动优化它。