A Multigrid Preconditioner for Tensor Product Spline Smoothing

TL;DR

This paper introduces a matrix-free multigrid preconditioner for tensor product spline smoothing that reduces memory usage and achieves scalable convergence in high-dimensional nonparametric regression.

Contribution

It presents a novel multigrid preconditioned conjugate gradient method for tensor product B-spline smoothing that is memory-efficient and scalable across fixed dimensions.

Findings

Achieves grid independent convergence in fixed dimensions.

Requires moderate memory, suitable for high-dimensional data.

Demonstrates scalability and efficiency in multivariate smoothing.

Abstract

Uni- and bivariate data smoothing with spline functions is a well established method in nonparametric regression analysis. The extension to multivariate data is straightforward, but suffers from exponentially increasing memory and computational complexity. Therefore, we consider a matrix-free implementation of a geometric multigrid preconditioned conjugate gradient method for the regularized least squares problem resulting from tensor product B-spline smoothing with multivariate and scattered data. The algorithm requires a moderate amount of memory and is therefore applicable also for high-dimensional data. Moreover, for arbitrary but fixed dimension, we achieve grid independent convergence which is fundamental to achieve algorithmic scalability.

Tables

Table 1. Table 1: Number of B-spline coefficients for varying spatial dimension P 𝑃 P .

	$P=2$	$P=3$	$P=4$	$P=5$
$K$	$1.225$	$42.875$	$\mathrm{1.500.625}$	$\mathrm{52.521.875}$

Table 2. Table 2: Required number of iterations of the considered methods for various dimensions with G = 5 𝐺 5 G=5 grids and computational time in seconds.

	CG		MGCG_JAC		MGCG_SSOR
P=1	22	(1.05)	2	(0.91)	1	(<0.01)
P=2	73	(4.41)	4	(4.59)	2	(0.37)
P=3	367	(103.31)	14	(79.9)	-
P=4	747	(2745.98)	19	(1927.60)	-

Table 3. Table 3: Condition number of the coefficient matrices for G = 5 𝐺 5 G=5 grids and P = 2 𝑃 2 P=2 dimensions.

	CG	MGCG_JAC	MGCG_SSOR
condition	$1933.27$	$1.82$	$1.03$

Equations

$$y_i=s(x_i)+{\epsilon }_i,i=1,\dots ,n.$$

$$\mathcal{K}:=\{a={\kappa }_0<\dots <{\kappa }_{m+1}=b\}.$$

$${\mathcal{S}}_q(\mathcal{K}):=\{s\in {\mathcal{C}}^{q-1}(\Omega ):s{|}_{\left[{\kappa }_{j-1},{\kappa }_j\right]}\in {\mathcal{P}}_q\left(\left[{\kappa }_{j-1},{\kappa }_j\right]\right),j=1,\dots ,m+1\}$$

$$\{{\phi }_{j,q}:j=1,\dots ,J\}$$

$$\{{\phi }_{j_p,q_p}^p:j_p=1,\dots ,J_p\}$$

$${\phi }_{j,q}:\Omega :={\Omega }_1\times \dots \times {\Omega }_P\to \mathbb{R},{\phi }_{j,q}(x)=\prod _{p=1}^P{\phi }_{j_p,q_p}^p(x^p),$$

$${\mathcal{S}}_q(\mathcal{K}):=\text{span}\{{\phi }_{j,q}:1\le j\le J:=(J_1,\dots ,J_P){}^{\prime }\},$$

$$s=\sum _{1\le j\le J}{\alpha }_j{\phi }_{j,q}$$

$$s=\sum _{1\le j\le J}{\alpha }_j{\phi }_{j,q}=\sum _{k=1}^K{\alpha }_k{\phi }_{k,q}.$$

$$\underset{s\in {\mathcal{S}}_q(\mathcal{K})}{\min }\sum _{i=1}^n{\left(s(x_i)-y_i\right)}^2+\lambda \underset{{\mathbb{R}}^P}{\int }{\left(\sum _{p_1=1}^P\sum _{p_2=1}^P\frac{{\partial }^2}{\partial x_{p_1}\partial x_{p_2}}s(x)\right)}^2\mathrm{d}x.$$

$$\mathcal{LS}(s):=\sum _{i=1}^n{\left(s(x_i)-y_i\right)}^2$$

$$\mathcal{R}(s):=\underset{{\mathbb{R}}^P}{\int }{\left(\sum _{p_1=1}^P\sum _{p_2=1}^P\frac{{\partial }^2}{\partial x_{p_1}\partial x_{p_2}}s(x)\right)}^2\mathrm{d}x,$$

$$\mathcal{LS}(s)=\Vert \Phi \alpha -y{\Vert }_2^2,$$

