A further look at time-and-band limiting for matrix orthogonal polynomials

TL;DR

This paper extends classical scalar results on time-and-band limiting to matrix-valued orthogonal polynomials, involving new integral and differential operators acting on matrix functions.

Contribution

It introduces a novel extension of scalar time-and-band limiting results to matrix orthogonal polynomials with associated integral and differential operators.

Findings

Extension of scalar time-band limiting results to matrix case

Development of new integral and differential operators for matrix functions

Foundation for further research in matrix-valued signal processing

Abstract

We extend to a situation involving matrix valued orthogonal polynomials a scalar result that originates in work of Claude Shannon and a ground-breaking series of papers by D. Slepian, H. Landau and H. Pollak at Bell Labs in the 1960's. While these papers feature integral and differential operators acting on scalar valued functions, we are dealing here with integral and differential operators acting on matrix valued functions.

Equations

$$S\tilde{D}=\tilde{D}S.$$

$$(A_w{)}_{w=0}^{\infty }⟼\sum _{w=0}^{\infty }{A}_w{Q}_w(x)$$

$$f⟼A_w={\int }_a^{b}f(x)W(x)Q_w^{\ast }(x)dx.$$

$${\chi }_{\Omega }\mathcal{F}f=g=\text{ known },{\chi }_{N}f=f.$$

$$Ef={\chi }_{\Omega }\mathcal{F}{\chi }_{N}f=g.$$

$$E^{\ast }E={\chi }_N{\mathcal{F}}^{-1}{\chi }_{\Omega }\mathcal{F}{\chi }_N\text{ and }{S}_2=E{E}^{\ast }={\chi }_{\Omega }\mathcal{F}{\chi }_N{\mathcal{F}}^{-1}{\chi }_{\Omega }.$$

$$(M{)}_{m,n}=(E^{\ast }E{)}_{m,n}={\int }_a^{\Omega }{Q}_m(x)W(x)Q{{}^{\ast }}_n(x)dx,0\le m,n\le N.$$

$$k(x,y)=\sum _{n=0}^N{Q}_n^{\ast }(x)Q_n(y).$$

$$w_{\alpha ,\beta }(x)=(1-x{)}^{\alpha }(1+x{)}^{\beta },$$

$$W(x)=W_{(\alpha,\beta)}=\frac{1}{2}\left(\begin{array}[]{cc}w_{\alpha,\beta}+w_{\beta,\alpha}&-w_{\alpha,\beta}+w_{\beta,\alpha}\\ -w_{\alpha,\beta}+w_{\beta,\alpha}&w_{\alpha,\beta}+w_{\beta,\alpha}\end{array}\right),\quad x\in[-1,1].$$

$$W_{p,m}(x)=(1-x^2{)}^{\frac{m}{2}-1}\left(\begin{array}{cc}p{x}^2+m-p& -mx\\ -mx& (m-p)x^2+p\end{array}\right),x\in [-1,1].$$

$$W_{\lambda}(x)=[x(2-x)]^{\lambda-3/2}\left(\begin{array}[]{cc}1&x-1\\ x-1&1\end{array}\right),\quad x\in[0,2],$$

$$P_n(x)=P_n^{(\alpha ,\beta )}(x)=\frac{1}{2}\left(\begin{array}{cc}{p}_n^{(\alpha ,\beta )}(x)+p_n^{(\beta ,\alpha )}(x)& -p_n^{(\alpha ,\beta )}(x)+p_n^{(\beta ,\alpha )}(x)\\ -p_n^{(\alpha ,\beta )}(x)+p_n^{(\beta ,\alpha )}(x)& p_n^{(\alpha ,\beta )}(x)+p_n^{(\beta ,\alpha )}(x)\end{array}\right),$$

$$p_n^{(\alpha ,\beta )}=\frac{(\alpha +1{)}_n}{n!}\text{mbox}{}_2{F}_1\left(\begin{array}{c}-n,n+\alpha +\beta \\ \alpha +1\end{array};\frac{1-x}{2}\right),$$

$$⟨P_n,P_n⟩=||P_n(x)|{|}^2={\int }_{-1}^1{P}_n(x)W(x)P_n(x{)}^{\ast }dx=h_n\mathrm{Id},$$

$$h_n=\frac{2^{\alpha +\beta +1}\Gamma (\alpha +n+1)\Gamma (\beta +n+1)}{(\alpha +\beta +2n+1)n!\Gamma (\alpha +\beta +n+1)}.$$

