在python中创建稀疏循环矩阵答案

【问题标题】：Create sparse circulant matrix in python在python中创建稀疏循环矩阵
【发布时间】：2017-01-28 14:06:58
【问题描述】：

我想在 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 sparse matrices 文档，但我很困惑（我是 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 和索引模式 - 这里有类似的东西吗？

【问题讨论】：

标签： python numpy matrix scipy sparse-matrix

【解决方案1】：

要创建一个密集循环矩阵，您可以使用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>

【讨论】：

正是我想要的！我知道scipy.linalg.circulant，但我不知道scipy.sparse.diags。非常感谢！