$$\mathcal{R}(s)=\alpha {}^{\prime }\left(\sum _{r\in {\mathbb{N}}_0^P;|r|=2}\frac{2}{r!}{\Psi }_r\right)\alpha ,$$

$$\Lambda :=\sum _{r\in {\mathbb{N}}_0^P;|r|=2}\frac{2}{r!}{\Psi }_r\in {\mathbb{R}}^{K\times K}$$

$$\underset{\alpha \in {\mathbb{R}}^K}{\min }\Vert \Phi \alpha -y{\Vert }_2^2+\lambda \alpha {}^{\prime }\Lambda \alpha .$$

$$\left(\Phi {}^{\prime }\Phi +\lambda \Lambda \right)\alpha \stackrel{!}{=}\Phi {}^{\prime }y.$$

$$K=\prod _{p=1}^P{J}_p=35^P$$

$${\Phi }_p\in {\mathbb{R}}^{n\times J_p},{\Phi }_p[i,j_p]:={\phi }_{j_p,q_p}^p(x_i^p),$$

$$\Phi {}^{\prime }=\underset{p=1}{\overset{P}{⨀}}{\Phi }_p{}^{\prime },$$

$${\Psi }_{r_p}^p\in {\mathbb{R}}^{J_p\times J_p},{\Psi }_{r_p}^p[j_p,{\ell }_p]={⟨{\partial }^{r_p}{\phi }_{j_p,q_p}^p,{\partial }^{r_p}{\phi }_{{\ell }_p,q_p}^p⟩}_{L^2({\Omega }_p)}.$$

$${\Psi }_r=\underset{p=1}{\overset{P}{⨂}}{\Psi }_{r_p}^p$$

$${\phi }_{j,q}^{2h}=\frac{1}{2^q}\sum _{i=0}^{q+1}\left(\genfrac{}{}{0ex}{}{q+1}{i}\right){\phi }_{2j-(q+1)+i,q}^h$$

$$s=\sum _{j=1}^{J^{2h}}{\alpha }_j^{2h}{\phi }_{j,q}^{2h}=\sum _{j=1}^{J^{2h}}{\alpha }_j^{2h}\frac{1}{2^q}\sum _{i=0}^{q+1}\left(\genfrac{}{}{0ex}{}{q+1}{i}\right){\phi }_{2j-(q+1)+i,q}^h.$$

$$s=\sum _{j=1}^{J^h}{\alpha }_j^h{\phi }_{j,q}^h$$

$$I_{2h}^h[i,j]:=\frac{1}{2^q}\left(\genfrac{}{}{0ex}{}{q+1}{i-2j+q+1}\right),i=1,\dots ,J^h,j=1,\dots ,J^{2h}.$$

$$I_{2h}^h=\underset{p=1}{\overset{P}{⨂}}{I}_{2h_p}^{h_p}\in {\mathbb{R}}^{K^h\times K^{2h}},$$

$$A:=\underset{p=1}{\overset{P}{⨂}}{A}_p\in {\mathbb{R}}^{m\times n},m:=\prod _{p=1}^P{m}_p,n:=\prod _{p=1}^P{n}_p$$

$$A:=\underset{p=1}{\overset{P}{⨀}}{A}_p\in {\mathbb{R}}^{m\times n},m:=\prod _{p=1}^P{m}_p$$

$$A=\underset{p=1}{\overset{P}{⨀}}{A}_p=\left[\begin{array}{c}\underset{p=1}{\overset{P}{⨂}}{A}_p[\cdot ,1],\dots ,\underset{p=1}{\overset{P}{⨂}}{A}_p[\cdot ,n]\end{array}\right]$$

Full text

Keywords: multidimensional smoothing $\mathbf{\cdot}$ tensor product B-splines $\mathbf{\cdot}$ memory efficiency $\mathbf{\cdot}$ multigrid methods

1 Introduction

In many statistical applications it is fundamental to investigate the relationship between explanatory variables and a variable of interest, which is generally modeled by a function of the covariates. This function attempts to capture the important patterns in the data while leaving out noise and other insignificant structures. This concept is known under several names like smoothing, (nonparametric) regression, or surface fitting. The input data is often multivariate and scattered, i.e. there is no inherent structure. To represent complex, e.g. multivariate and highly nonlinear data, the modeled function has to be very flexible in order to allow data-driven estimation of the complex effects. In one and two dimensions, there exist efficient smoothing methods, based on spline functions. We refer the reader to, e.g., Eilers and Marx, (1996), Ruppert et al., (2003), Wand and Ormerod, (2008), Fahrmeir et al., (2013). Unfortunately, the straightforward extension of these spline-based methods to multiple input variables suffers from an exponential growth of the number of parameters to be estimated within the spatial dimension. This issue is often referred to as the curse of dimensionality (cf. Bellman, (1957)). Therefore, the computational and especially memory complexity of the related estimation procedure becomes unjustifiably large, even for moderate spatial dimensions, such that smoothing methods are rarely applied for covariates of dimension larger than two (cf. Fahrmeir et al., (2013), p. 531). For gridded covariates an efficient algorithm, that extends the spline smoothing method of Eilers and Marx, (1996), is implemented in Eilers et al., (2006). However, for scattered covariates there is no satisfactory approach. The challenge is to deal with the large-scale linear systems which arise from the spline smoothing for scattered data sets with increasing covariate dimensions.