$$x{P}_n(x)=A_n{P}_{n+1}(x)+B_n{P}_n(x)+C_n{P}_{n-1}(x),n\ge 1,$$

$$A_n=$$

$$C_n=$$

$$T=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right).$$

$$(1-x^2)\frac{d}{dx}{P}_n(x)=-nx{P}_n(x)-\frac{n(\alpha -\beta )}{\alpha +\beta +2n}T{P}_n(x)+{\gamma }_{n-1}{P}_{n-1}(x),$$

$${\gamma }_{n-1}=\frac{2(n+\alpha )(n+\beta )}{\alpha +\beta +2n}.$$

$$\frac{{\kappa }_{n-1}}{{\kappa }_n{h}_{n-1}}\left(P_{n-1}^{\ast }(y)P_n(x)-P_n^{\ast }(y)P_{n-1}(x)\right)=(x-y)\sum _{k=0}^{n-1}\frac{P_k^{\ast }(y)P_k(x)}{h_k},$$

$${\kappa }_n=\frac{\Gamma (\alpha +\beta +2n+1)}{2^{n}n!\Gamma (\alpha +\beta +n+1)}.$$

$$\frac{{\kappa }_{n-1}}{{\kappa }_n}=\frac{2n(n+\alpha +\beta )}{(2n+\alpha +\beta )(2n+\alpha +\beta -1)}.$$

$$D=\frac{d^2}{d{x}^2}(1-x^2)+\frac{d}{dx}\left(-x(\alpha +\beta +2)+(\alpha -\beta )T\right),$$

$$D=\frac{d}{dx}\left(\frac{d}{dx}(1-x^2)W(x)\right)W(x{)}^{-1},$$

$$\frac{d}{dx}\left(\frac{d{P}_n}{dx}(x)(1-x^2)W(x)\right)W(x{)}^{-1}={\Lambda }_n{P}_n(x).$$

$$k(x,y)=\sum _{w=0}^N{Q}_w^{\ast }(x)Q_w(y).$$

$$(fS)(x)=\int_{-1}^{\Omega}f(y)W(y)\big{(}k(x,y)\big{)}^{*}dy.$$

$$\tilde{D}=\frac{d}{dx}\left((x-\Omega )\frac{d}{dx}p(x)\right)W(x{)}^{-1}+xA,$$

Full text

The work of the first author was partially supported by MTM2015-6588-C4-1-P (Ministerio de Economía y Competitividad), FQM-262, FQM-7276 (Junta de Andalucía), Feder Funds (European Union) and the program Campus de Excelencia Internacional of the Ministerio de Educación, Cultura y Deporte.

The third and fourth authors were partially suported by SeCyT-UNC and by CONICET grant PIP 112-200801-01533.

The work of the fourth author was also supported by FONDECYT 3160646.

1. Introduction

Claude Shannon, [28], posed the question of how to best use the values of the Fourier transform $\mathcal{F}f(k)$ of $f$ for values of $k$ in the band $[-W,W]$ when $f(t)$ is a time limited signal.

A detailed account of how this led to the series of papers by three workers at Bell labs in the 1960’s: David Slepian, Henry Landau and Henry Pollak, see [35, 20, 21, 31, 33] is given in [5, 14, 15] and need not be repeated here. Readers unfamiliar with these Bell Lab papers may want to look at these last references.

With this motivation at hand, we can give an account of what we do in this paper: we start with a (matrix valued) version of a second order differential operator. Here Shannon would have started with the (scalar valued) second derivative.

We then build the analog of the “time-and-band limiting” integral operator, which we will denote by $S$. We then show that the same “lucky accident” that the workers at Bell labs found holds here too: we can exhibit a second order differential operator , denoted by $\widetilde{D}$ such that

[TABLE]

This has, as in the original case of Shannon, very important numerical consequences: it gives a reliable way to compute the eigenvectors of $S$, something that cannot be done otherwise.

The eigenfunctions of $S$ and $\widetilde{D}$ are the same, but using the differential operator instead of the integral one, we have a manageable numerical problem: while the integral operator has a spectrum with eigenvalues that are extremely close together, the differential one has a very spread out spectrum, resulting in a stable numerical computation.

Previous explorations of the commutativity property above in the matrix valued case can be seen in [14, 5], dealing with a full matrix and a narrow banded one, and in [15], dealing with an integral operator. In the present paper we extend the work started in the previous references.

For more details on computational issues see [2, 18, 22]. For applications involving (sometimes) vector-valued quantities on the sphere, see [17, 26, 29, 30].

2. Preliminaries

