Computation of Induced Orthogonal Polynomial Distributions

27 Pages·2017·0.76 MB·English

by  Akil Narayan

ElectronicTransactionsonNumericalAnalysis.
ETNA
Volume50,pp.71–97,2018.
KentStateUniversityand
Copyright c 2018,KentStateUniversity. JohannRadonInstitute(RICAM)
(cid:13)
ISSN1068–9613.
DOI:10.1553/etna_vol50s71
COMPUTATIONOFINDUCEDORTHOGONALPOLYNOMIAL
DISTRIBUTIONS˚
AKILNARAYAN:
DedicatedtoWalterGautschiontheoccasionofhis90thbirthday
Abstract.Weprovidearobustandgeneralalgorithmforcomputingdistributionfunctionsassociatedtoinduced
orthogonalpolynomialmeasures. Weleverageseveraltoolsfororthogonalpolynomialstoprovideaspectrally-
accuratemethodforabroadclassofmeasures,encompassingthoseassociatedtoclassicalorthogonalpolynomial
families,whichisstableforpolynomialdegreesuptoatleast1000. Pairedwithotherstandardtoolssuchasa
numericalroot-findingalgorithmandinversetransformsampling,thisprovidesamethodologyforgeneratingrandom
samplesfromaninducedorthogonalpolynomialmeasure.Generatingsamplesfromthismeasureisoneingredientin
optimalnumericalmethodsforcertaintypesofmultivariatepolynomialapproximation.Forexample,samplingfrom
induceddistributionsforweighteddiscreteleast-squaresapproximationhasrecentlybeenshowntoyieldconvergence
guaranteeswithaminimalnumberofsamples.Wealsoprovidepublicly-availablecodethatimplementsthealgorithms
inthispaperforsamplingfrominduceddistributions.
Keywords.orthogonalpolynomials,induceddistributions,sampling
AMSsubjectclassifications.33C45,65D15
1. Introduction. Let µ be a probability measure on R such that a family of L2 -
dµ
orthonormalpolynomialstp u8 canbedefined.1 Thenon-decreasingfunction
n n“0
ż
x
F pxq“ p2ptqdµptq, xPR,
n n
´8
isaprobabilitydistributionfunctiononRsincep2 hasunitµ-integraloverR. Thispaper
n
ischieflyconcernedwithdevelopingalgorithmsfordrawingrandomsamplesofarandom
variable whose cumulative distribution function is F . The high-level algorithmic idea is
n
straightforward: developrobustalgorithmsforevaluatingF ,andsubsequentlyuseastandard
n
root-findingapproachtocomputeF´1pUqwhereU isacontinuousuniformrandomvariable
n
onr0,1s. (Thisiscolloquiallycalled“inversetransformsampling".) Thechallengethatthis
paperaddressesisinthecomputationalevaluationofF pxqforanynPN andforrelatively
n 0
generalµ. Borrowingterminologyfrom[8],wecallF theorder-ndistributioninducedbyµ.
n
In our algorithmic development, we focus on three classes of continuous measures µ
fromwhichinduceddistributionsspring: (1)Jacobidistributionsonr´1,1s,(2)Freud(i.e.,
exponential) distributions on R, and (3) “Half-line" Freud distributions on r0,8q. These
measuresencompassarelativelybroadselectionofcontinuousmeasuresµonR,including
allmeasuresassociatedwithclassicalorthogonalpolynomialfamilies. Ouralgorithmapplies
directlytothemeasureclassesidentifiedinTable1.1.
Theutilityofsamplingfromunivariateinduceddistributionshasrecentlycomeintolight:
Theauthorsinvariouspapers[11,19,3]notethatadditivemixturesofinduceddistributions
areoptimalsamplingdistributionsforconstructingmultivariatepolynomialapproximationsof
functionsusingweighteddiscreteleast-squaresfromindependentandidentically-distributed
randomsamples. “Optimal"meansthatthesedistributionsdefineasamplingstrategywhich
˚ReceivedApril27,2017.AcceptedSeptember25,2018.PublishedonlineonDecember4,2018.Recommended
byTomBella.A.NarayanispartiallysupportedbyNSFDMS-1720416,AFOSRFA9550-15-1-0467,andDARPA
EQUiPSN660011524053.
:DepartmentofMathematics,andScientificComputingandImaging(SCI)Institute,UniversityofUtah,USA
([email protected]).
1SeeSection2.1fortechnicalconditionsimplyingthis.
71
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
72 A.NARAYAN
TABLE1.1
Classesofmeasuresconsideredinthispaper.Γp¨qistheEulerGammafunction,andBp¨qistheBetafunction.
αą´1
Jacobi dµpα,βqpxq“ 1 p1´xqαp1`xqβ xPr´1,1s cpα,βq“2α`β`1Bpβ`1,α`1q
J cpJα,βq βą´1 J
αą0 ` ˘
Half-lineFreud dµpα,ρqpxq“ 1 xρexpp´xαq xě0 cpα,ρq“ 1Γ ρ`1
HF cpHαF,ρq ρą´1 HF α α
αą0 ` ˘
Freud dµpα,ρqpxq“ 1 |x|ρexpp´|x|αq xPR cpα,ρq“ 2Γ ρ`1
F cpFα,ρq ρą´1 F α α
providesstabilityandaccuracyguaranteeswithasamplecomplexitythatiscurrentlythought
tobethebest(smallest). Thisdistributionalsoarisesinrelatedsettings[12]. Theabilityto
samplefromaninduceddistribution, whichthispaperaddresses, thereforehassignificant
importanceformultivariateappliedapproximationproblems.
Induceddistributionscanalsohelpprovideinsightformoretheoreticalproblems. The
weighted pluripotential equilibrium measure is a multivariate probability measure that de-
scribesasymptoticdistributionsofoptimalsamplingpoints[1,19]. However,anexplicitform
for this measure is not known in general. The authors in [19] make conjectures about the
Lebesgue density associated to equilibrium measure in one case when its explicit form is
currentlyunknown. Whiletheseconjecturesremainunproven,univariateinduceddistributions
canbeusedtosimulatesamplesfromtheequilibriummeasure. Hence,samplingfrominduced
distributionscanbeusedtoprovidesupportingevidenceforthetheoreticalconjecturesin[19].
Theoutlineofthispaperisasfollows: inSection2wereviewmanystandardproperties
ofgeneralorthogonalpolynomialsystemsthatareexploitedforcomputinginduceddistribu-
tions. Section3containsadetaileddiscussionofournovelapproachforcomputingF pxq
n
for three classes of measures; this section also utilizes potential theory results in order to
approximateF´1p0.5q. Section4usestheprevioussection’salgorithmsinordertoformu-
n
lateanalgorithmicstrategyforcomputingF´1puq,u P r0,1s. Finally,Section5discusses
n
theabove-mentionedapplicationsofmultivariatepolynomialapproximationusingdiscrete
least-squares,andinvestigatingconjecturesforaweightedequilibriummeasure.
Code that reproduces many of the plots in this paper is available for download [17].
The code contains routines for accomplishing almost all of the procedures in this paper
includingevaluationandinversionofinduceddistributions(formanyofthedistributionsin
Table1.1),inversetransformsamplingformultivariatesamplingfromadditivemixturesof
induceddistributions,andfastversionsofallcodesthatutilizeapproximatemonotonespline
interpolantsforfastevaluationandinversionofdistributionfunctions. Thecodealsocontains
routinesthatreproduceFigure1.1(left,center),Figure2.1(right),Figure3.1(center,right),
Figure3.2(left),Figure3.3(left),andFigure5.1(left).
1.1. Asimpleexample. Themainalgorithmicnoveltiesofthispaperrevolvearound
evaluation of F pxq. In Figure 1.1 we show one example of the integrand p2pxqdµpxq
n n
and the associated F . One suspects that packaged integration routines should be able to
n
performrelativelywellinordertocomputeintegralsforsuchaproblem. Inourexperience
thisisfrequentlytrue,butcomesatapriceofincreasedcomputationaleffortandtime,and
decreasedrobustness. Theright-handpaneinFigure1.1showstimingsforMatlab’sbuilt-in
integralroutineversusthealgorithmsdevelopedinthispaper. Weseethatthealgorithms
inthispaperaremuchfaster,usuallyresultinginaroundanorderofmagnitudesavingswhen
compared against integral. We caution that such tests compare an adaptive numerical
integrationmethod,Matlab’sintegral,againstatailoredandoptimizedintegrationmethod,
thealgorithmofthispaper. Therefore,anadaptiveversionofouralgorithmislikelytobe
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
INDUCEDORTHOGONALPOLYNOMIALDISTRIBUTIONS 73
4 1 Matlab’sintegral
40 Induceddistributionalgorithm
3 0.8
2p(x)dµ(x)n2 F(x)n00..46 Time(ms)2300
1 0.2 10
0−1−0.8−0.6−0.4−0.2 x0 0.2 0.4 0.6 0.8 1 0−1−0.8−0.6−0.4−0.2 x0 0.2 0.4 0.6 0.8 1 00 20 40 60 n80 100 120 140
FIG.1.1.Left:Theinduced?distributionintegrandp2npxqdµpxqfortheJacobimeasureµpJα,βqinTable1.1
withparametersα“´0.8,β“ 101,andn“13.Center:TheassociateddistributionfunctionFnpxq.Right:
ComputationaltimingsforevaluationofFnp¨q(averagedover100runs)usingMatlab’sbuilt-inintegralroutine
comparedtothealgorithmsdevelopedinthispaper.TimingswereperformedinMatlab(2015b)onasingle-core1.7
GHzInteli5processorwith4GBofRAM.
lessefficientthanreportedhere,andthenon-adaptiveversioninthispaperisnaturallyvery
specializedandthereforelessusefulforgeneralintegrationproblems.
Ourexperimentation(usingMatlab’sintegral)alsorevealsthefollowingadvantages
ofusingthespecializedalgorithminthispaper:
‚ Theresultsfromouralgorithmappeartobemoreaccuratecomparedtostandard
routines, and use significantly less computation. This conclusion is based on our
testing with integral, and is true even if one modifies algorithm tolerances in
integral.
‚ Wehaveoccasionallyobservedintegralreturnnon-monotonicevaluationsfor
thedistributionF . Thistypicallyhappenswhenmostofthemassoftheintegrandis
n
concentratedfarfromtheboundaryofthesupportofµ. Non-monotonicbehavior
causesproblemsinperforminginversetransformsampling.
‚ Whendµpxqhasasingularityatboundariesofthesupportofµ,thenintegral
frequentlycomplainsaboutsingularitiesandfailuretoachieveerrortolerances.
SimilarresultsandconclusionsholdwhenthealgorithmofthispaperistestedagainstMat-
lab’squadgkroutine,anotherbuilt-inintegrationfunction. Wealsonotethattheintegrand
p2ptqdµptq can be highly oscillatory for large n, which suggests that asymptotics for this
n
integrandalongwithtechniquesfromoscillatoryquadraturecanbeeffectiveincomputing
theseintegrals. Thispaperdoesnotinvestigatesuchapproachesandinsteadfocusesonex-
ploitingpropertiesoforthogonalpolynomials. However,webelievethattoolsforoscillatory
quadraturecanbeeffectiveforthisproblem.
1.2. Historical discussion. The distribution function of the arcsine or “Chebyshev"
measureis
ż
1 x 1 1 1
Fpxq“ ? dt“ ` arcsinpxq, xPr´1,1s.
π 1´t2 2 π
´1
Itwasshownin[20]thatifµbelongstotheNevaiclassofmeasures,thenF pxq Ñ Fpxq
n
pointwiseforallxPr´1,1s. Furtherrefinementsonthisstatementweremadein[23]which
generalizedtheclassofmeasuresforwhichconvergenceholds. Generatingpolynomialsor-
thogonalwithrespecttothemeasureassociatedtoF isconsideredin[8]withageneralization
n
givenin[13].
Theauthorsin[11]proposedsamplingfromanadditivemixtureofinduceddistributions
using a Markov Chain Monte Carlo method for the purposes of computing polynomial
approximationsoffunctionsviadiscreteleast-squares;theworkin[19]investigatessampling
fromthen-asymptoticlimitoftheseadditivemixtures. Theauthorsin[3]leveragetheadditive
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
74 A.NARAYAN
mixturepropertytosamplefromthisdistributionusingamonotonesplineinterpolant. Atthe
veryleastthislattermethodrequiresmultipleevaluationsofF pxq. Toourknowledgethere
n
hasbeenessentiallynoinvestigationintorobustalgorithmsfortheevaluationofF forbroad
n
classesofmeasures,whichisthesubjectofthispaper.
2. Background.
2.1. Orthogonalpolynomials. Thissectioncontainsclassicalknowledge,mostofwhich
isavailablefromanyseminalreferenceonorthogonalpolynomials[24,4,20,6].
LetµbeaLebesgue-StieltjesprobabilitymeasureonR,i.e.,thedistributionfunction
ż
x
Fpxq“ dµptq,
´8
is non-decreasing and right-continuous on R, with Fp´8q “ 0 and Fp8q “ 1. For any
distributionF weusethenotationFcpxq:“1´Fpxqforitscomplementaryfunction.
Weassumethatµhasaninfinitenumberofpointsofincrease,andhasfinitepolynomial
momentsofallorders,i.e.,
ˇż ˇ
ˇ ˇ
ˇ ˇ
ˇ xndµpxqˇă8, n“0,1,....
R
Under these assumptions, a sequence of orthonormal polynomials tp u8 exists, with
n n“0
degp “j,satisfying
j
ż
p pxqp pxqdµpxq“δ ,
j k k,j
R
whereδ istheKroneckerdeltafunction. Wewillwritep p¨q“p p¨;µqtodenoteexplicit
k,j n n
dependence of p on µ when necessary. Such a family can be mechanically generated by
n
iterativeapplicationofathree-termrecurrencerelation:
a a
(2.1) xp pxq“ b p pxq`a p pxq` b p pxq,
n n n´1 n n n`1 n`1
wheretherecurrencecoefficientsa andb arefunctionsofthemomentsofµ. Theinitial
n n
conditionsp ”0andp ”1areusedtoseedtherecurrence. Withp definedinthisway,
´1 0 n
the(positive)leadingcoefficientofp hasvalue
n
źn 1
γ :“ a , p pxq“γ xn`¨¨¨ .
n n n
b
j“0 j
The polynomial p has n real-valued, distinct roots; when µ has a Lebesgue density with
n
supportgivenbyafiniteunionofintervals,thentherootslieinsidetheconvexhullofthe
supportofµ. Theserootstx un arenodesfortheGaussianquadraturerule,
k k“1
ż
ÿn (cid:32) (
fpxqdµpxq“ w fpx q, f Pspan 1,x,x2,...,x2n´1 ,
k k
R
k“1
wheretheweightsw areuniqueandpositive. Thesenodesandweightscanbecomputed
k
havingknowledgeonlyoftherecurrencecoefficientsa andb ;numerousmodernalgorithms
k k
accomplishthis,withahistoricallysignificantproceduregivenin[10].
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
INDUCEDORTHOGONALPOLYNOMIALDISTRIBUTIONS 75
2.2. Inducedorthogonalpolynomialsandmeasures. Withtp u8 theorthonormal
n n“0
polynomialfamilywithrespecttoµ,thecollectionofpolynomialsorthogonalwithrespectto
theweighteddistributionp2pxqdµpxqwithnfixedarecalledinducedorthogonalpolynomials.
n
Weadoptthisterminologyfrom[8].
DefinetheLebesgue-Stieltjesmeasureµ anditsassociateddistributionfunctionF by
n n
ż ż
x x
(2.2) F px;µq“F pxq:“ dµ ptq:“ p2ptqdµptq,
n n n n
´8 ´8
withtp p¨qu “tp p¨;µqu theorthonormalfamilyforµ. Notethatµ !µ,andF p8q“
n n n n n n
µ pRq “ 1 so that µ is also a probability measure. The measure µ has its own three-
n n n
termrecurrencecoefficientsa andb forj “ 0,...,thatdefineanewsetofL2 pRq-
j,n j,n µn
orthonormalpolynomials,whichcanbegeneratedthroughthecorrespondingversionofthe
mechanicalprocedure(2.1). Onesuchprocedureforgeneratingthesecoefficientsisgiven
in[8].
Wewillcallthemeasureµ the(order-n)inducedmeasureforµ,andF thecorrespond-
n n
ing(order-n)induceddistributionfunction.
Ourmaincomputationalgoalis,givenuPr0,1s,theevaluationofF´1puqforvarious
n
measures µ. The overall algorithm for accomplishing this is a root-finding method, e.g.,
bisectionorNewton’smethod. Thus,thegoaloffindingF´1puqalsoinvolvestheevaluation
n
of F pxq, which is the focus of Section 3. A good root-finding algorithm also requires a
n
reasonableinitialguessforthesolution. Thisinitialguessisprovidedbythemethodologyin
Section4.1.
2.3. Measuremodifications. Ouralgorithmsrelyontheabilitytocomputepolynomial
measuremodifications. Thatis,giventhethree-termcoefficientsa andb forµ,tocompute
n n
thecoefficientsra andrb forµrdefinedas
n n
dµrpxq“ppxqdµpxq,
whereppxqisapolynomial,non-negativeonthesupportofµ.Thisisawell-studiedproblem[9,
5,6,18]. Inparticular,onemayreducetheproblemtoiteratingovermodificationsbylinear
and quadratic polynomials. We describe in detail how to accomplish linear and quadratic
modificationsintheappendix,withtheparticulargoalofstructuringcomputationstoavoid
numericalunder-andover-flowwhennislarge.
ThefollowingcomputationaltasksdescribedintheAppendixareusedtoaccomplish
measuremodifications.
1. (AppendixA.1)Evaluationofr ,theratioofsuccessivepolynomialsintheorthogo-
n
nalsequence:
p pxq
r pxq:“ j .
j
p pxq
j´1
Above,werequirextolieoutsidethezerosetofp .
j´1
2. (AppendixA.2)Evaluationofanormalizedorweighteddegree-npolynomial:
p pxq
Cnpxq:“ bř n , ną0, xPR.
n´1p2pxq
j“0 j
NotethatC pxq{r pxq„1forlargeenough|x|.
n n
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
76 A.NARAYAN
50
4
Fn(x1)(x1≤x0)
Fbnc(x2)(x2>x0) 40
3
b
x) 30
(
µ
d 2 n
)
x
2(n 20
p
1
10 x0(n),α=385,β=π
x0(n),α=−1/π,β=100π
x0(n),α=0.5,β=2
ExactFn−1(0.5)
0
x1 x0≈Fn−1(0.5) x2 −1 −0.8−0.6−0.4−0.2 x0 0.2 0.4 0.6 0.8 1
FIG.2.1.Left:Avisualdescriptionofthehigh-levelalgorithmgivenbyequation(3.1).Shadedregionsindicate
definiteintegralvaluesFpnandFpnccomputedbythealgorithmsinthispaper.Right:Illustrationofexactvaluesof
Fn´1p0.5qversustheapproximatemediansx0pnqderivedin(3.4)and(3.13b)forJacobiandHalf-Freudweights,
respectively.Thevaluesofx0pnqaredesignedtobeproximaltoFn´1p0.5q.
3. (AppendixB)Polynomialmeasuremodifications: givenµanditsassociatedthree-
term recurrencecoefficients a and b , computation of thethree-term recurrence
n n
coefficientsassociatedwiththemeasuresµrandµrr,definedas
dµrpxq“˘px´y qdµpxq, y Rsuppµ,
0 0
dµrrpxq“px´z q2dµpxq, z PR.
0 0
The˘signinµrischosensothatµrisapositivemeasure.
3. EvaluationofF . Thissectiondevelopscomputationalalgorithmsfortheevaluation
n
oftheinduceddistributionF definedin(2.2). Thesealgorithmsdependfairlyheavilyon
n
theformofaLebesguedensitydµpxq(i.e.,apositiveweightfunction)forthemeasureµon
R,buttheideascanbegeneralizedtovariousmeasures. Weconsidertheclassesofweights
enumeratedinTable1.1:
‚ (Jacobi weights) dµ pxq “ p1 ´ xqαp1 ` xqβ for x P r´1,1s with parameters
J
α,β ą´1.
‚ (Freud weights) dµ pxq “ |x|ρexpp´|x|αq for x P R with parameters α ą 0,
F
ρą´1.
‚ (half-lineFreudweights)dµ pxq“xρexpp´xαqforxPr0,8qwithparameters
HF
αą0,ρą´1.
TheinduceddistributionF forFreudweightscanactuallybewrittenexplicitlyintermsof
n
thecorrespondinginduceddistributionforhalf-lineFreudweights,somostofthealgorithm
development concentrates on the Jacobi and half-line Freud cases. The strategy for these
twolattercasesisessentiallythesame: withµinoneoftheclassesaboveandnfixed,we
dividethecomputationintooneoftwoapproximations,dependingonthevalueofx. Each
approximationisaccurateforitscorrespondingvaluesofx. Formally,thealgorithmis
#
p
F pxq, xďx pµ,nq,
(3.1) Fnpxq“ 1´Fpncpxq, xąx0pµ,nq,
n 0
p p
whereF pxqandFcpxqrepresentcomputationalapproximationstoF pxqandFcpxq,respec-
n n n n
tively,andaretheoutputsfromthealgorithmsthatwewilldevelop. Apictorialdescriptionof
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
INDUCEDORTHOGONALPOLYNOMIALDISTRIBUTIONS 77
thisisgiveninFigure2.1(left). Wechoosetoseekavalueofx suchthatF´1p0.5q. Thisisa
0 n
heuristicchoicethatweobserveworkswellinpracticeacrossallthemeasuresweconsider. It
islikelythatotherchoicesforx canalsoproduceeffectivealgorithms,butwehavefoundour
0
choicetobesimpleandamenabletodirectcomputation. However,wecannotknowthisvalue
F´1p0.5qapriori,soweusepotential-theoreticargumentstocomputeavaluex “x pµ,nq
n 0 0
approximatingthemedianofF ,
n
1
(3.2) F px q« .
n 0
2
Forourchoiceofx ,wecanprovidenoestimatesforthisapproximateequality,butempirical
0
evidence in Figure 2.1 (right) shows that our choices are very close to the real median,
uniformlyinn. Thecomingsectionsconcentrateon,foreachclassofµmentionedabove,
p p
specifyingx anddetailingalgorithmsforF andFc.
0 n n
3.1. Jacobiweights. WeconsidercomputinginduceddistributionfunctionsforJacobi
measuresµpα,βqasdefinedinTable1.1. Whencircumstancesareclear,wewillwriteµpα,βq “
J J
µtoavoidnotationalclutter. Weseekthedistributionfunctionoftheinducedmeasure,
ż
x
(3.3) µpα,βqpr´1,xsq“µ pr´1,xsq“F pxq:“ p2ptqdµptq.
J,n n n n
´1
TocomputeF ,wespecifyanapproximatemedianx “x pµ,nqsatisfying(3.2),and
n 0 0
p p
construct algorithmic procedures to evaluate F pxq « F pxq (for x ď x ) and Fcpxq «
n n 0 n
Fcpxq(forx ą x ). Havingspecifiedthese,weuse(3.1)tocomputeourapproximationto
n 0
F .
n
3.1.1. Computationofx pnq. Withα,β ą´1andnPNfixed,considerthemeasure
0
µpα{2n,β{2nq. NotethatthismeasureisstillaJacobimeasuresince α ą ´1and β ą ´1.
J 2n 2n
Since
” ı ´ ¯
2n 2n
dµpα,βqptq9p1´tqαp1`tqβ “ p1´tqα{2np1`tqβ{2n 9 dµpα{2n,β{2nqptq ,
J J
theintegrandin(3.3)satisfies
” ´ ¯ ı
n 2
p2ptqdµpα,βqptq9 p ptq dµpα{2n,β{2nqptq .
n J n J
Thequantityunderthesquarebracketsontheright-handsideis,inthelanguageofpotential
theory,aweightedpolynomialofdegreen. Oneresultinpotentialtheorycharacterizesthe
“essential"supportofthisweightedpolynomial;inparticular,theweightedpolynomialdecays
quicklyoutsidethissupport. Theessentialsupportisaninterval,andwetakethemedianx
0
oftheinducedmeasureµ tobethecentroidofthisinterval.
n
TheessentialsupportoftheweightedpolynomialaboveisdemarcatedbytheMhaskar-
Rakhmanov-Saffnumbersfortheasymmetricweightdµpα{2n,β{2nqonr´1,1s.Thesenumbers
J
for this weight are computed explicitly in [22, pp 206-207]. When α,β ě 0, the support
intervalisrθ´∆,θ`∆s,with
a
β2´α2 4 npn`α`βqpn`αqpn`βq
θ “ , ∆“ .
p2n`α`βq2 p2n`α`βq2
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
78 A.NARAYAN
Sincethisistheintervalwheremostofthe“mass"oftheintegralin(3.3)lies,wesetx tobe
0
thecentroidθofthisinterval:
´ ¯
β2´α2
(3.4) x µpα,βq,n “ , α,β ą´1, ną0.
0 p2n`α`βq2
Notethatthedefinitionsofθand∆requireαandβtobenon-negative. Withoutmathematical
justification, weextendusageoftheformula(3.4)tovalidnegativevaluesofα,β aswell.
Figure2.1comparesx andF´1p0.5qforcertainchoicesofαandβ.
0 n
p
3.1.2. ComputationofF pxq. Firstassumethatxďx ,withx definedin(3.4),and
n 0 0
defineAPN as
0
A:“t|α|u ùñ α´APp´1,1q,
wheret¨uisthefloorfunction. Wetransformtheintegral(3.3)overr´1,xsontothestandard
intervalr´1,1sviathesubstitutionu“ 2 pt`1q´1:
x`1
ż
1 x
F pxq“ p2 ptqp1´tqαp1`tqβdt
n cpα,βq ´1 n
ˆJ ˙ ż ˆ ˙
x`1 β`1 cp0,βq 1 1 α´A
“ J 2´ pu`1qpx`1q U puqdµp0,βqpuq,
2 cpα,βq ´1 2 2n`A
J
whereU puqisadegree-p2n`Aqpolynomialgivenby
2n`A
(3.5)
ˆ ˙ˆ ˙
A
1 1
U puq“p2 pu`1qpx`1q´1 2´ pu`1qpx`1q
2n`A n 2 2
ˆ ˙ „ˆ ˙  ˆ ˆ ˙˙
x`1 2n`A źA 3´x źn 2 2
“γ2 ´u u´ px `1q´1 ,
n 2 1`x x`1 j,n
lk“oo1ooooooooomooooooooooonjl“oo1ooooooooooooooooooooomooooooooooooooooooooooon
paq pbq
where tx un are the n zeros of p p¨q. Since we explicitly know the polynomial roots
j,n j“1 n
of U , we can absorb the factor marked paq into the measure dµp0,βqpuq via A linear
2n`A
modifications,andwecanlikewiseabsorbthefactorpbqvianquadraticmodifications. (See
AppendixB.)Thus,define
(3.6) dµr puq“U puqdµp0,βqpuq,
n 2n`A
whichisamodifiedmeasurewhoserecurrencecoefficientscanbecomputedviasuccessive
applicationofthelinearandquadraticmodificationmethodsinAppendixB.Thus,
ˆ ˙ ż ˆ ˙
x`1 β`1 cp0,βq 1 1 α´A
(3.7) F pxq“ J 2´ pu`1qpx`1q dµr puq.
n 2 cpα,βq ´1 2 n
J
Theintegrandabovehasaroot(αą0)orsingularity(αă0)atu“ 3´x “1`21´x ě1`
x`1 1`x
21´x0;thisrootisfaroutsidetheintervalr´1,1sunlessβ isverylargeandbothnandαare
1`x0
small. Theintegrandisthereforeapositive,monotonic,smoothfunctiononr´1,1s,taking
valuesbetween1´xand2;weuseanorder-M µr -Gaussianquadraturetoefficientlyevaluate
n
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
INDUCEDORTHOGONALPOLYNOMIALDISTRIBUTIONS 79
it. Withpur ,wr qM denotingthenodesandweights,respectively,ofthisquadraturerule,
m m m“1
wecompute
ˆ ˙
ÿM 1 α´A
(3.8) Ipxq:“ wr 2´ pur `1qpx`1q ,
m m
2
ˆm“1 ˙
β`1
p x`1 Ipxq
(3.9) F pxq“ .
n 2 2apβ`1qBpβ`1,α`1q
TheentireprocedureissummarizedinAlgorithm1.
p
InordertocomputeFcpxqforxąx ,weusesymmetry. Sincexдxinterchanges
n 0
theparametersαandβ,then
ż ” ı ż ” ı
p 1 2 ´x 2 p
Fpα,βq,cpxq“ ppα,βqptq dµpα,βqptq“ ppβ,αqptq dµpβ,αqptq“Fpβ,αqp´xq.
n n n n
x ´1
` ˘ ` ˘
p
Notethatifx ą x µpα,βq,n ,then´x ă x µpβ,αq,n . Thus,Fc canbecomputedvia
0 0 n
p
thesamealgorithmforF ,butwithdifferentvaluesforα,β,andx. Thisisalsoshownin
n
Algorithm1.
p
Algorithm1:ComputationofF pxq,approximatingF pxqforµcorrespondingtoaJacobi
n n
polynomialmeasure.
input :α,β ą´1: Jacobipolynomialparameters.
input :nPN andxPr´1,1s: Orderofinducedpolynomialandmeasureµ and
0 n
valuex.
input :M `PN: Quad˘ratureorderforapproximatecomputationofFnpxq.
p
output:F x;µpα,βq .
n
p
1 Ifxąx0pµpα,βq,nq,return1´Fnp´x`;µpβ,αqq˘.
2 Computenzeros,txj,nunj“1ofpn “pn ¨;µpα,βq ,andleadingcoefficientγnofpn.
3 Computerecurrencecoefficientsaj andbj associatedtoµp0,βqfor
0ďj ďM `A`2n.
4 forj “1,...,ndo
5 Quadraticmeasuremodification(B.1b): updatean,bnfor
n“0,...,M `A`2pn´jqwithmodificationfactor
´ ´ ¯¯
2
u´ 2 px `1q´1 .
x`1 j,n
6 end
7 fork “1,...,Ado
8 Uselinearmodificatio´n(B.´1a)to¯u¯pdatean,bnforn“0,...,M `pA´kqwith
modificationfactor u´ 3´x .
1`x
9 end ` ˘
10 b0 Ðb0 x`21 2n`Aγn2 (Cf. theleadingcoefficientofthepolynomialU2n`Ain(3.5).)
11 ComputeM-pointGaussquadraturepurm,wrmqMm“1associatedwithmeasureµrnvia
tpa ,b quM .
j j j“0 ` ˘
p
12 ComputetheintegralI in(3.8),andreturnFn x;µpα,βq givenby(3.9).
THEOREM3.1. Withµpα,βq,nPN,andxPr´1,1sallgiven,assumethatxďx0with
p
x asin(3.4). ThentheoutputF pxqfromAlgorithm1usinganM-pointquadraturerule
0 n
ETNA
KentStateUniversityand
JohannRadonInstitute(RICAM)
80 A.NARAYAN
satisfies
(3.10) ˇˇˇF pxq´Fp pxqˇˇˇďCpα,β,n,MqźM b pµr q,
n n j n
j“0
where b pµr q are the b three-term recurrence coefficients associated to the x-dependent
j n j
measureµr definedin(3.6). TheconstantC is
n
ˆ ˙
2β`1´A x pnq`1 2M`β`1
0
Cpα,β,n,Mq“ .
pβ`1qBpβ`1,α`1q 4
NotethatCpα,β,n,Mqontheright-handsideof(3.10)isexplicitlycomputableonceα,
β,M,andnarefixed,andonlytheproductinvolvingb pµr qdependsonx. Furthermore,the
j n
b pµr qfactorsareexplicitlycomputedinAlgorithm1,sothatarigorouserrorestimatefor
j n
thealgorithmcanbecomputedbeforeitstermination.
Since x0`1 ă 1,thentheestimate(3.10)alsohintsatexponentialconvergenceofthe
4 2
quadraturerulestrategy,assumingthatthefactorsb pµr qcanbeboundedorcontrolled. We
j n
cannotprovidethesebounds, althoughwedoknowtheasymptoticbehaviorb pµr q Ñ 1
j n 4
asj Ñ 8[20,Remark3.1.10]. Also,sincex P r´1,1sforanyn ą 0then,uniformlyin
0
n ą 0, Cpα,β,n,Mq ď C1pα,βq4´M. Thus, forafixedxweexpectthattheestimatein
Theorem3.1behaveslike
ˇˇˇF pxq´Fp pxqˇˇˇďCpα,β,n,MqźM b pµr qÀC2pα,β,xq4´2M,
n n j n
j“0
showingexponentialconvergencewithrespecttoM. However,wecannotprovethislatter
statement.
ProofofTheorem3.1. Theresultisarelativelystraightforwardapplicationofknownerror
estimatesforGaussianquadraturewithrespecttonon-classicalweights. Weusethenotation
ofAlgorithm1: ur andwr denotetheM-pointµr -Gaussianquadraturenodesandweights,
m m n
respectively. WestartwiththeCorollarytoTheorem1.48in[6],statingthatiffp¨qisinfinitely
differentiableonr´1,1s,then
ˇ ˇ
ˇˇż1 ÿM ˇˇ fp2Mqpτqż1 źM
ˇ fpuqdµr puq´ wr fpur qˇ“ pu´ur q2dµr puq
ˇ n m m ˇ p2Mq! m n
´1 j“1 ´1j“1
forsomeτ P p´1,1q. Notingthatsinceur arethezerosofthedegree-M µr -orthogonal
m n
polynomial,then
ż ż
1 źM źM 1 źM
pu´ur q2dµr puq“ b pµr q p2 pu;µr qdµr puq“ b pµr q.
m n j n M n n j n
´1j“1 j“0 ´1 j“0
From (3.7), the integral that we would like to approximate has the integrand
` ˘
fpuq“ 2´ 1pu`1qpx`1q α´A. Thenforanyτ Pp´1,1qandanyxďx ,
2 0
ˇˇ ˇˇ ˆx`1˙2Mˆ 1 ˙α´A´2M2Mź´1
ˇfp2Mqpτqˇ“ 2´ px`1qpτ `1q |α´A´j|
2 2
j“0
ˆ ˙ ˆ ˙
x `1 2M 2Mź´1 x `1 2M2Mź´1
ď 0 2α´A´2M |α´A´j|“2α´A 0 |α´A´j|.
2 4
j“0 j“0