In order to overcome this issue, we apply two techniques in the solution process and investigate their performance. First, the matrix corresponding to the linear system is never explicitly formed and stored in memory. When using an iterative solution algorithm like the conjugate gradient (CG) method only the matrix-vector product is required, but not matrix entries explicitly. A second important ingredient is a suitable preconditioner for the linear system. For this purpose, we introduce multigrid techniques to the spline smoothing problem. With increasing space dimension and grid resolution the performance of the CG method usually deteriorates. We thus investigate a geometric multigrid preconditioned and matrix-free CG iteration, which significantly reduces memory and computational complexity compared to common estimation methods. Furthermore, the grid-independent convergence of the geometric multigrid is mandatory in order to achieve algorithmic scalability in the sense that simultaneously doubling the degrees of freedom and the number of processors leads to constant algorithmic running times.

The remainder of this paper is organized as follows. In Section 2, we formulate the spline smoothing problem in multiple dimensions based on tensor product B-splines and formulated the underlying large-scale linear system. In Section 3, we develop a matrix-free implementation of a multigrid preconditioned conjugate gradient method to solve the large-scale system with comparatively small memory and computational complexity. In Section 4, we apply the proposed method, also compared to traditional methods, on a test data set. Finally, Section 5 gives a conclusion of the paper.

2 Spline Smoothing in Multiple Dimensions

In statistics, smoothing a given data set describes the process of constructing a function that captures the important patterns in the data while leaving out noise and other fine scaled structures. More precisely, for given data $\{(x_{i},y_{i})\in\mathbb{R}^{P}\times\mathbb{R}\,|\,i=1,\ldots,n\}$, where the $y_{i}\in\mathbb{R}$ are observations of a continuous response variable and the $x_{i}\in\mathbb{R}^{P}$ represent the corresponding value of a continuous covariate, we seek a smooth, but further unspecified function $s\colon\Omega\subset\mathbb{R}^{P}\rightarrow\mathbb{R}$ such that

[TABLE]

A common assumption is that $\varepsilon_{1},\ldots,\varepsilon_{n}$ are independent and identically distributed (i.i.d.) random errors with zero mean, common variance $\sigma_{\varepsilon}^{2}$ and assumed to be independent of the covariates. We begin by defining a spline basis which forms the underlying, linear space for the representation of $s$.

2.1 Tensor Product Splines

Let $\Omega\coloneqq[a,b]$ be a bounded and closed interval partitioned by the knots

[TABLE]

Let $\mathcal{C}^{q}(\Omega)$ denote the space of $q$-times continuously differentiable functions and let $\mathcal{P}_{q}(\Omega)$ denote the space of polynomials of degree $q$. We call the function space

[TABLE]

the space of spline functions of degree $q\in\mathbb{N}_{0}$ with knots $\mathcal{K}$. It is a finite dimensional, linear space of dimension $J\coloneqq dim\left(\mathcal{S}_{q}(\mathcal{K})\right)=m+q+1$. With

[TABLE]

we denote its B-spline basis (cf. de Boor, (1978)), which is well suited for numerical applications. To extend the spline concept to $P$-dimensional covariates, a tensor product approach is common. Let $\mathcal{S}_{q_{p}}(\mathcal{K}_{p})$ be the spline space for the $p$-th covariate, $p=1,\ldots,P$, and let

[TABLE]

denote its B-spline basis. The function

[TABLE]

where $j\coloneqq(j_{1},\ldots,j_{P}){{}^{\prime}}$ and $q\coloneqq(q_{1},\ldots,q_{P}){{}^{\prime}}$ are multiindices, is called tensor product B-spline. We define the space of tensor product splines as their linear combination

[TABLE]

which is then a $K\coloneqq\prod_{p=1}^{P}J_{p}$ dimensional linear space. Note, that we use the same symbol for the univariable and the tensor product spline space as well as for B-splines and tensor product B-splines. The difference is that we make use of the multiindex notation for the tensor products. Every tensor product spline $s\in\mathcal{S}_{q}(\mathcal{K})$ therefore has a unique representation

