Fast Cross-validation in Harmonic Approximation

TL;DR

This paper introduces a fast, FFT-based method for computing cross-validation scores in harmonic approximation problems, significantly reducing computational costs across various geometries and node distributions.

Contribution

The authors develop a novel FFT-based approach to efficiently compute cross-validation scores for Tikhonov regularization, applicable to scattered and structured nodes on multiple manifolds.

Findings

Algorithm reduces computational cost to that of the regularization problem

Applicable to nodes on the torus, interval, and sphere

Numerical experiments confirm efficiency and accuracy

Abstract

Finding a good regularization parameter for Tikhonov regularization problems is a though yet often asked question. One approach is to use leave-one-out cross-validation scores to indicate the goodness of fit. This utilizes only the noisy function values but, on the downside, comes with a high computational cost. In this paper we present a general approach to shift the main computations from the function in question to the node distribution and, making use of FFT and FFT-like algorithms, even reduce this cost tremendously to the cost of the Tikhonov regularization problem itself. We apply this technique in different settings on the torus, the unit interval, and the two-dimensional sphere. Given that the sampling points satisfy a quadrature rule our algorithm computes the cross-validations scores in floating-point precision. In the cases of arbitrarily scattered nodes we propose an approximating algorithm with the same complexity. Numerical experiments indicate the applicability of our algorithms.

Figures

Equations

$$\mathbf{F}={\mathbf{F}}_{\mathcal{X},\mathcal{I}}=({\phi }_{\mathbf{n}}(x){)}_{x\in \mathcal{X},\mathbf{n}\in \mathcal{I}},$$

$$\mathbf{f}=(f_x{)}_{x\in \mathcal{X}}=\mathbf{F}\hat{\mathbf{f}}+\mathbf{\epsilon }\in {\mathbb{C}}^{|\mathcal{X}|},$$

$$J_{\lambda }\left(\tilde{\mathbf{f}}\right)={\Vert \mathbf{F}\tilde{\mathbf{f}}-\mathbf{f}\Vert }_{\mathbf{W}}^2+\lambda {\Vert \tilde{\mathbf{f}}\Vert }_{\hat{\mathbf{W}}}^2$$

$$P(\lambda )=\sum _{x\in \mathcal{X}}{\left|{\left[\mathbf{F}{\tilde{\mathbf{f}}}_{\lambda ,(x)}\right]}_x-f_x\right|}^2$$

$${\tilde{\mathbf{f}}}_{\lambda }={\left({\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{F}+\lambda \hat{\mathbf{W}}\right)}^{-1}{\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{f}.$$

$$\nabla J_{\lambda }\left({\tilde{\mathbf{f}}}_{\lambda }\right)=2{\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{F}{\tilde{\mathbf{f}}}_{\lambda }-2{\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{f}+2\lambda \hat{\mathbf{W}}{\tilde{\mathbf{f}}}_{\lambda }\stackrel{!}{=}0.$$

$${\mathbf{f}}_{(x)}=(f_y{)}_{y\in \mathcal{X}\setminus \{x\}}\in {\mathbb{C}}^{|\mathcal{X}|-1}$$

$${\mathbf{F}}_{(x)}=({\phi }_n(y){)}_{y\in \mathcal{X}\setminus \{x\},n\in \mathcal{I}}\in {\mathbb{C}}^{(|\mathcal{X}|-1)\times |\mathcal{I}|}$$

$${\mathbf{W}}_{(x)}=(W_{y,y^{\prime }}{)}_{y,y^{\prime }\in \mathcal{X}\setminus \{x\}}\in {\mathbb{C}}^{(|\mathcal{X}|-1)\times (|\mathcal{X}|-1)}$$

$${\tilde{\mathbf{f}}}_{(x)}={\left({\mathbf{F}}_{(x)}^{\mathsf{H}}{\mathbf{W}}_{(x)}{\mathbf{F}}_{(x)}+\lambda \hat{\mathbf{W}}\right)}^{-1}{\mathbf{F}}_{(x)}^{\mathsf{H}}{\mathbf{W}}_{(x)}{\mathbf{f}}_{(x)}\in {\mathbb{C}}^{|\mathcal{I}|}.$$

$$P(\lambda )=\sum _{x\in \mathcal{X}}{\left|{\left[\mathbf{F}{\tilde{\mathbf{f}}}_{(x)}\right]}_x-f_x\right|}^2$$

$$\frac{\Vert \mathbf{f}-\mathbf{F}{\tilde{\mathbf{f}}}_{{\lambda }^{\star }}{\Vert }_2^2}{{\inf }_{\lambda \ge 0}\Vert \mathbf{f}-\mathbf{F}{\tilde{\mathbf{f}}}_{\lambda }{\Vert }_2^2}\to 1$$

$$\begin{array}{rl}\mathbf{A}& ={\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{F}+\lambda \hat{\mathbf{W}},\\ {\mathbf{A}}_{(x)}& ={\mathbf{F}}_{(x)}^{\mathsf{H}}{\mathbf{W}}_{(x)}{\mathbf{F}}_{(x)}+\lambda \hat{\mathbf{W}}\end{array}$$

$${\mathbf{F}}_{x,:}=({\phi }_n(x){)}_{n\in \mathcal{I}}\in {\mathbb{C}}^{1\times |\mathcal{I}|}$$

$${\mathbf{A}}_{(x)}^{-1}={\mathbf{A}}^{-1}+\frac{{\mathbf{A}}^{-1}{w}_x{\mathbf{F}}_{x,:}^{\mathsf{H}}{\mathbf{F}}_{x,:}{\mathbf{A}}^{-1}}{1-w_x{\mathbf{F}}_{x,:}{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}}.$$

$${\mathbf{A}}_{(x)}=\mathbf{A}-w_x{\mathbf{F}}_{x,:}^{\mathsf{H}}{\mathbf{F}}_{x,:}.$$

$$\mathbf{H}:=\mathbf{F}{\mathbf{A}}^{-1}{\mathbf{F}}^{\mathsf{H}}\mathbf{W}=\mathbf{F}{\left({\mathbf{F}}^{\mathsf{H}}\mathbf{W}\mathbf{F}+\lambda \hat{\mathbf{W}}\right)}^{-1}{\mathbf{F}}^{\mathsf{H}}\mathbf{W}$$

$$P(\lambda )=\sum _{x\in \mathcal{X}}\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x^2}{(1-h_{x,x}{)}^2}$$

$${\tilde{\mathbf{f}}}_{(x)}={\mathbf{A}}_{(x)}^{-1}{\mathbf{F}}_{(x)}^{\mathsf{H}}{\mathbf{W}}_{(x)}{\mathbf{f}}_{(x)}={\mathbf{A}}_{(x)}^{-1}\left(\mathbf{b}-{\mathbf{F}}_{x,:}^{\mathsf{H}}{w}_x{f}_x\right).$$

$${\tilde{\mathbf{f}}}_{(x)}$$

$$=\left({\mathbf{A}}^{-1}+\frac{w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}{\mathbf{F}}_{x,:}{\mathbf{A}}^{-1}}{1-h_{x,x}}\right)\left(\mathbf{b}-{\mathbf{F}}_{x,:}^{\mathsf{H}}{w}_x{f}_x\right)$$

$$=\tilde{\mathbf{f}}+\frac{w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}{f}_x(h_{x,x}-1)+w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}{\mathbf{F}}_{x,:}\tilde{\mathbf{f}}-w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}{h}_{x,x}{f}_x}{1-h_{x,x}}$$

$$=\tilde{\mathbf{f}}+w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}\frac{{\mathbf{F}}_{x,:}\tilde{\mathbf{f}}-f_x}{1-h_{x,x}}$$

$$=\tilde{\mathbf{f}}+w_x{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x}{1-h_{x,x}}.$$

$${\left[\mathbf{F}{\tilde{\mathbf{f}}}_{(x)}\right]}_x-f_x$$

