Gaussian quadrature for matrix valued functions on the unit circle (Q1920176)

In this interesting paper, the author investigates analogues of the Gauss quadrature formula on the unit circle in the setting of \(p\times p\) complex valued matrices. Let \(\rho\) be a \(p\times p\) matrix valued distribution function on the unit circle, and introduce left and right inner products: \[ \langle P,Q\rangle_L= {1\over 2\pi} \int^{2\pi}_0 P(e^{i\theta})d\rho(\theta)Q(e^{i\theta})^*; \langle P,Q\rangle_R={1\over 2\pi} \int^{2\pi}_0 P(e^{i\theta})^*d\rho(\theta)Q(e^{i\theta}). \] Here, \(P,Q\in\mathbb{C}^{p\times p}[z]\), the set of \(p\times p\) matrix polynomials in a complex variable and * denotes the Hermitian conjugate. Corresponding to these left and right inner products, we can define left and right orthonormal matrix polynomials \(\{\phi^L_n(z)\}^\infty_{n=0}\), \(\{\phi^R_n(z)\}^\infty_{n=0}\). The author seeks a quadrature formula of the form \[ \langle F,G\rangle_L={1\over 2\pi} \int^{2\pi}_0 F(e^{i\theta})d\rho(\theta)G(e^{i\theta})^*\approx \sum^k_{j=1} F(z_j)\Lambda_jG(z_j)^*,\tag{\(*\)} \] where \(F\) and \(G\) are matrix valued functions defined on the unit circle. The \(\{\Lambda_j\}\) are \(p\times p\) matrices, and the \(\{z_j\}\) should lie on the unit circle for the quadrature to be generally applicable. How does one choose the points \(\{z_j\}\)? The zeros of the matrix orthogonal polynomials \(\{\phi^L_n(z)\}^\infty_{n=0}\), \(\{\phi^R_n(z)\}^\infty_{n=0}\) are unsuitable, as they lie in the open unit disk. The author introduces para-orthogonal polynomials \[ B_n(z,W_n)= \phi^L_n(z)+ W_n\widetilde\phi^R_n(z), \] where \(W_n\) is unitary, and \(\widetilde\phi^R_n(z)= z^n\phi^R_n(1/\overline z)^*\). The author proves that the zeros of \(B_n\) are eigenvalues of a unitary block lower Hessenberg matrix and so are on the unit circle. The \(k\) distinct zeros \(\{z_j\}\) are then chosen as the quadrature points, and the weights \(\{\Lambda_j\}\) are chosen to ensure that the quadrature formula \((*)\) is exact for matrix Laurent polynomials \(F\in\Lambda_{-s,t}\); \(G\in\Lambda_{-(n-1-t),(n-1-s)}\) for any \(0\leq s,t\leq n-1\). Here \(\Lambda_{-m,n}\) is the set of matrix Laurent polynomials whose powers of \(z\) range between \(-m\) and \(n\). The author provides explicit formulae for the weights, analogous to those in Gauss quadrature. Then several identities are proved, and a ``divide and conquer'' algorithm is provided to implement the quadrature rule. Finally, numerical examples are given.