Быстрый взвешенный расчет матрицы рассеяния

В этом вопросе шесть месяцев назад , 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, где каждый внешний продукт взвешен матрицей A, которая имеет информацию о местонахождении пары, поэтому новая строка:

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

К сожалению, быстрый способ рассчитать невзвешенную матрицу рассеяния, представленную в предыдущем ответе, не работает для этого, и, насколько я могу судить, сделать быстрые изменения непросто.

Линейная алгебра, безусловно, не мой сильный костюм, я не умею придумывать подобные вещи. Что такое быстрый способ вычисления взвешенной суммы парных внешних разностных продуктов?

    One Solution collect form web for “Быстрый взвешенный расчет матрицы рассеяния”

    Вот способ векторизации указанного вами расчета. Если вы так много делаете, то, возможно, стоит изучить, как использовать «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 действительно быстрые и работают параллельно.

    Спасибо за вопрос. Это было отличное упражнение и напомнило мне, насколько важно избегать создания больших временных массивов.

    Python - лучший язык программирования в мире.