$$={\mathbf{F}}_{x,:}\tilde{\mathbf{f}}+w_x{\mathbf{F}}_{x,:}{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x}{1-h_{x,x}}-f_x$$

$$={\left[\mathbf{F}\tilde{\mathbf{f}}\right]}_x+h_{x,x}\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x}{1-h_{x,x}}-[\mathbf{f}{]}_x$$

$$=[\mathbf{H}\mathbf{f}{]}_x+\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x}{1-h_{x,x}}+[\mathbf{f}-\mathbf{H}\mathbf{f}{]}_x-[\mathbf{f}{]}_x$$

$$=\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x}{1-h_{x,x}}$$

$$\mathbf{r}=\mathbf{F}\tilde{\mathbf{f}}-\mathbf{f}=\mathbf{H}\mathbf{f}-\mathbf{f}.$$

$$h_{x,x}=w_x{\mathbf{F}}_{x,:}{\mathbf{A}}^{-1}{\mathbf{F}}_{x,:}^{\mathsf{H}}$$

$$h=\frac{1}{|\mathcal{X}|}\sum _{x\in \mathcal{X}}{h}_{x,x}=\frac{1}{|\mathcal{X}|}\mathrm{trace}\mathbf{H}$$

$$V(\lambda )=\sum _{x\in \mathcal{X}}\frac{[\mathbf{H}\mathbf{f}-\mathbf{f}{]}_x^2}{(1-h{)}^2}={\left(\frac{|\mathcal{X}|{\Vert \mathbf{H}\mathbf{f}-\mathbf{f}\Vert }_2}{\mathrm{trace}(\mathbf{I}-\mathbf{H})}\right)}^2.$$

$$h_{x,x}<1$$

Full text

Keywords and phrases: cross-validation, regularization, fast evaluation of cross-validation score, discrete Fourier transforms, spherical Fourier transform, NFFT

*2010 AMS Mathematics Subject Classification:

65Txx, 65F22, 42A10. *

1 Introduction

Estimating a good regularization parameter is a frequent problem in approximation, statistics, and inverse problems. In this paper we restrict ourselves to the basic example of approximating a function from discrete function values. To make the setting concrete we fix for a finite index set $\mathcal{I}$ a family of basis functions $\varphi_{\bm{n}}\colon X\to\mathbb{C}$, $\bm{n}\in\mathcal{I}$ on a domain $X\subset\mathbb{R}^{d}$. Given a finite set of nodes $\mathcal{X}\subset X$ and the corresponding Fourier matrix

[TABLE]

we consider the problem of recovering Fourier coefficients $\bm{\hat{f}}=(\hat{f}_{\bm{n}})_{\bm{n}\in\mathcal{I}}\in\mathbb{C}^{|\mathcal{I}|}$ from noisy data

[TABLE]

where $\bm{\varepsilon}$ is zero mean Gaussian noise. More specifically, we look for minimizers of the Tikhonov functional

[TABLE]

and ask for the optimal regularization parameter $\lambda>0$, where $\|\bm{f}\|_{\bm{W}}^{2}={\bm{f}}^{\mathsf{H}}\bm{W}\bm{f}$ and $\|\bm{\hat{f}}\|_{\bm{\hat{W}}}^{2}=\bm{\hat{f}}^{\mathsf{H}}\bm{\hat{W}}\bm{\hat{f}}$ for the strictly positive diagonal weight matrices $\bm{W}\in\mathbb{R}^{|\mathcal{X}|\times|\mathcal{X}|}$ and $\bm{\hat{W}}\in\mathbb{R}^{|\mathcal{I}|\times|\mathcal{I}|}$ in space, respectively frequency domain.

Because of its importance to many practical problems there is a vast literature on many different strategies to determine an optimal regularization parameter $\lambda$, e.g. [8, 6, 15, 38]. The idea of so called cross-validation methods is to divide the set of nodes into a subset used to compute the approximation and a subset used for validating the goodness of fit. This procedure might be repeated for different splittings and eventually results for a fixed regularization parameter $\lambda$ in a cross-validation score. By minimizing this score an “optimal” regularization parameter is found. This approach was initially proposed by Golub, Heath and Wahba [15] in the setting of smoothing splines and since then has been applied to a wide range of parameter estimating problems.

In this paper we consider “leave-one-out” cross-validation, i.e., for fixed regularization parameter $\lambda$ and any node $x\in\mathcal{X}$ we compute the minimizer $\bm{\tilde{f}}_{\lambda,(x)}$ of the functional (1.1) restricted to the set of nodes $\mathcal{X}\setminus\{x\}$ and use

[TABLE]

as cross-validation score. A drawback with a purly data-driven regularization method like this is, that there is no guarantee for a good approximation of the solution of the regularized problem as formulated in the Bakushinskiĭ veto, [3]. Another difficulty of this approach is its numerical complexity. Indeed, computing $P(\lambda)$ for a single value of $\lambda$ requires solving the minimization problem (1.1) $\left|\mathcal{X}\right|$ times, which is too expensive for most applications. For spline interpolation on the interval or in higher dimensional domains different algorithms have been proposed to lower the computational costs. Those include Monte Carlo approximations [9], matrix decomposition methods [50, 43] and Krylow space methods [35].

For the specific setting of Fourier approximation on the torus ${\mathbb{T}^{d}}$ at regular lattice points a fast algorithm has been proposed by Tasche and Weyrich [45] which requires to solve the minimization problem only once for each regularization parameter $\lambda$. The idea of this paper is to generalize the approach in [45] to arbitrary sampling nodes and other domains like the unit interval or the two-dimensional sphere.

The paper is organized as follows. In the second chapter we introduce the necessary notations and prove in Theorem 2.5 a representation of the cross-validation score that depends only on one solution of (1.1), but includes the diagonal entries $h_{x,x}$ of the so called hat matrix $\bm{H}=\bm{F}(\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}})^{-1}\bm{F}^{\mathsf{H}}\bm{W}$. The efficient approximate computation of those diagonal entries for different settings is subject to the remaining chapters.

Our essential requirement for the exact fast computation of the diagonal entries $h_{x,x}$ is that the nodes $\mathcal{X}$ together with the weights $\bm{W}$ form an exact quadrature rule. This requirement is satisfied for regular tensor product grids and rank-1 lattices on the $d$-dimensional torus, Chebyshev nodes on the interval $[-1,1]$, and, e.g., Gauss Legendre nodes on the two-dimensional unit sphere $\mathbb{S}^{2}$. The corresponding formulae for the diagonal entries $h_{x,x}$ are given in the Theorems 3.2, 4.2, and 5.1. Together with fast Fourier algorithms on the torus [29, 27], for rank-1 lattices [24, 25], on the interval [11, 27], and on the sphere [28, 27] this allows the efficient evaluation of the cross-validation score $P(\lambda)$ with a numerical complexity close to $\mathcal{O}(|\mathcal{I}|+|\mathcal{X}|)$. Numerical examples for all these settings illustrate our findings.

In the case that no exact quadrature rule is known for the given interpolation nodes we suggest approximating them by the volume of the corresponding Voronoi cells. Our numerical tests in Section 3.4, 4.3, and 5.2 indicate a good approximation of the true cross-validation score, which is much more expensive to compute. The Matlab code of our algorithm as well as for all numerical experiments can be found on the GitHub repository https://github.com/felixbartel/fcv.

2 Cross-validation

Lets start this section by reminding that the minimizer of the Tikhonov functional (1.1) can be given explicitly.

Lemma 2.1**.**

The unique Tikhonov minimizer of (1.1) is

[TABLE]

Proof.

We look for stationary points by calculating the roots of the gradient of $J_{\lambda}$

[TABLE]