[TABLE]

and for computational reasons we uniquely identify the set of multiindices $\{j\in\mathbb{N}^{P}:1\leq j\leq J\}$ in descending lexicographical order as $\{1,\ldots,K\coloneqq\prod_{p=1}^{P}J_{p}\}$ such that

[TABLE]

We are now prepared to formulate the so called smoothing spline problem.

2.2 Tensor Product Smoothing Spline

Spline smoothing, also known as regularized or penalized spline regression, is a popular method in data smoothing. We define (tensor product) smoothing splines in $\mathcal{S}_{q}(\mathcal{K})$ as solution of

[TABLE]

The objective function consists of several parts. One the one hand, we have a least squares fitting term

[TABLE]

that measures the goodness-of-fit of the smoothing spline to the given observations. On the other hand, we have a regularization term

[TABLE]

that penalizes the roughness of the spline function. The term $\mathcal{R}(s)$ is weighted by a regularization or smoothing parameter $\lambda>0$, which balances the two competitive terms $\mathcal{LS}(s)$ and $\mathcal{R}(s)$. For $\lambda\rightarrow 0$ the smoothing spline tends to be the common least squares spline, which heavily overfits, or even interpolates, the observations. Conversely, for $\lambda\rightarrow\infty$, to much impact is given to the regularization term such that the smoothing spline tends to be the least squares hyperplane (cf. Green and Silverman, (1993), p. 159). More details can be found in the monographs Eubank, (1988), Wahba, (1990), and Green and Silverman, (1993).

Using the unique B-spline representation (2.2) we obtain

[TABLE]

where $\Phi\in\mathbb{R}^{n\times K}$ is element wise defined as $\Phi[i,k]\coloneqq\varphi_{k,q}(x_{i})$ and

[TABLE]

where each $\Psi_{r}\in\mathbb{R}^{K\times K}$ is element wise defined as $\Psi_{r}[k,\ell]=\left\langle\partial^{r}\varphi_{k,q},\partial^{r}\varphi_{\ell,q}\right\rangle_{L^{2}(\Omega)}$. Defining

[TABLE]

the smoothing spline (2.3) is equivalently formulated in terms of the B-spline coefficients as

[TABLE]

This is a regularized least squares problem and its solution is given by the solution of the linear system

[TABLE]

The main focus of this paper is a memory efficient evaluation and solution of (2.4).

2.3 Curse of Dimensionality

Spline smoothing is known to suffer from the so called curse of dimensionality, which describes an exponential growth of the number of B-spline coefficients $K$ within the dimension of the covariates $P$. Choosing $J_{p}=35$ B-spline functions for each direction $p=1,\ldots,P$, which is a quite realistic number, the resulting tensor product B-spline is given by

[TABLE]

parameters. Table 1 gives an overview on the increasing number of degrees of freedom under increasing dimensionality with which we are dealing with in this work.

Obviously, the method of smoothing splines becomes impracticable for dimensions $P\geq 3$ (cf. Fahrmeir et al., (2013), p. 531), since solving the large-scale linear system (2.4) requires an unjustified computational effort. One reason for this is that, even when the coefficient matrix $\Phi{{}^{\prime}}\Phi+\lambda\Lambda$ is stored in a sparse format, the exponential growth ensures that the memory capacity of a common digital computer is exceeded already for moderate spatial dimensions $P$.

Therefore, to make tensor product spline smoothing practicable for increasing spatial dimensions $P$, we need to develop computational and memory efficient methods to solve the large-scale linear system (2.4).

2.4 Tensor Product Properties

Before focusing on efficient solution methods in the next section, we state some important properties on the occurring matrices that are fundamental for the proposed algorithm and are based on the tensor product nature of the underlying splines. First, we define for each spatial direction $p=1,\ldots,P$ the matrix

[TABLE]

which corresponds to the matrix $\Phi$ with a one-dimensional covariate in direction $p$. Then, because of the scattered data structure, it holds

[TABLE]

where $\odot$ denotes the Khatri-Rao product. For each Gramian matrix $\Psi_{r}$ we analogously define a Gramian matrix for each spatial direction of the respective derivatives, that is

[TABLE]

Then it holds

[TABLE]

due to the tensor property.

A further important property of uniform B-splines in one variable is the subdivision formula. Let $\mathcal{S}_{q}(\mathcal{K}^{2h})$ and $\mathcal{S}_{q}(\mathcal{K}^{h})$ denote univariable spline spaces with uniform knot set $\mathcal{K}^{2h}$ and $\mathcal{K}^{h}$ with mesh sizes $2h$ and $h$, respectively. Then it holds (cf. Höllig, (2003), p. 32)

[TABLE]

