Create sparse circulant matrix in python
<p>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 <code>[i,i+1], [i,i+2], [i,i+N-2], [i,i+N-1]</code>, where I have assumed periodic boundary conditions for the indices (i.e. <code>[10^5,10^5]=[0,0], [10^5+1,10^5+1]=[1,1]</code> and so on). I looked at the scipy sparse matrices documentation but I am quite confused (I am new to Python). </p>

<p>I can create the matrix with numpy</p>

<pre><code>import numpy as np

def Bc(i, boundary):
    """(int, int) -&gt; int

    Checks boundary conditions on index
    """
    if i &gt; boundary - 1:
        return i - boundary
    elif i &lt; 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] 
</code></pre>

<p>However, this is quite slow and for large <code>N</code> 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.</p>

<p>I know how to do it in Mathematica, where one can use SparseArray and index patterns - is there something similar here?</p>

<p>To create a <em>dense</em> circulant matrix, you can use <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.circulant.html" rel="nofollow"><code>scipy.linalg.circulant</code></a>.  For example,</p>

<pre><code>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.    ]])
</code></pre>

<p>As you point out, for large <code>N</code>, that requires a lot of memory and most of the values are zero.  To create a scipy sparse matrix, you can use <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.diags.html" rel="nofollow"><code>scipy.sparse.diags</code></a>.  We have to create offsets (and corresponding values) for the diagonals above and below the main diagonal:</p>

<pre><code>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.    ]])
</code></pre>

<p>The matrix is stored in the "diagonal" format:</p>

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

<p>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:</p>

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

python numpy matrix scipy sparse-matrix