Быстрый способ преобразования верхней треугольной матрицы в симметричную матрицу

Question

У меня есть верхнетреугольная матрица значений np.float64, например:

array([[ 1.,  2.,  3.,  4.],
       [ 0.,  5.,  6.,  7.],
       [ 0.,  0.,  8.,  9.],
       [ 0.,  0.,  0., 10.]])

Я хотел бы преобразовать это в соответствующую симметричную матрицу, например:

array([[ 1.,  2.,  3.,  4.],
       [ 2.,  5.,  6.,  7.],
       [ 3.,  6.,  8.,  9.],
       [ 4.,  7.,  9., 10.]])

Преобразование может быть сделано на месте или в виде новой матрицы. Я бы хотел, чтобы это было как можно быстрее. Как я могу сделать это быстро?

Answer 1

np.where кажется довольно быстрым в неуместном сценарии без кэширования:

np.where(ut,ut,ut.T)

На моем ноутбуке:

timeit(lambda:np.where(ut,ut,ut.T))
# 1.909718865994364

Если у вас установлен Pythran, выможет ускорить это в 3 раза с почти нулевым усилием. Но учтите, что, насколько мне известно, pythran (в настоящее время) понимает только непрерывные массивы.

file <upp2sym.py>, компилируется с pythran -O3 upp2sym.py

import numpy as np

#pythran export upp2sym(float[:,:])

def upp2sym(a):
    return np.where(a,a,a.T)

Время:

from upp2sym import *

timeit(lambda:upp2sym(ut))
# 0.5760842661838979

Это почти так же быстро, как цикл:

#pythran export upp2sym_loop(float[:,:])

def upp2sym_loop(a):
    out = np.empty_like(a)
    for i in range(len(a)):
        out[i,i] = a[i,i]
        for j in range(i):
            out[i,j] = out[j,i] = a[j,i]
    return out

Время:

timeit(lambda:upp2sym_loop(ut))
# 0.4794591029640287

Мы также можем сделать это на месте:

#pythran export upp2sym_inplace(float[:,:])

def upp2sym_inplace(a):
    for i in range(len(a)):
        for j in range(i):
            a[i,j] = a[j,i]

Время

timeit(lambda:upp2sym_inplace(ut))
# 0.28711927914991975

Answer 2

Это самая быстрая процедура, которую я нашел до сих пор, которая не использует Cython или JIT, как Numba. На моей машине у меня уходит около 1,6 мкс для обработки массива 4x4 (среднее время по списку массивов 4K4 размером 100 КБ):

inds_cache = {}

def upper_triangular_to_symmetric(ut):
    n = ut.shape[0]
    try:
        inds = inds_cache[n]
    except KeyError:
        inds = np.tri(n, k=-1, dtype=np.bool)
        inds_cache[n] = inds
    ut[inds] = ut.T[inds]

Вот еще несколько вещей, которые я пробовал, но не так быстро:

Вышеприведенный код, но без кеша. Принимает около 8,3 мкс на массив 4x4:

def upper_triangular_to_symmetric(ut):
    n = ut.shape[0]
    inds = np.tri(n, k=-1, dtype=np.bool)
    ut[inds] = ut.T[inds]

Простой вложенный цикл Python. Принимает около 2,5 мкс на массив 4x4:

def upper_triangular_to_symmetric(ut):
    n = ut.shape[0]
    for r in range(1, n):
        for c in range(r):
            ut[r, c] = ut[c, r]

С плавающей запятой с использованием np.triu. Принимает около 11,9 мкс на массив 4x4:

def upper_triangular_to_symmetric(ut):
    ut += np.triu(ut, k=1).T

Cython-версия вложенного цикла Python. Я новичок в Cython, так что это не может быть полностью оптимизировано. Поскольку Cython добавляет эксплуатационные накладные расходы, мне интересно услышать ответы как от Cython, так и от чистого Numpy. Принимает около 0,6 мкс на массив 4x4:

cimport numpy as np
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def upper_triangular_to_symmetric(np.ndarray[np.float64_t, ndim=2] ut):
    cdef int n, r, c
    n = ut.shape[0]
    for r in range(1, n):
        for c in range(r):
            ut[r, c] = ut[c, r]

Answer 3

Вы в основном измеряете накладные расходы на вызов функций при таких крошечных проблемах

Еще один способ сделать это - использовать Numba. Давайте начнем с реализации только для одного (4x4) массива.

Только один массив 4x4

import numpy as np
import numba as nb

@nb.njit()
def sym(A):
    for i in range(A.shape[0]):
        for j in range(A.shape[1]):
            A[j,i]=A[i,j]
    return A


A=np.array([[ 1.,  2.,  3.,  4.],
       [ 0.,  5.,  6.,  7.],
       [ 0.,  0.,  8.,  9.],
       [ 0.,  0.,  0., 10.]])

%timeit sym(A)
#277 ns ± 5.21 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Более крупный пример

@nb.njit(parallel=False)
def sym_3d(A):
    for i in nb.prange(A.shape[0]):
        for j in range(A.shape[1]):
            for k in range(A.shape[2]):
                A[i,k,j]=A[i,j,k]
    return A

A=np.random.rand(1_000_000,4,4)

%timeit sym_3d(A)
#13.8 ms ± 49.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
#13.8 ns per 4x4 submatrix