Minimax Robust Function Reconstruction in Reproducing Kernel Hilbert Spaces
1 Abstract — In this paper, the problem of approximating a function belonging to an arbitrary real-valued reproducing kernel Hilbert space (RKHS) on d based on approximate observations of the function is considered. The observations are approximate in the sense that the actual observations are known only to belong to a convex set of admissible observations. A minimax optimal approximation for the function is sought that minimizes the supremum of the RKHS norm on the error between the true function and the chosen approximation subject only to the conditions that the true function belongs to a uniformly bounded uncertainty set of functions that satisfy the constraints on the observations and that the approximation is a member of the RKHS. Such a solution is referred to as a minimax robust reconstruction . The solution to the minimax robust reconstruction problem is characterized and it is shown to be equivalent to solving a straightforward convex optimization problem. The stability properties of the minimax robust reconstruction are investigated and the approach is motivated by characterizing the minimax robust reconstruction for several specific convex observational models and discussing relationships with other approaches to function approximation. Index Terms — RKHS, Smoothing Splines, Scattered Data Interpolation, Optimal Approximation, Chebyshev Center I. I NTRODUCTION he results presented in this paper are related to the general problems of scattered data interpolation and approximation, and in particular to the large body of results on scattered data interpolation and approximation in a reproducing kernel Hilbert space (RKHS) setting [1-11]. We consider the problem of reconstructing an unknown function : d f based on approximate observations of the function on a subset of points in d . By approximate observations, we mean that the actual observations (that is, the true values of the function f on the set ) are known only to belong to a convex set of admissible observations. In particular, we let f represent the observations (equivalently, the restriction of f to the set ), and we assume that the observations are known only to satisfy f , where is an arbitrary convex subset of . We further assume that f , where is an RKHS with norm , inner product Here and elsewhere, the notation represents the set of functions : g . , , and reproducing kernel : d d K , and that arbitrary functions in are not determined uniquely by their values on the set of points d . We seek a minimax robust reconstruction ˆ f of f that satisfies ˆsup inf sup gf f f f g f , (1) where : , f f M f , (2) represents a non-empty uncertainty set of admissible reconstructions of f , and M is any real number satisfying inf : f M f f . Note that putting some sort of bound on the norm for the set of admissible reconstructions is generally necessary in order to generate a well posed problem. That is, assuming that the set , f f is not empty and that a function f is not uniquely determined by its values on the set , it is straightforward to show that for any function g , we have sup : f g f f . Hence, without the bound on the norm, all approximations can be regarded as equally good (or equally bad) in terms of maximum possible deviation from the true function. On the other hand, we show below that the solution to Problem (1) does not actually depend on the value of M so that ˆ f solving (1) satisfies ˆ ˆlim MM f f , where ˆsup : ,inf sup : , . Mf g f f f f M fg f f M f
Since That is, we are only interested in the case when there are some degrees of freedom in the function reconstruction that are not determined by the set of observations.
Minimax Robust Function Reconstruction in Reproducing Kernel Hilbert Spaces
Richard J. Barton T :lim : , , M f f f f M f it follows that ˆ f solving (1) can be regarded as the unique meaningful solution to the problem ˆsup :inf sup : . f g f f f fg f f Robust function reconstruction in contexts similar to the one considered here has been studied by many different investigators, and an excellent survey of those results is presented in [12, 13]. The particular approach taken here is of interest because of the very general form of the set and because of the simple, intuitive form of the corresponding solution. Note that in the simplest case, the convex set of observations in (2) takes the form y , where N y is the nominal value of the function observed at N points in d . In this case, consists of only a single point, and the uncertainty set is just the set of all possible functions in that interpolate between the observed values and satisfy an (essentially arbitrary) upper bound on the norm of f . As we discuss in the next section, the solution of (1) in this case is equivalent to the minimum-norm interpolation for f . We show that with proper selection of the convex set , the solution of (1) can be identified with other well-known function approximation techniques such as smoothing splines and Moore-Penrose pseudo-inverses. Hence, part of the significance of this paper is that it presents a unified approach to function approximation that establishes for the first time the minimax optimality of function approximation techniques such as minimum-norm interpolation, function smoothing, and pseudo-inverses. These approaches have previously been justified and studied primarily on heuristic grounds. Furthermore, this approach provides considerable insight into stability and error analysis by stating and solving the approximation problem in a very geometrical context. Finally, from a practical point of view, the method easily accommodates a broad range of different observational models, accounts for the inevitable presence of noise in observations while still producing solutions that fall within a prescribed range of deviation from nominal behavior, and is computationally tractable. The remainder of this paper is organized as follows. The basic results on minimax robust reconstruction are stated and proven in Section II. A discussion of computational complexity, stability characteristics and connections with other function approximation techniques is given in Section III. The method is illustrated by solving an approximation problem involving two-dimensional thin-plate splines in Section IV, and some concluding remarks are given in Section V. II. R ESULTS
Recall that in an RKHS over d with kernel K , we have , K v for any d v . Furthermore, the reproducing property of the kernel K implies that , , f v f K v , for any f and any d v . Recall further that the set sp , , K u u , (3) consists of limits of functions of the form , n i ii s v a K v u , for arbitrary , , , n a a a and , , , n u u u , and that , , , , , f K u f K u u , where f represents the projection of f onto . It follows from these relationships that for an uncertainty set of the form (2), we have f if and only if f M and f . Since we assume that a function f is not uniquely determined by its values on the set , it also follows that ; that is is a proper subset of and so is . Our main result is as follows. Theorem 1 . Let be an RKHS over d with reproducing kernel K , let d be a (possibly infinite) set of observation points for an unknown function f , and assume that functions in are not uniquely determined by their values on the set . Let be given by (2), where is convex, and let be given by (3). Then a solution to Problem (1) always exists and is given by the unique element ˆ f satisfying ˆ min min f f f f f . (4) That is, the minimax robust reconstruction is the minimum norm element of . It should be noted that the solution to (1) does not depend on the choice of the upper bound inf : f M f f . To prove Theorem 1, we need the following two lemmas. Lemma 1 . ˆ f solves Problem (1) if and only if ˆ f and ˆ f satisfies Here and elsewhere, the over-bar notation indicates the closure of the set . ˆsup inf sup gf f f f g f ; that is, to find the optimal minimax reconstruction of the unknown function f , it is sufficient to consider functions ˆ f . Proof . Recall [14] that for any closed convex subset and any g , there is a unique element g , called the projection of g onto , such that
22 0 , g f g g f . Furthermore, g satisfies , 0, f g g g f . Letting and g , it follows that
22 0 020 0 020 2 20 020 2 sup sup 2 ,supsupsupinf sup . f ffffg f g f g g g fg g f g g gg fg g g fg fg f
Hence, inf sup inf sup g gf f g f g f , and since inf sup inf sup gg f f g f g f , we have inf sup inf sup inf sup gg gf f f g f g f g f . It follows that inf sup inf sup g gf f g f g f . Now, let g and f . Then g and g f g g g f f fg g f fg g g fg g f f If , 0 g g f f , then g f g f . If , 0 g g f f , let f f f f . Then f f , so f and g f g g g f f fg g f fg g g fg g f fg f It follows that,
22 2 max , g f g f g f . Hence, for all f , we have
22 2 2 sup max , f g f g f g f g f , which implies that sup sup inf sup gf f f g f g f g f , which in turn implies that inf sup inf supinf sup . g gf fg f g f g fg f Hence, inf sup inf supinf sup , g gf fg f g f g fg f which proves the lemma. ■ Lemma 2 . ˆ f solves Problem (1) if and only if ˆ ˆsup 2 ,inf sup 2 , . f g f M f f fM f g g (5)
Proof . Let g . Since, , the orthogonal complement of in contains non-zero elements, and for any f , h , we have M ff h fh , with f M and f f . Hence, 4 sup sup 2 ,sup 2 ,sup 2 ,sup 2 , . f ffff g f f f g gf f g gM f g gM f g g It follows that, inf supinf sup 2 , . g fg f g fM f g g ■ Proof of Theorem 1 . The results of Lemma 1 and Lemma 2 imply that Theorem 1 will be proven if we can show that a solution to Problem (5) always exists and is given by the unique element ˆ f satisfying Equation (4). Toward that end, we define the loss functional : as , 2 , f g M f g g . To characterize the solutions to (5), we seek a saddle point for the game , , ; that is, we seek , L R f g such that , ,, , , .
R L RL f g f gf g f g (6) Since is continuous, a solution to (6) satisfies , sup , inf sup , L R R gf f f g f g f g . Similarly, for any g , , , sup , L R L f f g f g f g . Hence, , inf sup ,
L R g f f g f g , and it follows that , inf sup ,inf sup 2 , . L R g fg f f g f gM f g g
Substituting ˆ R f g gives the desired solution to (5). To find necessary and sufficient conditions for , L R f g to solve (6), we note that the Cauchy-Schwarz inequality implies , , L L L L f g M f f f for all g , with equality if and only if L g f . Hence, , , , L R L f g f g g if and only if
R L g f . On the other hand, we have , , , R L R f g f g f if and only if L R L RR L RR L R f f g f gf g f gf g f g where the notation “ ” indicates that two statements are equivalent. Therefore, , L R f g solves (6) if and only if
R L g f and , , 0,, , 0,, 0,min , R R RR R RR RR g f g g g fg g g g gg g g gg g where the last statement follows from the definition of the projection of the origin onto ( i.e. , the minimum norm element of ). This establishes the theorem. ■ III. D ISCUSSION
It follows from Theorem 1 that solving Problem (1) is equivalent to finding the minimum-norm element of a convex set. From a practical point of view, one must also invert a possibly very large and ill-conditioned matrix, but the minimum-norm problem is still the most computationally complex part of the solution. Fortunately, finding the minimum-norm element of a convex set is a straightforward convex programming problem, and well known polynomial-time algorithms can be applied to find the solution [15-18]. To illustrate what is involved in solving a minimax robust reconstruction problem and to study the properties of the solutions, we characterize in this section the solutions for three particular convex uncertainty classes. Throughout this section, we assume that the vector , , , TN y y y y represents a nominal vector of function observations at N points , , , dN v v v , and we let the gram matrix for the problem be given by , , ,, , ,, , , NNN N N N
K v v K v v K v vK v v K v v K v vK v v K v v K v v K . For simplicity, we assume that K is invertible. We focus on the three uncertainty sets 5 , , 1, 2,3 i i i f f M i x , where , , , , , , TTN N x x x f v f v f v x , and
11 101 2 22 203 0,1, , 1 , ,, ,max , .
NN n nnNN n nnN n nn N x yx y x y x yx yx y Here and elsewhere, the notation , p y indicates a ball of radius around the nominal observation vector y with respect to the p -norm in N for
1, 2, p . To avoid the trivial solution ˆ 0 f , we assume that in each case, has been chosen such that i does not contain the origin. Then the subspace consists of functions of the form , , N dn nn f v h K v v v , where , , , TN h h h h , and i f if and only if f x Kh for some i x , in which case we have T T T f h Kh x K x h x . Furthermore, if i f with i f x Kh and i g with i g x Kh , then
11 2 1 2 1 2 , T T T f g h Kh x K x h x . It follows from Theorem 1 that the minimax robust reconstruction ˆ i f takes the form ˆ ˆ , , N dn nn f v h K v v v , with ˆˆ x Kh , where ˆ i x satisfies ˆˆ ˆ , T i f g g h x . That is, ˆ f is the minimum-norm element of i , or equivalently, ˆ f is the projection of the origin onto the closed convex set i . Hence, for all i g of the form , , N dn nn g v h K v v v , with x Kh , we must have ˆ ˆ, f f g , or equivalently, ˆ ˆˆ , T T i h x h x x . (7) Hence, solving Problem (1) for each of the uncertainty sets i ,
1, 2,3 i , is equivalent to identifying a pair ˆ ˆ, h x that satisfies inequality (7) with ˆˆ x Kh . These three problems have been solved by Verdu and Poor [19] in the context of robust matched filtering . The solutions are discussed below. : The pair ˆ ˆ, h x with ˆˆ x Kh satisfies Equation (7) if and only if ˆˆ sgn , 0,1, , 1 n n n n x y h n N , where ˆ ˆ0 if max ,. n n nn NN nn h h In the special case diag , , , N K , the robust matched filter ˆ h is a clipped version of the nominal matched filter h K y ; that is, ˆ , if ,sgn , if ,, n n nn nN n nn h h hh hh where max 0, x x . : The pair ˆ ˆ, h x with ˆˆ x Kh satisfies Equation (7) if and only if ˆˆ ,ˆ . x y hh where the symbol “ ” indicates the Euclidean norm of a vector. Alternatively, ˆ , . h K I yK I y : The pair ˆ ˆ, h x with ˆˆ x Kh satisfies Equation (7) if and only if ˆ ˆ, if 0,ˆ, if 0. n n nn n x y hy h If h K y and minmax iiN jnj n h K , then 6 ˆ , if 0,, if 0. n n nn n x y hy h A. Stability and Connections with Other Approaches
To get some insight into the stability of minimax robust reconstruction, we consider the solution for the uncertainty set discussed above. In this case, the convex set of admissible observations is given by , NN n nn x y x y . Interestingly, the solution for Problem (1) in this case is identical to the solution of a classical RKHS smoothing problem [11, 20-25] . That is, the function ˆ f given by ˆ ˆ , , N dn nn f v h K v v v , (8) with ˆ h K I y , (9) and
11 12 2
K I y K I y , (10) also satisfies
12 21 1 22 1 ˆ ˆinf .
N n nn N n nf n f f v yf f v y
Hence, ˆ f can be regarded as a generalized smoothing spline , where the smoothing parameter and the uncertainty parameter are chosen jointly to satisfy (10). It follows that minimax robust reconstruction can be regarded as a generalization of function smoothing. Note that in the classical smoothing problem, one generally starts with a fixed value of , so that the uncertainty parameter becomes a function of the selected value of . In fact, Equation (10) implies that as , y as , and that, in general, the region of implied uncertainty increases with the norm of the nominal observed vector y . Note also that the equivalence established here between minimax robust reconstruction and function smoothing reveals the previously unrecognized minimax optimality of smoothing, which is generally justified and developed in a much more heuristic context. As a result of the equivalence with smoothing in this particular case, one might conjecture that minimax robust reconstruction in general has the same stability properties as smoothing. That is, the general solution to the RKHS smoothing problem, given by (8) and (9) where is a fixed parameter for all nominal observation vectors y , is This is a generalization of a similar well-known relationship for univariate function approximation with polynomial splines [24, 25]. unconditionally stable relative to variations in y since K I will not be ill-conditioned even if K is. Unfortunately, the stability properties of minimax robust reconstruction in general are not so well defined, and although the approach is much less sensitive to ill-conditioned gram matrices than techniques such as minimum-norm interpolation, it is not unconditionally stable even in the case illustrated here. To see what is going on in some generality, let us assume that the set N that determines the uncertainty set is a closed convex set of small diameter with a fixed but arbitrary shape relative to its centroid y , which represents a nominal observation vector. In this case, a new nominal observation vector y changes the location but not the shape or the orientation of the set . While this is certainly not completely general, it covers a large class of useful uncertainty sets, and it is straightforward to visualize the behavior as y changes. The stability of minimax robust reconstruction in this case is determined by the rate of change in the solution ˆ f with respect to the norm relative to arbitrarily small changes in the coordinates of y . Furthermore, the quantity ˆ f is given by ˆ ˆ ˆ T f x K x , where ˆ x is the minimum-norm element of the set with respect to the norm T K x x K x . To illustrate, we adopt the unit eigenvectors of K as our coordinate system, where N are the eigenvalues with associated orthonormal (with respect to ) eigenvectors , , , N v v v . The set can be visualized (roughly) with respect to as an N -dimensional sphere centered at y . On the other hand, with respect to K , which is really what matters, can be visualized (again roughly) as an N -dimensional ellipsoid with major axis oriented along N v and minor axis oriented along v . Clearly, for visualization purposes, we have reverted to the earlier example with . This is illustrated in Figure 1 for the two-dimensional subspace spanned by v and N v . 7 (a) (b) Fig. 1. Illustration of minimax robust reconstruction geometry – (a) with respect to ; (b) with respect to K . It is now easy to see that the stability of the solution relative to y depends greatly on the location of y and the shape of . For example, in the situation illustrated in Figure 1(b), a small change in the coordinate of the nominal observation vector in the positive N v direction will produce a relatively large change in ˆ x , particularly if K is ill-conditioned. This is illustrated in Figure 2(a), which indicates a large change in ˆ x from the value illustrated in Figure 1(b). On the other hand, a small change of the coordinate in any other direction, for example, even in the negative N v direction, will have much less impact on ˆ x . This is illustrated in Figure 2(b), which now indicates a much smaller change in ˆ x from the value illustrated in Figure 1(b). Hence, even if K is ill-conditioned, the minimax robust reconstruction will be much less sensitive to small changes in observation on the average than interpolation. (a) (b) Fig. 2. Illustration of stability properties of minimax robust reconstruction – (a) small change in the positive N v direction; (b) small change in the negative N v direction. In fact, minimax robust reconstruction in general can be made unconditionally stable if the set is chosen appropriately. For example, suppose s u ,where s is a convex subset of the subspace spanned by the “stable” eigenvectors of K, and u is a convex subset of the subspace spanned by the “unstable” eigenvectors of K . If the origin is contained in the interior of u but is not contained in s , then Problem (1) will have a non-trivial ( ˆ 0 x ) but unconditionally stable solution. A simple example of this is illustrated in Figure 3. For this figure, we have again illustrated the two-dimensional subspace spanned by v and N v , but in this case we let s u , where s and u are arbitrary convex sets, which in one dimension just correspond to real-valued intervals. Hence, for this example, the uncertainty region becomes a simple rectangle that is much more elongated with respect to the K norm than with respect to the norm. (a) Here, the terms stable and unstable refer to the relative magnitudes of the eigenvalues. That is, the stable eigenvalues are those with relative magnitude larger than some threshold and the unstable eigenvalues are those with relative magnitude smaller than the same threshold. This is analogous to the manner in which the singular values would be chosen in the singular value decomposition of K . v v N y v v N y ˆ x g v v N y ˆ x g v v N y ˆ x g v v N y (b) Fig. 3. Illustration of unconditionally stable minimax robust reconstruction geometry – (a) with respect to ; (b) with respect to K (solid line indicates original location of , dashed line indicates position after small change in coordinate in the N v direction). As Figure 3 indicates, since the origin is contained in the interval representing u , a small change in the coordinates of the nominal observed vector y in the N v direction will produce no change in the coordinate of the minimax reconstruction ˆ x in that direction, which is zero in both cases. Interestingly, constructing in this fashion is analogous to using the Moore-Penrose pseudo-inverse [26] of K rather than the true inverse of K when finding a solution in the case of minimum-norm interpolation. That is, in both cases, the projection of the solution onto the unstable eigenvectors is zero. In fact, if the set contains only the point representing the projection of the nominal observed vector y onto the span of the stable eigenvectors, then this approach is identical to using the pseudo-inverse of K to solve the minimum-norm interpolation problem. It follows that minimax robust reconstruction can also be regarded as a generalization of regularization applied to minimum-norm interpolation. Here again, the equivalence established between minimax robust reconstruction and matrix regularization reveals the previously unrecognized minimax optimality property of a function approximation technique that is generally developed primarily in a heuristic context. IV. C ONCLUSION
In this paper, we have introduced and studied a multidimensional RKHS function approximation technique that can be viewed as a generalization of scattered data interpolation. The desired approximation is derived by solving a minimax problem with respect to a uniformly bounded uncertainty set of admissible functions that are required to satisfy a known set of convex constraints on the observations. We refer to this approach as minimax robust reconstruction. We have demonstrated that this approach to function approximation has several significant advantages. From a practical point of view, a major advantage is that either empirical or a-priori information about the statistical distribution of noise in the observations can be easily and explicitly incorporated in the selection of the uncertainty set . In addition, this approach places the function approximation problem in a very geometrical context that facilitates the choice of the uncertainty set and provides insight into sensitivity and error analysis. From a theoretical perspective, minimax robust reconstruction can be regarded as a unified approach to function approximation that includes many other popular techniques as special cases and establishes the previously unrecognized minimax optimality associated with many of these techniques. The work presented here can be extended easily to handle both approximation in a general Hilbert space and vector-valued functions in a very similar fashion. In addition, the approach can be extended to identify minimax robust reconstructions in situations where the appropriate norm on the function space is itself not completely identified. This might be useful, for example, in situations where the observed functions are regarded as realizations of a Gaussian random process for which the covariance function is known only to belong to a convex set of possible covariance functions. Extensions such as these, as well as other questions of interest regarding the proposed approach, will be studied in future work. R EFERENCES [1] A. Bouhamidi and A. Le Mehaute, "Multivariate Interpolating (m,l,s)-Splines,"
Advances in Computational Mathematics, vol. 11, pp. 287-314, 1999. [2] J. Duchon, "Splines Minimizing Rotation-Invariant Semi-Norms in Sobolev Spaces," in
Lecture Notes in Mathematics . vol. 571 Berlin: Springer, 1977, pp. 85-100. [3] W. R. Madych and S. A. Nelson, "Multivariate Interpolation and Conditionally Positive Definite Functions,"
Approximation Theory & its Applications, vol. 4, pp. 77-89, Dec. 1988. [4] W. R. Madych and S. A. Nelson, "Multivariate Interpolation and Conditionally Positive Definite Functions, II,"
Mathematics of Computation, vol. 54, pp. 211-230, Jan. 1990. [5] L. E. Mansfield, "On the Optimal Approximation of Linear Functionals in Spaces of Bivariate Functions,"
SIAM Journal on Numerical Analysis, vol. 8, pp. 115-126, March 1971. [6] J. Meinguet, "Multivariate Interpolation at Arbitrary Points Made Simple,"
Journal of Applied Mathematics and Physics, vol. 30, pp. 292-304, 1979. [7] F. J. Narcowich, "Recent Developments in Error Estimates for Scattered-Data Interpolation Via Radial Basis Functions,"
Numerical Algorithms, vol. 39, pp. 307-315, 2005. [8] F. J. Narcowich and J. D. Ward, "Scattered-Data Interpolation on R n : Error Estimates for Radial Basis and Band-Limited Functions," SIAM Journal on Mathematical Analysis, vol. 36, pp. 284-300, 2004. [9] F. J. Narcowich and J. D. Ward, "Norms of Inverses and Condition Numbers for Matrices Associated with Scattered Data,"
Journal of Approximation Theory, vol. 64, pp. 69-94, 1991. [10] L. L. Schumaker, "Fitting Surfaces to Scattered Data," in
International Symposium on Approximation Theory , University of Texas, Austin, TX, 1976, pp. 203-268. [11] H. L. Weinert, "Statistical Methods in Optimal Curve Fitting,"
Communications in Statistics: Simulation and Computation, vol. B7, pp. 417-435, 1978. [12] C. A. Micchelli and T. J. Rivlin,
Optimal Estimation in Approximation Theory . New York: Plenum Press, 1977. [13] C. A. Micchelli and T. J. Rivlin, "A Survey of Optimal Recovery," in
Optimal Estimation in Approximation Theory , C. A. Micchelli and T. J. Rivlin, Eds. New York: Plenum Press, 1977, pp. 1-54. [14] D. G. Luenberger,
Optimization by Vector Space Methods . New York: John Wiley & Sons, Inc., 1969. [15] S. Boyd and L. Vandenberghe,
Convex Optimization . Cambridge: Cambridge University Press, 2004. v v N y g ˆ x y g [16] H. H. Bauschke, P. L. Combettes, and D. R. Luke, "Finding Best Approximations Pairs Relative to Two Closed Convex Sets in Hilbert Spaces," Journal of Approximation Theory, vol. 127, pp. 178-192, 2004. [17] S. G. Nash and A. Sofer, "On the Complexity of a Practical Interior-Point Method,"
SIAM Journal on Optimization, vol. 8, pp. 833-849, August 1998. [18] R. Schaback and J. Werner, "Linearly Constrained Reconstruction of Functions by Kernels with Applications to Machine Learning,"
Advances in Computational Mathematics, vol. 25, pp. 237-258, 2006. [19] S. Verdu and H. V. Poor, "Minimax Robust Discrete-Time Matched Filters,"
IEEE Transactions on Communications, vol. COM-31, pp. 208-215, 1983. [20] P. Craven and G. Wahba, "Smoothing Noisy Data with Spline Functions: Estimating the Correct Degree of Smoothing by the Method of Generalized Cross-Validation,"
Numerische Mathematik, vol. 31, pp. 377-403, 1979. [21] G. M. Nielson, "Multivariate Smoothing and Interpolating Splines,"
SIAM Journal on Numerical Analysis, vol. 11, pp. 435-446, April 1974. [22] G. Wahba, "Smoothing Noisy Data with Spline Functions,"
Numerische Mathematik, vol. 24, pp. 383-393, 1975. [23] H. Wendland and C. Rieger, "Approximate Interpolation with Applications to Selecting Smoothing Parameters,"
Numerische Mathematik, vol. 101, pp. 729-748, 2005. [24] C. H. Reinsch, "Smoothing by Spline Functions,"
Numerische Mathematik, vol. 10, pp. 177-183, 1967. [25] C. H. Reinsch, "Smoothing by Spline Functions, II,"
Numerische Mathematik, vol. 16, pp. 451-454, 1971. [26] R. A. Horn and C. R. Johnson,
Matrix Analysis . Cambridge: Cambridge University Press, 1988.