and consequently for a uniform spline $s\in\mathcal{S}_{q}(\mathcal{K}^{2h})$ we obtain

[TABLE]

Since $s\in\mathcal{S}_{q}(\mathcal{K}^{h})$, it holds

[TABLE]

and the B-spline coefficients for the different meshes are therefore related through $\alpha^{(h)}=I_{2h}^{h}\alpha^{2h}$, where $I_{2h}^{h}\in\mathbb{R}^{J^{h}\times J^{2h}}$ is element wise defined as

[TABLE]

Because of the tensor property of B-splines, the formula carries over to multivariable B-splines and the corresponding B-spline coefficients are related through the matrix

[TABLE]

where $h\coloneqq(h_{1},\ldots,h_{P}){{}^{\prime}}$ denotes the mesh vector.

3 Geometric Multigrid Preconditioner

In this section we focus on the numerical behavior of the normal equation (2.4), which is equivalent to solving the smoothing spline problem (2.3). Since the coefficient matrix $A\coloneqq\Phi{{}^{\prime}}\Phi+\lambda\Lambda$ is symmetric and positive definite, the conjugated gradient (CG) method is the straightforward choice. For our application it is particularly important, that for the CG algorithm only matrix-vector products are required but not explicit entries of the matrix. From a computational point of view assembling $A$ in a sparse format is infeasible for the targeted problem sizes. This limits the choice for linear solvers and preconditioners. It is common practice to apply a preconditioner to the CG method in order to speed up convergence. Since the convergence rate depends on the condition of the matrix we encounter a worsening while improving the approximation quality by refining the spline space (2.1). This behavior is more dramatic, the higher the spatial dimension is, which is reflected in the numerical test cases in Section 4. For the choice of a preconditioner it is again important, that only matrix-vector products are involved. This cancels out the well-established, incomplete Cholesky factorizations of $A$ as a preconditioner and we thus concentrate on geometric multigrid methods in this paper.

The origins of multigrid methods date back into the late 70th (see for instance Brandt, (1977), Hackbusch, (1978), or Trottenberg et al., (2000) for an overview). Since then, they are successfully applied in the field of partial differential equations. The main idea is to build a hierarchy of - in a geometrical sense - increasingly fine discretizations of the problem. One then uses a splitting iteration like Jacobi or symmetric successive overrelaxation (SSOR) to smooth different frequencies of the error $e=x^{\ast}-x$ on different grids, where $x^{\ast}$ solves $Ax^{\ast}=b$ and $x$ is an approximation. Due to the decreasing computational complexity the problem can usually be solved explicitly on the coarsest grid, e.g. by a factorization of the system matrix. Between the grids, interpolation and restriction operations are used in order to transport the information back and forth. This works for hierarchical grids, i.e. that all nodes of one grid are also contained in the next finer grid. We achieve this by successively halving the grid spacing from $2h$ to $h$. A typical choice for the restriction is the transposed interpolation given by $I_{h}^{2h}=\left(I_{2h}^{h}\right){{}^{\prime}}$.

The outstanding feature of the multigrid method is that under certain circumstances it is possible to solve sparse, linear systems in optimal $\mathcal{O}(m)$ complexity, where $m$ is the number of discretization points on the finest grid. Algorithm 1 shows the basic structure of one V-cycle that serves as a preconditioner in the CG method. Here $g\in\{1,\dots,G\}$ denotes the grid levels from coarse $g=1$ to fine $g=G$, which coincides with the original problem. The main ingredients of Algorithm 1 are:

• •

the system matrices $A_{g}$,

• •

interpolation matrices $I_{g-1}^{g}$ and restriction matrices $I_{g}^{g-1}$,

• •

smoothing iteration method,

• •

the coarse grid solver $A_{1}^{-1}$.

Note that the smoother (line 6 and 11 in Algorithm 1) does not necessarily have to be convergent on its own. Typical choices are Jacobi or SSOR iterations, for which the convergence depends on the spectral radius of the iteration matrix. Convergence is not necessary for the smoothing property. It is thus recommendable only to apply a small number of pre/post-smoothing steps $\nu_{1}$ and $\nu_{2}$, such that a diverging smoother does not affect the overall convergence. This is further addressed in the discussion of Section 4.

It might be tempting to apply an algebraic (AMG) instead of a geometric multigrid preconditioner. The attractivity stems from the fact that AMG typically works as a black-box solver on a given matrix without relying a geometric description of grid levels and transfer operators. Yet, common AMG implementations explicitly require access to the matrix, which is prohibitively memory consuming for the tensor-product smoothing splines.

In the following we concentrate on a memory efficient realization of the matrix-vector product in line 7 of Algorithm 1, the interpolation/restriction in lines 8 and 10 and the smoothing operation in lines 6 and 11.

