Cálculo de matriz de dispersión ponderada rápida

En esta pregunta hace seis meses , jez fue lo suficientemente agradable como para ayudarme a encontrar una aproximación rápida para el producto externo de las diferencias de filas, 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 un Análisis Discriminante de Fisher Local, donde cada producto externo está 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, por lo que puedo decir, no es fácil hacer el cambio rápido.

El álgebra lineal definitivamente no es mi fuerte, no soy bueno para idear este tipo de cosas. ¿Cuál es una forma rápida de calcular la sum ponderada de productos externos de diferencia de filas por pares?

Aquí hay una manera de vectorizar el cálculo que ha especificado. Si haces un montón de este tipo de cosas, entonces puede valer la pena aprender a usar “numpy.tensordot”. Multiplica todos los elementos de acuerdo con la emisión de números estándar, y luego los sum a 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 bash (fbetter) hizo una gran matriz temporal de tamaño NxNxN. El segundo bash (fbest) nunca hace nada más grande que NxN. Esto funciona bastante bien hasta N ~ 1000.

Una prueba de tiempo

Además, el código se ejecuta más rápido cuando la matriz de salida es más pequeña. introduzca la descripción de la imagen aquí

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

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