CODEWITHSUNDEEP
pythonnumpyscipysparse-matrix
Ohm
Ohm

Reputation: 2442

Buildin a sparse 2D laplacian matrix using SciPy modules

I need to construct the 2D laplacian which looks like this:

2D_Laplacian

, where

enter image description here

, and I is the identity matrix. So far, I have done it using the diags method of scipy, but I wonder whether there is a smarter way to do it using the block_diag method. Has anyone tried to build the 2D laplacian with this method?

My current method to create this is by this function:

from scipy.sparse import diags
# Defining the size of the matrix
nx = 3
ny = 3
N  = nx*ny
main_diag = [-4.0 for i in xrange(N)]
side_diag = []
for i in xrange(1,N):
    if i%4 == 0:
        side_diag.append(0)
    else:
        side_diag.append(1)
up_down_diag = [1.0 for i in xrange(N-4)]
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

Upvotes: 3

Views: 6105

Answers (2)

Shiv
Shiv

Reputation: 31

N-dimensional laplacians can be expressed as Kronecker product of 1D laplacians:

import scipy.sparse as sp
def laplacian2D(N):
    diag=np.ones([N*N])
    mat=sp.spdiags([diag,-2*diag,diag],[-1,0,1],N,N)
    I=sp.eye(N)
    return sp.kron(I,mat,format='csr')+sp.kron(mat,I)

Connectivity pattern of 2D Laplacian - img

Upvotes: 3

hpaulj
hpaulj

Reputation: 231385

I replaced your use of lists with arrays:

import numpy as np
from scipy import sparse

nx, ny = 3, 3
N  = nx*ny
main_diag = np.ones(N)*-4.0
side_diag = np.ones(N-1)
side_diag[np.arange(1,N)%4==0] = 0
up_down_diag = np.ones(N-3)
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = sparse.diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

producing

[[-4.  1.  0.  1.  0.  0.  0.  0.  0.]
 [ 1. -4.  1.  0.  1.  0.  0.  0.  0.]
 [ 0.  1. -4.  1.  0.  1.  0.  0.  0.]
 [ 1.  0.  1. -4.  0.  0.  1.  0.  0.]
 [ 0.  1.  0.  0. -4.  1.  0.  1.  0.]
 [ 0.  0.  1.  0.  1. -4.  1.  0.  1.]
 [ 0.  0.  0.  1.  0.  1. -4.  1.  0.]
 [ 0.  0.  0.  0.  1.  0.  1. -4.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0. -4.]]

Is the [1 1 1 0 1 1 1 0] pattern for the side diagonal right?

For a small example like this it may run the same speed, but with large dimensions using arrays rather than lists should be faster - plus it is more consistent with numpy code underlying sparse.

With uniform diagonals like this diags looks pretty good.

I've only played with the block format for one other SO question. https://stackoverflow.com/a/34124377/901925

coo is good for matricies composed of overlapping smaller ones, such as in finite element stiffness. But recasting your diagonals into coo is tedious.

For what it's worth, sparse.diags uses dia_matrix, having converted the list of diagonals into the dia data matrix. You could look at that, but its layout isn't as obvious. To make a csr matrix, diags converts this dia format to coo, and then on to csr. But normally you shouldn't worry about all these conversions. Use the format that makes most sense in your problem, and let sparse handle the conversion details.

If you want to explore block format more, you need to outline how your problem can be viewed as blocks.

Upvotes: 4

Related Questions