Because $\bm{F}^{\mathsf{H}}\bm{F}$ is positive semidefinite, $\bm{W},\bm{\hat{W}}$, and $\lambda$ are strictly positive we find that $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}}$ is positive definite. In particular it is invertible such that the stationary point can be written as $\bm{\tilde{f}}_{\lambda}=(\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}})^{-1}\bm{F}^{\mathsf{H}}\bm{W}\bm{f}$. Using the positive definiteness we see that $\bm{\tilde{f}}_{\lambda}$ fulfills the required minimizing property. ∎

Remark 2.2**.**

Since for random nodes $\mathcal{X}$ and many important function systems the matrix $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}$ is invertible with probability one, see e.g. [4], we may relax the assumption on the frequency weights $\hat{w}_{n}$ to be only non negative. Especially, the zero order frequency weight is often set to zero to avoid penalizing constant functions.

As the leave-one-out cross-validation score depends on solving (2.1) for sets of nodes of the form $\mathcal{X}\setminus\{x\}$ we introduce the following notations for omitting a single node $x\in\mathcal{X}$. For $x\in\mathcal{X}$ and $\bm{f}\in\mathbb{C}^{|\mathcal{X}|}$ we denote by

[TABLE]

the vector of function values $\bm{f}$ with one node $x\in\mathcal{X}$ omitted. Accordingly, we denote by

[TABLE]

the Fourier matrix $\bm{F}$ with the row corresponding $x\in\mathcal{X}$ omitted and by

[TABLE]

the restriction of the spatial weight matrix $\bm{W}$ to the set of nodes $\mathcal{X}\setminus\{x\}$. With these notations the minimizer of the Tikhonov functional (1.1) reduced to the nodes $\mathcal{X}\setminus\{x\}$ can be written as

[TABLE]

Definition 2.3**.**

The ordinary cross-validation score for the Tikhonov functional (1.1) is defined as

[TABLE]

where $\bm{\tilde{f}}_{(x)}$ is defined by (2.2) and $[\bm{F}\bm{\tilde{f}}_{(x)}]_{x}$ denotes the entry of $\bm{F}\bm{\tilde{f}}_{(x)}$ corresponding to the node $x\in\mathcal{X}$.

Interpreting (2.3), we are comparing the predicted or smoothened value $[\bm{F}\bm{\tilde{f}}_{(x)}]_{\bm{x}}$ with the noisy data $f_{\bm{x}}$ for each node. They intuitively differ more in the case of under- or oversmoothing. So its minimum is a candidate for the smoothing parameter $\lambda$. Indeed, probabilistic optimality results of the form

[TABLE]

in probability for $\lambda^{\star}$ being the minimum of the cross-validation score have been shown for $|\mathcal{X}|\to\infty$ in [34] under the assumption of homogenious noise, i.e., $\varepsilon_{x}$ has the same variance for all $x\in\mathcal{X}$. In Theorem 3.5 of [20], weights were incoporated and the condition was generalized to the case where all $\sqrt{w_{x}\varepsilon_{x}}$ have a common variance for $x\in\mathcal{X}$. These optimality results are with respect to a LOSS function similar to the data-fitting term in the Tikhonov functional (1.1). However, in our numerical experiments, we will have a look at the unweighted $L_{2}$-error as we use the weights $w_{x}$ for numerical purposes and want to analyze how this disstorts the results. Unflattering is the fact that, for a single regularization parameter $\lambda$, the direct computation of the ordinary cross-validation score requires to solve $|\mathcal{X}|$ times the normal equation (2.1).

Our first goal is to relate the solution of the reduced problem (2.2) to the solution of the full problem (2.1). To this end we define the matrices

[TABLE]

which are decisive for the computation of $\bm{\tilde{f}}$ and $\bm{\tilde{f}}_{(x)}$, respectively, and show the following relationship between their inverse, cf. [15].

Lemma 2.4**.**

Let $\bm{A}$ and $\bm{A}_{(x)}$ be defined as in (2.4) and

[TABLE]

denote the row of the matrix $\bm{F}$ which corresponds to the node $x\in\mathcal{X}$. Then we have

[TABLE]

Proof.

The assertion of the lemma follows immediately by applying the Sherman-Morrison formula to

[TABLE]

∎

Our next goal is to make the repetitive solving of the normal equation (2.2) in (2.3) independent of the right-hand side $\bm{f}$. To this end we define the so called hat matrix

[TABLE]

which when applied to a data vector $\bm{f}$ solves the normal equation (2.1) and evaluates the resulting function at the nodes $\mathcal{X}$. The next lemma is a generalization of [15, equation 2.2] and [49, equation 4.2.9], and shows that for the computation of (2.3), with given diagonal entries $h_{x,x}$ of the hat matrix $\bm{H}$, it is sufficient to solve the normal equation (2.1) with respect to the data vector $\bm{f}$ only once.

Theorem 2.5**.**

The ordinary cross-validation score (2.3) can be written as

[TABLE]

with $h_{x,x}$, $x\in\mathcal{X}$ being the diagonal entries of the hat matrix $\bm{H}$ defined in (2.6).

Proof.

Let $\bm{b}=\bm{F}^{\mathsf{H}}\bm{W}\bm{f}$. Then

[TABLE]

Next we apply Lemma 2.4 and observe that the denominator in (2.5) can be expressed in terms of the diagonal entries $h_{x,x}$ of the hat matrix $\bm{H}$:

[TABLE]

Multiplying with $\bm{F}$ from the left-hand side and subtracting $f_{x}$ results in

[TABLE]

and hence each summand in (2.3) is equal to the corresponding summand in (2.7). ∎

Remark 2.6**.**

According to Theorem 2.5 the ordinary cross-validation score is nothing more than the weighted norm of the residue

[TABLE]

Although this means that the normal equation (2.1) has to be solved only once with respect to the data vector $\bm{f}$ the most expensive part remains, namely the computation of the diagonal entries

[TABLE]

for $x\in\mathcal{X}$, which again requires repetitive solving of the normal equation.

Replacing the diagonal entries $h_{x,x}$ with their mean value

[TABLE]

we obtain the so called generalized cross-validation score, cf. [49, section 4.3].

Definition 2.7**.**

The generalized cross-validation score is defined as

[TABLE]

Obviously, if all diagonal entries $h_{x,x}$ of $\bm{H}$ coincide we have $P(\lambda)=V(\lambda)$.

Lemma 2.8**.**

The diagonal elements $h_{x,x}$ of the hat matrix $\bm{H}$ satisfy

[TABLE]

for all $\lambda>0$ and $x\in\mathcal{X}$.

Proof.

Since $\bm{F}_{(x)}^{\mathsf{H}}\bm{W}_{(x)}\bm{F}_{(x)}$ is positive semidefinite and $\lambda\bm{\hat{W}}$ is strictly positive definite we see that $\bm{A}_{(x)}=\bm{F}_{(x)}^{\mathsf{H}}\bm{W}_{(x)}\bm{F}_{(x)}+\lambda\bm{\hat{W}}$ is invertible. On the other hand we know by the Sherman-Morrison formula that $\bm{A}_{(x)}=\bm{A}-w_{x}\bm{F}_{x,:}^{\mathsf{H}}\bm{F}_{x,:}$ is invertible if and only if $w_{x}\bm{F}_{x,:}A^{-1}\bm{F}_{x,:}^{\mathsf{H}}\neq 1$ . Therefore

[TABLE]

Since the minimizer $\bm{\tilde{f}}_{\lambda}$ of (1.1) converges to the zero vector as $\lambda\to\infty$, we obtain for $\bm{f}=\bm{e}_{x}$ and $\lambda\to\infty$

[TABLE]

Together with the fact that the diagonal entries $h_{x,x}$ depend continuously on $\lambda$ this proves the assertion. ∎

2.1 Algorithm to compute the ordinary or generalized cross-validation score

Concluding the previous statements we end up with a sheme to compute the cross-validation scores.

Remark 2.9**.**

