Как использовать диагональные матрицы матрицы/массива эффективности на основе остальной части столбца (избегая цикла?)

Question

Я пытаюсь закодировать решение, которое принимает функцию, и устанавливает центральные диагональные значения (то есть вдоль A [0,0] , A [1,1], ...., A [N, N]) на основе суммирования значений в столбце ниже диагональной ячейки.

Пример:

A = np.array([[0, 0, 0], 
       [3, 0, 0], 
       [4, 2, 0]]) 

B = function_to_be_built(A) 

B = np.array([[7, 0, 0], 
       [3, 2, 0], 
       [4, 2, 0]]) 

, где можно увидеть, что добавление операторов для ячейки [0,0] = 3 + 4 = 7 и клетки [1,1], 2 = 2.

Конечно, это можно сделать относительно легко, используя петлю for, но код предназначен для центрального ядра, поэтому эффективность имеет решающее значение. Я чувствую, что должен быть способ достичь этого эффективно с помощью numpy .... что-то вроде использования np.tril, поэтому выберите только значения ниже диагонали?

Может ли кто-нибудь помочь мне добраться до решения, пожалуйста?

Спасибо.

Answer 1

Если вы уверены, что диагональ и верхняя диагональ заполняются нулями, вы можете сделать:

In [43]: B = A + np.diag(A.sum(axis=0)) 

In [44]: B 
Out[44]: 
array([[7, 0, 0], 
     [3, 2, 0], 
     [4, 2, 0]]) 

Если вы не можете гарантировать, что эти области равны нулю, то вы можете сделать:

In [62]: B = A + np.diag(np.tril(A, -1).sum(axis=0)) 

In [63]: B 
Out[63]: 
array([[7, 0, 0], 
     [3, 2, 0], 
     [4, 2, 0]]) 

Answer 2

Вы можете взять нижнюю треугольную часть A, обобщать ее столбцами, преобразовать в диагональную матрицу, и добавить обратно в A:

A + np.eye(A.shape[0]) * np.tril(A).sum(axis=0) 

Answer 3

Вот решение с использованием np.add.at

y, x = np.tril_indices(len(A), -1) 
np.add.at(A, [x,x], A[y,x]) 

Это добавляет ниже диагональные суммы столбцов в месте на соответствующих диагональные элементы. Если они не гарантированно равны нулю, сделайте a

A[np.arange(len(A)), np.arange(len(A))] = 0 

до.

Обратите внимание, что в любом случае этот метод не требует, чтобы верхний треугольник был равен нулю.

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

# compute once 
flatinds = np.ravel_multi_index((y, x), A.shape) 
flatdiag = np.ravel_multi_index((x, x), A.shape) 
flshdiag = (len(A)+1) * np.arange(len(A)) 
# use every iteration 
Afl = A.ravel() 
Afl[flshdiag] = 0 # only if necessary 
np.add.at(Afl, flatdiag, Afl[flatinds]) 

Не смог устоять перед действительно быстрым. Он использует линейную индексацию и np.add.reduceat. Так как большинство сумм относятся к относительно длинным участкам столбца, он рассчитывается сделать полный (нелазный) транспонирование в начале и в конце вычисления. При N> = 5 он последовательно превосходит всех остальных. Если вы можете держать матрицы, участвующие в столбцам (Fortran) порядка, два транспонированные даже может быть сохранен и в этом случае это быстрее, чем все остальные, начиная с N == 3.

# precompute this (A.shape == (N, N)) 
finds = np.c_[N * np.arange(N), 1 + (N+1) * np.arange(N)].ravel()[1:-2].copy() 

def PPt(A): 
    A = A.T.copy() 
    Af = A.ravel() 
    Af[:-1:N+1] = np.add.reduceat(Af[:-N], finds)[::2] 
    Af[-1] = 0 
    return A.T.copy()