python - 从python中的矩阵向量化所有组合的角度计算
问题描述
请参阅下面的编辑
我有一个相当大的矩阵x
(3xn,n >> 1000),其中关于每列关系的信息有限。从这个矩阵x
中,我需要找到两列的最大角度和对应的索引。目前我正在使用两个for
循环,这需要很长时间,请参见下面的代码。
我知道我可以向量化向量的范数,例如 via np.linalg.norm(x, axis=0)
。另外,我发现了一些关于如何矢量化点积的信息[这里]。但是,这只会计算每列之间的点积,而不是列 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]
解决方案
原来,我正在寻找的实际上是[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
在某种程度上足够聪明,可以检测到冗余并自动优化它。
推荐阅读
- xml - 如何批量下载页面源
- python - 隐藏使用 TensorFlow 后端。TensorFlow 2.0 Alpha 中的消息
- c# - 如何为 azure function v1 编写单元测试
- c# - 遍历 DataSnapShot.Children 会停止代码的执行(Unity、Firebase)
- javascript - 使用 Apify 抓取多个页面
- python - 在我的数据框上使用 pivot_table 的问题
- list - Prolog 中的谓词“浅”或“深”是什么意思?
- javascript - React Native 数字输入掩码始终带有两位小数
- ruby-on-rails - 无法更新验证方法中的属性 - Rails
- go - 结构中的匿名映射列表的“复合文字中缺少类型”