将上三角矩阵转换为对称矩阵的快速方法

我有一个np.float64值的上三角矩阵,如下所示:

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

我想将其转换为相应的对称矩阵,如下所示:

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

转换可以就地完成,也可以作为新矩阵进行。我希望它尽快。我该如何快速做到这一点?

zhonglin0701 回答:将上三角矩阵转换为对称矩阵的快速方法

这是迄今为止我发现的最快的例程,它不使用Cython或Numba之类的JIT。我在计算机上花费约1.6μs的时间来处理4x4阵列(在100K 4x4阵列列表中的平均时间):

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]

以下是我尝试过的其他一些事情,它们不那么快:

以上代码,但没有缓存。每个4x4阵列大约需要8.3μs:

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

一个普通的Python嵌套循环。每个4x4阵列大约需要2.5μs:

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进行浮点加法。每个4x4阵列大约需要11.9μs:

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

Numba版本的Python嵌套循环。这是我发现的最快的结果(每个4x4阵列约0.4μs),并且最终在生产中使用,至少直到我开始遇到Numba的问题并不得不恢复为纯Python版本:

import numba

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

Cython版本的Python嵌套循环。我是Cython的新手,因此可能无法完全优化。由于Cython会增加运营开销,因此我想听听Cython和纯Numpy答案。每个4x4阵列大约需要0.6μs:

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,r]
,

np.where在无缓存,无缓存的情况下显得相当快:

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

在我的笔记本电脑上:

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

如果您安装了pythran,则可以以几乎零的努力将其加速3倍。但是请注意,据我所知,pythran(目前)仅了解连续数组。

文件<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
,

您主要是在测量此类微小问题上的函数调用开销

另一种方法是使用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.,8.,9.],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
本文链接:https://www.f2er.com/3156842.html

大家都在问