快速加权散射矩阵计算 [英] Fast weighted scatter matrix calculation

问题描述

六个月前的问题，jez是足够好，可以帮助我快速估算出行差异的外部乘积，即:

In this question six months ago, jez was nice enough to help me come up with a fast approximation for the outer product of row differences, i.e.:

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具有有关该对的局部性的信息，因此新行是:

That worked for finding a scatter matrix calculation for a form of Fisher Discriminant Analysis. But now I am trying to do Local Fisher Discriminant Analysis, where each outer product is weighted by a matrix A that has information about the pair's locality, so the new line is:

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

不幸的是，计算上一个答案中未加权的散布矩阵的快速方法对此不起作用，据我所知，进行快速更改并不容易.

Unfortunately, the quick way to calculate the unweighted scatter matrix presented in the previous answer doesn't work for this, and as far as I can tell it's not easy to make the quick change.

线性代数绝对不是我的强项，我也不擅长提出这类问题.计算成对的行差外积加权和的快速方法是什么?

Linear Algebra is definitely not my strong suit, I'm not good at coming up with these kinds of things. What is a fast way to calculate the weighted sum of pairwise row-difference outer products?

推荐答案

这里是一种矢量化您指定的计算的方法.如果您做很多这样的事情，那么值得学习如何使用"numpy.tensordot".它将根据标准的numpy广播将所有元素相乘，然后将其与通过kwrd轴"给出的两对轴求和.

Here is a way to vectorize the calculation you specified. If you do a lot of this kind of thing, then it may be worth learning how to use, "numpy.tensordot". It multiplies all elements according to standard numpy broadcasting, and then it sums over the pairs of axes given with the 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))

我的第一次尝试(更好)制作了一个大小为NxNxN的大型临时数组.第二次尝试(最好)永远不会使结果大于NxN.直到N〜1000都可以正常工作.

My first attempt (fbetter) made a large temporary array of size NxNxN. The second attempt (fbest) never makes anything bigger than NxN. This works pretty good up to N~1000.

此外，当输出数组较小时，代码将运行得更快.

Also, the code runs faster when the output array is smaller.

我安装了MKL，因此对tensordot的调用非常快并且可以并行运行.

I have MKL installed so the calls to tensordot are really fast and run in parallel.

感谢您的提问.这是一个很好的练习，并提醒我避免制作大型临时数组有多么重要.

Thanks for the question. This was a nice exercise and reminded me how important it is to avoid making large temporary arrays.

这篇关于快速加权散射矩阵计算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！