在python中创建稀疏循环矩阵 [英] Create sparse circulant matrix in python
问题描述
我想在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]等).我查看了稀疏的稀疏矩阵文档，但是我很困惑(我是Python的新手).
I want to create a large (say 10^5 x 10^5) sparse circulant matrix in Python. It has 4 elements per row at positions [i,i+1], [i,i+2], [i,i+N-2], [i,i+N-1], where I have assumed periodic boundary conditions for the indices (i.e. [10^5,10^5]=[0,0], [10^5+1,10^5+1]=[1,1] and so on). I looked at the scipy sparse matrices documentation but I am quite confused (I am new to Python).
我可以用numpy创建矩阵
I can create the matrix with 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创建和将其转换为稀疏矩阵并直接转到后者.
However, this is quite slow and for large N uses a lot of memory, so I want to avoid the creation with numpy and the converting to a sparse matrix and go directly to the latter.
我知道如何在Mathematica中做到这一点，在那里人们可以使用SparseArray和索引模式-这里有类似的东西吗?
I know how to do it in Mathematica, where one can use SparseArray and index patterns - is there something similar here?
推荐答案
To create a dense circulant matrix, you can use scipy.linalg.circulant. For example,
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.sparse.diags .我们必须为主要对角线上方和下方的对角线创建偏移量(和相应的值):
As you point out, for large N, that requires a lot of memory and most of the values are zero. To create a scipy sparse matrix, you can use scipy.sparse.diags. We have to create offsets (and corresponding values) for the diagonals above and below the main diagonal:
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. ]])
矩阵以对角线"格式存储:
The matrix is stored in the "diagonal" format:
In [226]: a
Out[226]:
<7x7 sparse matrix of type '<class 'numpy.float64'>'
with 28 stored elements (8 diagonals) in DIAgonal format>
您可以使用稀疏矩阵的转换方法将其转换为其他稀疏格式.例如，以下结果将得出CSR格式的矩阵:
You can use the conversion methods of the sparse matrix to convert it to a different sparse format. For example, the following results in a matrix in CSR format:
In [227]: a.tocsr()
Out[227]:
<7x7 sparse matrix of type '<class 'numpy.float64'>'
with 28 stored elements in Compressed Sparse Row format>
这篇关于在python中创建稀疏循环矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
