5

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$?

Natasha
  • 421
  • 4
  • 18
  • 2
    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. – GoHokies Nov 14 '18 at 20:34
  • 3
    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). – rchilton1980 Nov 14 '18 at 20:49
  • @rchilton1980 yes, and that's exactly what the LDL decomposition produces. – GoHokies Nov 14 '18 at 21:17
  • 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. – Wolfgang Bangerth Nov 15 '18 at 04:26
  • The whole point is that you can't find such a factorization because $D$ is sparse and so your factorization is sparse. But we know that the inverse of the Laplace matrix is not sparse, and so it's factorization can also not be sparse. – Wolfgang Bangerth Nov 15 '18 at 04:27
  • @WolfgangBangerth Could you please explain the last point? The Laplacian matrix in my example is not a full rank matrix and the inverse doesn't exist. May I ask for another clarification? Should the Laplacian always be a symmetric matrix? Can we factorize a Laplacian matrix that has complex eigenvalues? – Natasha Nov 15 '18 at 04:50
  • @GoHokies Many thanks for the reply. A naive question, when the positive semidefinite matrix is factored as $LDL^T$, will $L^T$ be the grad operator? – Natasha Nov 15 '18 at 05:30
  • 1
    Even if the inverse does not exist, take a pseudo-inverse. It will also not be sparse. – Wolfgang Bangerth Nov 15 '18 at 14:17
  • @WolfgangBangerth Would it be correct to write Laplacian,L, in terms of incidence matrix,I? $-L=I*I^T $? – Natasha Nov 16 '18 at 04:19
  • What is the incidence matrix? – Wolfgang Bangerth Nov 17 '18 at 00:28
  • @WolfgangBangerth If we consider the convection-dispersion that occurs in a cylinder of length L, suppose L is discretized into 4 segments(edge) with 5 nodes(vertex). Incidence matrix(I) can be defined as an edge x vertex matrix. When the fluid flow is in the forward direction, the edges are directed forward. – Natasha Nov 17 '18 at 02:36
  • @WolfgangBangerth The first row of I will be filled as follows, the first row represents edge 1 .columns represent 5 vertices.Edge 1 leaves the first vertex(-1) and enters second vertex . Row 1 of Incidence matrix = [ -1 1 0 0 0].This would represent the discretized form,forward difference of the first derivative – Natasha Nov 17 '18 at 02:37
  • @rchilton1980 Did you mean "using the (signed) edge-to-vertex adjacency matrix" , signed edge-to-vertex incidence matrix in place of adjacency matrix? Also ,will the last row (5th row)of D D=[-1 +1 0 0 0; 0 -1 +1 0 0; 0 0 -1 +1 0; 0 0 0 -1 +1]be zeros? – Natasha Nov 17 '18 at 02:52

2 Answers2

5

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.

lightxbulb
  • 2,082
  • 8
  • 14
3

@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.

Federico Poloni
  • 11,344
  • 1
  • 31
  • 59
  • 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. – lightxbulb Aug 06 '23 at 11:31
  • 1
    @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. – Federico Poloni Aug 06 '23 at 11:38