Let $W=W(x)$ be a weight matrix of size $R$ in the open interval $(a,b)$. By this we mean a complex $R\times R$-matrix valued integrable function $W$ on the interval $(a,b)$ such that $W(x)$ is positive definitive almost everywhere and with finite moments of all orders. Let $Q_{w}(x),w=0,1,2,...$, be a sequence of real valued matrix orthonormal polynomials with respect to the weight $W(x)$. Consider the following two Hilbert spaces: the space $L^{2}((a,b),W(t)dt)$, denoted here by $L^{2}(W)$, of all matrix valued measurable matrix valued functions $f(x)$, $x\in(a,b)$, satisfying $\int_{a}^{b}\operatorname{tr}\left(f(x)W(x)f^{*}(x)\right)dx<\infty$ and the space $\ell^{2}(M_{R},\mathbb{N}_{0})$ of all real valued $R\times R$ matrix sequences $(C_{w})_{w\in\mathbb{N}_{0}}$ such that $\sum_{w=0}^{\infty}\operatorname{tr}\left(C_{w}\,C_{w}^{*}\right)<\infty$.

The map $\mathcal{F}:\ell^{2}(M_{R},\mathbb{N}_{0})\longrightarrow L^{2}(W)$ given by

[TABLE]

is an isometry. If the polynomials are dense in $L^{2}(W)$, this map is unitary with the inverse $\mathcal{F}^{-1}:L^{2}(W)\longrightarrow\ell^{2}(M_{R},\mathbb{N}_{0})$ given by

[TABLE]

We denote our map by $\mathcal{F}$ to remind ourselves of the usual Fourier transform. Here $\mathbb{N}_{0}$ takes up the role of “physical space” and the interval $(a,b)$ the role of “frequency space”. This is, clearly, a noncommutative extension of the problem raised by C. Shannon since he was concerned with scalar valued functions and we are dealing with matrix valued ones.

The time limiting operator, at level $N$, acts on $\ell^{2}(M_{R},\mathbb{N}_{0})$ by simply setting equal to zero all the components with index larger than $N$. We denote it by $\chi_{N}$. The band limiting operator, at level $\Omega$, acts on $L^{2}(W)$ by multiplication by the characteristic function of the interval $(a,\Omega)$, $\Omega\leq b$. This operator will be denoted by $\chi_{\Omega}$. One could consider restricting the band to an arbitrary subinterval $(a_{1},b_{1})$. However, the algebraic properties exhibited here, see Section 4 and beyond, hold only with this restriction. A similar situation arises in the classical case going all the way back to Shannon.

Consider the problem of determining a function $f$, from the following data: $f$ has support on the finite set $\{0,\dots,N\}$ and its Fourier transform $\mathcal{F}f$ is known on the compact set $[a,\Omega]$. This can be formalized as follows

[TABLE]

We can combine the two equations into

[TABLE]

To analyze this problem we need to compute the singular vectors (and values) of the operator $E:\ell^{2}(M_{R},\mathbb{N}_{0})\longrightarrow L^{2}(W)$. These are given by the eigenvectors of the operators

[TABLE]

The operator $E^{*}E$, acting in $\ell^{2}(M_{R},\mathbb{N}_{0})$ is just a finite dimensional block-matrix $M$, and each block is given by

[TABLE]

The second operator $S=EE^{*}$ acts in $L^{2}((a,\Omega),W(t)dt)$ by means of the integral kernel

[TABLE]

Consider now the problem of finding the eigenfunctions of $E^{*}E$ and $EE^{*}$. For arbitrary $N$ and $\Omega$ there is no hope of doing this analytically, and one has to resort to numerical methods and this is not an easy problem. Of all the strategies one can dream of for solving this problem, none sounds so appealing as that of finding an operator with simple spectrum which would have the same eigenfunctions as the original operators. This is exactly what Slepian, Landau and Pollak did in the scalar case, when dealing with the real line and the actual Fourier transform. They discovered (the analog of) the following properties:

• •

For each $N$, $\Omega$ there exists a symmetric tridiagonal matrix $L$, with simple spectrum, commuting with $M$.

• •

For each $N$, $\Omega$ there exists a selfadjoint differential operator $D$, with simple spectrum, commuting with the integral operator $S=EE^{*}$.

To this day nobody has a simple explanation for these miracles, and this paper displays more instances where this holds. Indeed, there has been a systematic effort to see if the “bispectral property” first considered in [7], guarantees the commutativity of these two operators, a global and a local one. A few papers where this question has been taken up, include [8, 9, 10, 11, 12, 16, 24, 25].

