首页 文章

在python中创建稀疏循环矩阵

提问于
浏览
2

我想在Python中创建一个大的(比如10 ^ 5 x 10 ^ 5)稀疏循环矩阵 . 它在位置 [i,i+1], [i,i+2], [i,i+N-2], [i,i+N-1] 处每行有4个元素,其中我假设了索引的周期性边界条件(即 [10^5,10^5]=[0,0], [10^5+1,10^5+1]=[1,1] 等等) . 我查看了scipy稀疏矩阵文档,但我很困惑(我是Python的新手) .

我可以用numpy创建矩阵

import numpy as np

def Bc(i, boundary):
    """(int, int) -> int

    Checks boundary conditions on index
    """
    if i > boundary - 1:
        return i - boundary
    elif i < 0:
        return boundary + i
    else:
        return i

N = 100
diffMat = np.zeros([N, N])
for i in np.arange(0, N, 1):
    diffMat[i, [Bc(i+1, N), Bc(i+2, N), Bc(i+2+(N-5)+1, N), Bc(i+2+(N-5)+2, N)]] = [2.0/3, -1.0/12, 1.0/12, -2.0/3]

然而,这是非常缓慢的,并且对于大 N 使用大量内存,所以我想避免使用numpy创建并转换为稀疏矩阵并直接转到后者 .

我知道如何在Mathematica中做到这一点,在那里可以使用SparseArray和索引模式 - 这里有类似的东西吗?

1 回答

  • 4

    要创建密集的循环矩阵,可以使用scipy.linalg.circulant . 例如,

    In [210]: from scipy.linalg import circulant
    
    In [211]: N = 7
    
    In [212]: vals = np.array([2.0/3, -1.0/12, 1.0/12, -2.0/3])
    
    In [213]: offsets = np.array([1, 2, N-2, N-1])
    
    In [214]: col0 = np.zeros(N)
    
    In [215]: col0[offsets] = -vals
    
    In [216]: c = circulant(col0)
    
    In [217]: c
    Out[217]: 
    array([[ 0.    ,  0.6667, -0.0833,  0.    ,  0.    ,  0.0833, -0.6667],
           [-0.6667,  0.    ,  0.6667, -0.0833,  0.    ,  0.    ,  0.0833],
           [ 0.0833, -0.6667,  0.    ,  0.6667, -0.0833,  0.    ,  0.    ],
           [ 0.    ,  0.0833, -0.6667,  0.    ,  0.6667, -0.0833,  0.    ],
           [ 0.    ,  0.    ,  0.0833, -0.6667,  0.    ,  0.6667, -0.0833],
           [-0.0833,  0.    ,  0.    ,  0.0833, -0.6667,  0.    ,  0.6667],
           [ 0.6667, -0.0833,  0.    ,  0.    ,  0.0833, -0.6667,  0.    ]])
    

    正如您所指出的,对于大 N ,这需要大量内存,并且大多数值为零 . 要创建scipy稀疏矩阵,可以使用scipy.sparse.diags . 我们必须为主对角线上方和下方的对角线创建偏移量(和相应的值):

    In [218]: from scipy import sparse
    
    In [219]: N = 7
    
    In [220]: vals = np.array([2.0/3, -1.0/12, 1.0/12, -2.0/3])
    
    In [221]: offsets = np.array([1, 2, N-2, N-1])
    
    In [222]: dupvals = np.concatenate((vals, vals[::-1]))
    
    In [223]: dupoffsets = np.concatenate((offsets, -offsets))
    
    In [224]: a = sparse.diags(dupvals, dupoffsets, shape=(N, N))
    
    In [225]: a.toarray()
    Out[225]: 
    array([[ 0.    ,  0.6667, -0.0833,  0.    ,  0.    ,  0.0833, -0.6667],
           [-0.6667,  0.    ,  0.6667, -0.0833,  0.    ,  0.    ,  0.0833],
           [ 0.0833, -0.6667,  0.    ,  0.6667, -0.0833,  0.    ,  0.    ],
           [ 0.    ,  0.0833, -0.6667,  0.    ,  0.6667, -0.0833,  0.    ],
           [ 0.    ,  0.    ,  0.0833, -0.6667,  0.    ,  0.6667, -0.0833],
           [-0.0833,  0.    ,  0.    ,  0.0833, -0.6667,  0.    ,  0.6667],
           [ 0.6667, -0.0833,  0.    ,  0.    ,  0.0833, -0.6667,  0.    ]])
    

    矩阵以“对角线”格式存储:

    In [226]: a
    Out[226]: 
    <7x7 sparse matrix of type '<class 'numpy.float64'>'
        with 28 stored elements (8 diagonals) in DIAgonal format>
    

    您可以使用稀疏矩阵的转换方法将其转换为不同的稀疏格式 . 例如,以下结果以CSR格式生成矩阵:

    In [227]: a.tocsr()
    Out[227]: 
    <7x7 sparse matrix of type '<class 'numpy.float64'>'
        with 28 stored elements in Compressed Sparse Row format>
    

相关问题