For computing the Tikhonov-minimizer of (1.1) one can use the LSQR method for numerical stability. This can be accomplished with the coefficient matrix

[TABLE]

and the right-hand side

[TABLE]

where $\bm{0}$ is a column vector containing $|\mathcal{I}|$ zeros. The resulting system of equations $(\bm{M}^{\mathsf{H}}\bm{M})^{-1}\bm{M}^{\mathsf{H}}\bm{b}$, which the LSQR method solves, is equivalent to (2.1).

The computationally most expensive part of Algorithm 1 is the computation of the values $\bm{\tilde{f}}$ and $h_{x,x}$ for all $x\in\mathcal{X}$. In the subsequent sections we discuss some specific settings to speed up the process and propose an approximation of the ordinary and the generalized cross-validation score in more general cases.

3 Cross-validation on the torus

In this section, we seek for fast algorithms to compute the cross-validation score on the $d$-dimensional torus ${\mathbb{T}^{d}}$ with respect to the Fourier basis $\{\mathrm{e}^{2\pi\mathrm{i}\bm{n}\cdot\bm{x}}\}_{\bm{n}\in\mathbb{Z}^{d}}$ in $L_{2}({\mathbb{T}^{d}})$. With this setting the Fourier matrix $\bm{F}$ becomes

[TABLE]

for a finite node set $\mathcal{X}\subset{\mathbb{T}^{d}}$, a finite multi-index set $\mathcal{I}\subset\mathbb{Z}^{d}$ and $\bm{n}\cdot\bm{x}$ the Euclidean inner product. So $\mathcal{I}$ determines all possible frequencies and $\mathcal{X}$ the nodes of the transform. For the specific case of equispaced nodes $\mathcal{X}$ fast algorithms have been reported in [45]. In fact, our approach in this section can be seen as a generalization of [45] to more general sampling sets and leads to the same algorithm for equispaced data.

Our central goal is to find a simpler expression for the diagonal entries of the hat matrix $\bm{H}=\bm{F}(\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}})^{-1}\bm{F}^{\mathsf{H}}\bm{W}$ that allows us to apply Theorem 2.5 efficiently. The idea is to choose $\bm{W}$ such that $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}$ has diagonal form because the inverse of $\bm{A}=\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}}$ could then be calculated entry-wise.

3.1 Exact Quadrature

The first approach is to use quadrature rules. Because they are not limited to the torus we define them for general measure spaces so we can make use of them in subsequent sections.

Definition 3.1**.**

Let $(\mathcal{M},\mu)$ be a measure space and $\mathcal{P}\subset L_{1}(\mathcal{M})$ a set of integrable functions. We call a set of nodes $\mathcal{X}\subset\mathcal{M}$ and a list of weights $\bm{W}=\operatorname*{diag}(w_{\bm{x}})_{\bm{x}\in\mathcal{X}}$ an exact quadrature rule for $\mathcal{P}$, if for all $f\in\mathcal{P}$ we have

[TABLE]

For the tours we obtain the following

Theorem 3.2**.**

Let $\mathcal{I}\subset\mathbb{Z}^{d}$ be a finite multi-index set with Fourier weights $\bm{\hat{W}}=\operatorname*{diag}(\hat{w}_{\bm{n}})_{\bm{n}\in\mathcal{I}}$, $\mathcal{X}\subset{\mathbb{T}^{d}}$ a set of nodes with $\bm{W}$ their corresponding quadrature weights such that $(\mathcal{X},\bm{W})$ forms a quadrature rule which is exact for all trigonometric polynomials $\mathrm{e}^{2\pi\mathrm{i}\bm{n}\cdot}$ with frequencies $\bm{n}$ in the difference set $\mathcal{D}(\mathcal{I})\coloneqq\{\bm{n}_{1}-\bm{n}_{2}:\bm{n}_{1},\bm{n}_{2}\in\mathcal{I}\}$. Then

1. (i)

the inverse of $\bm{A}$, given in (2.4), where $\bm{F}$ is the Fourier matrix (3.1) on ${\mathbb{T}^{d}}$ is

[TABLE] 2. (ii)

the diagonal entry corresponding to the node $\bm{x}\in\mathcal{X}$ of the hat matrix $\bm{H}=\bm{F}\bm{A}^{-1}\bm{F}^{\mathsf{H}}\bm{W}$ becomes

[TABLE]

Proof.

Since the product of two exponential functions supported on the frequency set $\mathcal{I}$ has only frequencies in $\mathcal{D}(\mathcal{I})$, where the quadrature nodes and weights are exact, we have

[TABLE]

and, hence,

[TABLE]

This implies (i). For (ii) we compute the diagonal entries of $\bm{H}$ as

[TABLE]

∎

Corollary 3.3**.**

With the prerequisites of Theorem 3.2 we can compute $P(\lambda)$ and $V(\lambda)$ by Algorithm 1 in the same complexity as multiplying a vector with $\bm{F}$ or $\bm{F}^{\mathsf{H}}$ for a fixed regularization parameter $\lambda$.

Remark 3.4**.**

Theorem 3.2 requires $(\mathcal{X},\bm{W})$ to form an exact quadrature rule which directly alters the weighting of the data-fitting term in the underlying Tikhonov functional (1.1). Preferebly, one wants to have the weight $w_{x}$ proportional to the level of noise $\varepsilon_{x}$ for every node $x\in\mathcal{X}$, cf. [40]. In most of the following examples the weights as well as the Gaussian noise are uniform, thus fulfulling this property.

Case studies for specific exact quadrature rules on the torus are discussed in the following two subsections.

3.2 Equispaced Nodes

The simplest example of quadrature on the torus $\mathbb{T}^{d}$ is Gauss quadrature which consists of $N^{d}$ equispaced nodes

[TABLE]

with uniform weights $w_{\bm{x}}=N^{-d}$. The resulting quadrature formula is exact for all trigonometric polynomials supported on $\mathcal{I}_{2N}\coloneqq\mathbb{Z}^{d}\cap\prod_{t=1}^{d}[-N,N)=\mathcal{D}(\mathcal{I}_{N})$. Thus we can apply Theorem 3.2 for $\mathcal{X}$ and $\mathcal{I}=\mathcal{I}_{N}$. The corresponding Fourier matrix $\bm{F}=\bm{F}_{\mathcal{X},\mathcal{I}_{N}}$ describes the ordinary discrete Fourier transform for which the matrix-vector product can be computed in $\mathcal{O}(N^{d}\log N)$ with the fast Fourier transform.

For $d=1,2$ our algorithm coincides with the algorithm proposed in [45] with the only difference that the authors evaluated the data fitting term in the frequency domain and used specific Fourier weights $w_{n}=n^{4}$ as regularization term.

In order to illustrate our approach we chose as the test function $f$ the peaks function from Matlab, which is a sum of translated and scaled Gaussian bells. We evaluated this function on a grid of $1024\times 1024$ equispaced nodes $\mathcal{X}$ and corrupted the data by 10% Gaussian noise $\varepsilon_{\bm{x}}$, i.e., we set

[TABLE]

for all $\bm{x}\in\mathcal{X}$ as depicted in Figure 3.1, (a).

As regularization term we fixed isotropic Sobolev weights $\hat{w}_{\bm{n}}=1+\|\bm{n}\|_{2}^{s}$ for $\bm{n}\in\mathcal{I}_{N}$ and $s=3$ in Fourier space, which correspond to a function with $3$ derivatives in $L_{2}({\mathbb{T}^{d}})$. Varying the regularization parameter $\lambda\in[2^{-18},2^{-8}]$ we computed the Tikhonov minimizers $\bm{\tilde{f}}_{\lambda}$ according to (2.1). We then applied Parseval to the original $\bm{\hat{f}}$ and $\bm{\tilde{f}}_{\lambda}$ which is a byproduct from Algorithm 1 to compute the $L_{2}({\mathbb{T}^{d}})$-approximation errors as a function of the regularization parameter $\lambda$. The resulting curve is depicted in Figure 3.1, (b). Note that according to (3.2) all diagonal entries of the hat matrix are equal and, hence, the ordinary cross-validation score coincides with the generalized cross-validation score. The reconstruction $\bm{\tilde{f}}_{\lambda}$ with respect to the minimizer of the cross-validation score $P(\lambda)$ is depicted in Figure 3.1, (a).