3.1 Memory Efficient Matrix Operations

Matrix-free methods for the solution of the large-scale linear system (2.4) require the efficient computation of matrix-vector products of the occurring matrices. For our particular application, we exploit their special structure, namely Kronecker and Khatri-Rao product structure.

3.1.1 Kronecker Matrices

For given matrices $A_{p}\in\mathbb{R}^{m_{p}\times n_{p}}$ we consider the Kronecker matrix

[TABLE]

and aim for a matrix-free computation of matrix-vector products with $A$ by only accessing the Kronecker factors $A_{1},\ldots,A_{P}$. In the case of $m_{p}=n_{p}$ for all $p=1,\ldots,P$ an implementation is provided by Benoit et al., (2001) and we extend their idea to arbitrary factors by Algorithm 2 to form a matrix-vector product with a Kronecker matrix only depending on the Kronecker factors.

3.1.2 Khatri-Rao Matrices

For given matrices $A_{p}\in\mathbb{R}^{m_{p}\times n}$ we consider the Khatri-Rao matrix

[TABLE]

and aim for a matrix-free computation of matrix-vector products with $A$ and $A{{}^{\prime}}$ by only accessing the Khatri-Rao factors $A_{1},\ldots,A_{P}$. By definition of the Khatri-Rao product it holds

[TABLE]

such that

[TABLE]

for all $x\in\mathbb{R}^{n}$. This yields Algorithm 3 to form a memory-efficient matrix-vector product with a Khatri-Rao matrix only requiring the Khatri-Rao factors.

Similarly, it holds

[TABLE]

for all $y\in\mathbb{R}^{m}$. We can thus formulate Algorithm 4 to form a matrix-vector product with a transposed Khatri-Rao matrix by only accessing its Khatri-Rao factors.

For the implementation of the Jacobi smoother in the multigrid algorithm we need access to the diagonal of the matrix $AA^{{}^{\prime}}\in\mathbb{R}^{m\times m}$ without computing the entire matrix but only accessing the Khatri-Rao factors $A_{1},\ldots,A_{P}$. For $j=1,\ldots,m$ the $j$-th diagonal element of $AA{{}^{\prime}}$ is given by $e_{j}{{}^{\prime}}AA{{}^{\prime}}e_{j}=\|A{{}^{\prime}}e_{j}\|_{2}^{2}$, where $e_{j}$ denotes the $j$-th unit vector, and it holds

[TABLE]

This yields Algorithm 5 to extract the diagonal of $AA{{}^{\prime}}$ by only accessing its factors.

3.2 MGCG Algorithm for Spline Smoothing

We are now prepared to apply the multigrid V-cycle of Algorithm 1 as a preconditioner in the CG method to solve the linear system (2.4).

3.2.1 Hierarchy

For the multigrid method we require hierarchical grids denoted by $g=1,\ldots,G$, where $g=1$ is the coarsest and $g=G$ the finest. Since we do not assume initial knowledge on the data we propose to base the underlying spline space on $m_{p}=2^{g}-1$ equidistant knots in each space dimension $p=1,\ldots,P$. For $g=1,\ldots,G$ we then define the coefficient matrix

[TABLE]

With

[TABLE]

we denote the dimension of the spline space $\mathcal{S}_{q}(\mathcal{K}_{g})$ and hence the size of the coefficient matrix $A_{g}\in\mathbb{R}^{K_{g}\times K_{g}}$ on grid level $g$. Thus, we obtain a hierarchy of linear systems for which we can use the subdivision properties given in Section 2.

3.2.2 Smoothing Iteration

To apply the Jacobi method we additionally require the diagonal of the coefficient matrix, which is a vector of length $K_{g}$. In principle, this is not prohibitively memory consuming yet, since the coefficient matrix is not explicitly accessible, we cannot simply extract its diagonal. However, Algorithm 5 allows the memory efficient computation of the diagonal of $\Phi{{}^{\prime}}_{g}\Phi_{g}$, $g=1,\ldots,G$, and the diagonal of each matrix $\Psi_{r,g}$ is directly given by the diagonal property of Kronecker matrices. In contrast to the Jacobi smoother the SSOR method additionally requires explicit access to all elements in one of the triangular parts of the coefficient matrices. It would in principal be possible to compute the desired elements entry-wise within each iteration. Yet, this would be computationally very expensive, since these elements have to be computed repeatedly in each iteration and for each grid level. Therefore, we utilize the Jacobi method 6 as smoothing iteration.

We also test a SSOR smoother in the multigrid algorithm in Section 4, yet we apply this only to low dimensional tests, where we can explicitly assemble and keep the coefficient matrix in memory.

3.2.3 Grid Transfer

