5
$\begingroup$

I am trying to factorize the following Laplacian matrix in terms of $ D^TD$, D is the first derivative matrix. The tridiagonal form of the secon derivative matrix using Neumann boundary condition is given by,

$\begin{bmatrix} -1 & 1 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 &0\\ 0 & 1 & -2 &1&0\\0&0&1&-2&1\\0&0&0&1&-1\end{bmatrix}$.

In order to write the above in terms of the gradient( first derivative) operator, I tried factorizing the Laplacian using Cholesky factorization. But, I couldn't succeed in factorizing.

Any suggestions on how to find a decomposition of the form $D^TD$?

$\endgroup$
13
  • 2
    $\begingroup$ naïve Cholesky won't work, because your Laplacian matrix is not positive definite. there are many square roots (=choices of $D$) to select from. to recover the "upwind" discretization of the first order derivative operator, you'll want to use the LDL decomposition. see here for a more thorough discussion. $\endgroup$
    – GoHokies
    Nov 14, 2018 at 20:34
  • 3
    $\begingroup$ Isn't it just D=[-1 +1 0 0 0; 0 -1 +1 0 0; 0 0 -1 +1 0; 0 0 0 -1 +1]? I think you can also extend to higher dimensions, using the (signed) edge-to-vertex adjacency matrix. (Essentially a discrete version of the gradient operator). $\endgroup$ Nov 14, 2018 at 20:49
  • $\begingroup$ @rchilton1980 yes, and that's exactly what the LDL decomposition produces. $\endgroup$
    – GoHokies
    Nov 14, 2018 at 21:17
  • $\begingroup$ But the Laplace matrix is not just $D^T D$. That's because the Laplacian operator is the divergence of the gradient, so the two operators are transposes of each other. But the Laplace matrix is not of this form. $\endgroup$ Nov 15, 2018 at 4:26
  • 1
    $\begingroup$ Even if the inverse does not exist, take a pseudo-inverse. It will also not be sparse. $\endgroup$ Nov 15, 2018 at 14:17

2 Answers 2

5
$\begingroup$

It is sufficient if you consider a $D$ that uses forward or backward differences with reflecting boundaries: \begin{equation} D_f = \frac{1}{h}\begin{bmatrix} -1 & 1 & & \\ & \ddots &\ddots & \\ && -1 & 1 \\ &&&0 \end{bmatrix}, \quad D_b = \frac{1}{h}\begin{bmatrix} 0 & & & \\ -1 & 1 & & \\ & \ddots &\ddots & \\ && -1 & 1 \end{bmatrix}. \end{equation} Then you can factorise the 3-point stencil matrix with reflecting boundaries as \begin{equation} \frac{1}{h^2}\begin{bmatrix} -1 & 1 & &&\\ 1 & -2 & 1 && \\ &\ddots& \ddots & \ddots &\\ && 1&-2&1 \\ &&&1 & -1 \end{bmatrix} = -W = -D_{b}^TD_b = -D_f^TD_f. \end{equation} Note that neither $D^2_b$ nor $D^2_f$ will give you the desired matrix. The above discretisation is consistent with the following continuous reformulation $\partial_{xx} = -\partial_{x}^*\partial_{x}$, where $\partial_x^*$ is the adjoint of $\partial_x$, i.e. $D\approx\partial_x$ and $-D^T\approx\partial_{x}^*$.

The above works for the finite difference method as you may note. In the finite element method $\Delta \approx -W = -\int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j$ which is also reminiscent of the above, except summation now becomes an integral. For a grid triangulation of a rectangular domain however the FEM stiffness matrix is the same as the one resulting from the finite difference method. You can get the same by using the two-point flux approximation scheme from the finite volume method too.

$\endgroup$
3
$\begingroup$

@lightxbulb's answer gives the correct factorization already, but since you mention failed attempts with the Cholesky factorization, let me describe a method to discover the factorization numerically, in the "teach a man to fish" spirit.

As the comments note, the factorization fails because $A$ is not positive definite; in fact it is negative semidefinite, as you can check with eig(A). But even chol(-A) fails, because the matrix is singular.

However, you can compute a Cholesky factorization of a small perturbation $-A + \varepsilon I$, designed to make it positive definite; for instance in Octave

octave:1> A = toeplitz([-2 1 0 0 0]); A(1,1) = -1; A(end,end) = -1
A =

  -1   1   0   0   0
   1  -2   1   0   0
   0   1  -2   1   0
   0   0   1  -2   1
   0   0   0   1  -1

octave:2> chol(-A + sqrt(eps)*eye(size(A)))
ans =

   1.0000  -1.0000        0        0        0
        0   1.0000  -1.0000        0        0
        0        0   1.0000  -1.0000        0
        0        0        0   1.0000  -1.0000
        0        0        0        0   0.0003

This factor strongly suggests the result for the unperturbed problem.

$\endgroup$
2
  • $\begingroup$ While this works I wouldn't suggest perturbing the problem, because this essential makes it a reaction-diffusion problem. The consequent thing is to compute the QR decomposition and get the Cholesky decomposition from that for singular matrices. Fortunately the Laplacian matrix is well-known, and even its eigenvectors and eigenvalues are known. Although in the general case I would even question the idea of factorising this singular matrix, in practice this would suggest that one is trying to solve an ill-posed boundary value problem. $\endgroup$
    – lightxbulb
    Aug 6, 2023 at 11:31
  • 1
    $\begingroup$ @lightxbulb Sure, I suggested this method just as a heuristic to find out a formula to later be proved theoretically. I would also discourage using the perturbed computed factors for an actual numerical computation. $\endgroup$ Aug 6, 2023 at 11:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.