janmr blog

The 1D Laplace Matrix Using NumPy and SciPy

Let us consider the n×nn \times n matrix representation of the Laplace operator in one dimension, as introduced in the post Finite Difference Discretization of the 1D Laplace Operator:

[2100121001200002]\begin{bmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \cdots & 0 \\ 0 & 1 & -2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -2 \end{bmatrix}

(we ignore the 1/h21/h^2 factor for simplicity).

If we want to create this matrix using NumPy, we can use:

L = np.eye(n, k=-1) - 2 * np.eye(n) + np.eye(n, k=1)

where the k argument specifies the diagonal offset.

Note that this matrix representation is dense, meaning that each of the n2n^2 elements is stored explicitly, and in case of a matrix-vector multiplication, all n2n^2 elements are accessed. In the case of the Laplace matrix, this is highly inefficient, as roughly 3n3n elements are non-zero.

We therefore turn to SciPy's sparse matrix representations and create the Laplace matrix using scipy.sparse.eye:

L =     scipy.sparse.eye(n, k=-1) \
  - 2 * scipy.sparse.eye(n) \
  +     scipy.sparse.eye(n, k=1)

On my machine this creates a sparse matrix in the Compressed Sparse Row (CSR) format (stores an array of column indeces and value items for each row of the matrix). Note that SciPy supports several other sparse matrix formats and the output format is implementation-dependent, unless explicitly specified.

Another way to create the Laplace matrix is to use scipy.sparse.diags_array:

L = scipy.sparse.diags_array(
        [1, -2, 1],
        offsets=[-1, 0, 1],
        shape=(n, n)

On my machine this creates a sparse matrix in the diagonal storage format (stores one or more arrays of values relative to the main diagonal).

Which matrix format should be used for the Laplace matrix? As mentioned above, the dense format contains mostly zeros and is therefore inefficient, at least for large nn. The sparse formats are much more (space) efficient, but sparse matrices do come with some extra overhead. Furthermore, the most useful sparse format depends on the specific operation(s) to be performed.

Figure 1 below shows times spent on matrix-matrix multiplication, where the left-hand operand was a dense or sparse Laplace matrix and the right-hand operand was a dense matrix containing 100 columns of random values.

Dense vs sparse laplace matrix-matrix multiplication for small n
Figure 1. Dense vs sparse matrix-matrix multiplication for small n.

The figure shows that dense matrix-matrix multiplication is faster than sparse matrix-matrix multiplication for small nn. Both sparse formats will eventually be faster, but the break-even points depend on the specific machine and library versions used.

Figure 2 compares the two sparse matrix representations for larger nn. It is here clear that the CSR format is significantly faster than the diagonal storage format for matrix-matrix multiplication.

Laplace matrix-matrix multiplication using two different sparse representations
Figure 2. Laplace matrix-matrix multiplication using two different sparse representations.
Feel free to leave any question, correction or comment in this Mastodon thread.