For the transfer of a data vector from grid $g$ to grid $g+1$ we require a prolongation matrix $I_{g}^{g+1}\in\mathbb{R}^{K_{g+1}\times K_{g}}$. Due to the subdivision formula and the tensor product nature of the splines we can use (2.6) for this purpose. For the restriction matrix $I_{g+1}^{g}\in\mathbb{R}^{K_{g}\times K_{g+1}}$ we choose $I_{g+1}^{g}\coloneqq\left(I_{g}^{g+1}\right){{}^{\prime}}$, which yields the Garlerkin property

[TABLE]

The restriction and prolongation matrices do not fit into memory as well, but to apply the V-cycle Algorithm 1 only matrix-vector products with $I_{g}^{g+1}$ and $I_{g+1}^{g}$ are required. Since

[TABLE]

these product can be memory efficiently computed by Algorithm 2.

3.2.4 Coarse Grid Solver

On the coarsest grid $g=1$ the V-cycle algorithm requires the exact solution of a linear system with coefficient matrix $A_{1}\in\mathbb{R}^{K_{1}\times K_{1}}$. Since $K_{1}\ll K_{G}$ we assume an explicitly assembled coarse grid coefficient matrix in a sparse matrix format. This allows a factorization of $A_{1}$, that can precomputed and stored in memory. Keeping the factorization in memory is important, since each call of the multigrid preconditioner requires a solution of a linear system given by $A_{1}$. Due to the symmetry of the matrix we apply a sparse Cholesky factorization. Yet, in higher dimensional spaces even the coarse grid operators might be prohibitively memory consuming. Recall that the number of variables $K_{1}$ grows exponentionally with the space dimension $P$. We then apply the matrix-free CG algorithm also as a coarse grid solver. Note that the overall scalability is not affected since the computational costs of the coarse grid solver are fixed, even when the fine grid $G$ is further subdivided as $G\leftarrow G+1$.

3.2.5 V-Cycle

Putting everything together we obtain Algorithm 7, which performs one V-cycle of the multigrid method for the large-scale linear system (2.4) with negligible memory requirement.

The V-cycle can be interpreted as linear iteration with iteration matrix

[TABLE]

where $C_{\text{smooth}}$ denotes the iteration matrix of the utilized smoothing iteration (cf. Saad,, 2003, p. 446). Note that the iteration matrix is only for theoretical investigations and is never assembled in practical applications.

3.2.6 MGCG Method

Finally, applying the multigrid V-cycle as preconditioner for the CG method yields Algorithm 8 as memory efficient multigrid preconditioned conjugated gradient (MGCG) method to solve the large-scale linear system (2.4).

4 Numerical Results

In this section we demonstrate the performance of Algorithm 8 in different numerical test cases. On the one hand, we consider various spatial dimensions $P=1,\dots,4$ in order to investigate the computational complexity of the algorithm. On the other hand, we successively apply uniform grid refinements for a maximum grid level of $G=4,\dots,7$ in each dimension in order to inspect the scalability of the multigrid preconditioner. The underlying codes are programmed within R (version 3.4.4). The algorithmic building blocks, which are critical to computational performance, are accelerated by using the RCPP extension library (cf. Eddelbuettel and François,, 2011; Eddelbuettel,, 2013) and programmed in C++. These are in particular the matrix-vector products with the coefficient matrix (Algorithm 2, 3, 4), interpolation and restriction (Algorithm 2), and the Jacobi smoother (Algorithm 6). Throughout this section, we consider the following test data set

[TABLE]

obtained by uniformly random sampling of the normalized, multivariable sigmoid function

[TABLE]

enriched by normally distributed noise $\varepsilon\sim\mathcal{N}(0;0.1^{2})$, that is

[TABLE]

For $P=1$ and $P=2$ the related sigmoid functions are graphed in Figure 1. We consider these functions since they are irrational functions that show a similar behavior for varying covariate dimensions $P$. However, since the performance of the solution algorithms is in focus, the exact form of the generating function is of minor importance.

In a first test case, we fix the spatial dimension to $P=2$ and successively refine the B-spline basis on $G=1,\dots,7$ grids. We then consider four different tests. In each of them, $G\in\{4,\dots,7\}$ is chosen to be the maximum grid level for the V-cycle of the multigrid preconditioner. The coarsest grid $g=1$ is used for the coarse grid solver in each case. This leads to problem dimensions of $K_{1}=(2^{1}+3)^{2}=25$ on the coarse grid and $K_{G}=(2^{G}+3)^{2}$ on the finest grid. Note that the unpreconditioned CG method simply uses the discretization matrix $A_{G}$ on the $G$-th grid. The results are shown in Figure 2, where a logarithmic scale is used. We observe that the number of unpreconditioned conjugate gradient iterations significantly increases under grid refinements, whereas for both, the Jacobi and SSOR preconditioned MGCG algorithms, the number of iterations is almost constant (i.e. 1-2 for SSOR and 4 for Jacobi). This result illustrates that the multigrid preconditioner enables a scalable solver for the regularized least squares problem (2.4) determining the smoothing spline. Clearly, the computational times are increasing since we run only a single core code. A standard approach here is to distribute the matrix between an increasing number of processors in a cluster computer and one obtains almost constant running times in the sense of weak scalability. Note that this is not achievable by the plain CG solver due to the increasing number of iterations.