In Figure 3.1, (b) the actual $L_{2}({\mathbb{T}^{d}})$-approximation error is compared to the cross-validation score $P(\lambda)$ computed according to Algorithm 1 with use of the fast Fourier transform. We observe that the minimizers of both functionals coincide surprisingly well. For this numerical experiment the average running time for the evaluation of $P(\lambda)$ for a single value of $\lambda$ was $0.06$ seconds with the fast algorithm and more than $14$ hours for a direct implementation of (2.3).

3.3 Rank-1 Lattices

The approximation of high-dimensional multivariate periodic functions by trigonometric polynomials using particular finite index sets $\mathcal{I}$ in frequency domain is possible using special index sets [46, 10] on the domain $\mathcal{X}$. The most efficient method uses samples along rank-1 lattices and is based on a simple univariate FFT [26]. Rank-1 lattices are defined by

[TABLE]

where $M\bm{1}=(M,\dots,M)^{\mathsf{T}}\in\mathbb{Z}^{d}$. They are fully characterized by the generating vector $\bm{z}\in\mathbb{Z}^{d}$ and the lattice size $M$. There exist algorithms which, given a frequency index set $\mathcal{I}$ and $M$, compute a generating vector $\bm{z}$ such that $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}$ equals the identity matrix for $\bm{W}=\operatorname*{diag}(1/M)_{\bm{x}\in\mathcal{X}}$, cf. [26, 24, 39]. The advantage of rank-1 lattices is the variable index set $\mathcal{I}$ instead of the tensor-product approach like in Section 3.2. So depending on the function we can adapt to different decay properties of the Fourier coefficients. Furthermore there exist fast algorithms which evaluate the matrix-vector product with $\bm{F}$ or $\bm{F}^{\mathsf{H}}$ in $\mathcal{O}(M\log M+d|\mathcal{I}|)$ using only one one-dimensional fast Fourier transform.

To exemplify these ideas we looked at a sample function consisting of a tensor product of $L_{2}({\mathbb{T}^{d}})$-normed B-splines of order two in seven dimensions, i.e., $d=7$,

[TABLE]

where $\mathcal{X}_{A}$ denotes the indicator function. The Fourier coefficients of $f$ are

[TABLE]

Therefore the Fourier coefficients decay like $\mathcal{O}(\prod_{j=1}^{d}n_{j}^{-2})$ and a candidate for an index set $\mathcal{I}$ would be a $d$-dimensional hyperbolic cross

[TABLE]

In particular we used a radius of $N=16$, a reconstructing rank-1 lattice $\mathcal{X}$ from [26, table 6.2] with $M=1\,105\,193$ nodes and set the weights in Fourier space to $\hat{w}_{\bm{n}}=\prod_{j=1}^{d}\max(|n_{j}|^{2},1)$.

Applying Algorithm 1 to $f_{\bm{x}}=f(\bm{x})+\varepsilon$, $\bm{x}\in\mathcal{X}$, where $\varepsilon$ denotes 5% Gaussian noise we computed the cross-validation scores $P(\lambda)=V(\lambda)$ for $\lambda\in[2^{-9},2^{0}]$. Again both scores coincide since the diagonal entries (3.2) of the hat matrix are multiples of the constant weights $w_{\bm{x}}$ in spatial domain. For the multiplications with $\bm{F}$ and $\bm{F}^{\mathsf{H}}$ we made use of fast rank-1 lattices Fourier transforms. A comparison of the actual $L_{2}({\mathbb{T}^{d}})$-approximation error with the cross-validation score is plotted in Figure 3.2. We observe that the optimal $\lambda$ with respect to the $L_{2}({\mathbb{T}^{d}})$-error and the $\lambda$ chosen by cross-validation are very close in this example.

3.4 Approximative quadrature

In the case of scattered data approximation exact quadrature rules are typically not available. In principle, one can compute exact quadrature rules by determining the weighs $\bm{W}=\operatorname*{diag}(w_{\bm{x}})_{\bm{x}\in\mathcal{X}}$ as a solution of the linear system

[TABLE]

where $\bm{e}_{\bm{0}}$ is the vector which only contains zeros, except in the $\bm{0}$-th position where it is equal to one. These weights can be guaranteed to be non negative under certain conditions on the frequency index set $\mathcal{I}$ and the mesh norm

[TABLE]

of the nodes $\mathcal{X}$, cf. [12]. However, those conditions strongly restrict the polynomial degree and do not guaranty the quadrature weights to be non-oscillating. This may be problematic, since the weights directly alter the problem (1.1) we want to solve.

Our idea is to replace the exact quadrature weights $\bm{W}$ by approximative weights derived from the Voronoi tessellation of the node set $\mathcal{X}$.

Definition 3.5**.**

Let $\mathcal{M}$ be a Riemannian manifold with a distance function $\operatorname*{dist}(\cdot,\cdot)$. For a set of nodes $\mathcal{X}\subset\mathcal{M}$ we define the Voronoi cell $V_{\bm{x}}$ corresponding to $\bm{x}\in\mathcal{X}$ by

[TABLE]

The Voronoi weight $w_{\bm{x}}$ is the area of the Voronoi cell $V_{\bm{x}}$

[TABLE]

To emphasize the choice of the Voronoi weights as approximative quadrature weights we make a rough error estimate for the approximated quadrature using the Voronoi weights.

Theorem 3.6**.**

Let $f\colon\mathcal{M}\to\mathbb{C}$ be Lipschitz continuous with the Lipschitz constant $L$. Let further $\mathcal{X}$ be a set of nodes with mesh norm

[TABLE]

and Voronoi weights $w_{\bm{x}}$. Then

[TABLE]

Proof.

Since the disjoint union of all Voronoi cells $V_{\bm{x}}$ is $\mathcal{M}$ itself we can decompose the integral as follows

[TABLE]

Noting that the maximal distance of $\bm{x}$ to any other node of the corresponding Voronoi cell $V_{\bm{x}}$ cannot exceed $\delta_{\mathcal{X}}$, we use the Lipschitz continuity to estimate the leftover integrand

[TABLE]

∎

Remark 3.7**.**

1. (i)

Theorem 3.6 states that the error of the quadrature formula gets small for smooth functions in the sense of a small Lipschitz constant and for small mesh norms, like for approximately equidistributed nodes. 2. (ii)

Deterministic and probabilistic error estimates are available from **[19]** and **[4]**, respectively. 3. (iii)

The Voronoi decomposition is dual to the Delaunay triangulation and thus can be computed in $\mathcal{O}(|\mathcal{X}|\log|\mathcal{X}|)$ for the Euclidean distance in $\operatorname*{dist}(\cdot,\cdot)$.

Given that the error of the Voronoi quadrature is small we obtain approximately

[TABLE]

where $\bm{I}$ denotes the identity matrix. Inserting this into the definition of the hat matrix $\bm{H}$ we have formally

[TABLE]

Analogously to Theorem 3.2 we obtain for the diagonal entries $\tilde{h}_{\bm{x},\bm{x}}$ of the approximated hat matrix $\bm{\tilde{H}}$,

[TABLE]

Together with Theorem 2.5 and Definition 2.7 this motivates the following definition of approximated cross-validation scores.

Definition 3.8**.**

The approximated cross-validation score $\tilde{P}(\lambda)$ and the approximated generalized cross-validation score $\tilde{V}(\lambda)$ are defined by

