numpy中更快的矩阵计算答案

【问题标题】：Faster matrix calculation in numpynumpy中更快的矩阵计算
【发布时间】：2020-11-26 18:03:50
【问题描述】：

在给定一个 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”方式来代替循环......

编辑：虽然所有答案都有效并且比我的幼稚实现要快得多，但我选择了我基准测试最快的答案。谢谢！

【问题讨论】：

X的形状是什么？是向量还是二维数组？
@Jalo X 是一个向量，M 是一个二维数组
为什么不对所有G[i, j] 使用公式？似乎对角线被简单地分配给 X 并且没有计算下三角形。
@swag2198 公式是对称的，所以我只需要矩阵的上三角部分
你考虑过 numba 还是 cython？

标签： python numpy math

【解决方案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 的对角元素无关紧要并且它们包含一些零，只需添加一些偏移量以避免divide by zero 错误，例如：

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)

【讨论】：

这看起来很棒，而且速度更快！但是如果 M 的值为零，它会给出一个错误，这对我来说是一些对角线的情况（这就是我使用 if-else 子句的原因，很抱歉没有在问题中写出来......）
好的，如果某些M[i, j] 为零，那么规则是什么？
这是个好问题；我相信这并不重要，因为这些值永远不会被使用。我只是在我当前的代码中使用 numpy 的 nan_to_num 将它们设置为零。
我在回答中添加了一些注释。您可以偏移对角线以进行相同的计算，然后如果您愿意，稍后将G 的对角线分配给X。这基本上为您提供了原始代码的if i == j 部分。

【解决方案2】：

首先，你的功能不是你的方程式。作为这一行

 mid = (xi2 - xj2)/mij

应该是

 mid = (xi - xj)/mij

其次，我用 numpy 生成你的方程。

生成测试数据

test_m = np.array(
    [
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
    ]
)
test_x = np.array([5, 6, 7, 8, 9])

构建函数

def solve(m, x):
    x_size = x.shape[0]
    x = x.reshape(1, -1)
    reshaped_x = x.reshape(-1, 1)
    result = np.sqrt(
        m ** 2
        + ((reshaped_x - x) / m) ** 2
        + 2 * np.repeat(reshaped_x, x_size, axis=1) ** 2
        + 2 * np.repeat(x, x_size, axis=0) ** 2
    ) / 2
    return result

运行

print(solve(test_m, test_x))

实际上，结果部分可以像这样简单：

    result = np.sqrt(
        m ** 2
        + ((reshaped_x - x) / m) ** 2
        + 2 * reshaped_x ** 2
        + 2 * x ** 2
    ) / 2

【讨论】：

【解决方案3】：

用 googles colab 测试：

导入数字将 numpy 导入为 np 导入数学

# your implementation

def bench_1(n):
    #M, X are given as numpy arrays
    G = np.zeros((n,n))
    M = np.random.rand(n, n)
    X = np.random.rand(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
    return G

  %%timeit
  n = 1000
  bench_1(n)

  1 loop, best of 3: 1.61 s per loop

使用Numba编译函数：

@numba.jit(nopython=True, parallel=True)
def bench_2(n):
  #M, X are given as numpy arrays
  G = np.zeros((n,n))
  M = np.random.rand(n, n)
  X = np.random.rand(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
  return G

%%timeit
n = 1000
bench_2(n)

The slowest run took 88.13 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 9.8 ms per loop

【讨论】：