Efficient Solution Techniques for a Finite Element Thin Plate Spline Formulation

We present a new technique for solving the saddle point problem arising from a finite element based thin plate spline formulation. The solver uses the Sherman–Morrison–Woodbury formula to divide the domain into different regions depending on the properties of the data projection matrix. We analyse the conditioning of the resulting system on certain data distributions and use the results to develop effective preconditioners. We show our approach is efficient for a wide range of parameters by testing it on a number of different examples. Numerical results are given in one, two and three dimensions.

1. http://www.scilab.org/

2. www.python.org/

The authour would like to thank the anonymous referees for their valuable comments. The time and care they have taken in reviewing the paper is much appreciated.

Authors and Affiliations

1. Mathematical Sciences Institute, Australian National University, Canberra, ACT, 0200, Australia

   Linda Stals

Corresponding author

Correspondence to Linda Stals.

Appendix: Eigenvalue Calculations

The following calculations find the eigenvalues referenced in Sect. 6.

1.1 Eigenvalues for One Dimensional Grid

Let \(\varvec{w}^{\mu }\) be the vector with entries \(w^{\mu }_j = \sin \left( 2\pi \mu j h\right) \) and \(\varvec{z}^{\mu }\) be the vector with entries \(z^{\mu }_j = \cos \left( 2\pi \mu j h\right) \). In all of the theorems stated below \(m\) is taken to be odd.

Suppose the matrices \(\tilde{A}\), \(L\) and \(G_1\) are defined on a uniform grid in one dimension as in Sect. 5.1. Then the following theorem gives, most of, the eigenvalues of \(G_1L^{-1}G_2\).

Theorem 1

Given the eigenvectors \(\varvec{w}^{\mu }\), \(m-1\) eigenvalues of \(G_1L^{-1}G_1^T\) are

for \(1\le \mu \le m\) and \(\mu \ne (m+1)/2\).

Proof

After applying the stencils in a similar way as is done in Hackbusch [11] we find that for \(1 \le j \le m\),

For \(1<j<m\),

The eigenvalues of \(L\) are listed in Sect. 5.1. The boundary conditions need to be treated a little more carefully. If \(j = 1\), \(\left( L\varvec{z}^{\mu }\right) _j = \lambda ^{2\mu }\left( L\right) z^{\mu }_j+\frac{1}{h}.\) If \(j = m\), \(\left( L\varvec{z}^{\mu }\right) _j = \lambda ^{2\mu }\left( L\right) z^{\mu }_j+\frac{1}{h}.\) So, \(L\varvec{z}^{\mu } = \lambda ^{2\mu }\left( L\right) \varvec{z}^{\mu } + \frac{1}{h}\varvec{e}_1+ \frac{1}{h}\varvec{e}_m = \lambda ^{2\mu }\left( L\right) \varvec{z}^{\mu }+ L\varvec{z}^0\), where \(\varvec{e}_1\) and \(\varvec{e}_m\) are the first and \(m\)th unit vectors. Hence

For \(1 < j < m\),

To get \(\left( G_1\left( \varvec{z}^{\mu }-\varvec{z}^0\right) \right) _j\) = \( \sin (2\pi \mu h) w_j^{\mu }\) for the end points we use the fact that if \(j = 1\), \(\cos \left( 2\pi \mu h\left( j-1\right) \right) = 1\) and if \(j=m\), \(\cos \left( 2\pi \mu h\left( j+1\right) \right) =1\).

Finally, combining the above results gives

\(\square \)

Observe that if \(\mu = (m+1)/2\), \(\varvec{w}^{\mu } = 0\) so \(\varvec{w}^{(m+1)/2}\) is not an eigenvector. We have tried to find the one remaining eigenvector, but have not been able to. Fortunately, by using the theory in [1] and [13] we are able to continue our analysis even without knowing all of the eigenvalues.

Finding that \(G_1L^{-1}G_1^T\varvec{w}^{\mu } = \frac{h}{4} \sin ^{-2}\left( \pi \mu h \right) \sin ^2\left( 2\pi \mu h \right) \varvec{w}^{\mu }\) makes sense as we can view \(G_1G_1^T\) as the Laplacian defined on a grid of twice the spacing.