The next test under consideration is the comparison of computational complexity with respect to the spatial dimension. Table 2 shows the number of required iterations of the utilized methods in the dimensions $P=1,\dots,4$ together with the CPU time in seconds. The underlying hardware is one core of a Xeon W-2155 processor. One can observe that each of the preconditioned iterations is more expensive with respect to computational time. Yet, the MGCG is significantly faster due to the small number of iterations. Note that the computational time of the SSOR multigrid preconditioned iteration is not comparable since it does not rely on the matrix-free approach. This is only for the comparison of MGCG iterations. Since the SSOR smoother needs access to all matrix entries, we explicitly assemble the coefficient matrix $\left(\Phi{{}^{\prime}}\Phi+\lambda\Lambda\right)$ in R’s sparse format and use it for the algorithm.

Crucial for the convergence speed of the CG method in the unpreconditioned as well as in the preconditioned case is the condition of the respective coefficient matrix, that is

[TABLE]

for the pure CG method and

[TABLE]

for the MGCG method. An explicit form of the iteration matrix of the multigrid method is given in (3.1), i.e. one call of Algorithm 1 is equivalent to a matrix-vector product with $C_{\text{MG},G}$. Note that for the the Jacobi smoother and for the SSOR smoother different preconditioning matrices $C_{\text{MG},G,\text{JAC}}$ and $C_{\text{MG},G,\text{SSOR}}$ arise. Figure 3 visualizes the distribution of the eigenvalues of the corresponding coefficient matrices for the $P=2$ and $G=5$ case. Additionally, Table 3 shows the condition number of the (un)preconditioned iteration matrix for the CG algorithm. Note that we only compute the eigenvalues in the $P=2$ test case due to the computational complexity, since we have to assemble a matrix representation of the preconditiond system matrix. That is, we have to run the V-cycle on each of the $K$-dimensional unit vectors to obtain a matrix on which we then perform an eigenvalue decomposition. Here it can be seen that the multigrid preconditioner pushes the eigenvalues towards 1, which explains the significant lower number of required MGCG iterations. Although it seems that the multigrid preconditioner with SSOR smoother outperforms the Jacobi smoother, it is prohibitively expensive with respect to memory requirement. For the SSOR iteration the entire upper triangular part of the system matrix is required, which is not efficiently accessible in our matrix-free approach. In $P=3$ dimensions, the MGCG with SSOR smoother for the considered test problem requires approximately 30 GB of RAM, which is at the limit of the utilized computer system. In contrast, the matrix-free methods require approximately 16 MB (CG) and 78 MB (MGCG_JAC) of RAM. Due to the fact, that in the matrix free approach we have computational redundancy, the computational times of the full approach are significantly better. Yet, this is bounded to the low dimensional case because of the exponential growth of required memory.

Finally, Figure 4 shows the results of the smoothing spline approximation in two dimensions on the left-hand side and its corresponding residuals to the $100{.}000$ noisy data points on the right-hand side. From this we can see, that the resulting smoothing spline recovers the underlying test function with adequate precision.

5 Conclusions

In this paper, we have presented a memory efficient algorithm in order to determine a smoothing spline in increasing spatial dimensions. An important feature of our approach is the possibility to handle also scattered data in contrast to the already existing methods for gridded data. The main challenge is to deal with memory and computational complexity originating from the large-scale linear system arising from spline smoothing with scattered data sets. In order to overcome the issue of memory requirements, we have initially implemented a matrix-free conjugate gradient method, which comes along without assembling and storing the occurring matrices, but relies solely on matrix-vector products. For this purpose, we especially have exploited the inherent tensor product structure of the multivariable spline functions. Moreover, we address the issue of computational complexity by applying a geometric multigrid preconditioner to the CG algorithm. This enables almost constant iteration numbers for fixed dimensions even under arbitrary grid refinements, which provides an important building block for algorithmic scalability. In a representative numerical case study, we show the applicability and performance of the proposed method.

Acknowledgment

This work has been partly supported by the German Research Foundation (DFG) within the research training group ALOP (GRK 2126).