janmr blog

Let us now consider the discretization of the 2D Laplace operator using finite differences. Recall that the Laplace operator is defined as

$$\Delta u = \nabla \cdot \nabla u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2},$$

for some sufficiently smooth function $u: \mathbb{R}^2 \mapsto \mathbb{R}$.

By following the procedure as in the previous post, we get:

$$\begin{aligned} \frac{\partial^2 u(x, y)}{\partial x^2} &= \frac{u(x - h, y) - 2 u(x, y) + u(x + h, y)}{h^2} + \cal{O}(h^2), \\ \frac{\partial^2 u(x, y)}{\partial y^2} &= \frac{u(x, y - h) - 2 u(x, y) + u(x, y + h)}{h^2} + \cal{O}(h^2), \end{aligned}$$

and thus

$$\Delta u(x, y) = \frac{u(x - h, y) + u(x, y - h) - 4 u(x, y) + u(x + h, y) + u(x, y + h)}{h^2} + \cal{O}(h^2).$$

Let us now consider a uniform $(m+2) \times (n+2)$ grid with spacing $h$ in both the $x$ and $y$ directions. We set $x_i = i h$ for $i = 0, 1, \ldots, m+1$ and $y_j = j h$ for $j = 0, 1, \ldots, n+1$. We will assume that $u$ is zero on the boundary of the domain, i.e., $u(x_0, \cdot) = u(x_{m+1}, \cdot) = u(\cdot, y_0) = u(\cdot, y_{n+1}) = 0$.

Now let $D \in \mathbb{R}^{m \times n}$ be a matrix where $D_{i, j}$ is the finite difference approximation of $\Delta u(x_i, y_j)$, i.e.,

$$D_{i, j} = \frac{u(x_{i-1}, y_j) + u(x_i, y_{j-1}) - 4 u(x_i, y_j) + u(x_{i+1}, y_j) + u(x_i, y_{j+1})}{h^2}.$$

We now let $U \in \mathbb{R}^{m \times n}$ be the matrix with entries $U_{i, j} = u(x_i, y_j)$ and we seek an efficient way to compute $D$ from $U$.

One way to do this is to see that $D$ can come from $U$ by sweeping through the grid and applying a five-point stencil to each point:

But how to do this efficiently using matrix operations? We definitely want to avoid looping over the grid ourselves. Instead, we can split the five-point stencil into two parts, one horizontal and one vertical:

This shows that we can compute $D$ from $U$ as

$$D = L_m U + U L_n^T,$$

where $L_k$ is the 1D Laplace matrix of order $k$:

$$L_k = \frac{1}{h^2} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & -2 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{bmatrix} \in \mathbb{R}^{k \times k}.$$

Note how the first term, $L_m U$, applies the three-point stencil to each row of $U$ and that the second term, $U L_n^T$, applies the stencil to each column of $U$ (transposing $L_n$ is strictly not necessary here as $L_n$ is symmetric).

Note also how $L$ contains 3 non-zero values per row (except for the first and last rows), so assuming that it is represented in a sparse (compressed row) format, the matrix-matrix product $L_m U$ requires roughly $6 m n$ floating point operations (3 multiplications and 3 additions per output entry) and likewise for $U L_n^T$. That is a total of $12 m n$ operations to compute $D$ from $U$.

Another approach is to use the identities related to the $\operatorname{vec}$ and Kronecker product operators as seen in an earlier post. Using these, we get

$$\begin{aligned} \operatorname{vec}(D) &= \operatorname{vec}(L_m U + U L_n^T) \\ &= \operatorname{vec}(L_m U) + \operatorname{vec}(U L_n^T) \\ &= (I_n \otimes L_m) \operatorname{vec}(U) + (L_n \otimes I_m) \operatorname{vec}(U) \\ &= (I_n \otimes L_m + L_n \otimes I_m) \operatorname{vec}(U), \end{aligned}$$

where $I_k$ is the $k \times k$ identity matrix.

So the matrix operator here is

$$L_{m \times n} = I_n \otimes L_m + L_n \otimes I_m = L_n \oplus L_m$$

where $\oplus$ denotes the Kronecker sum (available in SciPy as kronsum). It contains (mostly) 5 non-zero values for each of its $m n$ rows, so the matrix-vector product $L_{m \times n} \operatorname{vec}(U)$ requires roughly $10 m n$ floating point operations.

Note that the approach of splitting the five-point stencil into two parts can also be done using a two-point stencil $1 \diamond 0 \diamond 1$ in one direction and $1 \diamond -4 \diamond 1$ in the other. This will reduce the number of floating point operations down to $4 m n + 6 m n = 10 m n$, just like the Kronecker sum approach.