快速加权散点矩阵计算

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.tensordot"。它根据标准 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 的调用非常快并且运行是并行的。

感谢提问。这是一个很好的练习，提醒我避免制作大型临时数组的重要性。

快速加权散点矩阵计算

Fast weighted scatter matrix calculation

python

numpy

matrix

linear-algebra