[TABLE]

where $\tilde{h}=\frac{1}{|\mathcal{X}|}\sum_{\bm{x}\in\mathcal{X}}\tilde{h}_{\bm{x},\bm{x}}$.

Remark 3.9**.**

Algorithm 1 is easily modified to compute the approximated scores by replacing all occurrences of $h_{\bm{x},\bm{x}}$ by $\tilde{h}_{\bm{x},\bm{x}}$. The computationally most expensive part remains the computation of the Tikhonov minimizer $\bm{\tilde{f}}=\bm{H}\bm{f}=\bm{F}(\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}})^{-1}\bm{F}^{\mathsf{H}}\bm{W}\bm{f}$. Making use of the NFFT [29, 27] the matrix-vector multiplications with $\bm{F}$ and $\bm{F}^{\mathsf{H}}$ can be performed with $\mathcal{O}(|\mathcal{I}|\log|\mathcal{I}|+|\mathcal{X}|)$ numerical operations.

In the remainder of this sections we illustrate that the approximated cross-validation scores can be computed drastically faster while providing a good approximation to the minimizer of the actual cross-validation score. To this end we considered scattered sampling points on the one-dimensional torus $\mathbb{T}$ as well as on the two-dimensional torus $\mathbb{T}^{2}$. In order generate sufficiently nonuniform sampling points we drew random samples with respect to the uniform distribution on $[0,1]$ and $[0,1]^{2}$ and squared them. This leads to sampling sets that are more dense towards the point [math] and the edges $0\times[0,1]$ and $[0,1]\times 0$.

In the one-dimensional example we used $|\mathcal{X}|=128$ sampling points and the index set $\mathcal{I}_{64}^{1d}=\{-32,\dots,31\}$. In the two-dimensional example we chose $|\mathcal{X}|=8192$ and $\mathcal{I}_{64}^{2d}=\mathcal{I}_{64}^{1d}\times\mathcal{I}_{64}^{1d}$. In both cases this corresponds to an oversampling factor of two. As in Subsection 3.2 we chose as a test function the Matlab peaks function with fixed second argument zero in the one-dimensional case. Finally, we added $5\%$ Gaussian noise as depicted in Figure 3.3, (a) and Figure 3.4, (a).

As in both cases the weights $w_{\bm{x}}$ are far from uniform we may expect some difference between the ordinary cross-validation score $P(\lambda)$ and the generalized cross-validation score $V(\lambda)$. This can be observed in the one-dimensional example, cf. Figure 3.3. In the two-dimensional example both functionals coincide surprisingly well, cf. Figure 3.4. Furthermore, the weights are completely uncorrelated to the homogeneous noise level and, thus, contradict with Remark 3.4.

Judging the approximation of the exact cross-validation scores $P(\lambda)$ and $V(\lambda)$ by the approximated scores $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$ we observe that for small regularization parameters $\lambda$ the score $\tilde{P}(\lambda)$ contains several peaks for both examples. These peaks occur because we overestimate the diagonal entries $h_{x,x}$ such that they attain values around one and summands of the ordinary cross-validation score (2.7) diverge. In contrast Lemma 2.8 says that these diagonal entries are always smaller than one. Nevertheless, the minimizer of all four functionals $P,\tilde{P},V,\tilde{V}$ are very close together for the two-dimensional example while for the one-dimensional example the minimizer of $P$ and $\tilde{P}$ are closer to the $L_{2}({\mathbb{T}^{d}})$-optimal regularization parameter compared to the minimizer of $V$ and $\tilde{V}$. A natural idea to avoid the oscillatory regions of the functional $\tilde{P}$ would be to use the minimizer of $\tilde{V}$ as initial guess for minimizing $\tilde{P}$.

The central reason for preferring the functional $\tilde{P}$ and $\tilde{V}$ over the functionals $P$, $V$ is that they are faster to compute. Indeed, if we fix the number of iterations for computing the Tikhonov minimizer, $P(\lambda)$ and $V(\lambda)$ can be acquired in $\mathcal{O}(|\mathcal{I}||\mathcal{X}|\log|\mathcal{I}|+|\mathcal{X}|^{2})$ numerical operations for a single regularization parameter $\lambda$, which compares to $\mathcal{O}(|\mathcal{I}|\log|\mathcal{I}|+|\mathcal{X}|)$ numerical operations for the evaluation of $\tilde{P}$ and $\tilde{V}$. In our toy example the computation of $P$ and $V$ took approximately $1.61$ for the one-dimensional and $1278$ seconds for the two-dimensional case, while the computation of $\tilde{P}$ and $\tilde{V}$ was performed within $0.02$ and $0.16$ seconds averaged over all tested $\lambda$.

4 Cross-validation on the unit interval

In this section, we consider the cross-validation scores for nonperiodic functions on the unit interval $[-1,1]$ with respect to the Chebyshev polynomials

[TABLE]

up to polynomial degree $N\in\mathbb{N}$. In this setting the Fourier matrix becomes

[TABLE]

for a set of nodes $\mathcal{X}$.

4.1 Exact Quadrature

Similarly as for functions on the torus we consider the case of exact quadrature first. Therefore we remind that the Chebyshev polynomials are orthogonal with respect to the inner product

[TABLE]

and are normalized such that

[TABLE]

Assuming that the nodes $\mathcal{X}$ and weights $\bm{W}$ form a quadrature rule that is exact up to polynomial degree $2N-2$ the diagonal entries of the hat matrix $\bm{H}$ can be given expicitly using the following theorem.

Theorem 4.1**.**

Let the nodes $\mathcal{X}$ and the weights $\bm{W}=\operatorname*{diag}(w_{x})_{x\in\mathcal{X}}$ form a quadrature rule which is exact up to polynomial degree $2N-2$ and let $\bm{\hat{W}}=\operatorname*{diag}(\hat{w}_{0},\dots,\hat{w}_{N-1})$ be the weights in frequency domain. Then the diagonal entries $h_{x,x}$ of the hat matrix $\bm{H}=\bm{F}(\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}})^{-1}\bm{F}\bm{W}$ corresponding to the nodes $x\in\mathcal{X}$ satisfy

[TABLE]

Proof.

Similar to Theorem 3.2 we obtain $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}+\lambda\bm{\hat{W}}=\operatorname*{diag}(\pi+\lambda\hat{w}_{0},\pi/2+\lambda\hat{w}_{1},\dots,\pi/2+\lambda\hat{w}_{N-1})$. Putting this into the formula for the diagonal elements of the hat matrix obtain

[TABLE]

In combination with the addition theorem $\cos(2x)=2\cos^{2}x-1$ this proves the assertion. ∎

4.2 Chebyshev nodes

The most basic examples of an exact quadrature formula on the interval is probably quadrature at Chebyshev nodes. In order to restate Theorem 4.1 for this case we require the discrete cosine transforms from first up to third kind

[TABLE]

with $\bm{\gamma}(0)=\bm{\gamma}(N)\coloneqq\sqrt{2}/2$ and $\bm{\gamma}(n)\coloneqq 1$ for $n=1,\dots,N-1$. The corresponding matrix-vector products can be calculated using $\mathcal{O}(N\log N)$ arithmetic operations, cf. [39, Chapter 6].

Using the fact that $\bm{C}_{N}^{\mathrm{III}}$ is orthonormal, cf. [39, Sec. 3.5], we acquire

[TABLE]

If we multiply with $\operatorname*{diag}(\sqrt{\pi},\sqrt{\pi/2},\dots,\sqrt{\pi/2})$ from both sides we obtain

[TABLE]

Putting this into the form $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}$ we see that the Chebyshev nodes of first kind

[TABLE]

and the uniform weights $w_{x}=\pi/N$ form a quadrature rule which is exact up to degree $2N-2$.

For these specific nodes Theorem 4.1 simplifies to:

Theorem 4.2**.**

