Fast associated classical orthogonal polynomial transforms
FFast associated classical orthogonal polynomial transforms B ROCK K LIPPENSTEIN AND R ICHARD M IKAËL S LEVINSKY * Department of Mathematics, University of Manitoba, Winnipeg, CanadaFebruary 17, 2021
Abstract
We discuss a fast approximate solution to the associated classical – classical orthogonal polynomial connection problem.We first show that associated classical orthogonal polynomials are solutions to a fourth-order quadratic eigenvalue problemwith polynomial coefficients such that the differential operator is degree-preserving. Upon linearization, the discretization ofthis quadratic eigenvalue problem is block upper-triangular and banded. After a perfect shuffle, we extend a divide-and-conquerapproach to the upper-triangular and banded generalized eigenvalue problem to the blocked case, which may be acceleratedby one of a few different algorithms. Associated orthogonal polynomials arise from iterated Stieltjes transforms of orthogonalpolynomials; hence, fast approximate conversion to classical cases combined with fast discrete sine and cosine transformsprovides a modular mechanism for synthesis of singular integral transforms of classical orthogonal polynomial expansions.
Let L ( D, d µ ) denote the Hilbert space of square integrable functions on D = supp( µ ) ⊂ R with positive Borel measure d µ with inner-product (cid:104) f, g (cid:105) = (cid:90) D f g d µ . We shall denote by { p n ( x ) } ∞ n =0 an orthogonal polynomial sequence. Orthogonalpolynomials satisfy a three-term recurrence relation which may be cast in the following form: p n +1 ( x ) = ( A n x + B n ) p n ( x ) − C n p n − ( x ) , p ( x ) = 1 , p − ( x ) = 0 , (1)where A n − A n C n > for n ≥ .Given a family of orthogonal polynomials, the associated orthogonal polynomials are those polynomials with the same initialconditions that use the recurrence coefficients with indices offset by c ∈ N units: p n +1 ( x ; c ) = ( A n + c x + B n + c ) p n ( x ; c ) − C n + c p n − ( x ; c ) . (2)Thanks to Favard’s theorem [1], the associated polynomials are indeed orthogonal with respect to non-negative measures µ ( x ; c ) with D = supp( µ ( · , c )) ⊂ R .Associated orthogonal polynomials diagonalize the following integral transforms [2, vol. 2, p. 162 (6)] with a removablesingularity: p n − ( x ; c + 1) = 1 A c (cid:82) D d µ ( t ; c ) (cid:90) D p n ( x ; c ) − p n ( t ; c ) x − t d µ ( t ; c ) . (3)Classical orthogonal polynomials are characterized by Bochner [3] and Krall [4] as the polynomial solutions of the degree-preserving second-order linear differential equation: (cid:0) σ D + τ D (cid:1) p n = λ n p n , (4)where σ and τ are polynomials in x independent of n that satisfy deg( σ ) ≤ and deg( τ ) ≤ , and the eigenvalues λ n = n [( n − σ (cid:48)(cid:48) + 2 τ (cid:48) ] . We call this factorization degree-preserving because the degrees of the polynomial variable coefficientsdo not exceed the orders of the respective differential operators. Table 1 summarizes this classical characterization data. * Corresponding author. Email: [email protected] a r X i v : . [ m a t h . NA ] F e b Table 1: Classical Orthogonal Polynomial data for Eqs. (1) and (4).Name Jacobi Laguerre HermiteHilbert space L ([ − , , (1 − x ) α (1 + x ) β d x ) L ( R + , x α e − x d x ) L ( R , e − x d x ) Special notation P ( α,β ) n ( x ) L ( α ) n ( x ) H n ( x ) A n (2 n + α + β + 1)(2 n + α + β + 2)2( n + 1)( n + α + β + 1) − n + 1 2 B n ( α − β )(2 n + α + β + 1)2( n + 1)( n + α + β + 1)(2 n + α + β ) 2 n + α + 1 n + 1 0 C n ( n + α )( n + β )(2 n + α + β + 2)( n + 1)( n + α + β + 1)(2 n + α + β ) n + αn + 1 2 nh n = (cid:104) p n , p n (cid:105) α + β +1 Γ( n + α + 1)Γ( n + β + 1)(2 n + α + β + 1)Γ( n + α + β + 1) n ! Γ( n + α + 1) n ! √ π n n ! σ x − − x − τ α − β + ( α + β + 2) x x − α − xλ n n ( n + α + β + 1) n n There is an appreciable literature on the associated classical orthogonal polynomial differential equation. First, it is built fromspecial cases for Hermite and Laguerre polynomials [5]. This is followed upon by the Jacobi polynomial case [6] for a completeclassical picture. The differential equation is linear and fourth-order with polynomial variable coefficients that depend on both x and n ; hence the associated classical polynomials are members of the Laguerre–Hahn family [7, 8].Zarzo, Ronveaux, and Godoy [9] present the degree-preserving differential equation in terms of some of the classical orthogonalpolynomial data, notably σ and τ . It reads: (cid:8) − σ D − σσ (cid:48) D + (cid:2) τ − τ σ (cid:48) − σ (cid:48) + (2 n + 4 c ) στ (cid:48) − (4 + n − n + 4 c − nc − c ) σσ (cid:48)(cid:48) (cid:3) D + (cid:2) τ τ (cid:48) + (2 n − c ) σ (cid:48) τ (cid:48) − τ σ (cid:48)(cid:48) + ( n − n − c + 2 nc + 2 c ) σ (cid:48) σ (cid:48)(cid:48) (cid:3) D (cid:9) p n ( x ; c ) (5) = n ( n + 2) [( n + 2 c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] [( n + 2 c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] p n ( x ; c ) . Foupouagnigni, Koepf, and Ronveaux [10] refine this presentation by including the classical eigenvalues in the description: (cid:8) − σ D − σσ (cid:48) D + (cid:2) τ + 2 τ (cid:48) σ − τ σ (cid:48) − σσ (cid:48)(cid:48) + 2( λ n + c + λ c − ) σ − σ (cid:48) (cid:3) D +3 [ τ τ (cid:48) + ( λ n + c + λ c − ) σ (cid:48) − ( τ + σ (cid:48) ) σ (cid:48)(cid:48) ] D} p n ( x ; c ) (6) = (cid:2) ( λ n + c − λ c − ) − ( λ n + c + λ c − ) σ (cid:48)(cid:48) + τ (cid:48) ( σ (cid:48)(cid:48) − τ (cid:48) ) (cid:3) p n ( x ; c ) . In this way, Eq. (6) is the associated analogue of Eq. (4).There has also been significant interest in factorizing the fourth-order differential equation. For the first order of association, c = 1 , Ronveaux [11] presents a homogeneous degree-preserving factorization of Eq. (6): (cid:2) σ D + ( τ + σ (cid:48) ) D + τ (cid:48) − λ n +1 (cid:3) (cid:2) σ D + (2 σ (cid:48) − τ ) D + σ (cid:48)(cid:48) − τ (cid:48) − λ n +1 (cid:3) p n ( x ; 1) = 0 . (7)More generally, Lewanowicz [12] shows that the same differential operators in Eq. (7) factorize Eq. (6) incompletely: (cid:2) σ D + ( τ + σ (cid:48) ) D + τ (cid:48) − λ n +1 (cid:3) (cid:2) σ D + (2 σ (cid:48) − τ ) D + σ (cid:48)(cid:48) − τ (cid:48) − λ n +1 (cid:3) p n ( x ; c )= ( c −
1) [( n + c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] (cid:2) σ D + 3 σ (cid:48) D − n ( n + 2) σ (cid:48)(cid:48) (cid:3) p n ( x ; c ) . (8)Foupouagnigni, Koepf, and Ronveaux [10] also present a factorization of the fourth-order differential equation that is notdegree-preserving, as the variable coefficients are also described by the classical orthogonal polynomials: (cid:8) σp c − D + (cid:2) ( τ + σ (cid:48) ) p c − − σp (cid:48) c − (cid:3) D − (cid:2) ( λ n + c − λ c − − τ (cid:48) ) p c − − τ p (cid:48) c − (cid:3)(cid:9) × (cid:8) σp c − D − p c − (cid:2) ( τ − σ (cid:48) ) p c − + 2 σp (cid:48) c − (cid:3) D (9) + (cid:2) τ − σ (cid:48) ) p c − p (cid:48) c − + 2 σp (cid:48) c − + ( σ (cid:48)(cid:48) − τ (cid:48) − λ n + c − λ c − ) p c − (cid:3)(cid:9) p n ( x ; c ) = 0; We shall require a refinement of Eq. (6).
Theorem 1.1.
Let: A = − σ D − σσ (cid:48) D + (cid:2) τ + 2 τ (cid:48) σ − τ σ (cid:48) − σσ (cid:48)(cid:48) + 4 λ c − σ − σ (cid:48) (cid:3) D + 3 [ τ τ (cid:48) + 2 λ c − σ (cid:48) − ( τ + σ (cid:48) ) σ (cid:48)(cid:48) ] D + [2 λ c − σ (cid:48)(cid:48) − τ (cid:48) ( σ (cid:48)(cid:48) − τ (cid:48) )] I , and (10) B = 2 σ D + 3 σ (cid:48) D + σ (cid:48)(cid:48) I . (11) The quadratic eigenvalues { λ ± n } ∞ n =0 of the problem defined by: {A + λ ± n B} p ± n ( x ; c ) = ( λ ± n ) p ± n ( x ; c ) . (12) satisfy: λ ± n = σ (cid:48)(cid:48) ( n + 1) ± ( n + 1) (cid:112) ( σ (cid:48)(cid:48) − τ (cid:48) ) + 8 λ c − σ (cid:48)(cid:48) , (13) = ± ( n + 1) { [2 c − ± ( n + 1)] σ (cid:48)(cid:48) + 2 τ (cid:48) } . (14) In particular, it follows from λ + n = λ n + c − λ c − > that p + n ( x ; c ) = p n ( x ; c ) , the associated classical orthogonal polynomials,and Eq. (6) is alternatively formulated as a quadratic eigenvalue problem with known eigenvalues.Proof. As in the classical case of Eq. (4), the eigenvalues of Eq. (12) are found by applying A and B to x n and extracting theleading coefficient. We find: − σ (cid:48)(cid:48) n ( n + 1) ( n + 2)4 + (cid:0) τ (cid:48) − τ (cid:48) σ (cid:48)(cid:48) + 2 λ c − σ (cid:48)(cid:48) (cid:1) ( n + 1) + σ (cid:48)(cid:48) ( n + 1) λ ± n = ( λ ± n ) . Eq. (13) solves this quadratic equation, and when we write λ c − in terms of c , σ (cid:48)(cid:48) and τ (cid:48) , we find that the discriminant is aperfect square, concluding with the alternative form in Eq. (14). Finally, it is readily confirmed that: λ + n = ( n + 1)2 [( n + 2 c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] , = ( n + c − c + 1)2 [( n + 2 c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] , = ( n + c )2 [( n + c − c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] − ( c − n + c + c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] , = ( n + c )2 [( n + c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] − ( c − c − σ (cid:48)(cid:48) + 2 τ (cid:48) ] = λ n + c − λ c − . The full spectrum of Eq. (12) is helpful in determining if the positive family of polynomial eigenfunctions is linearly indepen-dent of the negative family. As degree-graded polynomials form a flag, it is also the case for each set of eigenfunctions, and if λ + n (cid:54) = λ − n then the solutions are distinct. Based on the eigenvalues for the particular classical cases, we can conclude that bothsets of eigenfunctions are linearly independent in the Laguerre and Hermite cases ( σ (cid:48)(cid:48) = 0 ) as: λ − n < λ − n − < · · · < λ − < < λ +0 < · · · < λ + n − < λ + n , (15)but Jacobi polynomials require further investigation. If we let γ = α + β + 2 c − , then a nice formula for the eigenvalues is: λ ± n = ( n + 1)( n ± γ + 1) , (16)and it follows that λ − n = λ + n − γ . In contrast to the Laguerre and Hermite cases, λ − n are asymptotically positive and if γ ∈ Z ,then there are countably infinite equal eigenvalues among the two families. The condition γ ∈ Z corresponds to a set of lines inthe αβ -Jacobi parameter plane. In particular, if γ = 0 , then the linear independence of both families degenerates as eigenvaluesof the polynomial solutions of the same degree are equal. This occurs only if c = 1 and α + β = − .The orthogonal polynomial connection problem between families { p n ( x ) } ∞ n =0 and { q (cid:96) ( x ) } ∞ (cid:96) =0 is to find connection coefficientsorganized in an upper-triangular matrix V such that: p n ( x ) = n (cid:88) (cid:96) =0 V (cid:96),n q (cid:96) ( x ) . (17) Fast orthogonal polynomial transforms have a rich history [13]. Notable methods include: the use of asymptotics with rigorouserror bounds to related synthesis and analysis to discrete sine and cosine transforms [14–17]; hierarchical and fast multipolemethods (FMMs) [18–20]; Hadamard matrix factorizations [21]; and, divide-and-conquer methods for structured eigenvalueproblems [22, 23].To develop fast algorithms for the associated classical – classical connection problem, it will be important to review the rathergeneral method of Olver, Slevinsky, and Townsend [23] for the classical cases. By using multiplication by x , differentiation,raising, and lowering operators for classical orthogonal polynomials, it is shown that Eq. (4) can be expressed as: (cid:0) σ D + τ D (cid:1) (cid:0) p p · · · (cid:1) = (cid:0) p p · · · (cid:1) Λ , (18) (cid:0) σ D + τ D (cid:1) (cid:0) q q · · · (cid:1) V = (cid:0) q q · · · (cid:1) V Λ , (19) (cid:0) r r · · · (cid:1) AV = (cid:0) r r · · · (cid:1) BV Λ , (20)where { r n ( x ) } ∞ n =0 is yet another classical family and A and B are upper-triangular and banded matrices. Dropping thepolynomials on both sides of Eq. (20), we obtain an upper-triangular and banded generalized eigenvalue problem for theconnection coefficients, AV = BV Λ .For example, describing the Jacobi–Jacobi connection problem in terms of the matrices in Appendix A, we find: (cid:16) σ D + τ ( α,β ) D (cid:17) (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) = (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) Λ , (21) (cid:16) σ D + τ ( α,β ) D (cid:17) (cid:16) P ( γ,δ )0 P ( γ,δ )1 · · · (cid:17) V = (cid:16) P ( γ,δ )0 P ( γ,δ )1 · · · (cid:17) V Λ , (22) (cid:16) P ( γ +1 ,δ +1)0 P ( γ +1 ,δ +1)1 · · · (cid:17) (cid:104) − L ( γ +1 ,δ +1)( γ +2 ,δ +2) D ( γ +2 ,δ +2)( γ,δ ) + τ ( α,β ) ( M ( γ +1 ,δ +1) ) D ( γ +1 ,δ +1)( γ,δ ) (cid:105) V = (cid:16) P ( γ +1 ,δ +1)0 P ( γ +1 ,δ +1)1 · · · (cid:17) R ( γ +1 ,δ +1)( γ,δ ) V Λ . (23)Appendix B compiles the analogous identities for Laguerre polynomials.The connection coefficients convert not only the polynomials but also their expansions via: (cid:0) p p · · · (cid:1) c c ... = (cid:0) q q · · · (cid:1) V c c ... = (cid:0) q q · · · (cid:1) d d ... , (24)and the implied matrix-vector product d = V c .We will extend this technique to any of the associated classical differential equations (5) – (12) to find structured forms for theassociated classical – classical connection problem. We have shown that the classical orthogonal polynomial connection coefficients are generalized eigenvectors of upper-triangularand banded eigenproblems of the form: AV = BV Λ . (25)It is clear that the generalized eigenvalues are the ratios of the diagonal entries of the matrices, Λ i,i = A i,i /B i,i , and we mayapproach the generalized eigenvectors by dividing and conquering.In what follows, we consider the n × n truncation of this system. That is, we consider A, B ∈ R n × n . Let s = (cid:98) n (cid:99) and let b bethe upper bandwidth of both A and B . We suppose that V has the form: V = (cid:18) V , V , (cid:19) (cid:18) I V , I (cid:19) , (26)where V , ∈ R s × s , V , ∈ R s × ( n − s ) , and V , ∈ R ( n − s ) × ( n − s ) are all roughly the same size. Block dividing A , B , and Λ conformably, we have: (cid:18) A , A , A , (cid:19) (cid:18) V , V , V , V , (cid:19) = (cid:18) B , B , B , (cid:19) (cid:18) V , Λ V , V , Λ V , Λ (cid:19) . This decomposes into the two diagonal subproblems A i,i V i,i = B i,i V i,i Λ i , for i = 1 , and the off-diagonal problem: A , V , V , + A , V , = B , V , V , Λ + B , V , Λ , a two-term matrix equation for the block V , .Since V − , B − , A , V , = Λ , the matrix equation can be cast into the following form: Λ V , − V , Λ = V − , B − , ( B , V , Λ − A , V , ) . Since A and B are banded with an upper bandwidth of b , it follows that A , and B , are matrices with only a small number ofnonzero entries, arranged in a lower-triangular fashion in their bottom left corners. Thus, the matrix-matrix product A , V , results in a rank- b matrix with nonzero entries only in the last b rows. Similarly, due to the zero pattern in the matrix A , V , − B , V , Λ only the last b columns of B − , are needed, which are calculated with back substitution using the last b columns ofthe s × s identity. We represent this block as the outer product of an s × b matrix X and a b × ( n − s ) matrix Y : V − , B − , ( A , V , − B , V , Λ ) =: XY (cid:62) . To conquer this eigenproblem, we must then solve the diagonal Sylvester matrix equation: Λ V , − V , Λ = − XY (cid:62) . (27)In [23, 24], the solution is examined component-wise where it is found to be a Hadamard product of the rank- b matrix − XY (cid:62) and a Cauchy matrix, which is approximated hierarchically: ( V , ) (cid:96),n = ( − XY (cid:62) ) (cid:96),n (Λ ) (cid:96),(cid:96) − (Λ ) n,n . (28)Due to the separation of spectra in Λ and Λ , this hierarchical approach offers an approximate matrix-vector product involving V , in O ( bn log n log(1 /(cid:15) )) flops and thereby V in O ( bn log n log(1 /(cid:15) )) flops. By the FMM, the approximate matrix-vector product with V , may be further accelerated to O ( bn log(1 /(cid:15) )) flops and thereby V in O ( bn log n log(1 /(cid:15) )) flops. As (Λ ) (cid:96),(cid:96) < (Λ ) n,n for all (cid:96) and n , the integral representation: V , = (cid:90) ∞ exp( t Λ ) XY (cid:62) exp( − t Λ ) d t, (29)may be approximated via quadrature rules in order to accelerate the linear algebra [24]. Eq. (29) is still valid if Λ and Λ arenot diagonal matrices, but the condition on the spectra remains.Diagonal matrix equations such as Eq. (27) may also be approximated by the factored alternating direction implicit (fADI)method [25–27]. For normal Sylvester matrix equations, the precise number of shifts and their coordinates are computablevalues related to the Jacobian elliptic function and elliptic integral solution of Zolotarev’s third problem [28, 29], that ofminimizing the maximum absolute value of a rational function on one set divided by its minimum absolute value on another— sets which enclose the disjoint spectra of Λ and Λ . For nonnormal Sylvester matrix equations, the relationship betweenthe shifts and the error after J steps of fADI is more complicated.The structured form of the generalized eigenvectors permits inversion and transposition with the same computational complex-ity. For example: V − = (cid:18) I − V , I (cid:19) (cid:18) V − , V − , (cid:19) ,V (cid:62) = (cid:18) IV (cid:62) , I (cid:19) (cid:18) V (cid:62) , V (cid:62) , (cid:19) ,V −(cid:62) = (cid:18) V −(cid:62) , V −(cid:62) , (cid:19) (cid:18) I − V (cid:62) , I (cid:19) . With the structured forms for the eigenvectors, any submultiplicative operator norm (cid:107)·(cid:107) can be estimated recursively via: (cid:107) V (cid:107) ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) V , V , (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) I V , I (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) , ≤ max {(cid:107) V , (cid:107) , (cid:107) V , (cid:107)} ( (cid:107) I n (cid:107) + (cid:107) V , (cid:107) ) . Estimates on the induced -norms of the leaves follow from Hölder’s inequality and by comparison with the Frobenius normso that singular values need not be computed [30]. Estimates for the norms of the inverse transforms (and thereby conditionnumber estimates) also follow naturally from the recursive argument above. We consider the quadratic eigenvalue problem in Eq. (12). Expanding both { p ± n ( x ; c ) } ∞ n =0 in the same classical orthogonalpolynomial basis, say { q (cid:96) ( x ) } ∞ (cid:96) =0 , the resulting discretizations are matrix equations of the form: AV − + BV − Λ − = CV − (Λ − ) , and (30) AV + + BV + Λ + = CV + (Λ + ) , (31)where A, B, C are upper-triangular and banded matrices, V ± are the upper-triangular connection coefficients, and Λ ± are therespective diagonal quadratic eigenvalue matrices. For example, with Jacobi polynomials, the same procedure we use for theclassical transforms results in: A = − L ( γ +2 ,δ +2)( γ +4 ,δ +4) D ( γ +4 ,δ +4)( γ,δ ) + 5 σ (cid:48) ( M ( γ +2 ,δ +2) ) L ( γ +2 ,δ +2)( γ +3 ,δ +3) D ( γ +3 ,δ +3)( γ,δ ) + (cid:2) τ + 2 τ (cid:48) σ − τ σ (cid:48) − σσ (cid:48)(cid:48) + 4 λ c − σ − σ (cid:48) (cid:3) ( M ( γ +2 ,δ +2) ) D ( γ +2 ,δ +2)( γ,δ ) (32) + 3 [ τ τ (cid:48) + 2 λ c − σ (cid:48) − ( τ + σ (cid:48) ) σ (cid:48)(cid:48) ] ( M ( γ +2 ,δ +2) ) R ( γ +2 ,δ +2)( γ +1 ,δ +1) D ( γ +1 ,δ +1)( γ,δ ) + [2 λ c − σ (cid:48)(cid:48) − τ (cid:48) ( σ (cid:48)(cid:48) − τ (cid:48) )] R ( γ +2 ,δ +2)( γ,δ ) ,B = 2 σ ( M ( γ +2 ,δ +2) ) D ( γ +2 ,δ +2)( γ,δ ) + 3 σ (cid:48) ( M ( γ +2 ,δ +2) ) R ( γ +2 ,δ +2)( γ +1 ,δ +1) D ( γ +1 ,δ +1)( γ,δ ) + σ (cid:48)(cid:48) R ( γ +2 ,δ +2)( γ,δ ) , (33) C = R ( γ +2 ,δ +2)( γ,δ ) . (34)Letting W ± = V ± Λ ± , the quadratic eigenvalue problem may be linearized, collecting both families of eigenvectors in thesame × block equation: (cid:18) A BI (cid:19) (cid:18) V − V + W − W + (cid:19) = (cid:18) CI (cid:19) (cid:18) V − V + W − W + (cid:19) (cid:18) Λ − Λ + (cid:19) . (35)Any × block equation with upper-triangular and banded structure may be permuted to an upper-triangular and bandedequation involving × blocks. This permutation, known as the perfect shuffle in reference to a deck of cards, can bedescribed as taking the odd columns before the even columns of the n × n identity, P = I n, [1:2:2 n, n ] . Then: P (cid:18) A BI (cid:19) P (cid:62) (cid:124) (cid:123)(cid:122) (cid:125) =: A P (cid:18) V − V + W − W + (cid:19) P (cid:62) (cid:124) (cid:123)(cid:122) (cid:125) =: V = P (cid:18) CI (cid:19) P (cid:62) (cid:124) (cid:123)(cid:122) (cid:125) =: B P (cid:18) V − V + W − W + (cid:19) P (cid:62) (cid:124) (cid:123)(cid:122) (cid:125) = V P (cid:18) Λ − Λ + (cid:19) P (cid:62) (cid:124) (cid:123)(cid:122) (cid:125) =: Λ , is in fact an upper-triangular and banded generalized eigenvalue problem with × blocks: AV = BV Λ . (36)Eq. (36) is almost in the form of Eq. (25): were we to relate the former to the latter, we would enable the divide-and-conquerstrategy that already accelerates the classical orthogonal polynomial connection problem. Since V is almost upper-triangular,we solve the n generalized eigenvalue problems of size × on the main diagonal and use the solutions to produce a sequenceof Givens rotations [30] to upper-triangularize V . If: V = Q V R V , and if: B Q V = Q B R B , then it follows that Q (cid:62) B A Q V is upper-triangular and banded. That is, the problem: (cid:0) Q (cid:62) B A Q V (cid:1) R V = R B R V Λ , is an upper-triangular and banded generalized eigenvalue problem of size n × n . The divide-and-conquer approach can beused to approximate R V , which together with Q V , provides a structured form for the generalized eigenvectors. Calculating Q V costs O ( n ) flops, while B Q V , Q B , R B , and Q (cid:62) B A Q V cost O ( nb ) flops, where b is the upper bandwidth of both A and B . The classical orthogonal polynomials all have well-separated spectra. From Theorem 1.1 and considering the associatedJacobi polynomial spectra in Eq. (16), if γ is an integer, there are eigenspaces with algebraic multiplicity greater than one. Forexample, for first associated Legendre polynomials, ( α, β, c ) = (0 , , and: Λ = λ − λ +0 λ − λ +1 λ − λ +2 . . . = . . . . We have added the partitioning lines to help illustrate the problem that arises in the divide-and-conquer approach wheneigenvalues are semi-simple: the component-wise formula for V , in Eq. (28) is invalid at any indices (cid:96) and n such that (Λ ) (cid:96),(cid:96) = (Λ ) n,n . Instead, given a tolerance (cid:15) > , we define ˆ V , by: ( ˆ V , ) (cid:96),n = ( − XY (cid:62) ) (cid:96),n (Λ ) (cid:96),(cid:96) − (Λ ) n,n if | (Λ ) (cid:96),(cid:96) − (Λ ) n,n | > (cid:15) max {| (Λ ) (cid:96),(cid:96) | , | (Λ ) n,n |} , , so that V , = ˆ V , + S , , where S , is a sparse matrix with O ( γ ) nonzero entries. We find the nonzero entries in S , bycomparing V , V , to “true” eigenvectors found by shifting-and-inverting in O ( n ) flops per eigenvector. If γ = α + β +2 c − so that λ − n ≡ λ + n , the linearization of the quadratic eigenvalue problem degenerates, a situation that occurs only if c = 1 along the line α + β = − . We provide a solution for this problem in § 3.2. n E x e c u t i o n T i m e ( s ) n log n ) 10 n - n o r m R e l a t i v e E rr o r n )( n ) Figure 1: Left: Computational timings for a fast approximate matrix-vector product with first associated Legendre–Legendreconnection coefficients. Right: Error growth given first associated Legendre coefficients c k = ( k + 1) − for k = 0 , . . . , n − .Figure 1 shows timings for the precomputation of the structured form of V and the execution of a matrix-vector productwith a set of first associated Legendre coefficients with modest decay in single and double precision. The timings scale with O ( n log n ) because our implementation [31] uses a hierarchical approximation of the Cauchy matrices; it may be improvedto O ( n log n ) by an implementation of the FMM. The figure also shows the low algebraic degree-dependent scaling of the -norm relative error.It has been observed that the accuracy of the divide-and-conquer approach depends on the conditioning of the eigenvectorsthemselves [22, 23]. For classical orthogonal polynomial connection problems within the same family, the -norm conditionnumber may be related to parameter differences: for ultraspherical–ultraspherical connection problems, the conditioning isrelated to an algebraic power of n on the order of the absolute difference of the ultraspherical parameters, say O ( n | λ − µ | ) . Thisgeneralizes to the Jacobi and Laguerre families as well. Hence, while it is theoretically possible to accelerate an out-of-family connection problem such as Hermite to Legendre, the conditioning of the problem may be a severe numerical impediment.Generally speaking, the ill-conditioning of a connection problem may be reduced (but not eliminated) by orthonormalization.For the associated classical – classical connection problem, the coupling of the positive and negative eigenspaces is a novelcause for concern. Heuristically, Eq. (4) promotes nonuniform oscillatory behaviour in the polynomial solutions on D =supp( µ ) . Eq. (12) promotes a similar oscillatory behaviour for the positive solutions p + n ( x ; c ) , but for solutions with a negativeeigenvalue, the polynomials p − n ( x ; c ) are permitted to develop an oscillation-free boundary layer. Such solutions are incon-sistent with any classical orthogonal polynomial family to which the associated polynomials are connected. This mismatchbetween the behaviours of p − n ( x ; c ) and q (cid:96) ( x ) may be the cause for a catastrophic growth in the conditioning of the coupledconnection problem of Eq. (36), eigenvalue separation notwithstanding. Figure 2 demonstrates that the normalized first as-sociated connection problems are all reasonably well-conditioned, but that the condition numbers of the coupled connectionproblem of Eq. (36) grows too rapidly to be useful in the Hermite, Laguerre, and high-parameter Jacobi connection problems.We believe the (block) triangular structure justifies the accuracy in the numerical estimates on the condition numbers even asthey may grow beyond /(cid:15) to astronomical figures. n ( V ) P ( x ; 1) P ( x ) P (1.23,4.56) ( x ; 1) P (1.23,4.56) ( x ) L (7.89) ( x ; 1) L (7.89) ( x ) H ( x ; 1) H ( x )( n )( n )(log n ) 10 n ( V ) P ( x ; 1) P ( x ) P (1.23,4.56) ( x ; 1) P (1.23,4.56) ( x ) L (7.89) ( x ; 1) L (7.89) ( x ) H ( x ; 1) H ( x ) Figure 2: Left: The degree-dependence of κ ( ˜ V ) for first associated Legendre, Jacobi, Laguerre, and Hermite connectionproblems. Right: The degree dependence of κ ( ˜ V ) for the same connection problems. In both plots, a tilde overtop thepolynomials indicates orthonormalization, related to the standard normalization via p n ( x ; c ) = (cid:112) h n + c ˜ p n ( x ; c ) , where h n isdefined in Table 1.Whereas the formulation of the associated classical – classical connection problem as a linearized quadratic eigenvalue problemformally solves our problem, the condition number of specific problems prohibits the method from delivering any reasonableresults. For these reasons, we describe starting points for two alternative methods in the next two subsections to provideavenues for future research. The divide-and-conquer strategy may be applied directly to either of Eqs. (30) and (31) without coupling the solutions. Bothsolutions of the quadratic eigenvalue problem are of the form: AV + BV Λ = CV Λ . We suppose that V has the factored form in Eq. (26). Block dividing all matrices conformably: (cid:18) A , A , A , (cid:19) (cid:18) V , V , V , V , (cid:19) + (cid:18) B , B , B , (cid:19) (cid:18) V , Λ V , V , Λ V , Λ (cid:19) = (cid:18) C , C , C , (cid:19) (cid:18) V , Λ V , V , Λ V , Λ (cid:19) . As the diagonal blocks have the same structure as the original problem, we assume they have been solved before seeking asolution to the off-diagonal block V , , that now satisfies: A , V , V , + A , V , + B , V , V , Λ + B , V , Λ = C , V , V , Λ + C , V , Λ . Since V − , C − , A , V , = Λ − V − , C − , B , V , Λ , the matrix equation can be cast in the following form: Λ V , − V , Λ − V − , C − , B , V , (Λ V , − V , Λ ) = V − , C − , (cid:0) C , V , Λ − A , V , − B , V , Λ (cid:1) , =: − XY (cid:62) . (37)As in the classical divide-and-conquer scheme, rank( XY (cid:62) ) depends only on the bandwidths of A , B , and C .So far, Eq. (37) does not utilize the quadratic aspect of the eigenvalue problem, and we would have obtained an apparently-similar equation for a two-parameter eigenvalue problem with known parameters, say Γ and Ω , had we started with Eq. (5).However, the quadratic character of the eigenvalue problem leads to a considerable simplification. Lemma 3.1.
For any two square matrices A and B and a conformable matrix Z : A Z − ZB = A ( AZ − ZB ) + ( AZ − ZB ) B. (38)Lemma 3.1 allows us to solve Eq. (37) by factoring Λ V , − V , Λ and using the following two step procedure. First, wewould solve for W in: (cid:0) Λ − V − , C − , B , V , (cid:1) W + W Λ = − XY (cid:62) . (39)Then, we would solve for V , in: Λ V , − V , Λ = W. (40)Although this scheme conserves the problem size and uncouples positive from negative eigenvectors, it is not our schemeof preference as Eq. (39) is a nonnormal matrix equation for W . Although it is likely that there exists a low-rank solution,none of the aforementioned algorithms provides computable bounds on the rank of W : the analogy to a Hadamard productwith a Cauchy matrix is lost when one of the terms in the matrix equation is not diagonal; the (optimal) choice of shifts infADI must either take into account any ill-conditioning in Eq. (39) that results from diagonalizing Λ − V − , C − , B , V , ormust be chosen according to its field of values; and, the quadrature-based approach based on the semi-infinite integral eitherrequires the fast approximate action of a nonnormal matrix exponential or an estimate of the departure from normality for thequadrature error. Theoretical properties on nonnormal matrix equations are sparse in the literature. Notably, Baker, Embree,and Sabino [32] discuss specific structured nonnormal matrix equations which are found to admit low-rank solutions despite(and even aided by) the nonnormality.Algorithmic considerations aside, it is important for invertibility to show that the spectra in the two terms in Eqs. (39) and (40)are disjoint. For Eq. (40), this follows from Theorem 1.1. For Eq. (39), the disjointedness of the spectra requires more attention.Since V , , B , , and C , are upper-triangular and Λ is diagonal: σ (Λ − V − , C − , B , V , ) = σ (Λ − C − , B , ) . (41)By Eq. (12), it follows that the spectrum of C − , B , is equal to the set of the first s eigenvalues of B . The three classicalfamilies behave differently.1. For Jacobi polynomials, σ = x − , σ (cid:48) = 2 x , and σ (cid:48)(cid:48) = 2 so that: σ D + 3 σ (cid:48) D + σ (cid:48)(cid:48) I = 2 (cid:16) σ D + τ ( , ) D + I (cid:17) , hence σ ( C − , B , ) = { λ ( , ) ν + 2 } s − ν =0 , where λ ( , ) n are the eigenvalues of the second kind Chebyshev differentialequation.2. For Laguerre polynomials, σ = − x , σ (cid:48) = − , and σ (cid:48)(cid:48) = 0 so that: σ D + 3 σ (cid:48) D + σ (cid:48)(cid:48) I = − x D − D , hence σ ( C − , B , ) = { } .3. Finally, for Hermite polynomials, σ = − and σ (cid:48) = σ (cid:48)(cid:48) = 0 so that: σ D + 3 σ (cid:48) D + σ (cid:48)(cid:48) I = − D , hence σ ( C − , B , ) = { } . For Laguerre and Hermite polynomials, it is clear that σ (Λ ) ∩ σ ( − Λ ) = ∅ .For Jacobi polynomials, a sufficient condition for σ (Λ − C − , B , ) ∩ σ ( − Λ ) = ∅ is that: λ ( α,β ) s + c + λ ( α,β ) s + c − > (cid:16) λ ( , ) s − + λ ( α,β ) c − + 1 (cid:17) . Canceling terms, this reduces to: (2 s + 1)( α + β + 2 c ) > . Since the order of association c is nontrivially at least one and the Jacobi parameters must satisfy α, β > − , there is always aspectral gap , guaranteeing that the two step procedure defined by Eqs. (39) and (40) indeed has a solution. We consider the degree-preserving factored differential equation, Eq. (7) or equivalently Eq. (8) with c = 1 , where we expand { p n ( x ; 1) } ∞ n =0 in the corresponding classical orthogonal polynomial basis, { p (cid:96) ( x ) } ∞ (cid:96) =0 . We separate this case for two reasons:firstly, it is reasonable to consider this connection based on the diagonalization of the integral transform in Eq. (3); secondly,the form of the matrix equation satisfied by the connection coefficients is simplified considerably.The latent truth underlying Eq. (7) is the inhomogeneous second-order problem [9–11]: (cid:2) σ D + (2 σ (cid:48) − τ ) D + ( σ (cid:48)(cid:48) − τ (cid:48) ) (cid:3) p n ( x ; 1) = λ n +1 p n ( x ; 1) + ( σ (cid:48)(cid:48) − τ (cid:48) ) A D p n +1 ( x ) . (42)The particular form of the inhomogeneity is relatively easy to discretize: expanding { p n ( x ; 1) } ∞ n =0 in { p (cid:96) ( x ) } ∞ (cid:96) =0 and con-verting this representation to the (also classical) basis of {D p (cid:96) +1 ( x ) } ∞ (cid:96) =0 results in a “forced” upper-triangular and bandedgeneralized eigenvalue problem: AV = BV Λ + Γ , (43)where Γ is a diagonal forcing matrix.Dividing and conquering, the resulting matrix equation that must be solved is: V − , B − , A , V , V , − V , Λ = V − , B − , ( B , V , Λ − A , V , ) =: − XY (cid:62) . (44)Eq. (44) is almost the same as Eq. (27), though we can no longer say that V − , B − , A , V , = Λ .It is an interesting observation that the differential operator on the left-hand side of Eq. (42) is the formal adjoint of the classicalone. That is, the diagonalization of B − , A , in Eq. (44) is related to finding the polynomial solutions to the eigenvalueproblem: (cid:2) σ D + (2 σ (cid:48) − τ ) D + ( σ (cid:48)(cid:48) − τ (cid:48) ) (cid:3) q n ( x ) = D [ σq n ] − D [ τ q n ] = ω n q n ( x ) . (45) Lemma 3.2.
The eigenvalues ω n of Eq. (45) are given by: ω n = n + 12 [( n + 2) σ (cid:48)(cid:48) − τ (cid:48) ] . (46) Proof.
The proof follows by comparing coefficients of the monomial x n .The implied diagonalization is an issue for all three classical families. For Jacobi polynomials, it is akin to connecting { P ( − α, − β ) n ( x ) } ∞ n =0 to { P ( α,β ) (cid:96) ( x ) } ∞ (cid:96) =0 , which is ill-conditioned if max { α, β } (cid:29) . For Hermite and Laguerre polynomi-als, the differential equations are not related to classical orthogonal polynomial eigenproblems with different parameters.If A , R , = B , R , Ω , then letting Z , = R − , V , V , in Eq. (44) we find the diagonalized form: Ω Z , − Z , Λ = − R − , V , XY (cid:62) , and for a solution to this matrix equation to exist, we must show that σ (Ω ) ∩ σ (Λ ) = ∅ . A sufficient condition is that: max { ω , ω s − } < λ s +1 , ormax (cid:110) σ (cid:48)(cid:48) − τ (cid:48) , s s + 1) σ (cid:48)(cid:48) − τ (cid:48) ] (cid:111) < s + 12 [ sσ (cid:48)(cid:48) + 2 τ (cid:48) ] . though that gap may shrink in the limit as α + β → − + when c = 1 ! Canceling terms, this is true since τ (cid:48) > for all classical families.Using the fact that A , V , = B , V , Λ + Γ , an alternative to solving Eq. (44) is to solve the nonnormal matrix equation: Λ V , − V , Λ + V − , B − , Γ V , = − XY (cid:62) . (47) Remark 3.3.
The quadratic eigenvalue problem approach fails when the positive and negative eigenfunctions degenerate, asituation that occurs only if c = 1 and α + β = − . We note that for Jacobi polynomials, σ (cid:48)(cid:48) − τ (cid:48) = − α + β + 1) .Therefore, in view of the right-hand side of Eq. (42), the forcing term is precisely , which allows us to solve the connectionproblem between { p n ( x ; 1) } ∞ n =0 and in fact any classical orthogonal polynomial basis { q (cid:96) ( x ) } ∞ (cid:96) =0 via an upper-triangular andbanded generalized eigenvalue problem. We collect the known elegant results on the associated classical – classical orthogonal polynomial connection problem. Generalformulæ for the Jacobi and Laguerre connection coefficients are given by Lewanowicz [12] as sums of generalized hyperge-ometric functions. The formulæ are likely too complicated to be useful in practice, though we refer the interested readerto the results in case they feel differently. Earlier, Lewanowicz [33] finds formulæ for two special cases in the associatedJacobi–Jacobi connection coefficients, apart from two typographical errors corrected here. The first associated ultraspherical–ultraspherical case is due to Watson [34, §3.15.2]; see also Paszkowski [35]. The first associated Legendre–Legendre caseis discovered independently by Temme [36, Eq. (8.30)]. For the associated Hermite–Hermite problem, the formulæ are dueto Askey and Wimp [5]. For the generalized Laguerre polynomials, L ( α ) n ( x ; c ) , the special cases of α = in [12] and also α = − are related to the Hermite problem.The connection coefficients are conveniently expressed in terms of the gamma function [37, §5], defined by: Γ( z ) := (cid:90) ∞ x z − e − x d x, for (cid:60) z > , and its analytic continuation to z ∈ C \ {− N } by the recurrence Γ( z + 1) = z Γ( z ) . If one or more happen to be singular inthe formulæ below, then we take limiting values to preserve continuity. Lemma 4.1 (Lewanowicz [33]) . The connection coefficients between the associated Jacobi polynomials { P ( α, ) n ( x ; c ) } ∞ n =0 and the Jacobi polynomials { P ( α, ) (cid:96) ( x ) } ∞ (cid:96) =0 are given by: V (cid:96),n = (cid:40) u ( α,c ) n t ( α,c ) (cid:96),n h ( α,c ) (cid:96),n v ( α ) (cid:96) for (cid:96) ≤ n, otherwise. where: u ( α,c ) n = Γ( n + 1)Γ( n + )Γ(2 n + 2 α + 2 c + 2)Γ( α + c + )Γ( c + 1)Γ(2 n + 2)Γ( α + 2 c + )Γ( n + α + c + )Γ( n + c + 1) ,t ( α,c ) (cid:96),n = Γ( n − (cid:96) + − α )Γ( α + 2 c + )Γ( n − (cid:96) + 2 c )Γ( − α )Γ( n − (cid:96) + α + 2 c + )Γ( n − (cid:96) + 1) ,h ( α,c ) (cid:96),n = Γ( n + (cid:96) + 2)Γ( n + (cid:96) + α + 2 c + )Γ( n + (cid:96) + α + )Γ( n + (cid:96) + 2 α + 2 c + 2) ,v ( α ) (cid:96) = (2 (cid:96) + α + ) Γ( (cid:96) + α + )Γ( (cid:96) + ) . Since P ( α,β ) n ( x ; c ) = ( − n P ( β,α ) n ( − x ; c ) , the connection coefficients between { P ( ,β ) n ( x ; c ) } ∞ n =0 and { P ( ,β ) (cid:96) ( x ) } ∞ (cid:96) =0 aresimilar. If both source and target families satisfy p n ( − x ) = ( − n p n ( x ) , then this symmetry imparts in the connection coefficients achessboard pattern of zeros. Lemma 4.2 (Lewanowicz [33]) . The connection coefficients between the associated Jacobi polynomials { P ( α,α ) n ( x ; c ) } ∞ n =0 and the Jacobi polynomials { P ( α,α ) (cid:96) ( x ) } ∞ (cid:96) =0 are given by: V (cid:96),n = (cid:40) u ( α,c ) n t ( α,c ) (cid:96),n h ( α,c ) (cid:96),n v ( α ) (cid:96) for (cid:96) ≤ n, (cid:96) + n even , otherwise. where: u ( α,c ) n = Γ( n + c + α + 1)Γ( c + 1)Γ( c + 2 α + 1)Γ( c + α + 1)Γ( n + c + 1)Γ( c ) ,t ( α,c ) (cid:96),n = Γ( + α )Γ( n − (cid:96) +12 − α )Γ( n − (cid:96) + c )Γ( n − (cid:96) +12 + c + α )Γ( − α )Γ( n − (cid:96) + 1) ,h ( α,c ) (cid:96),n = Γ( n + (cid:96) +22 )Γ( n + (cid:96) +12 + c + α )Γ( n + (cid:96) +32 + α )Γ( n + (cid:96) +22 + c + 2 α ) ,v ( α ) (cid:96) = ( (cid:96) + α + ) Γ( α + 1)Γ( (cid:96) + 2 α + 1)Γ( (cid:96) + α + 1)Γ(2 α + 1) . Lemma 4.3 (Askey and Wimp [5]) . The connection coefficients between the associated Hermite polynomials { H n ( x ; c ) } ∞ n =0 and the Hermite polynomials { H (cid:96) ( x ) } ∞ (cid:96) =0 are given by: V (cid:96),n = ( − n − (cid:96) Γ( n − (cid:96) + c )Γ( c )Γ( n − (cid:96) + 1) Γ( n + (cid:96) +22 )Γ( (cid:96) + 1) for (cid:96) ≤ n, (cid:96) + n even , otherwise. The dual purpose of this section is to suggest that other classes of fast transforms may be developed for these special cases.The lemmata above demonstrate that some special associated connection problems have a diagonally scaled Hadamard productstructure between a Toeplitz and a Hankel matrix. This theoretically permits the fast factorization approach in [21]; however, incertain parameter régimes, the two vectors defining the Toeplitz and Hankel parts grow and/or decay so rapidly that a numericalimplementation would exhibit overflow and/or underflow, respectively. Explicit formulæ for the connection coefficients mayalso enable adaptations of the FMM, extending the approach in [19, 20], though similar problems of growth and decay arepresent in the analysis of the off-diagonal numerical rank of the subblocks.
We have alluded to the fact that some associated connection problems may be ill-conditioned. It is of importance, then, toestablish modest bounds on the condition number in at least one scenario.By Lemma 4.2, we see that the first associated Legendre–Legendre connection coefficients are given by: V (cid:96),n = (cid:96) + 1)( n − (cid:96) + 1)( n + (cid:96) + 2) for (cid:96) ≤ n, (cid:96) + n even , otherwise. The -norm condition number, κ ( V ) , is equal to the ratio of the largest to the smallest singular values. It would have been tooeasy to estimate upper and lower bounds, respectively, for the largest and smallest singular values by the formulæ in [38–41].However, due to the slow off-diagonal decay, strictly upper-triangular absolute row and column sums are unbounded as n → ∞ .In consequence, the best lower bounds for the smallest singular value are eventually for n sufficiently large.We turn to the theory of M -matrices [42–47]. We wish to show that V is an inverse M -matrix; that is, its spectrum is in theclosed right-half complex plane and V − has non-positive off-diagonal entries. Given that V is a triangular matrix, it is clearthat the spectra of V and V − are positive. Combining the non-positive off-diagonal property of V − with triangular backsubstitution, it would follow from elementary row operations that: | ( V − ) (cid:96),n | ≤ V (cid:96),n V (cid:96),(cid:96) V n,n = (2 (cid:96) + 2)(2 n + 2)2(2 n + 1)( n − (cid:96) + 1)( n + (cid:96) + 2) for (cid:96) ≤ n, (cid:96) + n even , otherwise. If this inequality were true, we would use Hölder’s inequality [30] for the condition number, κ ( V ) ≤ (cid:112) κ ( V ) κ ∞ ( V ) , andestimate squared-logarithmic growth from each of κ ( V ) and κ ∞ ( V ) , proving Theorem 5.1. Theorem 5.1.
For the associated Legendre–Legendre connection problem, κ ( V ) = O (log n ) .Proof. Discretizing Eq. (42) results in the “forced” upper-triangular and banded generalized eigenvalue problem in Eq. (43)where: A (cid:96),n = (cid:96) ( (cid:96) + 1)( (cid:96) + 2)2(2 (cid:96) + 1) for (cid:96) = n, − ( (cid:96) + 2) ( (cid:96) + 3)2(2 (cid:96) + 5) for (cid:96) = n − , otherwise, B (cid:96),n = (cid:96) + 22(2 (cid:96) + 1) for (cid:96) = n, − (cid:96) + 22(2 (cid:96) + 5) for (cid:96) = n − , otherwise, and the diagonal matrices have entries Λ n,n = ( n + 1)( n + 2) and Γ n,n = − ( n + 2) .It is important to also note that A = B Ω , where Ω is a diagonal matrix with Ω n,n = n ( n + 1) . Multiplying Eq. (43) by B − ,we find: B − AV − V Λ = Ω V − V Λ = B − Γ . Subsequent multiplication by V − from the left and the right results in: Λ V − − V − Ω = − V − B − Γ V − = ( − V Γ − BV ) − . Now, if − V Γ − BV is an M -matrix, its inverse is non-negative. Given that the difference in diagonal scalings Λ and Ω changessign on the first super-diagonal, it would follow from the component-wise formula: ( V − ) (cid:96),n = [( − V Γ − BV ) − ] (cid:96),n Λ (cid:96),(cid:96) − Ω n,n , (cid:96) (cid:54) = n + 1 , that V − is an M -matrix.To begin, it is easy to show − V Γ − BV has positive entries on the main diagonal. Next, Let (cid:96) ≥ and m ≥ . Then: ( − V Γ − BV ) (cid:96),(cid:96) +2 m = m (cid:88) k =1 (cid:20) (cid:96) + 1)(2 k + 1)(2 (cid:96) + 2 k + 2) − (cid:96) + 1)(2 k − (cid:96) + 2 k ) (cid:21) m − k + 1)(2 (cid:96) + 2 k + 2 m + 2)+ 2(2 (cid:96) + 1)(2 (cid:96) + 2)(2 m + 1)(2 (cid:96) + 2 m + 2) . This sum is non-positive if and only if: m (cid:88) k =1 (cid:20) k + 1)(2 (cid:96) + 2 k + 2) − k − (cid:96) + 2 k ) (cid:21) m − k + 1)(2 (cid:96) + 2 k + 2 m + 2)+ 1(2 (cid:96) + 2)(2 m + 1)(2 (cid:96) + 2 m + 2) ≤ . Or equivalently: m (cid:88) k =1 (cid:20) (cid:96) + 1(2 k − k + (cid:96) ) − (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) (2 m + 1)( (cid:96) + m + 1)(2 m − k + 1)( k + (cid:96) + m + 1) ≥ . Consider the infinite telescoping series: ∞ (cid:88) k =1 (cid:20) (cid:96) + 1(2 k − k + (cid:96) ) − (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) = 1 . We use this clever form of unity to restate the inequality that we must prove as: m (cid:88) k =1 (cid:20) (cid:96) + 1(2 k − k + (cid:96) ) − (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) (cid:20) (2 m + 1)( (cid:96) + m + 1)(2 m − k + 1)( k + (cid:96) + m + 1) − (cid:21) ≥ ∞ (cid:88) k = m +1 (cid:20) (cid:96) + 1(2 k − k + (cid:96) ) − (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) . As the right-hand side is also telescoping, we must show: m (cid:88) k =1 (cid:20) (cid:96) + 1(2 k − k + (cid:96) ) − (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) (cid:20) (2 m + 1)( (cid:96) + m + 1)(2 m − k + 1)( k + (cid:96) + m + 1) − (cid:21) ≥ (cid:96) + 1(2 m + 1)( (cid:96) + m + 1) . Since: (2 m + 1)( (cid:96) + m + 1)(2 m − k + 1)( k + (cid:96) + m + 1) − k (2 k + 2 (cid:96) + 1)(2 m − k + 1)( k + (cid:96) + m + 1) , we simplify (canceling (cid:96) + 1 from both sides): m (cid:88) k =1 (cid:20) k + 2 (cid:96) + 1(2 k − k + (cid:96) ) − k + 2 (cid:96) + 1(2 k + 1)( k + (cid:96) + 1) (cid:21) k (2 m − k + 1)( k + (cid:96) + m + 1) ≥ m + 1)( (cid:96) + m + 1) . Since: k + 2 (cid:96) + 1 k + (cid:96) > , and 2 k + 2 (cid:96) + 1 k + (cid:96) + 1 < , and: k m − k + 1 ≥ m − > m , the left-hand side is greater than or equal to: m (cid:88) k =1 (cid:18) k − − k + 1 (cid:19) m ( k + (cid:96) + m + 1) . Finally, since: k + (cid:96) + m + 1 ≥ (cid:96) + 2 m + 1 , we again use the telescoping series to find that the left-hand side is greater than or equal to: (cid:18) − m + 1 (cid:19) m ( (cid:96) + 2 m + 1) = 2(2 m + 1)( (cid:96) + 2 m + 1) . The proof follows since: (cid:96) + 2 m + 2 (cid:96) + 2 m + 1 > . Numerical evidence in Table 2 suggests that the smallest singular value σ n ( V ) tends to a constant as n → ∞ , so that κ ( V ) = O (log n ) . Table 2: The least singular value of the first associated Legendre – Legendre connection coefficients. n σ n ( V ) n σ n ( V )4 0 . . . . . . . . . . Given a function f ∈ L ( D, d µ ) , we consider its (weighted) Hilbert transform: H D { f } ( x ) = 1 π − (cid:90) D f ( t ) t − x d µ ( t ) , for x ∈ D, where the dashed integral is interpreted as a Cauchy principal value. Applications of the Hilbert transform arise as a conse-quence of it being the solution operator to certain Riemann–Hilbert problems [48]. The obvious algorithm to compute theHilbert transform is to use singular integral quadrature rule. These are generally useful for evaluation at a single point. Butsuch schemes cannot rapidly evaluate the weighted Hilbert transform of a degree- ( n − polynomial at n points in O ( n log n ) flops.A common alternative strategy [49, 50] to compute Hilbert transforms is to cleave the singularity as: H D { f } ( x ) = 1 π (cid:90) D f ( t ) − f ( x ) t − x d µ ( t ) + f ( x ) π − (cid:90) D d µ ( t ) t − x . Thus, only the Hilbert transform of the measure requires principal value treatment.From the cleaved representation, the weighted Hilbert transform of any finite orthogonal polynomial expansion: f ( x ) = n − (cid:88) k =0 c k p k ( x ) , is given by: H D { f } ( x ) = A (cid:82) D d µ ( t ) π n − (cid:88) k =0 c k +1 p k ( x ; 1) + (cid:32) n − (cid:88) k =0 c k p k ( x ) (cid:33) π − (cid:90) D d µ ( t ) t − x (48)Solving the connection problem between associated orthogonal polynomials p k ( x ; 1) and the original polynomials p k ( x ) , wecan work with a common basis. With fast synthesis with p k ( x ) , this process enables the rapid computation of the Hilberttransform on the same grid, provided we compute the Hilbert transform of the measure. x y f ( x ) = e x {sech[4sin( x )]} exp( x )[ 1,1] { f }( x ) 10 n + 110 | c n | Figure 3: Left: a function and its Hilbert transform, ω = 80 and n = 8192 . Right: the approximate Legendre coefficients of f ( x ) .Figure 3 shows a particular function and its Hilbert transform on the unit interval with the uniform measure, d µ = d x , bysampling f on a Chebyshev grid, analyzing it in a Chebyshev series, converting the Chebyshev series to a Legendre series, andusing Eq. (48).There are other strategies for the unit interval. With a uniform measure, Olver [51] uses the Joukowsky transform to mapthe unit interval to the unit circle in the complex plane, and identifies a set of special functions that incorporate the inherent discontinuity in the uniform measure mapped to the unit circle. With a non-negatively weighted measure, d µ = w ( x ) d x ,Hasegawa and Torii [52] expand f ( x ) in Chebyshev polynomials of the first kind: f ( x ) = n − (cid:88) k =0 c T k T k ( x ) , and, using the formula due to Elliott [53, Appendix 1]: T k +1 ( x ) − T k +1 ( t ) x − t = U k ( t ) + 2 k (cid:88) j =1 U k − j ( t ) T j ( x ) , find: H D { f } ( x ) = 1 π n − (cid:88) k =0 c T k +1 (cid:90) − U k ( t ) + 2 k (cid:88) j =1 U k − j ( t ) T j ( x ) w ( t ) d t + (cid:32) n − (cid:88) k =0 c T k T k ( x ) (cid:33) π − (cid:90) − w ( t ) d tt − x . If: w j = (cid:90) − U j ( t ) w ( t ) d t, then by reversing the order of summation: H D { f } ( x ) = 1 π n − (cid:88) k =0 c T k +1 w k + 2 n − (cid:88) j =1 T j ( x ) n − (cid:88) k = j c T k +1 w k − j + (cid:32) n − (cid:88) k =0 c T k T k ( x ) (cid:33) π − (cid:90) − w ( t ) d tt − x . The discrete convolutions can be cast as an upper-triangular Toeplitz matrix-vector product, which can be applied in O ( n log n ) operations via the fast Fourier transform. We have developed fast approximate solutions to the associated classical – classical orthogonal polynomial connection problembased on the differential Eqs. (7) and (12). We have described when we anticipate these solutions to be successful and whenthe ill-conditioning of the problem warrants the development of alternative approaches. Promising alternatives require the fastapproximate solution of structured nonnormal matrix equations, a challenging area of active research.We have not fully explored the use of the differential Eqs. (5), (6), (8), and (9) for fast transforms. Of these, Eqs. (5) and (6)are essentially the same and Eq. (9) seems to be the least likely candidate for success: the polynomial variable coefficientsof the two factors in Eq. (9) have degrees at most c + 1 and c , respectively, the degree-expanding nature of which wouldcreate nontrivial lower bands in the discretizations. The lower bands can be dealt with by applying a sufficiently high-orderdifferential operator to ensure each factor is degree-preserving. However, the bandwidths of these differential discretizationsdepend on the order c of association, tying it to the complexity of any algorithm whose complexity depends on the bandwidth.There are other non-classical connection problems that may be accelerated by identifying similar structural relationships.Semi-classical orthogonal polynomials are those polynomials orthogonal with respect to a weight function that satisfies a first-order linear homogeneous differential equation with a rational coefficient. It has been shown [54] that the polynomials satisfya second-order linear homogeneous differential equation with all three coefficients variable in x and n . Rational measuremodifications do not necessarily satisfy a differential equation (unless they are modifying classical measures) but the identitiesthey satisfy [55] share enough properties to enable a structured solution of the connection problem. References [1] J. Favard. Sur les polynômes de Tchebicheff.
C. R. Acad. Sci. Paris , 200:2052–2053, 1935.[2] A. Erdélyi et al., editors.
Higher Transcendental Functions , volume 2. McGraw-Hill, New York, 1953.[3] S. Bochner. Über Sturm–Liouvillesche Polynomsysteme.
Math. Z. , 29:730–736, 1929. [4] H. L. Krall. Certain differential equations for Tchebycheff polynomials. Duke Math. J. , 4:705–718, 1938.[5] R. Askey and J. Wimp. Associated Laguerre and Hermite polynomials.
Proc. Roy. Soc. Edinburgh , 96:15–37, 1984.[6] J. Wimp. Explicit formulas for the associated Jacobi polynomials and some applications.
Can. J. Math. , 39:983–1000,1987.[7] E. Laguerre. Sur la réduction en fractions continues d’une fraction qui satisfait à une équation différentielle linéaire dupremier order dont les coefficients sont rationnels.
J. de Math. , 1:135–165, 1885.[8] W. Hahn. On differential equations for orthogonal polynomials.
Funkcialaj Ekvacioj , 21:1–9, 1978.[9] A. Zarzo, A. Ronveaux, and E. Godoy. Fourth-order differential equation satisfied by the associated of any order of allclassical orthogonal polynomials. A study of their distribution of zeros.
J. Comp. Appl. Math. , 49:349–359, 1993.[10] M. Foupouagnigni, W. Koepf, and A. Ronveaux. Factorization of fourth-order differential equations for perturbed classi-cal orthogonal polynomials.
J. Comp. Appl. Math. , 162:299–326, 2004.[11] A. Ronveaux. Fourth-order differential equation for numerator polynomials.
J. Phys. A: Math. Gen. , 21:L749–L753,1988.[12] S. Lewanowicz. Results on the associated classical orthogonal polynomials.
J. Comp. Appl. Math. , 65:215–231, 1995.[13] J. Keiner.
Fast Polynomial Transforms . Logos Verlag, Berlin, 2011.[14] S. A. Orszag. Fast eigenfunction transforms. In
Science and Computers , pages 13–30. Academic Press, New York, 1986.[15] A. Mori, R. Suda, and M. Sugihara. An improvement on Orszag’s fast algorithm for Legendre polynomial transform.
Trans. Info. Process. Soc. Japan , 40:3612–3615, 1999.[16] N. Hale and A. Townsend. A fast, simple, and stable Chebyshev–Legendre transform using an asymptotic formula.
SIAMJ. Sci. Comput. , 36:A148–A167, 2014.[17] R. M. Slevinsky. On the use of Hahn’s asymptotic formula and stabilized recurrence for a fast, simple, and stableChebyshev–Jacobi transform.
IMA J. Numer. Anal. , 38:102–124, 2018.[18] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.
J. Comp. Phys. , 73:325–348, 1987.[19] B. K. Alpert and V. Rokhlin. A fast algorithm for the evaluation of Legendre expansions.
SIAM J. Sci. Stat. Comput. ,12:158–179, 1991.[20] J. Keiner. Computing with expansions in Gegenbauer polynomials.
SIAM J. Sci. Comput. , 31:2151–2171, 2009.[21] A. Townsend, M. Webb, and S. Olver. Fast polynomial transforms based on Toeplitz and Hankel matrices.
Math. Comp. ,87:1913–1934, 2018.[22] J. Keiner. Gegenbauer polynomials and semiseparable matrices.
Elec. Trans. Numer. Anal. , 30:26–53, 2008.[23] S. Olver, R. M. Slevinsky, and A. Townsend. Fast algorithms using orthogonal polynomials.
Acta Numerica , 29:573–699,2020.[24] L. Grasedyck. Singular value bounds for the Cauchy matrix and solutions of Sylvester equations. Technical Report 13,University of Kiel, 2001.[25] D. W. Peaceman and Jr. H. H. Rachford. The numerical solution of parabolic and elliptic differential equations.
J. SIAM ,3:28–41, 1955.[26] E. L. Wachspress.
Iterative Solution of Elliptic Systems: And Applications to the Neutron Diffusion Equations of ReactorPhysics . Prentice-Hall, 1966.[27] P. Benner, R.-C. Li, and N. Truhar. On the ADI method for Sylvester equations.
J. Comp. Appl. Math. , 233:1035–1045,2009.[28] D. I. Zolotarev. Application of elliptic functions to questions of functions deviating least and most from zero.
Zap. Imp.Akad. Nauk St. Petersburg , 30:1–59, 1877. [29] B. Beckermann and A. Townsend. Bounds on the singular values of matrices with displacement structure. SIAM Rev. ,61:319–344, 2019.[30] G. H. Golub and C. F. Van Loan.
Matrix Computations . The Johns Hopkins University Press, fourth edition, 2013.[31] R. M. Slevinsky. https://github.com/MikaelSlevinsky/FastTransforms . GitHub , 2018.[32] J. Baker, M. Embree, and J. Sabino. Fast singular value decay for Lyapunov solutions with nonnormal coefficients.
SIAMJ. Matrix Anal. Appl. , 36:656–668, 2015.[33] S. Lewanowicz. Results on the associated Jacobi and Gegenbauer orthogonal polynomials.
J. Comp. Appl. Math. ,49:137–143, 1993.[34] A. Erdélyi et al., editors.
Higher Transcendental Functions , volume 1. McGraw-Hill, New York, 1953.[35] S. Paszkowski. Polynômes associés aux polynômes orthogonaux classiques. Technical report, Publications ANO-136,Univ. Sci. Techn. Lille, 1984.[36] N. Temme.
Special Functions: An Introduction to the Classical Functions of Mathematical Physics . Wiley Interscience,1996.[37] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.
NIST Handbook of Mathematical Functions .Cambridge U. P., Cambridge, UK, 2010.[38] L. Qi. Some simple estimates for singular values of a matrix.
Linear Algebra Appl. , 56:105–119, 1984.[39] C. R. Johnson. A Gersgorin-type lower bound for the smallest singular value.
Linear Algebra Appl. , 112:1–7, 1989.[40] Y. Yi-Sheng and G. Dun-he. A note on a lower bound for the smallest singular value.
Linear Algebra Appl. , 253:25–38,1997.[41] C. R. Johnson and T. Szulc. Further lower bounds for the smallest singular value.
Linear Algebra Appl. , 272:169–179,1998.[42] R. A. Willoughby. The inverse M -matrix problem. Linear Algebra Appl. , 18:75–94, 1977.[43] R. J. Plemmons. M -matrix characterizations. I—nonsingular M -matrices. Linear Algebra Appl. , 18:175–188, 1977.[44] C. R. Johnson. Inverse M -matrices. Linear Algebra Appl. , 47:195–216, 1982.[45] I. N. Imam. Tridiagonal and upper triangular inverse M -matrices. Linear Algebra Appl. , 55:93–104, 1983.[46] M. Lewin. On inverse M -matrices. Linear Algebra Appl. , 118:83–94, 1989.[47] C. R. Johnson and R. L. Smith. Inverse M -matrices, II. Linear Algebra Appl. , 435:953–983, 2011.[48] T. Trogdon and S. Olver.
Riemann–Hilbert Problems, Their Numerical Solution and the Computation of NonlinearSpecial Functions . SIAM, 2016.[49] F. W. King.
Hilbert Transforms , volume 1. Cambridge University Press, 2009.[50] F. W. King.
Hilbert Transforms , volume 2. Cambridge University Press, 2009.[51] S. Olver. Computing the Hilbert transform and its inverse.
Math. Comp. , 80:1745–1767, 2011.[52] T. Hasegawa and T. Torii. Hilbert and Hadamard transforms by generalized Chebyshev expansion.
J. Comp. Appl. Math. ,51:71–83, 1994.[53] D. Elliott. Truncation errors in two Chebyshev series approximations.
Math. Comp. , 19:234–248, 1965.[54] A. P. Magnus. Painlevé-type differential equations for the recurrence coefficients of semi-classical orthogonal polynomi-als.
J. Comp. Appl. Math. , 57:215–237, 1995.[55] V. B. Uvarov. The connection between systems of polynomials that are orthogonal with respect to different distributionfunctions.
Zh. Vychisl. Mat. Mat. Fiz. , 9:1253–1262, 1969. A Jacobi Recurrence Relations
It will be helpful to set γ = α + β + 1 to simplify the following formulæ. Proposition A.1. D (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) = (cid:16) P ( α +1 ,β +1)0 P ( α +1 ,β +1)1 · · · (cid:17) D ( α +1 ,β +1)( α,β ) (49) where: D ( α +1 ,β +1)( α,β ) = 12 γ + 1 γ + 2 γ + 3 . . . . (50) Proposition A.2. x (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) = (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) M ( α,β ) (51) where: M ( α,β ) = β − αγ +1 2( α +1)( β +1)( γ +1)( γ +2)2 γ +1 β − α ( γ +1)( γ +3) 2( α +2)( β +2)( γ +3)( γ +4)4( γ +1)( γ +2)( γ +3) β − α ( γ +3)( γ +5) 2( α +3)( β +3)( γ +5)( γ +6)6( γ +2)( γ +4)( γ +5) β − α ( γ +5)( γ +7) . . .. . . . . . . (52) Proposition A.3. (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) = (cid:16) P ( α +1 ,β +1)0 P ( α +1 ,β +1)1 · · · (cid:17) R ( α +1 ,β +1)( α,β ) (53) where: R ( α +1 ,β +1)( α,β ) = ( α − β )( γ +1)( γ +1)( γ +3) − ( α +2)( β +2)( γ +3)( γ +4)( γ +1)( γ +2)( γ +2)( γ +3) ( α − β )( γ +2)( γ +3)( γ +5) − ( α +3)( β +3)( γ +5)( γ +6)( γ +2)( γ +3)( γ +4)( γ +5) ( α − β )( γ +3)( γ +5)( γ +7) − ( α +4)( β +4)( γ +7)( γ +8) . . . . . . . . . . (54) Proposition A.4. (1 − x ) (cid:16) P ( α +1 ,β +1)0 P ( α +1 ,β +1)1 · · · (cid:17) = (cid:16) P ( α,β )0 P ( α,β )1 · · · (cid:17) L ( α,β )( α +1 ,β +1) (55) where: L ( α,β )( α +1 ,β +1) = α +1)( β +1)( γ +1)( γ +2)4( α − β )( γ +1)( γ +3) 4( α +2)( β +2)( γ +3)( γ +4) − γ +2)( γ +3) 8( α − β )( γ +3)( γ +5) 4( α +3)( β +3)( γ +5)( γ +6) − γ +4)( γ +5) 12( α − β )( γ +5)( γ +7) . . . − γ +6)( γ +7) . . .. . . . (56)The four operators can be composed naturally. For example, the second derivative results in D ( α +2 ,β +2)( α,β ) = D ( α +2 ,β +2)( α +1 ,β +1) D ( α +1 ,β +1)( α,β ) . B Laguerre Recurrence Relations
Proposition B.1. D (cid:16) L ( α )0 L ( α )1 · · · (cid:17) = (cid:16) L ( α +1)0 L ( α +1)1 · · · (cid:17) D ( α +1)( α ) (57) where: D ( α +1)( α ) = − − . . . . (58) Proposition B.2. x (cid:16) L ( α )0 L ( α )1 · · · (cid:17) = (cid:16) L ( α )0 L ( α )1 · · · (cid:17) M ( α ) (59) where: M ( α ) = α + 1 − α − − α + 3 − α − − α + 5 − α − − α + 7 . . .. . . . . . . (60) Proposition B.3. (cid:16) L ( α )0 L ( α )1 · · · (cid:17) = (cid:16) L ( α +1)0 L ( α +1)1 · · · (cid:17) R ( α +1)( α ) (61) where: R ( α +1)( α ) = − − . . . . . . . (62) Proposition B.4. x (cid:16) L ( α +1)0 L ( α +1)1 · · · (cid:17) = (cid:16) L ( α )0 L ( α )1 · · · (cid:17) L ( α )( α +1) (63) where: L ( α )( α +1) = α + 1 − α + 2 − α + 3 − . . .. . . . (64)The four operators can be composed naturally. For example, the second derivative results in D ( α +2)( α ) = D ( α +2)( α +1) D ( α +1)( α ) . It isalso true that M ( α ) = L ( α )( α +1) R ( α +1)( α ))