Given the particular form of the eigenvectors of \(G_1L^{-1}G_1^T\), we can also find the eigenvectors of \(\alpha \tilde{A}^{-1}+L^{-1}\left( G_1L^{-1}G_1^T\right) L^{-1}\).

Theorem 2

\(m-1\) eigenvalues of \(\tilde{E}_{\alpha } = \alpha \tilde{A}^{-1}+L^{-1}\left( G_1L^{-1}G_1^T\right) L^{-1}\) are

with corresponding eigenvectors \(\varvec{w}^{\mu }\), for \(1\le \mu \le m\) and \(\mu \ne (m+1)/2\).

Proof

From the analysis in Sect. 5.1 we know that \(\varvec{w}^{\mu }\) is an eigenvector of \(\tilde{A}\) with eigenvalue \(\lambda ^{2\mu }\left( \tilde{A}\right) \), and similarly \(\varvec{w}^{\mu }\) is an eigenvector of \(L\) with eigenvalue \(\lambda ^{2\mu }\left( L\right) \). So, using the results from Theorem 2 we get

\(\square \)

Using the arguments given in Theorem 2 we can readily find the eigenvalues of some of the preconditioned systems discussed in Sect. 6.1.1.

Theorem 3

\(m-1\) eigenvalues of \(\tilde{A}\tilde{E}_{\alpha } = \alpha I +\tilde{A} L^{-1}\left( G_1L^{-1}G_1^T\right) L^{-1}\) are

with corresponding eigenvectors \(\varvec{w}^{\mu }\) for \(1\le \mu \le m\) and \(\mu \ne (m+1)/2\).

Theorem 4

\(m-1\) eigenvalues of \(L\tilde{E}_{\alpha }L = \alpha L\tilde{A}^{-1}L + \left( G_1L^{-1}G_1^T\right) \) are

with corresponding eigenvectors \(\varvec{w}^{\mu }\) for \(1\le \mu \le m\) and \(\mu \ne (m+1)/2\).

Theorem 5

\(m-1\) eigenvalues of

are

with corresponding eigenvectors \(\varvec{w}^{\mu }\) for \(1\le \mu \le m\) and \(\mu \ne (m+1)/2\).

Proof

As \(\varvec{u}^{\mu }\) defined in Sect. 5.1 are eigenvectors for \(\tilde{C}_{\alpha }\), we see that \(\varvec{w}^{\mu }\) are also eigenvectors of \(\tilde{C}_{\alpha }\) with corresponding eigenvalues \(\lambda ^{2\mu }\left( \tilde{C}_{\alpha }\right) \). So, to find the eigenvalues of \(\tilde{C}^{-1}_{\alpha } \tilde{E}_{\alpha }\), we divide \(\lambda ^{\mu }\left( \tilde{E}_{\alpha }\right) \) by \(\lambda ^{2\mu }\left( \tilde{C}_{\alpha }\right) \). \(\square \)

1.2 Eigenvalues for Two Dimensional Grids

To find eigenvalues for a two dimensional grid, we use an approach similar to the one used for one dimensional grids in Sect. 9.1.

Let \(\varvec{w}^{\mu ,\nu }\) be the vector with entries \(w^{\mu , \nu }_{j, k} = \sin \left( 2\pi (\mu j + \nu k) h\right) \) and \(\varvec{z}^{\mu }\) be the vector with entries \(z^{\mu ,\nu }_{j,k} = \cos \left( 2\pi (\mu j +\nu k)h\right) \). In all of the theorems stated below, \(m\) is taken to be odd.

Suppose the matrices \(\tilde{A}\), \(G_1\) and \(G_2\) are defined on a uniform grid in two dimensions as described in Sect. 5.2, then the following theorem gives some of eigenvalues of the resulting system of equations.

Theorem 6

\(m-1\) eigenvalues of \(\tilde{A}\) are

with corresponding eigenvectors \(\varvec{w}^{\mu , \nu }\), for \(1\le \mu , \nu \le m\) and both \(\mu \) and \(\nu \) not equal to \((m+1)/2\).

Proof

Firstly, \(\varvec{w}^{\mu , \nu } = \varvec{0}\) if \(\mu j+\nu k\) is a multiple of \((m+1)/2\) for all \(1 \le j,k \le m\). Consequently, \(\varvec{w}^{\mu , \nu } \ne \varvec{0}\) unless \(\mu \) = \(\nu \) = \((m+1)/2\).

