python - numpy 中更快的矩阵计算

Question

给定一个 nxn 矩阵 M 和一个 n 向量 X: this paper ，是否有一些更快的变体来计算以下矩阵(来自 ) ？

我目前的计算如下:

#M, X are given as numpy arrays
G = np.zeros((n,n))
for i in range(0,n):
    for j in range(i,n):
        xi = X[i]
        if i == j:
            G[i,j] = abs(xi)
        else:
            xi2 = xi*xi
            xj = X[j]
            xj2 = xj*xj
            mij = M[i,j]
            mid = (xi2 - xj2)/mij
            top =  mij*mij + mid*mid + 2*xi2 + 2*xj2
            G[i,j] = math.sqrt(top)/2

这非常慢，但我怀疑有一种更好的“numpythonic”方法可以代替循环......

编辑:虽然所有答案都有效并且比我天真的实现要快得多，但我选择了我基准测试中最快的一个。谢谢!

Answer 1

实际上非常简单。

import math
import numpy as np

n = 5
M = np.random.rand(n, n)
X = np.random.rand(n)

您的代码和结果:

G = np.zeros((n,n))
for i in range(0,n):
    for j in range(i,n):
        xi = X[i]
        if i == j:
            G[i,j] = abs(xi)
        else:
            xi2 = xi*xi
            xj = X[j]
            xj2 = xj*xj
            mij = M[i,j]
            mid = (xi2 - xj2)/mij
            top =  mij*mij + mid*mid + 2*xi2 + 2*xj2
            G[i,j] = math.sqrt(top)/2

array([[0.77847813, 5.26334534, 0.8794082 , 0.7785694 , 0.95799072],
       [0.        , 0.15662266, 0.88085031, 0.47955479, 0.99219171],
       [0.        , 0.        , 0.87699707, 8.92340836, 1.50053712],
       [0.        , 0.        , 0.        , 0.45608367, 0.95902308],
       [0.        , 0.        , 0.        , 0.        , 0.95774452]])

使用广播:

temp = M**2 + ((X[:, None]**2 - X[None, :]**2) / M)**2 + 2 * (X[:, None]**2) + 2 * (X[None, :]**2)
G = np.sqrt(temp) / 2

array([[0.8284724 , 5.26334534, 0.8794082 , 0.7785694 , 0.95799072],
       [0.89251217, 0.25682736, 0.88085031, 0.47955479, 0.99219171],
       [0.90047282, 1.10306597, 0.95176428, 8.92340836, 1.50053712],
       [0.85131766, 0.47379576, 0.87723514, 0.55013345, 0.95902308],
       [0.9879939 , 1.46462011, 0.99516443, 0.95774481, 1.02135642]])

请注意，您没有直接对对角线元素使用公式，而仅针对 G 的上三角区域进行计算。我只是实现了计算所有 G[i, j] 的公式。

注意:如果M的对角线元素不重要并且它们包含一些零，只需添加一些偏移量以避免被零除 错误如下:

M[np.arange(n), np.arange(n)] += 1e-5

# Do calculation to get G

# Assign diagonal to X
G[np.arange(n), np.arange(n)] = abs(X)