We recall that while [5, 14] deal with the full matrix $E^{*}E$ alluded to above, here, as well as in [15], we deal with the integral operator $S=EE^{*}$ mentioned above.

3. An example of matrix valued orthogonal polynomials

For $\alpha,\beta>-1$, the scalar Jacobi weight is given by

[TABLE]

supported in the interval $[-1,1]$. In this paper we consider a Jacobi type weight matrix of dimension two

[TABLE]

Let us observe that a particular case of these weight matrices have been studied in [14, 15] and [5]. In fact, the weight matrix considered in [23, 14, 15] is

[TABLE]

Taking $\alpha=p-1$, $\beta=p+1$ in (2), we obtain a multiple of the previous weight for the special case $m=2p$.

On the other hand, taking $\beta=\alpha-1$ in (2) we obtain a linear translation $x\longrightarrow 1+x$ of the weight considered in [5],

[TABLE]

with $\lambda=\alpha+\displaystyle\frac{1}{2}$.

A sequence of matrix orthogonal polynomials with respect to the matrix valued inner product going with (2), is given by

[TABLE]

where $p_{n}^{(\alpha,\beta)}$ are the classical Jacobi polynomials

[TABLE]

which are orthogonal with respect to the weight $w_{\alpha,\beta}$ (see for instance [36, Chapter VI]).

The norm of the polynomials $P_{n}(x)$ is given by

[TABLE]

where $\mathrm{Id}$ is the matrix identity of size $2\times 2$ and

[TABLE]

The sequence of orthogonal polynomials $(P_{n})_{n}$ satisfies the three term recurrence relation

[TABLE]

where

[TABLE]

with

[TABLE]

The sequence of matrix orthogonal polynomials $(P_{n})_{n}$ satisfies the following differentiation formula

[TABLE]

where

[TABLE]

We also have the following Christoffel-Darboux formula for the sequence of orthogonal polynomials $(P_{n}(x))_{n}$, introduced for a general sequence of matrix orthogonal polynomials in [6].

[TABLE]

with

[TABLE]

Observe that $\kappa_{n}\textrm{Id}$ is the leading coefficient of the matrix polynomial $P_{n}(x)$ and we also have

[TABLE]

The matrix polynomial $P_{n}$, for each $n\geq 0$, is an eigenfunction of the second order differential operator

[TABLE]

with scalar eigenvalues $\Lambda_{n}=-n(n+\alpha+\beta+1)$. We have that the differential operator $D$ can be factorized as

[TABLE]

and therefore the sequence of matrix orthogonal polynomials $(P_{n}(x))_{n}$ satisfies

[TABLE]

4. Time and band limiting. Integral and differential operators

Given a sequence of matrix orthonormal polynomials $\{Q_{w}\}_{w\geq 0}$ with respect to the weight $W$, we fix a natural number $N$ and $\Omega\in(-1,1)$ and we consider the integral kernel

[TABLE]

It defines the integral operator $S$ acting on $L^{2}((-1,\Omega),W)$ “from the right-hand side”:

[TABLE]

The restriction to the interval $[1,\Omega]$ implements “band-limiting” while the restriction to the range $0,1,\ldots,N$ takes care of “time-limiting”. In the language of [12], where the authors were dealing with scalar valued functions defined on spheres, the first restriction gives a “spherical cap” while the second one amounts to truncating the expansion in spherical harmonics.

We search for a selfadjoint differential operator $\widetilde{D}$, defined in $[-1,\Omega]$, commuting with the integral operator $S$.

The main result of this section is the following

Theorem 4.1**.**

Let $p(x)=(1-x^{2})W(x)$. The symmetric second-order differential operator

[TABLE]

with $A=N(N+\alpha+\beta+2)\mathrm{Id}$, commutes with the integral operator $S$ given in (14).

Explicitly, we have

[TABLE]

where the coefficients $\widetilde{E}_{j},\ j=0,1,2$, are given by

[TABLE]

and $T$ is the permutation matrix given in (7).

Let us observe that the differential operator $\widetilde{D}$ is somehow related to the differential operator $D$, given in (11) and (12). Explicitly, we have

[TABLE]

We first show that the operator in (15) is indeed symmetric with respect to the (matrix valued) inner product defined by (2).

Proposition 4.2**.**

The differential operator $\widetilde{D}$ is a symmetric operator with respect to

[TABLE]

Proof.

For an appropriate dense set of functions $f,g$, we have

[TABLE]

