快速加权散射矩阵计算

Question

In this question six months ago，jez很好地帮我想出了行差的外积的快速近似值，即：

K = np.zeros((len(X), len(X)))
for i, Xi in enumerate(X):
  for j, Xj in enumerate(X):
    dij = Xi - Xj
    K += np.outer(dij, dij)

这是为了寻找散度矩阵计算的一种形式的Fisher判别分析。但现在我要做局部Fisher判别分析，其中每个外积由矩阵a加权，矩阵a包含关于对的局部性的信息，所以新行是：

K += A[i][j] * np.outer(dij, dij)

不幸的是，快速计算前一个答案中给出的未加权散射矩阵的方法不适用于此，而且据我所知，快速更改并不容易。在

线性代数绝对不是我的强项，我不擅长想出这些东西。什么是一种快速计算成对行差外积加权和的方法？在

Answer 1

下面是一种将指定计算矢量化的方法。如果你做了很多这样的事情，那么学习如何使用它可能是值得的。”纽比·坦索多特". 它根据标准的numpy广播对所有元素进行乘法，然后用kwrd“axes”对给定的轴对求和。在

代码如下：

# Imports
import numpy as np
from numpy.random import random

# Original calculation for testing purposes 
def ftrue(A, X):
  ""
  K = np.zeros((len(X), len(X)))
  KA_true = np.zeros((len(X), len(X)))

  for i, Xi in enumerate(X):
    for j, Xj in enumerate(X):
      dij = Xi - Xj
      K += np.outer(dij, dij)
      KA_true += A[i, j] * np.outer(dij, dij) 
  return ftrue

# Better: No Python loops. But, makes a large temporary array.
def fbetter(A, X):
  ""
  c = X[:, None, :] - X[None, :, :]
  b = A[:, :, None] * c           # ! BAD ! temporary array size N**3
  KA_better = np.tensordot(b, c, axes = [(0,1),(0,1)])
  return KA_better

# Best way: No Python for loops. No large temporary arrays
def fbest(A, X):
  ""
  KA_best = np.tensordot(A.sum(1)[:,None] * X, X, axes=[(0,), (0,)])
  KA_best += np.tensordot(A.sum(0)[:,None] * X, X, axes=[(0,), (0,)])
  KA_best -= np.tensordot(np.dot(A, X), X, axes=[(0,), (0,)])
  KA_best -= np.tensordot(X, np.dot(A, X), axes=[(0,), (0,)])
  return KA_best


# Test script
if __name__ == "__main__":

  # Parameters for the computation 
  N = 250
  X = random((N, N))
  A = random((N, N))

  # Print the error
  KA_better = fbetter(A, X)
  KA_best = fbest(A, X)

  # Test against true if array size isn't too big
  if N<100:
    KA_true = ftrue(A, X)
    err = abs(KA_better - KA_true).mean()
    msg = "Mean absolute difference (better): {}."
    print(msg.format(err))

  # Test best against better
  err = abs(KA_best - KA_better).mean()
  msg = "Mean absolute difference (best): {}."
  print(msg.format(err))

我的第一次尝试（fbetter）创建了一个大小为NxNxN的大型临时数组。第二次尝试（fbest）不会产生比NxN更大的结果。在N~1000的范围内，这种方法非常有效。在

另外，当输出数组较小时，代码运行得更快。

我已经安装了MKL，所以对tensordot的调用非常快，并且并行运行。在

谢谢你的问题。这是一个很好的练习，它提醒我避免生成大型临时数组是多么重要。在