In this paper it is shown that many features from polynomial spline methods used in nonparametric regression and smoothing procedures carry over to the class of L-splines where L is a linear differential operator of order 4 with constant coefficients. Special attention is given to the question whether an analogue of the Reinsch algorithm is valid and criteria are given such that the associated matrix R is strictly diagonal dominant.Project KP-06-N32-8; Project KP-06N42-2 with Bulgarian NSF; Grant No BG05M2OP001-1.001-0003,financed by the Science and Education for Smart Growth Operational Program (2014-2020) and co-financed by the European Union through the European structural and Investment funds

Kounchev, Ognyan

Render, Hermann

Tsachev, Tsvetomir

Proceeding of the Bulgarian Academy of Sciences

Bulgarian Digital Mathematics Library at IMI-BAS

Доклади на Българската академия на наукитеComptes rendus de l’Académie bulgare des SciencesTome 75, No 1, 2022MATHEMATICSNumerical mathematicsL-SPLINE INTERPOLATION FOR DIFFERENTIALOPERATORS OF ORDER 4WITH CONSTANT COEFFICIENTSOgnyan Kounchev, Hermann Render∗, Tsvetomir TsachevReceived on April 20, 2021Presented by P. Kenderov, Member of BAS, on October 28, 2021AbstractIn this paper it is shown that many features from polynomial spline meth-ods used in nonparametric regression and smoothing procedures carry over tothe class of L-splines where L is a linear differential operator of order 4 withconstant coefficients. Special attention is given to the question whether ananalogue of the Reinsch algorithm is valid and criteria are given such that theassociated matrix R is strictly diagonal dominant.Key words: L-splines, interpolation, differential operators of order 42020 Mathematics Subject Classification: 65D05, 65D071. Introduction. Spline interpolation and smoothing provide an importanttechnique in nonparametric methods for data analysis, [1–4]. As outlined in [5] theextension of cubic spline smoothing to the setting of L-splines allows to includeadditional prior knowledge for estimating the data and to retain several aspectsof the data. For example, for the estimation of the gross domestic product datait is assumed that the model curve is a linear combination of the functions(1) 1, exp (γt) , sinωt, cosωt,The work of H. Render and Ts. Tsachev was funded under project KP-06-N32-8, while thework of O. Kounchev was funded under project KP-06N42-2 with Bulgarian NSF. The work ofthe first-named author has been partially supported by Grant No BG05M2OP001-1.001-0003,financed by the Science and Education for Smart Growth Operational Program (2014-2020) andco-financed by the European Union through the European structural and Investment funds.DOI:10.7546/CRABS.2022.01.0211while for melanoma data the model curve is linear combination of(2) 1, t, sinωt, cosωt.In both examples the model curve is a solution of a linear differential operator Lwith constant coefficients of order N + 1, given by(3) L = L(λ0,...,λN ) =N∏j=0(ddx− λj),where λ0, . . . , λN are complex numbers. Let us recall that a function g : [a, b]→ Cis an L-spline of orderN+1 with knots t1 < · · · < tn if g isN−1 times continuouslydifferentiable such that the restriction of g to each interval [tj , tj+1] is a solutionof the equation LN (g) = 0.Let us note that special classes of L-splines are an important component ofthe theory of polysplines and its applications, cf. [6–11].In many textbooks, like [12,13] an efficient algorithm for L-spline smoothingis described in detail and usually it is based on reproduction kernel techniquesassociated to the differential operator. In this paper we want to describe analgorithm for linear differential operators L of order 4 which mimics the orginalalgorithm of Reinsch from 1967 as closely as possible, [14]. An important featureof this algorithm is the fact that for given real numbers t1 < · · · < tn withn > 2 there exists a symmetric positive definite (n− 2)× (n− 2) matrix R and an× (n− 2) matrix Q and such that(4) QTg = Rγfor any natural cubic spline g with knots t1 < · · · < tn, wheregT = (g (t1) , g (t2) , . . . , g (tn))(5)γT =(g′′ (t2) , g′′ (t3) , . . . , g′′ (tn−1)).(6)The matrix R can be described explicitly: it has tridiagonal form and for thediagonal and off-diagonal we haveRjj =13(hj−1 + hj) and Rj,j+1 = Rj+1,j =16hj ,where hj = tj+1 − tj . Moreover a simple algorithm allows to compute the inter-polating natural spline in O (n) arithmetic operations.Given a general linear differential operator L(λ0,...,λ3) with constant coeffi-cients of order 4 we define a (n− 2)× (n− 2) matrix R and a n× (n− 2) matrixQ such that the relationQTg = Rγ12 O. Kounchev, H. Render, Ts. Tsachevholds for the vectors g and γ, defined bygT = (g (t1) , g (t2) , . . . , g (tn))(7)γT =(L(λ0,λ1)g (t2) , L(λ0,λ1)g (t3) , . . . , L(λ0,λ1)g (tn−1)).(8)The matrix R is tridiagonal and a simple algorithm allows to compute the inter-polating natural spline in O(n) arithmetic operations. To facilitate the expositionwe follow the notations given in [13].The matrices R and Q depend on the differential operator given by the val-ues λ0, λ1, λ2, λ3. In order to describe the matrices we need the notion of thefundamental function ΦΛN with respect to the vector ΛN = (λ0, . . . , λN ): thereexists a unique solution ΦΛn of the equation LNu = 0 such that ΦΛN (0) = · · · =Φ(N−1)ΛN(0) = 0 and Φ(N)ΛN (0) = 1, see [6,15]. In Section 2 we shall recall severalproperties of the fundamental function ΦΛN and describe its relationship with theGreen function. In Section 3 we will show that the matrix R is tridiagonal andthat the diagonal and off-diagonal entries are given byRj,j = ρ (tj − tj−1)− ρ (− (tj+1 − tj))(9)Rj,j+1 = σ (tj+1 − tj) and Rj+1,j = −σ (− (tj+1 − tj)) ,(10)where ρ and σ are functions defined byρ (x) =Φ′(λ0,...,λ3) (x) Φ(λ0,λ1) (x)− Φ(λ0,...,λ3) (x) Φ′(λ0,λ1)(x)Φ(λ0,λ1) (x) Φ(λ2,λ3) (x),(11)σ (x) =Φ(λ0,...,λ3) (x)Φ(λ0,λ1) (x) Φ(λ2,λ3) (x).(12)In a similar way the matrix Q is defined. In Section 4 we will show that theentries in (9) and (10) are positive provided that λ0, . . . , λ3 are real. We providein Section 5 a simple sufficient criterion such that the matrix R is strictly diagonaldominant. In the case that (λ0, λ1) and (λ2, λ3) are conjugation invariant we showthat for any ε > 0 there exists δ > 0 such that the estimate|Rj,j−1|+ |Rj,j+1| ≤(12+ ε)|Rjj |holds for all t1 < · · · < tn such that tj+1 − tj ≤ δ for j = 1, . . . , n − 1. For theexample in (2) we shall provide exact estimates. We show that for the special setΛ3 = (a,−a,−b, b) ∈ R4 with 0 ≤ a ≤ b the estimateRj,j−1 +Rj,j+1 ≤12Rjjholds. This estimate is important for the study and applications of polysplines(see [7–11,16]) which requires a deep analysis of L-splines.The present paper generalizes the results obtained in [17] for special classesof L-splines, see also [18].The detailed proofs will appear elsewhere.C. R. Acad. Bulg. Sci., 75, No 1, 2022 132. The fundamental function. We assume that we are given a vectorΛ ∈ CN+1, with the first coordinate having index 0, namelyΛN = (λ0, . . . , λN ) .For the differential operator LN defined in (3) we introduce the solution set:E (λ0, . . . , λN ) :={f ∈ CN (R,C) : LNf = 0},where CN (R,C) is the space of all N -times continuously differentiable complex-valued functions f : R → C. The elements in E(λ0,...,λN ) are called exponentialpolynomials or L-polynomials, and λ0, . . . , λN are called exponents or frequencies(see e.g. Chapter 3 in [19]).We say that the space E(λ0,...,λN ) is closed under complex conjugation, if forf ∈ E(λ0,...,λN ) the complex conjugate function f is again in E(λ0,...,λN ). It is easyto see that for complex numbers λ0, . . . , λN the space E(λ0,...,λN ) is closed undercomplex conjugation if and only if there exists a permutation σ of the indices{0, . . . , N} such that λj = λσ(j) for j = 0, . . . , N. In other words, E(λ0,...,λN ) isclosed under complex conjugation if and only if the vector ΛN = (λ0, . . . , λN ) isequal up to reordering to the conjugate vector ΛN . In this case we say that ΛNis conjugation invariant.It is well known that for ΛN = (λ0, . . . , λN ) ∈ CN+1 there exists a uniquesolution ΦΛN ∈ E(λ0,...,λN ) to the Cauchy problemΦΛN (0) = . . . = Φ(N−1)ΛN(0) = 0 and Φ(N)ΛN (0) = 1.We shall call ΦΛN the fundamental function in E(λ0,...,λN ) (see e.g. [6,15]). Anexplicit formula for ΦΛN is(13) ΦΛN (x) =12πi∫Γrexz(z − λ0) · · · (z − λN )dz,where Γr is the path in the complex plane defined by Γr (t) = reit, t ∈ [0, 2π],surrounding all the complex numbers λ0, . . . , λN . In particular, for N = 0 wehaveΦ(λ0) (x) = eλ0x.Note that (13) implies the useful formula(14)(ddx− λN+1)Φ(λ0,...,λN+1) (x) = Φ(λ0,...,λN ) (x) .If ΛN = ΛN , then we see that the conjugate of ΦΛN is again in E(λ0,...,λN ),and by uniquessness of the fundamental function we conclude:14 O. Kounchev, H. Render, Ts. TsachevProposition 1. If ΛN = ΛN , then ΦΛN (t) is a real-valued function.Proposition 2. If λ0, . . . , λN are real numbers, then Φ(λ0,...,λN ) (t) > 0 forall t > 0.Proposition 2 is based on the fact that Φ(λ0,...,λN ) (·) has at most N zeros(including multiplicities) and all of them are concentrated at t = 0.We may prove directly the following:Proposition 3. Let λ0, λ1 be complex numbers. Then(15) Φ(λ0,λ1) (−x) Φ(λ2,λ3) (−x) = e−(λ0+λ1+λ2+λ3)xΦ(λ0,λ1) (x) Φ(λ2,λ3) (x) .The function Φ(λ0,λ1) (x) Φ(λ2,λ3) (x) is even if and only if λ0 + · · ·+ λ3 = 0.3. Computation of the matrix R. We assume that t1 < · · · < tn are realnumbers such that(16) Φ(λ0,λ1) (tj+1 − tj) 6= 0 and Φ(λ2,λ3) (tj+1 − tj) 6= 0for j = 1, . . . , n−1. In the general case condition (16) holds if tj− tj is sufficientlysmall for j = 1, . . . , n − 1. For complex numbers λ0, . . . , λ3 we define the lineardifferential operatorsL1 =(ddt− λ0)(ddt− λ1)and L2 =(ddt− λ2)(ddt− λ3).The following is the main result of this section:Theorem 4. Let λ0, . . . , λ3 be complex numbers and t1 < · · · < tn such that(16) holds, and let g be a natural spline for Λ3 = (λ0, λ1, λ2, λ3) with interpolationdata g1, . . . , gn . Then there exist a n×(n− 2) matrix Q and a symmetric (n− 2)×(n− 2) matrix R such that(17) QTg = Rγ.Theorem 5. Let λ0, . . . , λ3 be complex numbers and t1 < · · · < tn such that(16) holds, and let g be a natural spline for the operator L1L2 with interpolationdata g1, . . . , gn . Then the tridiagonal matrix R is given byRj,j = ρ (tj − tj−1)− ρ (− (tj+1 − tj))Rj,j+1 = σ (tj+1 − tj) and Rj+1,j = −σ (− (tj+1 − tj))and ρ and σ are defined in (11) and (12).In the same manner one proves the following result:Theorem 6. The n × (n− 2) matrix Q is given for i = 1, . . . , n and j =2, 3, . . . , n− 1, by(18) qj−1,j = −Φ′(λ0,λ1) (t− tj)Φ(λ0,λ1) (tj−1 − tj), qj+1,j =Φ′(λ0,λ1) (t− tj)Φ(λ0,λ1) (tj+1 − tj)andqjj =Φ′(λ0,λ1) (t− tj+1)Φ(λ0,λ1) (tj − tj+1)−Φ′(λ0,λ1) (t− tj−1)Φ(λ0,λ1) (tj − tj−1).C. R. Acad. Bulg. Sci., 75, No 1, 2022 154. Positivity of the entries of the matrix R. In the following we wantto analyze under which conditions the entries of the matrix R are real or evenpositive. If the vector ΛN = (λ0, . . . , λN ) is conjugation invariant, the fundamen-tal function ΦΛN is real-valued. If (λ0, λ1) and (λ2, λ3) are conjugation invariant,then (λ0, λ1, λ2, λ3) is conjuation invariant and we infer thatΦ(λ0,λ1),Φ(λ2,λ3),Φ(λ0,λ1,λ2,λ3) and Φ′(λ0,λ1,λ2,λ3)are real-valued functions. Hence under the additional assumption thatΦ(λ0,λ1) (t) 6= 0 and Φ(λ2,λ3) (t) 6= 0 for t ∈ (0, δ)for some δ > 0, the functions ρ (x) and σ (x) are real-valued and well-defined andall entries of the matrices R and Q are real-valued for all t1 < · · · < tn such that|tj+1 − tj | < δ for j = 1, . . . , n− 1.Corollary 7. Suppose that λ0, . . . , λ3 are real numbers. Then Rj,j and Rj,j+1and Rj−1,j are positive numbers for all t1 < · · · < tn.5. Diagonal dominance of the matrix R. The notion of diagonal domi-nance of matrices is very important in numerical analysis, cf. [20]. We start withthe following general criterion:Theorem 8. Assume that ρ (x) and −ρ (−x) are positive on the interval(0, δ) , and let Mδ > 0 be such that(19) |σ (x)| ≤Mδ |ρ (−x)| for all x ∈ [−δ, δ] .Then the estimate|Rj,j−1|+ |Rj,j+1| ≤Mδ |Rj,j |holds for all partitions t1 < · · · < tn such that tj+1 − tj ≤ δ for j = 1, . . . , n.We illustrate the result by two examples:Proposition 9. If Λ3 = (0, 0,−b, b) for some real b 6= 0, then Rj,j−1 +Rj,j+1 ≤12Rj,j for all t1 < · · · < tn.If we allow complex frequencies, we obtain a weaker result:Theorem 10. Let 0 < δ < π. Then for Λ3 = (0, 0,−iβ, iβ) with real β > 0we haveRj,j−1 +Rj,j+1 ≤sin δ − δδ cos δ − sin δRj,j < Rj,jfor all t1 < · · · < tn such that β (tj+1 − tj) ≤ δ for j = 1, . . . , n− 1.Theorem 11. Suppose that (λ0, λ1) and (λ2, λ3) are conjugation invariantand ε > 0. Then there exists δ > 0 such that the estimate|Rj,j−1|+ |Rj,j+1| ≤(12+ ε)|Rj,j |holds for all partitions t1 < · · · < tn such that tj+1 − tj ≤ δ for j = 1, . . . , n.16 O. Kounchev, H. Render, Ts. TsachevWe call the set ΛN symmetric if there exists a permutation π of the set{0, . . . , N} such that −λj = λπ(j) for j = 0, . . . , N, or, symbolically, −ΛN = ΛN .If ΛN is symmetric, thenΦΛN (−x) = (−1)N ΦΛN (x) .This formula shows that for odd N the complex-valued function ΦΛN is odd, andfor even N the function ΦΛN is even.One may prove the following main result concerning symmetric sets ΛN .Theorem 12. The matrix R for Λ3 = (λ0, . . . , λ3) is symmetric for anychoice of interpolation points t1 < · · · < tn if and only if Λ3 is symmetric.For a special class of examples which occur naturally in the study of polysplineson parallel strips (cf. [6]), we obtain an explicit estimate.Theorem 13. Let 0 < a < b and Λ3 = (a,−a,−b, b). Then the followingestimate holdsRj,j−1 +Rj,j+1 ≤12Rj,j.REFERENCES[1] de Boor C. (2001) A Practical Guide to Splines, Revised edition, Applied Mathe-matical Sciences, 27, New York, Springer-Verlag.[2] Schumaker L. L. (1981) Spline Functions: Basic Theory, Pure and Applied Math-ematics, A Wiley-Interscience Publication, New York, John Wiley & Sons, Inc.[3] Hastie T., R. Tibshirani, J. Friedman (2009) The elements of statistical learn-ing: Data Mining, Inference, and Prediction, Second edition, Springer Series inStatistics, New York, Springer.[4] Gu C. (2013) Smoothing Spline ANOVA Models, Second edition, Springer Series inStatistics, 297, New York, Springer.[5] Ramsay J. O., N. Heckman, B. W. Silverman (1997) Spline smoothing withmodel-based penalties, Behavior Research Methods, Instruments & Computers,29(1), 99–106.[6] Kounchev O. I. (2001) Multivariate Polysplines: Applications to Numerical andWavelet Analysis, London–San Diego, Academic Press.[7] Kounchev O., H. Render (2002) Wavelet Analysis of cardinal L-splines andConstruction of multivariate Prewavelets, In: Approximation theory X. Wavelets,splines, and applications (Eds Chui, K. Charles et al.), Innovations in Applied Math-ematics, Nashville, TN, Vanderbilt University Press, 333–353.[8] Kounchev O., H. Render (2005) Cardinal interpolation with polysplines on an-nuli, J. Approx. Theory, 137(1), 89–107.[9] Kounchev O., H. Render (2004) The Approximation order of Polysplines, Proc.Amer. Math. Soc., 132(2), 455–461.[10] Kounchev O., H. Render (2013) Polyharmonic functions of infinite order onannular regions, Tôhoku Math. J., 65(2), 199–229.2 C. R. Acad. Bulg. Sci., 75, No 1, 2022 17[11] Kounchev O., H. Render (2005) Polyharmonic splines on grids Z×aZn and theirlimits, Math. Comp., 74(252), 1831–1841.[12] Ramsay J. O., B. W. Silverman (2005) Functional Data Analysis, Second edition,Springer Series in Statistics, New York, Springer.[13] Green P. J., B. Silverman (1994) Nonparametric regression and generalized lin-ear models. A roughness penalty approach. Monographs on Statistics and AppliedProbability, 58, London, Chapman & Hall.[14] Reinsch C. H. (1967) Smoothing by Spline Functions. I, II, Numer. Math., 10,177–183.[15] Micchelli Ch. A. (1976) Cardinal L-splines, In: Studies in Spline Functions andApproximation Theory (Eds S. Karlin et al.), Academic Press, New York, 203–250.[16] Kounchev O., M. Wilson (2003) Application of PDE methods to Visualizationof Heart Data (Eds Michael J. Wilson and Ralph R. Martin), In: Mathematics ofSurfaces, Lecture Notes in Computer Science 2768, Springer-Verlag, 377–391.[17] Kounchev O., H. Render, Ts. Tsachev (2020) On a class of L-splines of order4: fast algorithms for interpolation and smoothing, BIT Numerical Mathematics,60(4), 879–899.[18] Kounchev O., H. Render (2021) Error estimates for interpolation with piecewiseexponential splines of order two and four, J. Comp. Appl. Math., 391, paper No113464.[19] Berenstein C. A., R. Gay (1995) Complex analysis and special topics in harmonicanalysis, New York, Springer-Verlag.[20] Golub G., Charles F. Van Loan (2012) Matrix Computations, Johns HopkinsUniversity Press.Institute of Mathematics and InformaticsBulgarian Academy of SciencesAkad. G. Bonchev St, Bl. 81113 Sofia, Bulgariae-mail: kounchev@math.bas.bgtsachev@math.bas.bg∗School of Mathematics and StatisticsUniversity College DublinBelfield, Dublin 4, Irelande-mail: hermann.render@ucd.ie18 O. Kounchev, H. Render, Ts. Tsachev

L-Spline Interpolation for Differential Operators of Order 4 with Constant Coefficients

http://sci-gems.math.bas.bg/jspui/bitstream/10525/4428/1/2022-01-02-ft.pdf

L-Spline Interpolation for Differential Operators of Order 4 with Constant Coefficients

Abstract

Similar works

Full text

Available Versions

Bulgarian Digital Mathematics Library at IMI-BAS