For \(1 \le j, k \le m\),

In the equation above we applied the stencil to determine the action of the matrix on the given vector. See Hackbusch [11]. \(\square \)

Theorem 7

Consider \(\lambda ^{\mu ,\nu }(L)\) given in Sect. 5.2, \(\lambda ^{\mu ,\nu }(\tilde{A})\) given in Theorem 6 and let

and

Then, \(m-1\) eigenvalues of

are

with corresponding eigenvectors \(\varvec{w}^{\mu }\), for \(1\le \mu , \nu \le m\) and both \(\mu \) and \(\nu \) not equal to \((m+1)/2\).

Proof

For \(1 \le j, k \le m\),

Similarly \(G_2\varvec{w}^{\mu , \nu } = \rho ^{\mu , \nu }(G_2)\varvec{z}^{\mu ,\nu }\).

For \(1 < j, k < m\),

If \(j = 1\) and \(1 < k < m\) use \(\cos \left( 2\pi \left( \mu (j-1) + \nu k\right) h\right) =1\) to get

After taking the other boundary conditions into account we find that \(L^{-1}\varvec{z}^{\mu , \nu } = \left( \lambda ^{2\mu , 2\nu }(L)\right) ^{-1}\left[ \varvec{z}^{\mu , \nu }-\varvec{z}^{0, 0}\right] \).

For \(1 < j, k < m\),

Once again, after take care with the boundary conditions we find that \(G_1\left( \varvec{z}^{\mu , \nu }-\varvec{z}^{0,0}\right) = \rho ^{\mu , \nu }(G_1)\varvec{w}^{\mu ,\nu }\). Furthermore, we can use the same procedure to show \(G_2\left( \varvec{z}^{\mu , \nu }-\varvec{z}^{0,0}\right) = \rho ^{\mu , \nu }(G_2)\varvec{w}^{\mu ,\nu }\).

We now combine the above results to show that

with corresponding eigenvectors \(\varvec{w}^{\mu , \nu }\), and both \(\mu \) and \(\nu \) are not \((m+1)/2\).

Finally, as \(\tilde{A}\), \(L\) and \(G_1L^{-1}G_1^T+G_2L^{-1}G_2^T\) all have the same eigenvectors,

\(\square \)

Theorem 8

Consider \(\tilde{E}_{\alpha }\), \(\lambda ^{\mu ,\nu }(L)\), \(\lambda ^{\mu ,\nu }(\tilde{A})\), \(\rho ^{\mu , \nu }(G_1)\) and \(\rho ^{\mu , \nu }(G_2)\) given in Theorem 7, then \(m-1\) eigenvalues of \(\tilde{A}\tilde{E}_{\alpha }\) are

with corresponding eigenvectors \(\varvec{w}^{\mu }\), for \(1\le \mu , \nu \le m\) and both \(\mu \) and \(\nu \) not equal to \((m+1)/2\).

Proof

See Theorems 6 and 7. The results follow since \(\varvec{w}^{\mu }\) is an eigenvector for both \(\tilde{A}\) and \(\tilde{E}_{\alpha }\).

Theorem 9

Consider \(\tilde{E}_{\alpha }\), \(\lambda ^{\mu ,\nu }(L)\), \(\lambda ^{\mu ,\nu }(\tilde{A})\), \(\rho ^{\mu , \nu }(G_1)\) and \(\rho ^{\mu , \nu }(G_2)\) given in Theorem 7. Furthermore let \(\tilde{C}_{\alpha } = \alpha \tilde{A}^{-1}+h^2 L^{-2}\). Then \(m-1\) eigenvalues of \(\tilde{C}_{\alpha }^{-1}\tilde{E}_{\alpha }\) are

with corresponding eigenvectors \(\varvec{w}^{\mu }\), for \(1\le \mu , \nu \le m\) and both \(\mu \) and \(\nu \) not equal to \((m+1)/2\).

Proof

See Theorems 6 and 7. The results follows since \(\varvec{w}^{\mu }\) is an eigenvector of \(\tilde{A}\), \(\tilde{E}_{\alpha }\) and \(L\). \(\square \)

Reprints and permissions