On the Inversion of Polynomials of Discrete Laplace Matrices

The efficient inversion of matrix polynomials is a critical challenge in computational mathematics. We design a procedure to determine the inverse of matrices polynomial of multidimensional Laplace matrices. The method is based on eigenvector and eigenvalue expansions. The method is consistent with previous expressions of the inverse discretized Laplacian in one spatial dimension \citep{Vermolen_2022}. Several examples are given.

💡 Research Summary

The paper addresses the problem of explicitly inverting matrix polynomials that arise from discretizations of the Laplace operator in multiple dimensions. While the inversion of a single discrete Laplace matrix is a well‑studied task, many modern applications require the inverse of a polynomial in that matrix (for example, operators of the form −Δ + αI, or higher‑order combinations such as (−Δ)^2 + β(−Δ) + γI). Existing numerical approaches typically rely on iterative solvers (Krylov subspace methods, multigrid) or on analytic fundamental solutions, but these either incur substantial computational cost or suffer from truncation errors in higher dimensions.

The authors propose a unified, analytic framework based on eigenvalue–eigenvector expansions. For a symmetric positive‑definite matrix A (the finite‑difference representation of the Laplacian with Dirichlet boundary conditions), the Principal Axis Theorem guarantees an orthogonal set of eigenvectors {vₖ} and real eigenvalues {λₖ}. Any right‑hand side b can be expanded as b = ∑βₖ vₖ with βₖ = (b,vₖ)/(vₖ,vₖ). Solving Ax = b then yields the closed‑form solution
 x = ∑ (βₖ/λₖ) vₖ,
provided λₖ ≠ 0. This expression mirrors the separation‑of‑variables technique for continuous PDEs and does not require the eigenvectors to be normalized.

To obtain the inverse matrix itself, the paper treats each column gₖ of A⁻¹ as the solution of Agₖ = eₖ (the k‑th canonical basis vector). Substituting the previous solution formula gives
 gₖ = ∑ (vⱼₖ/λⱼ) vⱼ,
where vⱼₖ denotes the k‑th component of eigenvector vⱼ. Consequently, the (i,k) entry of A⁻¹ can be written as
 (A⁻¹)ᵢₖ = ∑ (vⱼₖ vⱼᵢ)/λⱼ.
Because A is symmetric positive‑definite, A⁻¹ is also symmetric, and the formula automatically respects that symmetry.

The central theoretical contribution is a generalization to matrix polynomials. For any polynomial P(x)=∑_{m=0}^k a_m x^m, the authors prove that the eigenvectors of P(A) are identical to those of A, while the eigenvalues are transformed as P(λₖ). The proof proceeds in two steps: first, by induction, they show A^m vₖ = λₖ^m vₖ; second, they linearly combine the powers to obtain P(A)vₖ = P(λₖ)vₖ. Hence, if P(λₖ)≠0 for all k, the polynomial matrix P(A) is invertible and its inverse is given by
 (P(A)⁻¹)ᵢₖ = ∑ (vⱼₖ vⱼᵢ)/P(λⱼ).
Thus, once the eigen‑decomposition of the original Laplacian matrix A is known, the inverse of any polynomial in A can be assembled without further matrix factorizations.

The authors illustrate the theory with a detailed one‑dimensional example. For the Dirichlet problem −u’’ = f on (0,1) discretized with step h = 1/(n+1), the resulting tridiagonal matrix A has eigenvectors vⱼ whose components are vⱼₖ = √2 sin(jπkh) (the sampled continuous eigenfunctions) and eigenvalues λⱼ = 2h⁻²(1 − cos(jπh)) = 4h⁻² sin²(jπh/2). Substituting these into the general formula reproduces the well‑known explicit inverse entries
 (A⁻¹)ᵢₖ = h²/(n+1)·min(i,k)·(n+1 − max(i,k)),
and also matches the alternative expression obtained via the eigen‑expansion:
 (A⁻¹)ᵢₖ = h² ∑_{j=1}^n sin(jπkh) sin²(jπh/2) sin(jπih).
The paper provides an algebraic proof of the equivalence, noting that the two forms are related through discrete Fourier transform identities.

For higher dimensions, the authors argue that fundamental‑solution based methods become inefficient because the continuous Green’s function is logarithmic (in 2‑D) or involves higher‑order singularities, leading to non‑vanishing truncation errors. In contrast, the eigen‑expansion approach remains exact as long as the eigenpairs of the discrete Laplacian are known (which can be obtained analytically for regular grids or numerically via standard eigensolvers). They also discuss the case of Neumann boundary conditions, where the Laplacian matrix becomes singular; while A⁻¹ does not exist, certain polynomial combinations P(A) can be nonsingular, allowing the same eigen‑based inversion technique to be applied.

In conclusion, the paper presents a conceptually simple yet powerful method for inverting matrix polynomials derived from discrete Laplace operators. By leveraging the spectral decomposition of the base matrix, the approach yields closed‑form expressions for both the solution of linear systems and the inverse of any polynomial in the matrix, with guaranteed symmetry and positive‑definiteness when appropriate. The authors suggest future work on extending the framework to non‑symmetric or non‑normal matrices, handling nonlinear operators, and developing high‑performance parallel implementations for large‑scale scientific computing.