Python: эффективно рассчитать среднее значение недиагональных элементов для фильтра Кадзова

Question

В настоящее время я использую фильтр Гадзоу в Python.

Чтобы поместить в некоторый контекст. Вы начинаете с одномерного массива (давайте возьмем диапазон (10) в качестве примера) и строим из него ганкелеподобную матрицу, такую ​​как

H= [[0, 1, 2, 3, 4, 5],
    [1, 2, 3, 4, 5, 6],
    [2, 3, 4, 5, 6, 7],
    [3, 4, 5, 6, 7, 8],
    [4, 5, 6, 7, 8, 9]])

После этого вы выполняете некоторую линейную алгебру с этой матрицей, которая не является проблемой Затем следует самый трудоемкий этап, который представляет собой проблему усреднения.

В новой матрице B вы усредняете элементы полученной матрицы. В первом ряду вы усредняете все элементы по пути, заданному точными значениями в H. Так что-то вроде выключенных диагоналей, но идущих сверху справа внизу слева. Во втором срезе вы игнорируете первую строку и т. Д.

Матрица $ H $ будет инвариантной на этом этапе анализа, но, например, Матрица

1 2 2 1
1 1 1 1
1 1 1 1

станет

1 1.5 1.33 1
1 1   1    1
1 1   1    1

Хорошо, я надеюсь, вы понимаете проблему. Мой (рабочий, но неэффективный) код

def av_diag(A,i,j):
    dim = A.shape
    # get the "borders" of A
    lim = min((dim[0]-i,j+1))
    # calculate the mean
    return np.mean([A[i+it,j-it] for it in range(lim)])

def avHankel(A):
    # get the mean for all elements by nested list comprehension
    return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))]) 

Это займет некоторое время для моих данных, содержащих 2048 точек данных, что приведет к матрице 1024x1023.

И я был бы рад за возможные трюки, чтобы ускорить это.

Спасибо

Answer 1

Вы можете свернуть входную матрицу с фильтром , чтобы ускорить ваш код. Матрица filter может быть определена таким образом, чтобы на каждом этапе свертки она извлекала только антидиагоналы в заданных координатах. По сути, ваша матрица filter является просто матрицей антиидентификации. Наконец, поскольку свертка будет суммировать только элементы антидиагоналей, вы должны разделить выходные данные на правильное количество выборок, чтобы получить среднее значение:

import numpy as np
from scipy.signal import fftconvolve
from time import time


def av_diag(A,i,j):
    dim = A.shape
    lim = min((dim[0]-i,j+1))
    return np.mean([A[i+it,j-it] for it in range(lim)])

def avHankel(A):
    return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))])


def fast_avHankel(A):
    m, n = A.shape
    filt = np.eye(m)[:,::-1]
    Apad = np.pad(A, ((0, m-1), (m-1, 0)), mode = "constant", constant_values = 0)
    Asum = fftconvolve(Apad, filt, mode = "valid")
    Adiv = np.array([ [ min(m-i, j+1) for j in range(n) ] for i in range(m) ])
    return Asum / Adiv


if __name__ == "__main__":
    A = np.random.rand(500, 500)

    starttime = time()
    Hold = avHankel(A)
    print(time() - starttime)    # 10.6 seconds on a laptop

    starttime = time()
    Hnew = fast_avHankel(A)
    print(time() - starttime)    # 0.26 seconds on a laptop