Let $\mathcal{X}=\{x_{1},\dots,x_{m}\}$ be the Chebyshev nodes of first kind and $w_{x_{m}}=\pi/N$. Then the diagonal entries $h_{x_{m},x_{m}}$ of the hat matrix (2.6) can be written as

[TABLE]

with

[TABLE]

Proof.

According to Theorem 4.1 we have

[TABLE]

Using the coefficients $\bm{b}$ the first sum can be expressed with twice the frequency

[TABLE]

which is the cosine transform of first kind. ∎

Corollary 4.3**.**

For fixed $\lambda$ the ordinary cross-validation score $P(\lambda)$ and the generalized cross-validation score $V(\lambda)$ on the unit intervall for Chebyshev nodes of first kind can be computed in $\mathcal{O}(N\log N)$ using Algorithm 1.

Proof.

Because multiplying with $\bm{F}$ and $\bm{F}^{\mathsf{H}}$ can be done using the fast discrete cosine transform, we see that applying the hat matrix can be done in $\mathcal{O}(N\log N)$ and Theorem 4.2 allows us to compute the diagonal entries of the hat matrix in $\mathcal{O}(N\log N)$. ∎

To emphasize our results numerically we chose the peaks sample function $f$ from Matlab and fixed the second argument to zero. We evaluated $f$ in $N=128$ Chebyshev nodes and added 5% Gaussian noise as one can see in Figure 4.1, (a). To choose $\bm{\hat{W}}$ we made use of the following statement from [48, Theorem 7.1] which relates the smoothness of $f$ to the decay of the Chebyshev coefficients $\bm{a}$: If for $\nu\geq 0$ the derivatives up to $f^{(\nu-1)}$ are absolute continuous and $f^{(\nu)}$ has bounded variation $V$ then $|a_{k}|\leq 2V/(\pi(k-\nu)^{\nu+1})$. Because in general we do not know anything about the smoothness of the function $f$ we chose $\hat{w}_{n}=n^{3}$ as weights which corresponds to a function with one absolute continuous derivative. We used Algorithm 1 to calculate the ordinary cross-validation score $P(\lambda)$ and the generalized cross-validation score $V(\lambda)$ for $\lambda\in[2^{-16},2^{-11}]$ and plotted the regularization for the $\lambda$ with the smallest corresponding ordinary cross-validation score as one can see in Figure 4.1, (b).

We observe that the ordinary cross-validation score and the generalized cross-validation score differ only slightly and their minima are close to the $L_{2}([-1,1],(1-x^{2})^{-1/2})$ optimal $\lambda$.

4.3 Approximative Quadrature

In this section, we consider arbitrary, ordered nodes $x_{m}\in\mathcal{X}\subset[-1,1]$, $m=0,\ldots,M$. The corresponding cosine transforms can be computed using the nonequispaced discrete cosine transform, cf. [11], in $\mathcal{O}(N\log N+|\mathcal{X}|)$ where $N$ is the bandwidth. As in Section 3.4 we determine approximate quadrature weights $w_{x_{m}}$ for $m=0,\ldots,M$ that allow us to efficiently estimate the diagonal entries of the hat matrix $\bm{H}$. Since we consider the unit interval with the non-uniform weight function $(1-x^{2})^{-\frac{1}{2}}$ it is not a good idea to compute Voronoi weights directly. Instead, we consider the corresponding periodic approximation problem on the unit circle with constant weight by substituting $y_{m}=\arccos x_{m}\in[0,\pi]$ and use Voronoi weights with respect to $y_{m}$, i.e.,

[TABLE]

Remark 4.4**.**

Let $x_{m}$ be the Chebyshev nodes of first kind. Then the quadrature weights (4.1) coincide with the exact quadrature weights from Section 4.1.

Analogously to Section 3.4 we use the approxmiated hat matrix $\bm{\tilde{H}}$ from (3.3) for ease of computation.

Remark 4.5**.**

Using Theorem 4.1 the diagonal entries can be calculated efficiently with one nonequispaced discrete cosine transform.

With the given tools we can modify Algorithm 1 to compute $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$ from Definition 3.8 in $\mathcal{O}(N\log N+|\mathcal{X}|)$ arithmetic operations given a fixed number of iterations to compute the Tikhonov minimizer.

To exemplify our results we chose $128$ uniformly distributed nodes on the unit interval which we perturbed by $5$% Gaussian noise. Note that uniformly distributed nodes are far from optimal in the setting of polynomial interpolation on the interval. As in case of exact quadrature we set the bandwith equal to the number of nodes, i.e., $N=|\mathcal{X}|$. As the Voronoi weights resamble quadrature weights the choice of the bandwidth $N$ is critical because in the case of $|\mathcal{X}|<N$ one can not expect to get an exact quadrature formula. As test function we used again the Matlab peaks function with fixed second argument. Then we computed $P(\lambda)$, $\tilde{P}(\lambda)$, $V(\lambda)$, and $\tilde{V}(\lambda)$ for $\lambda\in[2^{-18},2^{-11}]$. The results can be seen in Figure 4.2.

We observe that all cross-validation scores follow the shape of the $L_{2}([-1,1],(1-x^{2})^{-1/2})$-error and their minima are close to the optimal $\lambda$. Again, $\tilde{P}(\lambda)$ is affected by oscillations for small $\lambda$ which are caused by diagonal entries $h_{x_{m},x_{m}}$ close to $1$. The computation of the exact $P(\lambda)$ and $V(\lambda)$ averaged over all $\lambda$ takes $4.07$ seconds whereas the approximated $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$ outperform this with $0.04$ seconds.

5 Cross-validation on the two-dimensional sphere

Approximation on the two-dimensional sphere ${\mathbb{S}^{2}}\coloneqq\{\bm{x}\in\mathbb{R}^{3}:\|\bm{x}\|_{2}=1\}$ has been subject of mathematical research for a long time. The base for approximation from scattered data is formed by positive quadrature rules, Marcinkiewicz-Zygmund inequalities which are investigated in the papers [52, 36, 33, 5], and by bounds for best approximations [44, 51, 21]. Based on these result the relationship between the mesh norm, the separation distance of the sampling points, and optimal approximation rates has been analyzed in the papers [13, 31, 28, 47]. Approximation from noisy data has been considered in [1] and a priori and a posteriori estimates of the approximation error with respect to the regularization parameter have been proven in [38].

Following the approach of the previous sections we again consider the weighted Tikhonov functional (1.1). The analogue of the exponential functions become the spherical harmonics $\{Y_{n,k}\}_{n=0,\dots,\infty,k=-n,\dots,n}$, cf. [14, 2, 37, 7], which we assume to be normalized such that they form an orthonormal basis in $L_{2}({\mathbb{S}^{2}})$. For nodes $\mathcal{X}\subset{\mathbb{S}^{2}}$ and a maximum polynomial degree $N\in\mathbb{N}$ the Fourier matrix $\bm{F}$ becomes

[TABLE]

As for the weight matrix $\bm{\hat{W}}$ in Fourier space we consider isotropic weights $\bm{\hat{W}}=\operatorname*{diag}(\hat{w}_{n,k})_{h=0,\ldots,N,\,k=-n,\ldots,n}$ that depend only the polynomial degree, i.e., $\hat{w}_{n,k}=\hat{w}_{n}$.

5.1 Exact Quadrature

There are several approaches to exact quadrature on the two-dimensional sphere. The most direct approach is to consider tensor products of Gauss quadrature rules on the circle and the unit interval $[-1,1]$, cf. [39, Section 9.6]. A relaxation of this idea is to choose the points equally spaced at fixed latitudinal circles which also allows for an explicit computation of the quadrature weights, cf. [42].

