transpuesta - resta de matrices python

Cálculo matricial de dispersión ponderado rápido

Question

En esta pregunta hace seis meses , Jez fue lo suficientemente amable como para ayudarme a llegar a una aproximación rápida para el producto externo de las diferencias de fila, es decir:

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)

Eso funcionó para encontrar un cálculo de matriz de dispersión para una forma de análisis discriminante de Fisher. Pero ahora estoy tratando de hacer el Análisis Discriminante Local de Fisher, donde cada producto externo es ponderado por una matriz A que tiene información sobre la localidad del par, por lo que la nueva línea es:

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

Desafortunadamente, la forma rápida de calcular la matriz de dispersión no ponderada presentada en la respuesta anterior no funciona para esto, y hasta donde puedo decir, no es fácil hacer el cambio rápido.

Linear Algebra definitivamente no es mi fuerte, no soy bueno para hacer este tipo de cosas. ¿Cuál es una forma rápida de calcular la suma ponderada de los productos exteriores de diferencia de fila pairwise?

Answer 1

Aquí hay una forma de vectorizar el cálculo que ha especificado. Si haces mucho de este tipo de cosas, entonces valdría la pena aprender a usar, "numpy.tensordot". Multiplica todos los elementos de acuerdo con la difusión numpy estándar, y luego suma los pares de ejes dados con el kwrd, "ejes".

Aquí está el código:

# 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))

Mi primer intento (fbetter) hizo una gran matriz temporal de tamaño NxNxN. El segundo intento (fbest) nunca hace algo más grande que NxN. Esto funciona bastante bien hasta N ~ 1000.

Además, el código se ejecuta más rápido cuando la matriz de salida es más pequeña.

Tengo MKL instalado, por lo que las llamadas a tensordot son realmente rápidas y se ejecutan en paralelo.

Gracias por la pregunta. Este fue un ejercicio agradable y me recordó lo importante que es evitar hacer grandes arreglos temporales.