Since the factor $(x-\Omega)p(x)$ vanishes at $x=-1$ and $x=\Omega$ we get

[TABLE]

Therefore

[TABLE]

Completing the proof of the proposition. ∎

Proof of Theorem 4.1.

From [15, Proposition 3.1] we have that a symmetric differential operator $\widetilde{D}$ commutes with an integral operator $S$ with kernel $k$ if and only if

[TABLE]

(Here we use $D_{x}$ to stress that $D$ acts on the variable $x$).

Let $D$ be the differential operator introduced in (11) and let $(Q_{n})_{n}$ be the orthonormal sequence of matrix valued polynomials given by

[TABLE]

where $h_{n}=\langle P_{n},P_{n}\rangle$ is explicitly given in (6).

These polynomials are eigenfunctions of the differential operator $D$, i.e. $Q_{n}D=\Lambda_{n}Q_{n}$, with $\Lambda_{n}=-n(n+\alpha+\beta+1)$ and they satisfy

[TABLE]

with

[TABLE]

Then

[TABLE]

and similarly

[TABLE]

Thus

[TABLE]

From the Cristoffel-Darboux formula, given in (10), we obtain that

[TABLE]

It is easy to verify that $\displaystyle\tilde{\gamma}_{n-1}\frac{\kappa_{n}\,h_{n-1}^{1/2}}{\kappa_{n-1}h_{n}^{1/2}}=\alpha+\beta+2n+1.$ Therefore, by exchanging the order of summation we get

[TABLE]

Now from (18), and using that $\Lambda_{n}+A-n=-n(n+\alpha+\beta+2)+N(N+\alpha+\beta+2)$ we get

[TABLE]

Hence the operators $\widetilde{D}$ and $S$ commute. ∎

Remark 4.3*.*

It is worth noticing that if one takes the new operator $\widetilde{D}^{(1)}=\widetilde{D}T$, one obtains an operator commuting with the integral operator $S$ in (14), linearly independent with the operator $\widetilde{D}$. The space of the differential operators of order two commuting with $S$ is generated by $T$, $\widetilde{D}T$ and $\widetilde{D}$.

5. A Chebyshev type example

As a particular example, if we put $\alpha=\displaystyle\frac{1}{2}$, $\beta=-\displaystyle\frac{1}{2}$ in (2), we have the Chebyshev type weight

[TABLE]

which was introduced in [1, Page 586] for a different purpose.

The monic family of polynomials orthogonal with respect to this weight matrix, is given explicitly in terms of the Chebyshev polynomials of the second kind $U_{n}(x)$,

[TABLE]

and it satisfies a first order differential equation as pointed out in [3, 4].

Moreover, the polynomials $P_{n}(x)=P^{(\alpha,\alpha-1)}_{n}(x)$ in (5) satisfy, for any $\alpha>0$, the first order differential equation

[TABLE]

as shown in [5].

One considers here the integral operator $S$ in (14) defined by the integral kernel

[TABLE]

Particularly, the norm of the polynomials $\widetilde{P}_{n}$ is given by $\displaystyle||\widetilde{P}_{n}||=\frac{\sqrt{\pi}}{2^{n}}$.

Hence, for this particular example, the commuting operator $\widetilde{D}$ is given by

[TABLE]

where $T$ is the permutation matrix in (7).

6. Conclusion and outlook

The main result derived in the previous sections is the existence of an explicit differential operator $\widetilde{D}$, which, as we proved, commutes with $S$.

If one compares this result with the one in the celebrated series of papers by D.Slepian, H.Landau and H.Pollak one may say that we are at the stage of their first paper. What is needed now is an argument to conclude that the eigenfunctions of $\widetilde{D}$ will automatically be eigenfunctions of the integral operator $S$.

In the series of papers mentioned above the simplicity of the spectrum of $\widetilde{D}$ follows from classical Sturm-Liouville theory and this guarantees that they have found a good way to compute the eigenvectors of $S$. In our situation, things could eventually be reduced to that case, but in principle $\widetilde{D}$, as well as $S$, have “matrix valued eigenvalues”, and the appropriate notion of “simple spectrum” requires careful handling.

To be quite explicit: there are two technical points that we are not addressing in this paper. The first one is the issue of “simple spectrum” in the matrix valued context. A good place to look for this sort of question is [27]. A second issue is that of making precise the correct selfadjoint extension ( i.e. boundary conditions) for our second order differential operator $\widetilde{D}$. This issue has been implicit starting with the first papers of Slepian, Landau and Pollak, and has been considered quite explicitly in the very recent paper [19].