A second approach is to choose the quadrature nodes approximately uniform and determine the quadrature weights by solving a linear system of equations. Given that the quadrature nodes are sufficiently well separated and the oversampling factor is sufficiently high, the resulting quadrature weights can guarantied to be nonnegative, cf. [36]. The computation of these quadrature weights can be implemented efficiently using fast spherical Fourier techniques, cf. [32, 30, 27].

A third approach, called Chebyshev quadrature, consists of fixing the weights to be constant and seeking quadrature nodes with a high degree of exactness. The resulting nodes are known as spherical t-designs. Efficient algorithms for their computation are described in [18] with the resulting spherical designs being available at [16]. Finally, one can try to compute both quadrature nodes and weights in an optimization schema, cf. [17].

For this section it is sufficient that the nodes $\mathcal{X}$ and the weights $\bm{W}=\operatorname*{diag}(w_{\bm{x}})_{\bm{x}\in\mathcal{X}}$ form an exact quadrature rule of degree $2N$, i.e., $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}=\mathbf{I}$. Under this assumption the diagonal entries of the hat matrix

[TABLE]

can by computed efficiently as it is stated in the following theorem.

Theorem 5.1**.**

Let the nodes $\mathcal{X}$ and the weights $\bm{W}$ form a quadrature formula $Q_{\mathcal{X},\bm{W}}$ that is exact for all spherical harmonics up to polynomial degree $2N$ then the diagonal entry corresponding to $\bm{x}$ of $\bm{H}$ satisfies

[TABLE]

Proof.

Since $\bm{F}^{\mathsf{H}}\bm{W}\bm{F}=\bm{I}$ we obtain analogously to the proof of Theorem 3.2

[TABLE]

Looking into the diagonal entry corresponding to $x$ and using the addition theorem of spherical harmonics, cf. [37, Theorem 5.11], we obtain the formula

[TABLE]

∎

Corollary 5.2**.**

For fixed $\lambda$ the ordinary cross-validation score $P(\lambda)$ and the generalized cross-validation score $V(\lambda)$ on the two-dimensional sphere given quadrature nodes and weights can be computed in $\mathcal{O}(N^{2}\log N+|\mathcal{X}|)$ using Algorithm 1.

Proof.

Due to Theorem 5.1 we can compute $h_{\bm{x},\bm{x}}$ in linear time. Using equation (5.1) applying the hat matrix has the same computational cost as one multiplication with $\bm{F}$ and one with $\bm{F}^{\mathsf{H}}$. Using the nonequispaced fast spherical Fourier transform (NFSFT, cf. [32]) this can be done in $\mathcal{O}(N^{2}\log N+|\mathcal{X}|)$. ∎

In order to illustrate Theorem 5.1 we consider a quadrature rule consisting of $21\,000$ approximately equidistributed nodes and equal weights $w_{\bm{x}}=4\pi/21\,000$ that is exact up to polynomial degree $2N=200$, as reported in [16]. Since by Theorem 5.1 the diagonal entries $h_{\bm{x},\bm{x}}$ of the hat matrix are constant multiples of the constant spatial weights $w_{\bm{x}}$ the ordinary cross-validation score and the generalized cross-validation score coincide for this setting. For weights in frequency domain we have chosen $\hat{w}_{n}=(2n)^{2s}$ for $s=3$ which corresponds to a function with $3$ derivatives in $L_{2}({\mathbb{S}^{2}})$.

The test function consists of a sum of quadratic B-splines to which we added an error of 5% Gaussian noise for each node as one can see in Figure 5.1, (a). This function was suggested in [23]. We calculated $V(\lambda)$ and $P(\lambda)$ for $\lambda\in[2^{-38},2^{-25}]$ using Algorithm 1 with the help of the Matlab toolbox MTEX, cf. [22]. Furthermore we calculated the $L_{2}({\mathbb{S}^{2}})$-error using Parseval from the original $\bm{\hat{f}}$ and $\bm{\tilde{f}}$ which are a byproduct of Algorithm 1.

As it is illustrated in the Figure 5.1, (b) the minimum of the cross-validation score is very close to the minimum of the approximation error.

5.2 Approximative quadrature

In the case function values are provided at nodes not forming a suitable quadrature rule we follow the previous ideas of Section 3.4 and 4.3 and use the approximated hat matrix $\bm{\tilde{H}}$ from (3.3) instead of $\bm{H}$ itself. This way we acquire $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$ as in Definition 3.8. In place of quadrature weights we use a spherical Voronoi decomposition, cf. [41].

The only changes to Algorithm 1 are the prior computation of the Voronoi weights and the necessity of solving a linear system of equations for computing the Tikhonov minimizer $\bm{\tilde{f}}$. This system of linear equation was solved with a conjugate gradient method, as implemented in the Matlab lsqr command, limited to 20 iterations.

In order to illustrate the efficiency of approximative quadrature weights for estimating the cross-validation score we consider the same test function and $\bm{\hat{W}}$ as in Example 5.1 and apply Algorithm 1 with polynomial degree $N=30$ to $|\mathcal{X}|=2(N+1)^{2}=1922$ random nodes, which corresponds to an oversampling factor of two. Figure 5.2 compares the different cross-validation scores $P(\lambda)$, $V(\lambda)$, $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$ for $\lambda\in[2^{-38},2^{-25}]$. The interval has been choosen by hand in such a way that it reflects the characteristic features of the functionals. All scores have their minimum close to the minimum of the of the actual approximation error. On the downside, we again observe several peaks in the approximated ordinary cross-validation score for small values of $\lambda$. So it is important to start the minimization process with a large $\lambda$. We also want to note that the computation of the exact $P(\lambda)$ and $V(\lambda)$ took $227$ seconds averaged over all $\lambda$ in contrast to $0.12$ seconds for the approximated $\tilde{P}(\lambda)$ and $\tilde{V}(\lambda)$.

6 Conclusion

In this paper we presented a fast algorithm for the computation of the leave-one-out cross-validation score $P(\lambda)$ for the Tikhonov regularizer (1.1). While many approximations of the cross-validation score do not consider every node separately, but fewer and larger validation sets, c.f. [6], we considered the full leave-one-out cros-validation score. In contrast to other approaches we did not restrict ourselves to spline interpolation on the interval at equispaced nodes, but considered more general domains and samplings. The key points of Algorithm 1 are explicit formulas for the diagonal elements $h_{x,x}$ of the hat matrix $\bm{H}$ which we were able to derive in the Theorems 3.2, 4.2, and 5.1, for approximation on the torus, the interval, and the two-dimensional sphere, respectively. Generalizations to other domains, e.g., the rotation group $SO(3)$, are possible following the framework presented in this paper. For all these domains FFT-like algorithms can be applied to achieve quasilinear complexity with respect to the number of nodes for the computation of the Tikhonov minimizer as well as for the leave-one-out cross-validation score.

The efficiency of our approach has been illustrated in several numerical experiments with respect to the different domains. For the nodes we distinguished two settings. For nodes belonging to a quadrature rule, like equispaced nodes or rank-1 lattices on the torus, our Algorithm 1 computes the cross-validation score $P(\lambda)$ with floating point precision, cf. Corollaries 3.3, 4.3, and 5.2. For arbitrary nodes we accomplished in Remarks 3.9, 4.5 and Corollary 5.2 a good approximation using Voronoi weights in place of the quadrature weights. The numerical experiments confirm our theoretical results. In all test scenarios our algorithm was several orders of magnitude faster then the direct reference implementation.

In some cases the approximated leave-one-out cross-validation score $\tilde{P}(\lambda)$ suffered from peaks for $\lambda$ smaller than the optimal one, cf. Subsection 3.4. Anyway, in our test cases we had no problems finding the global minimum by initializing the line search algorithm with a sufficiently large $\lambda$ and thus avoiding the oscillatory region.

All relevant Matlab code, including the algorithm for the fast computation of the leave-out-one cross-validation score, its minimizer and all numerical examples of this paper can be found on the GitHub repository https://github.com/felixbartel/fcv.