# Continued Fractions and Orthogonal Polynomials in Several Variables

Tomas Sauer  Lehrstuhl für Mathematik mit Schwerpunkt Digitale
Bildverarbeitung & FORWISS, Universitẗ Passau, Innstr. 43, D-94053 Passau,
Germany  and  Yuan Xu  Department of Mathematics, University of Oregon,
Eugene, OR 97403–1222, USA<EMAIL_ADDRESS>
###### Abstract.

We extend the close interplay between continued fractions, orthogonal
polynomials, and Gaussian quadrature rules to several variables in a special
but natural setting which we characterize in terms of moment sequences. The
crucial condition for the characterization is the commutativity of the
multiplication operators on finite polynomial subspaces modulo an ideal.
Moreover, starting from the orthogonal polynomials or the three-term
recurrence, our method constructs a sequence of moment sequences that provides
an approximation of the maximal order for recovering the moment sequence that
defines the orthogonality.

###### Key words and phrases:

Orthogonal polynomials, cubature, continued fraction, several variables,
moment sequence

###### 2020 Mathematics Subject Classification:

41A21, 42C05, 65D15, 65D32

The second author was partially supported by Simons Foundation Grant #849676

##  1\. Introduction

It is well-known that continued fractions are closely related to orthogonal
polynomials and both, in turn, are related to the Gaussian quadrature rules.
The interrelation between the three is elegant and revealing. We are
interested in a possible extension of such relations in several variables. A
moment of reflection, however, shows several obvious obstacles to extending
from one variable to the multi-variable setting. As a start, it is not even
clear how to define a continued fraction expansion utilizing multivariate
orthogonal polynomials.

In the literature, there exist already several generalizations of multivariate
continued fractions, especially in terms of partitions and lattices; see, for
example, [9, 20]. Our approach will be fundamentally different and focus on
rational approximations of formal power series based on polynomials defined by
a three-term recurrence formula. This is in the spirit of Gauss’ original
approach when developing Gaussian quadrature formulas from [5]. Therefore, it
is natural that our extension will be based on the three-term recurrences in
$d$ variables, which then consist of $d$ relations, one for each coordinate
multiplication. They are known for multivariate orthogonal polynomials. The
_matrix_ coefficients of the three-term relations will be embedded into an
appropriate definition of continued fractions, which, however, can only give a
complete extension of the univariate case when there is an appropriate
extension of the Gaussian quadrature rules to several variables. The latter
relies on the existence of the maximal number of _common_ zeros of orthogonal
polynomials of the same degree, which may or may not exist. Thus, the ideal
generated by all such orthogonal polynomials of the same total degree must
play a fundamental role, which translates to special structures of the matrix
coefficients and to the moment sequences that define the orthogonality. This
approach is fundamentally different from the ones found in [9, 20] and
therefore yields yet another generalization of continued fractions to several
variables, a generalization that is natural in our context.

Our main result shows that there is a natural and complete extension to
several variables whose existence can be characterized in terms of the moment
sequence.

In one direction, we start with a moment functional ${\mathcal{L}}$, defined
via a moment sequence $\mu$ in ${\mathbb{R}}^{d}$, and assume that the moment
functional ${\mathcal{L}}$ is _definite_ so that orthogonal polynomials
defined via ${\mathcal{L}}$ exist. These polynomials then satisfy three-term
relations which are characteristic properties of orthogonality by Favard’s
theorem. The three-term relations define the Jacobi matrices, which are
representations of multiplication operators by coordinates. Our main
assumption is that the truncated Jacobi matrices commute, which leads to the
existence of a maximal number of common zeros of orthogonal polynomials and
then to Gaussian cubature rules. If ${\mathcal{L}}$ is positive definite, this
is how the Gaussian cubature rules are characterized in [4]. In the above
narrative. however, we do not require anymore that ${\mathcal{L}}$ is positive
definite. This is in accordance with the well-known and classical theory of
univariate continued fractions and the fast convergence of the associated
convergents, see [16], since the existence of a so-called _associated
continued fraction_ is equivalent to the definiteness. Moreover, the
positivity of ${\mathcal{L}}$ can even be expressed in terms of the signs of
the recurrence coefficients in the univariate case. In several variables, it
is well-known that the issue of positivity is significantly more complex,
which is essentially due to the fact that positive polynomials are no longer
sums of squares (cf. [19]). In the framework that we present here, many
results have to be reconsidered and modified. For example, without positivity,
the common zeros may not be real or simple anymore. As a result, it is
necessary to recast many of the results on the structure of OPs in several
variables for a functional ${\mathcal{L}}$ that is only definite but not
necessarily positive definite, which we shall carry out in some detail and
provide proof whenever necessary. In return, this framework enables us to give
a definition of continued fractions that can be regarded as an extension from
one to several variables.

In the other direction, we start with a polynomial sequence that satisfies
three-term relations and assume that their truncated multiplication operators
commute. By Favard’s theorem, these polynomials are orthogonal with respect to
a linear functional ${\mathcal{L}}$ which can be regarded as a moment
functional defined by $\mu_{\alpha}={\mathcal{L}}(x^{\alpha})$. Using
continued fractions, we can then define a sequence of moment sequences
$\mu^{n}$ and show that $\mu^{n}\to\mu$ with an optimal rate in the sense that
a maximal number of terms in $\mu^{n}$ coincides with those in $\mu$; in other
words, we can efficiently recover the moment sequence that determines the
orthogonality. The main requirement for this setup lies in the commutativity
of the truncated Jacobi matrices, which can in turn be reduced to a
restriction on the moments.

The paper is organized as follows. In the next section, we recall the
classical results of one variable, which serves as a roadmap for our extension
to several variables. In the third section, we provide a detailed account of
orthogonal polynomials in several variables and Gaussian cubature rules, where
the results are presented and proved when necessary, without assuming the
positive definiteness. In the fourth section, the _multivariate_ continued
fraction expansion is defined and used to recover the moment sequence and we
provide a characterization of the moments for which such an extension is
possible.

##  2\. Preliminary

We recall results in one variable in this section. To that end, let
${\mathcal{L}}$ be a linear moment functional, defined at least on the
polynomials, and let $\mu_{n}={\mathcal{L}}(x^{n})$ , $n=0,1,2,\ldots$, denote
the moments of ${\mathcal{L}}$. In particular, if $w$ is a (not necessarily
nonnegative) weight function on ${\mathbb{R}}$, then

${\mathcal{L}}(f)=\int_{{\mathbb{R}}}f(x)w(x)\mathrm{d}x$

is such a moment functional. We assume that all $\mu_{n}$ are finite and call
${\mathcal{L}}$, or $w$, definite if the Hankel matrices

$M_{n}:=\left[\mu_{k+j}:\begin{matrix}0\leq k\leq n\\\ 0\leq j\leq
n\end{matrix}\right]$

are nonsingular for all $n=0,1,\ldots$, and positive definite if, in addition,
the matrices are positive definite. If $w\neq 0$ is nonnegative, then the
Hankel matrices are positive definite for all $n$, so that ${\mathcal{L}}(f)$
is positive definite. If $w$ changes sign on its domain, however, the
determinant of $M_{n}$ could be zero for some $n$, in which case orthogonal
polynomials could degenerate in their degrees; see, for example, [11, 21]. Our
definition of definiteness avoids these degenerate cases.

A real polynomial $p_{n}$ of degree $n$ is called orthogonal with respect to
${\mathcal{L}}$ if ${\mathcal{L}}(p_{n}q)=0$ for all $q$ of degree at most
$n-1$. It is known that if ${\mathcal{L}}$ is definite, then the orthogonal
polynomials $\\{p_{n}\\}_{n\geq 0}$ exist and they satisfy a three-term
relation

(2.1) $p_{n+1}(x)=\left(a_{n}x+b_{n}\right)p_{n}(x)+c_{n}\,p_{n-1}(x),\qquad
a_{n}c_{n}\neq 0,\qquad n\geq 0.$

Moreover, this relation characterizes the orthogonality according to Favard’s
theorem. For ${\mathcal{L}}$ definite, the polynomial $p_{n}$ has $n$ zeros
counting multiplicity. If ${\mathcal{L}}$ is positive definite, the zeros are
all real and simple. If ${\mathcal{L}}$ is definite but not positive definite,
the zeros can be complex and appear in conjugate pairs; that is, if $z$ is a
zero of $p_{n}$ then so is $\bar{z}$.

When $w$ is positive definite, the zeros of $p_{n}$, $\\{x_{k}:1\leq k\leq
n\\}$, are nodes of the Gaussian quadrature rule,

(2.2)
$\int_{{\mathbb{R}}}f(x)w(x)\mathrm{d}x=\sum_{k=1}^{n}{\lambda}_{k}f(x_{k}),\qquad\deg
f\leq 2n-1,$

where ${\lambda}_{k}>0$. This remains true if ${\mathcal{L}}$ is definite but
not positive definite, and $p_{n}$ and $p_{n+1}$ do not have common zeros, for
which ${\lambda}_{k}\neq 0$ for $1\leq k\leq n$. The quadrature rule (2.2) has
been studied extensively. For signed weight function, it holds for all even
$n$ if $w$ is a weight function on $[-1,1]$ and
$\det\left[\mu_{2(k+j)-1}\right]_{k,j=1}^{n}\neq 0$, for example; see [6,
Section 3.1]. Another interesting case is the oscillatory weight functions.
Let $p_{n}(w)$ denote the orthogonal polynomial of degree $n$ with respect to
$w$. The oscillatory weight function $\widehat{w}$ is defined by

$\widehat{w}(t)=p_{m}(w;t)w(t),\qquad t\in{\mathbb{R}}.$

In this case, the moment $\mu_{k}(\widehat{w})=0$ if $k<m$ and
$\mu_{m}(\widehat{w})>0$, so that the $m+1$ Hankel matrix satisfies $\det
M_{m+1}>0$. In particular, the polynomial $p_{n}(\widehat{w})$, $n=m+1$, is
well defined. Moreover, it has all real and simple zeros in some cases; for
example, when $w(t)=(1-t^{2})^{{\lambda}-\frac{1}{2}}$ and $0<{\lambda}\leq 2$
(see, for example, [6, 13]).

We are mostly interested in the case when $w$ is definite and $p_{n}(w)$ has
distinct real simple zeros, which holds whenever ${\mathcal{L}}$ is positive
definite. Below is an example when ${\mathcal{L}}$ is definite but not
positive definite.

######  Example 2.1.

For ${\lambda}>-\frac{1}{2}$ and $\mu>0$, we define the weight function

$w_{{\lambda},\mu}(t)=t^{2\mu+1}(1-t)(1-t^{2})^{{\lambda}-\frac{1}{2}},\qquad-1<t<1.$

Then $w_{{\lambda},\mu}$ changes sign on $(-1,1)$. The orthogonal polynomial
$p_{n}(w_{{\lambda},\mu})$ has $n$ simple, real zeros in $[-1,1]$ and they are
nodes of the Gaussian quadrature of degree $2n-1$.

Indeed, let $C_{n}^{({\lambda},\mu)}$ be the generalized Gegenbauer
polynomials that are orthogonal with respect to the weight function
$|t|^{2\mu}(1-t^{2})^{{\lambda}-\frac{1}{2}}$ on $[-1,1]$. These polynomials
can be given explicitly in terms of the Jacobi polynomials
$P_{n}^{({\alpha},{\beta})}$ (cf. [4, Section 1.5.2]). In particular, up to a
multiple constant,

$\displaystyle
C_{2n}^{(\lambda,\mu)}(t)=b_{n}P_{n}^{(\lambda-\frac{1}{2},\mu-\frac{1}{2})}(2t^{2}-1).$

Then it can be easily verified (cf. [1]) that

$p_{2n}(w_{{\lambda},\mu};t)=C_{2n}^{({\lambda},\mu+\frac{1}{2})}(t)\quad\hbox{and}\quad
p_{2n+1}(w_{{\lambda},\mu};t)=(1+t)C_{2n}^{({\lambda}+1,\mu+\frac{1}{2})}(t)$

for $n=0,1,2,\ldots$. Moreover, it is easy to see that, if $\mu$ is an
integer, then

$\int_{-1}^{1}\left|p_{2n}(w_{{\lambda},\mu};t)\right|^{2}w_{{\lambda},\mu}(t)\mathrm{d}t=-\int_{-1}^{1}\left|C_{2n}^{({\lambda},\mu)}\right||t|^{2\mu+2}(1-t^{2})^{{\lambda}-\frac{1}{2}}\mathrm{d}t<0.$

The properties of the zeros of $p_{n}(w_{{\lambda},\mu})$ follow from those of
the zeros of the Jacobi polynomials. In particular, it follows that
$w_{{\lambda},\mu}$ admits Gaussian quadrature rules of degree $2n-1$ for all
$n$.

The recurrence relation relates orthogonal polynomials closely to continued
fractions, a fact used by Gauss to obtain the quadrature rules in [5]. Indeed,
defining the (infinite) continued fraction as

$r(x)=\frac{b_{1}|}{|a_{1}(x)}+\frac{b_{2}|}{|a_{2}(x)}+\frac{b_{3}|}{|a_{3}(x)}\cdots=\frac{b_{1}}{a_{1}(x)+\dfrac{b_{2}}{a_{2}(x)+\dfrac{b_{3}}{\ddots}}}.$

Its _convergents_ of order $n$, defined by

$r_{n}(x)=\frac{b_{1}|}{|a_{1}(x)}+\frac{b_{2}|}{|a_{2}(x)}+\dots+\frac{b_{n}|}{|a_{n}(x)}=\frac{p_{n}(x)}{q_{n}(x)},$

are rational functions, where the numerator $p_{n}$ and denominator $q_{n}$
are polynomials and they satisfy a recurrence relation

$\displaystyle p_{n}(x)$ $\displaystyle=a_{n}(x)p_{n-1}(x)+b_{n}p_{n-2},\qquad
p_{-1}=1,\,p_{0}=0,$ $\displaystyle q_{n}(x)$
$\displaystyle=a_{n}(x)p_{n-1}(x)+b_{n}p_{n-2},\qquad q_{-1}=0,\,q_{0}=1.$

If the polynomials $a_{n}$ are _affine polynomials_ , i.e.,
$a_{n}(x)=\alpha_{n}x+\beta_{n}$, then the above recurrence is the classical
_three-term recurrence relation_ for orthogonal polynomials. In fact, the
polynomials $q_{n}$ _are_ orthogonal polynomials provided that
$\alpha_{n}\,c_{n}<0$.

A classical result in continued fractions says that a formal Laurent series

$\mu(z)=\sum_{n=1}^{\infty}\mu_{n-1}z^{-n}$

has an _associated_ continued fraction expansion, i.e.,
$\mu(z)-r_{n}(z)=O\left(z^{2n}\right)$, if and only if

$\det\left[\mu_{j+k}:\begin{array}[]{c}j=0,\dots,n\\\
k=0,\dots,n\end{array}\right]\neq 0,\qquad n\in{\mathbb{N}}_{0}.$

If $\mu$ is the moment sequence associated with a definite linear functional,
then the denominator polynomials $q_{n}$ are the respective orthogonal
polynomials and the maximal exactness of Gaussian quadrature is eventually a
consequence of the fact that the Laurent series for the quadrature moments is
$p_{n}/q_{n}$. This is the central idea in Gauss’ original paper [5].

Continued fraction expansions are closely related to _Prony’s problem_ which
consists of recovering a multi-exponential function

$f(x)=\sum_{k=1}^{m}f_{k}\,\zeta_{k}^{x},\qquad\zeta_{k}\in{\mathbb{C}}\setminus\\{0\\},$

from integer samples $f(n)$, $n\in{\mathbb{N}}$. In fact, interpreting these
samples of $f$ as the moment sequence $\mu_{n}=f(n)$, the central connection
is that the infinite _Hankel operator_
$\begin{bmatrix}\mu_{j+k}:j,k\in{\mathbb{N}}_{0}\end{bmatrix}$ has finite rank
in this case and, by Kronecker’s Theorem, the associated Laurent series
$\mu(z)$ is a rational function, hence has a _finite_ continued fraction
expansion. Moreover, the denominator of the rational function is the so-called
_Prony polynomial_ whose zeros are exactly $\zeta_{1},\dots,\zeta_{m}$. The
monograph [16] provides more details on the connection between continued
fractions and Prony’s problem, and how the continued fraction expansion and
the construction of the associated recurrence coefficients can be used to
solve Prony’s problem even numerically.

##  3\. Orthogonal polynomials and cubature rules

We now consider OPs in several variables with respect to a linear functional
${\mathcal{L}}$ or, equivalently, a moment sequence $\mu$. As mentioned in the
introduction, we shall need results on the structures of these polynomials for
${\mathcal{L}}$ being definite, but not necessarily positive definite. This
requires some modifications of the existing theory that has been established
for the positive definite case; most of these modifications are
straightforward but by no means all of them. For example, the proof for the
existence of a maximal number of common zeros requires additional algebraic
tools. For the record and the reader’s convenience, we will point out the
modifications and provide proof whenever necessary.

Let ${\mathcal{L}}$ be a linear functional. For
${\alpha}\in{\mathbb{N}}_{0}^{d}$, let
$\mu_{\alpha}={\mathcal{L}}(x^{\alpha})$ denote the ${\alpha}$th moment of
${\mathcal{L}}$. We assume for simplicity that $\mu_{0}={\mathcal{L}}(1)\neq
0$ and that all moments $\mu_{\alpha}$ are real-valued. Let $\Pi_{n}^{d}$
denote the space of polynomials of total degree at most $n$ in $d$ variables
and let ${\mathcal{P}}_{n}^{d}$ denote the space of homogeneous polynomials of
degree $n$. Then

$r_{n}:=\dim\Pi_{n}^{d}=\binom{n+d}{d}\quad\hbox{and}\quad
r_{n}^{0}:=\dim{\mathcal{P}}_{n}^{d}=\binom{n+d-1}{d-1}.$

For $P\in\Pi_{n}^{d}$, we write $\widehat{P}\in{\mathbb{R}}^{r_{n}}$ for the
coefficient vector so that
$P(x)=\widehat{P}^{\mathsf{T}}{\boldsymbol{x}}^{n}=\sum_{\alpha}\widehat{P}_{\alpha}x^{\alpha}$,
where ${\boldsymbol{x}}^{n}:=\left[x^{\alpha}:|\alpha|\leq n\right]$. We
arrange $\widehat{P}$ in a graded fashion as

$\widehat{P}=\left[\begin{array}[]{c}\widehat{P}_{0}\\\ \vdots\\\
\widehat{P}_{n}\end{array}\right],\qquad\widehat{P}_{j}:=\left[\widehat{P}_{\alpha}:|\alpha|=j\right],\quad
j=0,\dots,n,$

and decompose ${\boldsymbol{x}}^{n}$ accordingly. The moment sequence $\mu$
defines a linear functional ${\mathcal{L}}$ on $\Pi$ by linearity as

${\mathcal{L}}(P)=\sum_{|\alpha|\leq n}p_{\alpha}\mu_{\alpha},\qquad
P(x)=\sum_{|\alpha|\leq n}p_{\alpha}x^{\alpha}\in\Pi_{n}^{d}.$

A polynomial $P\in\Pi_{n}^{d}$ is called orthogonal with respect to the linear
functional ${\mathcal{L}}$ if

${\mathcal{L}}(PQ)=0,\qquad\forall Q\in\Pi_{n-1}^{d}.$

In the usual setting of OPs in several variables, we assume that
${\mathcal{L}}$ is (square) positive definite, which requires
${\mathcal{L}}(p^{2})>0$ for all non-zero polynomials. This ensures the
existence of an orthonormal basis of polynomials. In this section, we state
the structure of OPs in the general setting, which are known in the positive
definite case, but require some modification if ${\mathcal{L}}$ is definite
but not necessarily positive definite.

###  3.1. Structure of orthogonal polynomials

For $n=0,1,2,\ldots$, let

$M_{n}=\left[\mu_{{\alpha}+{\beta}}:\begin{matrix}|{\alpha}|\leq n\\\
|{\beta}|\leq n\end{matrix}\right]$

be the moment matrix of size $r_{n}\times r_{n}$. Clearly,

${\mathcal{L}}(PQ)=\widehat{P}^{\mathsf{T}}M_{n}\widehat{Q},\qquad
P,Q\in\Pi_{n}^{d}.$

We assume that the matrices $M_{n}$, $n=0,1,2,\ldots$, are non-singular, which
ensures the existence of orthogonal polynomials. Indeed, the usual
construction of OPs from the moments (cf. [4, Section 3.3.2]) applies, which
gives a monic OP of the form

$P_{\alpha}(x)=x^{\alpha}+q_{\alpha},\qquad q_{\alpha}\in\Pi_{n-1}^{d},$

for each ${\alpha}\in{\mathbb{N}}_{0}^{d}$. Let ${\mathcal{V}}_{n}^{d}$ be the
space of OPs of degree $n$. Then $\dim{\mathcal{V}}_{n}^{d}=r_{n}^{0}$. Let
$\widehat{\mathbb{P}}_{n}$ be the column vector

$\widehat{\mathbb{P}}_{n}(x)=\left[P_{\alpha}:|{\alpha}|=n\right]=\widehat{\boldsymbol{P}}_{n}^{\mathsf{T}}{\boldsymbol{x}}^{n},\qquad\widehat{\boldsymbol{P}}_{n}:=\left[\hat{P}_{\alpha}:|{\alpha}|=n\right]\in{\mathbb{R}}^{r_{n}\times
r_{n}^{0}}.$

Then $\widehat{\mathbb{P}}_{n}$ consists of a basis of
${\mathcal{V}}_{n}^{d}$. Partitioning $M_{n}$ as

$M_{n}=\left[\begin{array}[]{cc}M_{n-1}&M_{n-1,n,}\\\
M_{n-1,n}^{\mathsf{T}}&M_{n,n}\end{array}\right],$

with

$M_{n-1,n}=\left[\mu_{\alpha+\beta}:\begin{array}[]{c}|\alpha|\leq{n-1}\\\
|\beta|=n\end{array}\right]\quad\hbox{and}\quad
M_{n,n}=\left[\mu_{\alpha+\beta}:\begin{array}[]{c}|\alpha|=n\\\
|\beta|=n\end{array}\right],$

it can easily be verified that

$\widehat{\boldsymbol{P}}_{n}=\left[\begin{array}[]{c}-M_{n-1}^{-1}M_{n-1,n}\\\
I_{r_{n}^{0}\times r_{n}^{0}}\end{array}\right],$

and that the matrix

${\mathcal{L}}(\widehat{\mathbb{P}}_{n}\widehat{\mathbb{P}}_{n}^{\mathsf{T}})=\widehat{\boldsymbol{P}}_{n}^{\mathsf{T}}M_{n}\widehat{\boldsymbol{P}}_{n}=M_{n,n}-M_{n-1,n}^{\mathsf{T}}M_{n-1}^{-{\mathsf{T}}}M_{n-1,n}$

is the symmetric and non-singular Schur complement of $M_{n,n}$ in $M_{n}$.
This matrix satisfies a decomposition

${\mathcal{L}}(\widehat{\mathbb{P}}_{n}\widehat{\mathbb{P}}_{n}^{\mathsf{T}})=Q^{\mathsf{T}}\Lambda
Q,\qquad\Lambda=|\Lambda|^{\frac{1}{2}}S_{n}|\Lambda|^{\frac{1}{2}},$

where $Q$ is an orthogonal matrix and
$\Lambda=\mathrm{diag}\\{{\lambda}_{1},\ldots,{\lambda}_{r_{n}^{0}})$ is a
non-singular diagonal matrix, which is further decomposed by using
$|\Lambda|^{\frac{1}{2}}=\mathrm{diag}\\{|{\lambda}_{1}|^{\frac{1}{2}},\ldots,|{\lambda}_{r_{n}^{d}}|^{\frac{1}{2}})$
and a _signature matrix_ $S_{n}\in{\mathbb{R}}^{r_{n}^{0}\times r_{n}^{0}}$,
the latter being a diagonal matrix with the diagonal elements being either $1$
or $-1$. We now define

${\mathbb{P}}_{n}=|\Lambda|^{-\frac{1}{2}}Q\widehat{\mathbb{P}}_{n},\qquad
n=0,1,2,\ldots.$

Then ${\mathbb{P}}_{n}$ consists of a basis of orthogonal polynomials of
degree $n$, so that

${\mathcal{L}}({\mathbb{P}}_{n}{\mathbb{P}}_{m}^{\mathsf{T}})=\delta_{n,m}S_{n},\qquad
n\neq m,$

which we call a _sign-orthonormal_ basis. If ${\mathcal{L}}$ is positive
definite, then $S_{n}$ is the identity matrix, and the sign-orthonormal basis
becomes an orthonormal basis.

###  3.2. Recurrences

Orthogonal polynomials are described by a three-term recurrence relation as
follows.

######  Proposition 3.1.

Let ${\mathbb{P}}_{n}$ be a sign-orthonormal basis of ${\mathcal{V}}_{n}^{d}$,
$n=0,1,2,\ldots$ and ${\mathbb{P}}_{-1}:=0$. Then there exist unique matrices
$A_{n,i}:r_{n}^{0}\times r^{0}_{n+1}$ and $B_{n,i}:r_{n}^{0}\times r_{n}^{0}$,
such that

(3.1)
$x_{i}S_{n}{\mathbb{P}}_{n}(x)=A_{n,i}{\mathbb{P}}_{n+1}(x)+B_{n,i}{\mathbb{P}}_{n}(x)+A_{n-1,i}^{\mathsf{T}}{\mathbb{P}}_{n-1}(x),\qquad
1\leq i\leq d,$

where the matrix $B_{n}$ is symmetric and

(3.2)
$\mathrm{rank}\,A_{n,i}=r_{n}^{0}\quad\hbox{and}\quad\mathrm{rank}\,(A_{n,1}^{\mathsf{T}},\ldots,A_{n,d}^{\mathsf{T}})^{\mathsf{T}}=r_{n+1}^{0},$

Indeed, since the vectors ${\mathbb{P}}_{n}$, $n\geq 0$, are a basis of
$\Pi^{d}$, there must exist matrices $A_{n,i},B_{n,i},C_{n,i}$ of proper size
such that

$x_{i}S_{n}{\mathbb{P}}_{n}(x)=A_{n,i}{\mathbb{P}}_{n+1}(x)+B_{n,i}{\mathbb{P}}_{n}(x)+C_{n,i}{\mathbb{P}}_{n-1}(x),$

and using the orthogonality, it is easy to see that

(3.3)
$A_{n,i}S_{n+1}=S_{n}{\mathcal{L}}(x_{i}{\mathbb{P}}_{n}{\mathbb{P}}_{n+1}^{\mathsf{T}}),\qquad
B_{n,i}S_{n}=S_{n}{\mathcal{L}}(x_{i}{\mathbb{P}}_{n}{\mathbb{P}}_{n}^{\mathsf{T}}),$

which shows that $B_{n}$ is symmetric since the signature matrix $S_{n}$
satisfies $S_{n}^{\mathsf{T}}=S_{n}=S_{n}^{-1}$; moreover, using (3.3) with
$n$ replaced by $n-1$ in the second equation, we obtain

$C_{n,i}S_{n-1}=S_{n}{\mathcal{L}}(x_{i}{\mathbb{P}}_{n}{\mathbb{P}}_{n-1}^{\mathsf{T}})=S_{n}(S_{n}A_{n-1,i}^{\mathsf{T}}S_{n-1})=A_{n-1,i}^{\mathsf{T}}S_{n-1},$

which implies that $C_{n,i}=A_{n-1,i}^{\mathsf{T}}$.

The rank conditions can be established as in the positive definite case: for
example, since ${\mathbb{P}}_{n}$ is sign-orthonormal, its leading coefficient
matrix $\widehat{\boldsymbol{P}}_{n}$, defined by

${\mathbb{P}}_{n}(x)=\widehat{\boldsymbol{P}}_{n}^{\mathsf{T}}{\boldsymbol{x}}^{n}_{0}+\ldots,\qquad{\boldsymbol{x}}^{n}_{0}:=\left[x^{\alpha}:|{\alpha}|=n\right],$

is a non-singular matrix of size $r_{n}^{d}$. Comparing the leading
coefficients of (3.1) gives

$S_{n}\widehat{\boldsymbol{P}}_{n}L_{n,i}=A_{n,i}\widehat{\boldsymbol{P}}_{n+1},$

where $L_{n,i}$ is the shift matrix of dimension $r_{n}^{0}\times
r_{n+1}^{0}$, defined by
$L_{n,i}{\boldsymbol{x}}^{n+1}_{0}=x_{i}{\boldsymbol{x}}^{n}_{0}$; this matrix
has full rank, thus $A_{n,i}$ has rank $r_{n}^{0}$.

The coefficient matrices of the three-term relation for sign-orthonormal
polynomials satisfy a set of commutativity relations: for $1\leq i,j\leq d$
and $k\geq 0$,

(3.4)
$\displaystyle\begin{split}A_{n,i}S_{n+1}A_{n+1,j}&=A_{n,j}S_{n+1}A_{n+1,i},\\\
A_{n,i}S_{n+1}B_{n+1,j}+B_{n,i}S_{n}A_{n,j}&=B_{n,j}S_{n}A_{n,i}+A_{n,j}S_{n}B_{n+1,i},\\\
A_{n-1,i}^{\mathsf{T}}S_{n-1}A_{n-1,j}+B_{n,i}S_{n}B_{n,j}+&A_{n,i}S_{n+1}A_{n,j}^{\mathsf{T}}\\\
=A_{n-1,j}^{\mathsf{T}}S_{n-1}A_{n-1,i}&+B_{n,j}S_{n}B_{n,i}+A_{n,j}S_{n+1}A_{n,i}^{\mathsf{T}},\end{split}$

where $A_{-1,i}:=0$. This follows from computing
${\mathcal{L}}(x_{i}x_{j}{\mathbb{P}}_{n}{\mathbb{P}}_{n\pm 1}^{\mathsf{T}})$
via the three-term relation in two different ways. For example,

$\displaystyle
S_{n}{\mathcal{L}}(x_{i}x_{j}{\mathbb{P}}_{n}{\mathbb{P}}_{n+1}^{\mathsf{T}})S_{n+1}^{\mathsf{T}}$
$\displaystyle={\mathcal{L}}(x_{i}S_{n}{\mathbb{P}}_{n}(x_{j}S_{n+1}{\mathbb{P}}_{n+1})^{\mathsf{T}}$
$\displaystyle=A_{n,i}{\mathcal{L}}({\mathbb{P}}_{n+1}{\mathbb{P}}_{n+1}^{\mathsf{T}})A_{n+1,j})=A_{n,i}S_{n+1}A_{n+1,j},$

and the left-hand side is unchanged if the order of $i$ and $j$ is exchanged,
which shows the first identity; the others are proved in exactly the same way.

In fact, for _any_ set of basis functions
${\mathbb{F}}_{n}\in\left(\Pi_{n}^{d}\right)^{r_{n}^{0}}$ and the
multiplication rule

(3.5)
$x_{i}{\mathbb{F}}_{n}(x)=A_{n,i}{\mathbb{F}}_{n+1}(x)+B_{n,i}{\mathbb{F}}_{n}(x)+C_{n,i}{\mathbb{F}}_{n-1},\qquad
n\in{\mathbb{N}}_{0},\,1\leq i\leq d,$

the consistency conditions
$x_{i}x_{j}{\mathbb{F}}_{n-1}=x_{j}x_{i}{\mathbb{F}}_{n-1}$ implies the
identities

(3.6) $\displaystyle\begin{split}A_{n-1,i}A_{n,j}&=A_{n-1,j}A_{n,i},\\\
A_{n-1,i}B_{n,j}+B_{n-1,i}A_{n-1,j}&=A_{n-1,j}B_{n,i}+B_{n-1,j}A_{n-1,i},\\\
A_{n-1,i}C_{n,j}+B_{n-1,i}B_{n-1,j}+&C_{n-1,i}A_{n-2,j}\\\
=A_{n-1,j}C_{n,i}&+B_{n-1,j}B_{n-1,i}+C_{n-1,j}A_{n-2,i}\\\
B_{n-1,i}C_{n-1,j}+C_{n-1,i}B_{n-2,j},&=B_{n-1,j}C_{n-1,i}+C_{n-1,j}B_{n-2,i},\\\
C_{n-1,i}C_{n-2,j}&=C_{n-1,j}C_{n-2,i},\end{split}$

which reduces to (3.4) for ${\mathbb{F}}_{n}={\mathbb{P}}_{n}$. Note that the
conditions (3.6) and (3.4), respectively, become trivial for $d=1$ and thus
are a purely multivariate phenomenon, i.e., are relevant for $d\geq 2$ only.

We can use the matrices occurring in the three-term recurrence (3.1) to define
a family of tridiagonal infinite matrices ${\mathsf{J}}_{i}$, $1\leq i\leq d$,
as follows:

(3.7) ${\mathsf{J}}_{i}:=\left[\begin{matrix}B_{0,i}&A_{0,i}&&\bigcirc\cr
A_{0,i}^{\mathsf{T}}&B_{1,i}&A_{1,i}&&\cr&A_{1,i}^{\mathsf{T}}&B_{2,i}&\ddots\cr\bigcirc&&\ddots&\ddots\end{matrix}\right],\qquad
1\leq i\leq d.$

The relations (3.4) are then equivalent to the formal commutativity of
${\mathsf{J}}_{i}$, that is,
${\mathsf{J}}_{i}{\mathsf{S}}{\mathsf{J}}_{j}={\mathsf{J}}_{j}{\mathsf{S}}{\mathsf{J}}_{i}$,
$1\leq i,j\leq d$, where
${\mathsf{S}}=\mathrm{diag}\\{S_{0},S_{1},S_{2},\ldots\\}$ is the infinite
block diagonal matrix formed by the sign matrices of the functional
${\mathcal{L}}$.

Since $A_{n,i}:r_{n}^{0}\times r_{n+1}^{0}$ has full rank, there exist
matrices $D_{n,i}:r_{n}^{0}\times r_{n+1}^{0}$, $i=1,\dots,d$, such that

$\sum_{i=1}^{d}D_{n,i}^{\mathsf{T}}A_{n,i}=I.$

The matrices $D_{n,i}$ are not unique and we can take, for example,
$(D_{n,1},\ldots D_{n,d})$ as the generalized inverses of
$(A_{n,1},\ldots,A_{n,d})$. The three-term relations (3.1) lead to a recursive
definition of ${\mathbb{P}}_{n}$, given by

(3.8)
$\displaystyle{\mathbb{P}}_{n+1}(x)=\sum_{i=1}^{d}D_{n,i}^{\mathsf{T}}x_{i}S_{n}{\mathbb{P}}_{n}(x)-E_{n}{\mathbb{P}}_{n}(x)-F_{n}{\mathbb{P}}_{n-1}(x),$

where the matrices $E_{n}$ and $F_{n}$ are given by

$E_{n}=\sum_{i=1}^{d}D_{n,i}^{\mathsf{T}}B_{n,i}\quad\hbox{and}\quad
F_{n}=\sum_{i=1}^{d}D_{n,i}^{\mathsf{T}}A_{n-1,i}^{\mathsf{T}},$

respectively. The relation (3.8) can be used for a recursive definition of the
sequence of ${\mathbb{P}}_{n}$ if $A_{n,i}$, $B_{n,i}$ are given. In one
variable, the sequence of polynomials so defined automatically satisfies the
three-term relation (3.1). In several variables, however, we have to require
the matrices $A_{n,i}$ and $B_{n,i}$ to satisfy some necessary conditions.
Indeed, a slightly modified version of the proof given in [4, Theorem 3.5.1]
for the positive definite setting shows that if $A_{n,i}$ and $B_{n,i}$
satisfy the commuting conditions (3.4) and the rank condition (3.2), then a
sequence of polynomials $\\{{\mathbb{P}}_{n}\\}_{n\geq 0}$ defined recursively
by (3.8) satisfies the three-term relations (3.1).

The three-term relations characterize the orthogonality by Favard’s theorem.
We call the linear functional ${\mathcal{L}}$ _quasi-determinate_ if there is
a basis ${\mathcal{B}}$ of $\Pi^{d}$ such that, for any $P,Q\in{\mathcal{B}}$,

${\mathcal{L}}(PQ)=0\quad\hbox{if $P\neq Q$ and ${\mathcal{L}}(P^{2})\neq
0$}.$

######  Theorem 3.2.

Let $\\{{\mathbb{P}}_{n}\\}_{n\geq
0}=\\{P_{\alpha}\in\Pi_{n}^{d}:|{\alpha}|=n\in{\mathbb{N}}_{0}\\}$ be an
arbitrary polynomial sequence. Then the following statements are equivalent.

  1. (i)

There exists a linear functional ${\mathcal{L}}$ which defines a quasi-
definite linear functional on $\Pi^{d}$ and which makes
$\\{P_{\alpha}^{n}:{\alpha}\in{\mathbb{N}}_{0}^{d}\\}$ an orthogonal basis of
polynomials.

  2. (ii)

For $n\geq 0$, $1\leq i\leq d$, there exist matrices $A_{n,i}$ and $B_{n,i}$
such that

    1. (a)

the polynomials ${\mathbb{P}}_{n}$ satisfy the three-term relation (3.1);

    2. (b)

the matrices in the relation satisfy the rank conditions (3.2).

For the OPs in the positive definite case, the above theorem holds with
$S_{n}=I$, $n\geq 0$. The proof of the theorem follows with obvious
modification of that in the positive definite case; see [4, Theorem 3.3.7].

As a consequence of the three-term relation, we can derive a Christoffel-
Darboux formula. Let

(3.9)
$K_{n}(x,y)=\sum_{k=0}^{n}{\mathbb{P}}_{k}(x)^{\mathsf{T}}{\mathsf{S}}_{k}\overline{{\mathbb{P}}_{k}(y)},\qquad
x,y\in{\mathbb{C}}^{d}.$

Then $K_{n}$ is the reproducing kernel of the space $\Pi_{n}^{d}$ since, for
$0\leq j\leq n$ and $x,y\in{\mathbb{R}}^{d}$,

${\mathcal{L}}(K_{n}(x,\cdot){\mathbb{P}}_{j}^{\mathsf{T}})=\sum_{k=0}^{n}{\mathbb{P}}_{k}(x)^{\mathsf{T}}S_{k}{\mathcal{L}}({\mathbb{P}}_{k}{\mathbb{P}}_{j}^{\mathsf{T}})={\mathbb{P}}_{j}(x)^{\mathsf{T}}.$

######  Theorem 3.3.

For $n=0,1,2,\ldots$, the following Christoffel-Darboux formula holds

(3.10)
$K_{n}(x,y)=\frac{\left[A_{n,j}{\mathbb{P}}_{n+1}(x)\right]^{\mathsf{T}}\overline{{\mathbb{P}}_{n}(y)}-{\mathbb{P}}_{n}(x)^{\mathsf{T}}\left[A_{n,j}\overline{{\mathbb{P}}_{n+1}(y)}\right]}{x_{j}-\overline{y}_{j}},\quad
x,y\in{\mathbb{C}}^{d},$

for $1\leq j\leq d$. In particular, taking the limit shows that

(3.11)
$K_{n}(x,x)=\left[A_{n,i}{\mathbb{P}}_{n+1}(x)\right]^{\mathsf{T}}\partial_{j}{\mathbb{P}}_{n}(x)-{\mathbb{P}}_{n}(x)^{\mathsf{T}}\left[A_{n,i}\partial_{j}{\mathbb{P}}_{n+1}(x)\right].$

###### Proof.

By the three-term relation (3.1),

$\displaystyle
x_{j}S_{k}{\mathbb{P}}_{k}(x)^{\mathsf{T}}\overline{{\mathbb{P}}_{k}(y)}-\overline{y}_{j}S_{k}{\mathbb{P}}_{k}(x)^{\mathsf{T}}\overline{{\mathbb{P}}_{k}(y)}$
$\displaystyle\quad=\left([A_{k,i}{\mathbb{P}}_{k+1}(x)]^{\mathsf{T}}\overline{{\mathbb{P}}_{k}(y)}-[A_{k,i}\overline{{\mathbb{P}}_{k+1}(y)}]^{\mathsf{T}}{\mathbb{P}}_{k}(x)\right)$
$\displaystyle\quad-\left([A_{k-1,i}{\mathbb{P}}_{k}(x)]^{\mathsf{T}}\overline{{\mathbb{P}}_{k-1}(y)}-[A_{k-1,i}\overline{{\mathbb{P}}_{k}(y)}]^{\mathsf{T}}{\mathbb{P}}_{k-1}(x)\right),$

and summing this identity over $k$ gives (3.10). Since ${\mathbb{P}}_{n}$ has
real coefficients, we see that
$\overline{\mathbb{P}}_{n}(x)={\mathbb{P}}_{n}(\bar{x})$. Hence, we can write

$\displaystyle\left[A_{n,j}{\mathbb{P}}_{n+1}(x)\right]^{\mathsf{T}}\overline{{\mathbb{P}}_{n}(y)}-{\mathbb{P}}_{n}(x)^{\mathsf{T}}\left[A_{n,j}\overline{{\mathbb{P}}_{n+1}(y)}\right]$
$\displaystyle=\left[A_{n,j}{\mathbb{P}}_{n+1}(x)\right]^{\mathsf{T}}({\mathbb{P}}_{n}(\bar{y})-{\mathbb{P}}_{n}(x))+{\mathbb{P}}_{n}(x)^{\mathsf{T}}A_{n,j}({\mathbb{P}}_{n+1}(x)-{\mathbb{P}}_{n+1}(\bar{y})).$

Dividing by $x_{j}-\overline{y}_{j}$ and taking the limit
$\overline{y}_{j}\mapsto x$ proves (3.11). ∎

###  3.3. Jacobi matrices and common zeros

We now recall the _truncated Jacobi matrices_ ${\mathsf{J}}_{n,i}$, defined
for each $n\in{\mathbb{N}}_{0}$ by

(3.12) ${\mathsf{J}}_{n,i}:=\left[\begin{matrix}B_{0,i}&A_{0,i}&&&\bigcirc\cr
A_{0,i}^{\mathsf{T}}&B_{1,i}&A_{1,i}&&\cr&\ddots&\ddots&\ddots&\cr&&A_{n-3,i}^{\mathsf{T}}&B_{n-2,i}&A_{n-2,i}\cr\bigcirc&&&A_{n-2,i}^{\mathsf{T}}&B_{n-1,i}\end{matrix}\right]\in{\mathbb{R}}^{r_{n-1}\times
r_{n-1}},\qquad 1\leq i\leq d.$

An element $\lambda\in{\mathbb{C}}^{d}$ is called a generalized joint
eigenvalue of ${\mathsf{J}}_{n,1},\ldots,{\mathsf{J}}_{n,d}$, if there is
$\xi\neq 0$, $\xi\in{\mathbb{C}}^{d}$, such that

${\mathsf{J}}_{n,i}\xi=\lambda_{i}{\mathsf{S}}_{n}\xi,\qquad
i=1,\ldots,d,\qquad{\mathsf{S}}:=\left[\begin{array}[]{ccc}S_{0}&&\\\
&\ddots&\\\ &&S_{n}\end{array}\right];$

the vector $\xi$ is called a joint eigenvector of the generalized eigenvalue
${\lambda}$.

######  Theorem 3.4.

A point $z\in{\mathbb{C}}^{d}$ is a zero of ${\mathbb{P}}_{n}$ if and only if
it is a generalized joint eigenvalue of
${\mathsf{J}}_{n,1},\ldots,{\mathsf{J}}_{n,d}$, the associated joint
eigenvector of $z$ is
$({\mathbb{P}}_{0}^{\mathsf{T}}(z),\ldots,{\mathbb{P}}_{n-1}^{\mathsf{T}}(z))^{\mathsf{T}}$.
Furthermore, if $z$ is a complex-valued generalized joint eigenvalue, then
$\bar{z}$ is also a generalized joint eigenvalue.

###### Proof.

The proof that a common zero is the same as a joint eigenvalue follows almost
verbatim the proof for the positive definite case. In particular,
${\mathbb{P}}_{n}$ can have at most $r_{n-1}=\dim\Pi_{n-1}^{d}$ zeros since
${\mathsf{J}}_{n,i}$ is of size $r_{n-1}\times r_{n-1}$. Now, since the
moments are real, so are the matrices ${\mathsf{J}}_{n,i}$. If $z$ is a joint
eigenvalue with joint eigenvector $\xi$, then
${\mathsf{J}}_{n,i}\xi=z{\mathsf{S}}_{n}\xi$; taking the conjugate of this
identity shows that $\bar{z}$ is also a joint eigenvalue with $\bar{\xi}$ as
its joint eigenvector. ∎

In the positive definite case, ${\mathsf{S}}_{n}$ is the identity matrix and
the generalized eigenvalues are just the ordinary eigenvalues, so that the
common zeros of ${\mathbb{P}}_{n}$ are exactly the eigenvalues of
${\mathsf{J}}_{n,i}$. Since ${\mathsf{J}}_{n,i}$ are symmetric, this implies
that all zeros are real when ${\mathcal{L}}$ is positive definite. We call
$x\in{\mathbb{C}}^{d}$ a simple zero of ${\mathbb{P}}_{n}$ if the Jacobian
$D{\mathbb{P}}_{n}$ of the polynomials in
${\mathbb{P}}_{n}:{\mathbb{R}}^{d}\to{\mathbb{R}}^{r_{n}^{0}}$ is nonsingular.

######  Proposition 3.5.

If $x$ is a common zero of ${\mathbb{P}}_{n}$ and $K_{n}(x,x)\neq 0$, then $x$
must be simple and $x$ cannot be a zero of ${\mathbb{P}}_{n+1}$.

This is a consequence of the identity (3.11). In the positive definite
setting, ${\mathsf{S}}_{n}$ is identity, so that $K_{n}(x,x)\geq
K_{0}(x,x)=1$. Hence, all zeros of ${\mathbb{P}}_{n}$ must be simple by (3.11)
and, furthermore, ${\mathbb{P}}_{n}$ and ${\mathbb{P}}_{n+1}$ do not have
common zeros. In the non-positive definite case, $K_{n}(x,x)$ may not be
positive because the sign changes in ${\mathsf{S}}_{n}$ and it can be zero.
However, since ${\mathsf{S}}_{n}$ is a signature matrix, it follows that

$K_{n}(x,x)=\sum_{k=0}^{n}\sum_{|{\alpha}|=k}{\varepsilon}_{\alpha}|P_{\alpha}(x)|^{2},\qquad{\varepsilon}_{\alpha}=\pm
1,$

so that $K_{n}(x,x)$ is real–valued even if $x\in{\mathbb{C}}^{d}$.

By Theorem 3.4, the polynomials in ${\mathbb{P}}_{n}$ can have at most
$r_{n-1}=\dim\Pi_{n-1}^{d}$ common zeros. The following theorem characterizes
when the number of common zeros is maximal. To prove it, we need some
algebraic preparations, which we shall tailor to our needs.

The _ideal_ generated by ${\mathbb{P}}_{n}$ is defined as

$\left\langle{\mathbb{P}}_{n}\right\rangle:={\mathbb{P}}_{n}^{\mathsf{T}}\left(\Pi^{d}\right)^{r_{n}^{0}}=\left\\{\sum_{|\alpha|=n}P_{\alpha}Q_{\alpha}:Q_{\alpha}\in\Pi^{d}\right\\},\qquad{\mathbb{P}}_{n}=\left[P_{\alpha}:|\alpha|=n\right],$

and ${\mathbb{P}}_{n}$ is called an _H-basis_ provided that

(3.13)
$\left\langle{\mathbb{P}}_{n}\right\rangle\cap\Pi_{m}^{d}={\mathbb{P}}_{n}^{\mathsf{T}}\left(\Pi_{m-n}^{d}\right)^{r_{n}^{0}},\qquad
m\geq n,$

which means that any polynomial of degree $m$ in the ideal
$\left\langle{\mathbb{P}}_{n}\right\rangle$ has a representation with respect
to ${\mathbb{P}}_{n}$ in which the degree of no term in the sum exceeds $m$.
Since ${\mathbb{P}}_{n}$ consists of $r_{n}^{0}$ linearly independent
polynomials of degree $n$ whose leading forms span ${\mathcal{P}}_{n}^{d}$, it
follows that $\Pi^{d}/\left\langle{\mathbb{P}}_{n}\right\rangle=\Pi_{n-1}^{d}$
if and only if ${\mathbb{P}}_{n}$ is an H-basis. For more information on
H-bases cf. [8, 17]. Moreover, ${\mathbb{P}}_{n}$ defines a _reduced
multiplication_ by the factors $Q\mapsto(x_{i}Q)_{{\mathbb{P}}_{n}}$ on
$\Pi_{n-1}^{d}$ by setting

$(x_{i}Q)(x)=x_{i}Q(x)-{\mathbb{P}}_{n}(x)^{\mathsf{T}}q,$

where the factor $q\in{\mathbb{R}}^{r_{n}^{0}}$ is chosen such that
$\deg\left(x_{i}Q(x)-{\mathbb{P}}_{n}(x)^{\mathsf{T}}q\right)\leq n-1$. In
other words, if $\deg x_{i}Q(x)=n$, then $q$ is chosen such that the
homogeneous leading form of ${\mathbb{P}}_{n}(x)^{\mathsf{T}}$ coincides with
the leading form of $x_{i}Q(x)$, otherwise we have $q=0$. Of course,
$(x_{i}\cdot)_{{\mathbb{P}}_{n}}$ is a linear map from $\Pi_{n-1}^{d}$ to
$\Pi_{n-1}^{d}$ and its matrix representation with respect to the basis
${\mathbb{P}}_{0},\dots,{\mathbb{P}}_{n-1}$ is
${\mathsf{M}}_{n,i}:={\mathsf{S}}_{n}{\mathsf{J}}_{n,i}$, $i=1,\dots,d$. More
precisely, if
$Q=\widehat{Q}_{0}^{T}{\mathbb{P}}_{0}+\cdots+\widehat{Q}_{n-1}^{T}{\mathbb{P}}_{n-1}$,
then

$(x_{i}Q)_{{\mathbb{P}}_{n}}=\left[\widehat{Q}_{0}^{T},\dots,\widehat{Q}_{n-1}^{T}\right]{\mathsf{S}}_{n}{\mathsf{J}}_{n,i}\left[\begin{array}[]{c}{\mathbb{P}}_{0}\\\
\vdots\\\ {\mathbb{P}}_{n-1}\end{array}\right].$

This is a a direct consequence of multiplying both sides of (3.1) by $S_{n}$
and the definition in (3.12). The matrices ${\mathsf{M}}_{n,i}$ are called the
_multiplication tables_ for the coordinate multiplication modulo
${\mathbb{P}}_{n}$. They can be used to characterize whether
${\mathbb{P}}_{n}$ is an H-basis; the following result is a special case of
characterization for border bases, given for example in [12].

######  Lemma 3.6.

${\mathbb{P}}_{n}$ is an H-basis if and only if the multiplication tables
${\mathsf{M}}_{n,i}$ commute.

Indeed, if ${\mathbb{P}}_{n}$ is an H-basis, then $M_{n,i}$ describes the
action of multiplication by $x_{i}$ on
$\Pi^{d}/\left\langle{\mathbb{P}}_{n}\right\rangle=\Pi_{n-1}^{d}$ and since
multiplication is commutative, so are the matrices that represent it. If, on
the other hand, the matrices commute, then for any $|\alpha|\leq n-1$ the
matrix
${\mathsf{M}}^{\alpha}={\mathsf{M}}_{n,1}^{\alpha_{1}}\cdots{\mathsf{M}}_{n,d}^{\alpha_{d}}$
is well defined and its first row contains the coefficients of the monomial
$x^{\alpha}$ with respect to ${\mathbb{P}}_{0},\dots,{\mathbb{P}}_{n-1}$.
Hence, there are $r_{n-1}$ linearly independent reduced polynomials, hence
$\dim\Pi^{d}/\left\langle{\mathbb{P}}_{n}\right\rangle\geq r_{n-1}$ and
therefore ${\mathbb{P}}_{n}$ is an H-basis and, in particular, has $r_{n-1}$
common zeros counting multiplicities.

######  Theorem 3.7.

The polynomials in ${\mathbb{P}}_{n}$ have $r_{n-1}$ common zeros if and only
if

(3.14) $\displaystyle
A_{n-1,i}S_{n}A_{n-1,j}^{\mathsf{T}}=A_{n-1,j}S_{n}A_{n-1,i}^{\mathsf{T}},\qquad
1\leq i,j\leq d.$

###### Proof.

According to Lemma 3.6, to have the maximal number of common zeros is
equivalent to the commuting of the ${\mathsf{S}}_{n}{\mathsf{J}}_{n,i}$ which
is in turn equivalent to
${\mathsf{J}}_{n,i}{\mathsf{S}}_{n}J_{n,j}={\mathsf{J}}_{n,j}{\mathsf{S}}_{n}J_{n,i}$,
$1\leq i\neq j\leq d$. Using the commutativity of the Jacobi matrices, it is
not difficult to see that we only need to consider the commutativity of
$Q_{i}S_{n}Q_{j}^{\mathsf{T}}$, where $Q_{i}$ and $Q_{j}$ denote the last row
of ${\mathsf{J}}_{n,i}$ and ${\mathsf{J}}_{n,j}$, respectively; that is,

$\displaystyle A$
${}_{n-2,i}^{\mathsf{T}}S_{n-2}A_{n-2,j}+B_{n-1,i}S_{n-1}B_{n-1,j}=A_{n-2,j}^{\mathsf{T}}S_{n-2}A_{n-2,i}+B_{n-1,j}S_{n-1}B_{n-1,i},$

which is equivalent to (3.14) by the third identity of (3.4) with $n$ replaced
by $n-1$. ∎

The same argument can also be applied to an arbitrary recurrence relation of
the form (3.5) with the multiplication table

(3.15)
${\mathsf{J}}_{n,i}^{\prime}:=\left[\begin{matrix}B_{0,i}&A_{0,i}&&&\bigcirc\cr
C_{1,i}&B_{1,i}&A_{1,i}&&\cr&\ddots&\ddots&\ddots&\cr&&C_{n-2,i}&B_{n-2,i}&A_{n-2,i}\cr\bigcirc&&&C_{n-1,i}&B_{n-1,i}\end{matrix}\right],\qquad
1\leq i\leq d,$

which leads to the following variant of Theorem 3.7, that in particular
includes the monic orthogonal polynomials $\widehat{\mathbb{P}}_{n}$.

######  Corollary 3.8.

The polynomials ${\mathbb{F}}_{n}$ from (3.5) have $r_{n-1}$ common zeros if
and only if

(3.16) $A_{n-2,i}\,C_{n-1,j}=A_{n-2,j}\,C_{n-1,i},\qquad 1\leq i,j\leq d.$

The conditions (3.14) and (3.16) are highly non-trivial. For example, the
condition does not hold if ${\mathcal{L}}$ is centrally symmetric, which means
that ${\mathcal{L}}$ satisfies

${\mathcal{L}}(x^{\alpha})=0,\qquad|\alpha|\in 2{\mathbb{N}}_{0}+1;$

if ${\mathcal{L}}f=\int_{\Omega}f(x)W(x)\mathrm{d}x$, then ${\mathcal{L}}$ is
centrally symmetric if $W(x)=W(-x)$ whenever $x\in\Omega$ implies
$-x\in\Omega$. It is known that if ${\mathcal{L}}$ is positive definite and
centrally symmetric, then ${\mathbb{P}}_{n}$ has no common zeros for $n$ even
and only one common zero (the origin) if $n$ is odd [4, Section 3.7].

When ${\mathbb{P}}_{n}$ has maximal common zeros, we can consider
interpolation and cubature rules as in one variable.

######  Proposition 3.9.

Assume (3.14) that holds. Let $Z_{n}=\\{(\xi_{k,1:}\ldots,\xi_{k,d}),1\leq
k\leq r_{n-1}\\}$ be the set of common zeros of ${\mathbb{P}}_{n}$. Assume
that $\xi\in Z_{n}$ are not zeros of ${\mathbb{P}}_{n+1}$. Then the unique
polynomial, $L_{n}f$, in $\Pi_{n}^{d}$ that satisfies

$L_{n}f(\zeta)=f(\zeta),\qquad\zeta\in Z_{n},$

for any generic function $f$ is given by the formula

(3.17) $L_{n}f(x)=\sum_{\zeta\in
Z_{n}}f(\zeta)\frac{K_{n}(x,\zeta)}{K_{n}(\zeta,\zeta)}.$

###### Proof.

By our assumption on the zeros, $K_{n}(\zeta,\zeta)\neq 0$ for $\zeta\in
Z_{n}$. By Theorem 3.4, if $\zeta$ is in $Z_{n}$, then so is $\bar{\zeta}$. If
$\zeta,\eta\in Z_{n}$ and $\zeta\neq\eta$, then there must be a $j$ such that
$\zeta_{j}\neq\overline{\eta_{j}}$. Hence, by (3.10), $K_{n}(\zeta,\eta)=0$.
Consequently, it follows that $L_{n}f$ in (3.17) satisfies the interpolation
conditions. ∎

Applying ${\mathcal{L}}$ to $L_{n}f$ leads to a cubature rule on $Z_{n}$,
which is a Gaussian cubature.

######  Theorem 3.10.

Let ${\mathcal{L}}(f)=\int_{\Omega}f(x)w(x)\mathrm{d}x$ for a weight $w$
defined on $\Omega$ and assume that $Z_{n}$ is defined as in Proposition 3.9.
Then

(3.18) $\int_{\Omega}f(x)w(x)\mathrm{d}x=\sum_{\zeta\in
Z_{n}}\frac{1}{K_{n}(\zeta,\zeta)}f(\zeta),\qquad f\in\Pi_{2n-1}^{d}.$

Moreover, the cubature rule is real-valued for all polynomials.

###### Proof.

That the cubature rule is of degree $2n-1$ follows as in the positive definite
case. If $\zeta\in Z_{n}$ is a complex number, then $\bar{\zeta}\in Z_{n}$.
Since $f$ is a real-valued polynomial, $\overline{f}(\zeta)=f(\bar{\zeta})$.
Moreover, $K_{n}(\zeta,\zeta)=K_{n}(\bar{\zeta},\bar{\zeta})$ is real, so that

$\frac{1}{K_{n}(\zeta,\zeta)}f(\zeta)+\frac{1}{K_{n}\left(\bar{\zeta},\bar{\zeta}\right)}\overline{f(\zeta)}=\frac{2}{K_{n}(\zeta,\zeta)}\Re
f(\zeta)$

is real-valued for each complex-valued $\zeta\in Z_{n}$. ∎

Since the Gaussian cubature rule requires ${\mathbb{P}}_{n}$ to have $r_{n-1}$
simple common zeros, they exist rarely and do not exist, for example, if
${\mathcal{L}}$ is centrally symmetric. Nevertheless, they do exist in some
cases. In the next subsection, we recall a class of weight functions for which
they exist.

###  3.4. Example of Gaussian cubature rules

We start with a weight function $w$ on ${\mathbb{R}}$, which is definite and
may change signs in its support set $[a,b]$. Let $p_{n}(w;x)$ be the
orthogonal polynomial with respect to $w$ that satisfies

$\int_{{\mathbb{R}}}p_{n}(w;t)p_{m}(w;t)w(t)\mathrm{d}t={\varepsilon}_{n}\delta_{n,m},\qquad
n\neq m,\quad{\varepsilon}=\pm 1.$

We assume that $p_{n}(w)$ has $n$ simple, real zeros in $[a,b]$ so that $w$
admits the Gaussian quadrature rule (2.2) of degree $2n-1$,

$\int_{\mathbb{R}}f(t)w(t)\mathrm{d}t=\sum_{k=1}^{n}{\lambda}_{k,n}f(t_{k,n}),\qquad\deg
f\leq 2n-1.$

We define two families of symmetric polynomials, indexed by partitions in

$\Lambda^{d}=\\{{\lambda}\in{\mathbb{N}}_{0}^{d+1}:{\lambda}_{1}\geq{\lambda}_{2}\geq\cdots\geq{\lambda}_{d}\\}.$

Let $\mathcal{S}_{d}$ be the symmetric group of $d$ elements. For
${\alpha}\in\Lambda^{d}$, define

$Q_{\alpha}^{-\frac{1}{2}}(x):=\sum_{\beta\in\mathcal{S}_{d}}p_{\alpha_{1}}(x_{\beta_{1}})\ldots
p_{\alpha_{d}}(x_{\beta_{d}}).$

Then $\\{Q_{\alpha}^{-\frac{1}{2}}:{\alpha}\in\Lambda^{d}\\}$ is a family of
symmetric polynomials and satisfies

(3.19)
$\frac{1}{d!}\int_{{\mathbb{R}}^{d}}Q_{\alpha}^{-\frac{1}{2}}(x)\,Q_{\beta}^{-\frac{1}{2}}(x)\,\prod_{i=1}^{d}w(x_{i})\mathrm{d}x=m_{1}!\ldots
m_{d^{\prime}}!\,{\varepsilon}_{{\alpha}_{1}}\ldots{\varepsilon}_{{\alpha}_{d}}\delta_{{\alpha},{\beta}},$

where $d^{\prime}$ is the number of distinct elements in ${\alpha}$ and
$m_{i}$ is the number of occurrences of the $i$-th distinct element in
$\alpha$. Moreover, define

$J_{\alpha}(x):=\det\,[p_{\alpha_{i}+d-i}(x_{j})]_{i,j=1}^{d}\quad\hbox{and}\quad
J(x)=\prod_{1\leq i<j\leq d}(x_{i}-x_{j}).$

Then the family of polynomials
$\\{Q_{\alpha}^{\frac{1}{2}}:{\alpha}\in\Lambda^{d}\\}$, defined by

$Q_{\alpha}^{\frac{1}{2}}(x)=\frac{J_{\alpha}(x)}{J(x)},\qquad{\alpha}\in\Lambda^{d},$

consists of symmetric orthogonal polynomials that satisfy

(3.20)
$\frac{1}{d!}\int_{{\mathbb{R}}^{d}}Q_{\alpha}^{\frac{1}{2}}(x)\,Q_{\beta}^{\frac{1}{2}}(x)\,\left[J(x)\right]^{2}\prod_{i=1}^{d}w(x_{i})\mathrm{d}x={\varepsilon}_{{\alpha}_{1}}\ldots{\varepsilon}_{{\alpha}_{d}}\delta_{{\alpha},{\beta}};$

see, for example, [4, Sect. 5.4]. These symmetric polynomials can be mapped to
a family of ordinary orthogonal polynomials. Let

$\mathcal{R}=\\{x\in{\mathbb{R}}^{d}\mid\,x_{1}<x_{2}<\cdots<x_{d},\,x_{j}\in\mathrm{supp}(w)\\}.$

The mapping $x\in\mathcal{R}\mapsto u\in\Omega$, defined by

(3.21) $u_{k}=e_{k}(x_{1},\ldots x_{d}):=\sum_{1\leq i_{1}<\ldots<i_{k}\leq
d}x_{i_{1}}\cdots x_{i_{d}},\quad k=1,2,\ldots,d,$

sends $\mathcal{R}$ onto a domain $\Omega$ and the Jacobian of the map is
$J(x)$. Under the mapping (3.21), the polynomials
$Q_{\alpha}^{\pm\frac{1}{2}}(x)$ become polynomials in the variable(s) $u$,
and these polynomials are orthogonal with respect to the weight function
$W_{\pm\frac{1}{2}}$ on $\Omega$, where

$W_{-\frac{1}{2}}(u)=\prod_{i=1}^{d}w(x_{i})[\Delta(u)]^{-\frac{1}{2}}\quad\hbox{and}\quad
W_{\frac{1}{2}}(u)=\prod_{i=1}^{d}w(x_{i})[\Delta(u)]^{\frac{1}{2}}.$

######  Proposition 3.11.

Let $Q_{\alpha}^{\pm\frac{1}{2}}$ be the symmetric orthogonal polynomials in
(3.19) and (3.20). For ${\alpha}\in\Lambda^{d}$ and
$n={\alpha}_{1}\geq\ldots\geq{\alpha}_{d}$, define

$P_{\alpha}^{n,-\frac{1}{2}}(u)=Q_{\alpha}^{-\frac{1}{2}}(x)\quad\hbox{and}\quad
P_{\alpha}^{n,\frac{1}{2}}(u)=Q_{\alpha}^{\frac{1}{2}}(x),$

where $u=u(x)$ as in (3.21). Then
${\mathbb{P}}_{n}^{\pm\frac{1}{2}}=\\{P_{\alpha}^{n,\pm\frac{1}{2}}:n={\alpha}_{1}\geq\ldots\geq{\alpha}_{d}\geq
0\\}$ is an orthogonal basis of ${\mathcal{V}}_{n}(W_{\pm\frac{1}{2}})$.

The polynomials in $\\{P_{\alpha}^{n,\pm\frac{1}{2}}:|{\alpha}|=n\\}$ have
$\dim\Pi_{n-1}^{d}$ common zeros and they are nodes of the Gaussian cubature
rules. More precisely, we have the following [3]:

######  Theorem 3.12.

The weight function $W_{\pm\frac{1}{2}}$ admits Gaussian cubature rules of
degree $2n-1$ for $n=1,2,\ldots$. More precisely,

(3.22)
$\int_{\Omega}f(u)W_{-\frac{1}{2}}(u)\mathrm{d}u=\sum_{{\gamma}_{1}=1}^{n}\sum_{{\gamma}_{2}=1}^{{\gamma}_{1}}\cdots\sum_{{\gamma}_{d}=1}^{{\gamma}_{d-1}}\Lambda_{{\gamma},n}^{(-\frac{1}{2})}f(u_{{\gamma},n}),$

where $u_{{\gamma},n}$ is the image of
$x_{{\gamma},n}=(t_{{\gamma}_{1},n},t_{{\gamma}_{2},n},\ldots,t_{{\gamma}_{d},n})$
under the mapping (3.21), and

(3.23)
$\int_{\Omega}f(u)W_{\frac{1}{2}}(u)\mathrm{d}u=\sum_{{\gamma}_{1}=1}^{n+d-1}\sum_{{\gamma}_{2}=1}^{{\gamma}_{1}-1}\cdots\sum_{{\gamma}_{d}=1}^{{\gamma}_{d-1}-1}\Lambda_{{\gamma},n}^{(\frac{1}{2})}f(u_{{\gamma},n+d-1}),$

where $u_{{\gamma},n+d-1}$ is the image of
$x_{{\gamma},n+d-1}=(t_{{\gamma}_{1},n+d-1},t_{{\gamma}_{2},n+d-1},\ldots,t_{{\gamma}_{d},n+d-1})$
under the mapping (3.21) and, furthermore,

(3.24)
$\Lambda_{{\gamma},n}^{(-\frac{1}{2})}=\frac{{\lambda}_{\hat{\gamma}_{1},n}^{m_{1}}\dots\lambda_{\hat{\gamma}_{d^{\prime}},n}^{m_{d^{\prime}}}}{m_{1}!\cdots
m_{d^{\prime}}!}\quad\hbox{and}\quad\Lambda_{{\gamma},n}^{(\frac{1}{2})}=J(t_{{\gamma},n})^{2}{\lambda}_{{\gamma}_{1},n}\dots\lambda_{{\gamma}_{d},n},$

where $\hat{\gamma}_{1}\ldots,\hat{\gamma}_{d^{\prime}}$ denote the distinct
elements in ${\gamma}$ and $m_{1},\ldots,m_{d^{\prime}}$ denote their
respective multiplicities.

The proof of this theorem is purely algebraic and holds also for the case that
$w$ may not be positive. For another example of the Gaussian cubature formula,
see [10], where the domain is the image of the simplex under the mapping of
elementary symmetric functions. For $d=2$, the domain is bounded by a
hypocycloid.

##  4\. Continued fraction expansions of moment sequences

To recover an analogy for the continued fractions in several variables, we let
$\mu=\left(\mu_{\alpha}:\alpha\in{\mathbb{N}}_{0}\right)$ be a moment sequence
such that $\det M_{n}\neq 0$, $n\geq 0$, and consider the associated formal
Laurent series

$\mu(z)=\sum_{{\alpha}\in{\mathbb{N}}_{0}^{d}}\mu_{\alpha}z^{-\alpha},\qquad
z\in z\in({\mathbb{C}}\setminus 0)^{d}.$

If there is a rational function $r_{n}(z)$ of degrees $\leq n$ for $n\geq 0$,
such that

$\mu(z)-r_{n}(z)=O\left(z^{-2n}\right),\qquad\text{i.e..}\qquad\mu(z)-r_{n}(z)=\sum_{|\alpha|\geq
2n}\lambda_{\alpha}z^{-\alpha},$

then $r_{n}(z)$ can be regarded as an analogy for the continued fraction in
several variables.

Our goals in this section are twofold. The first one is to construct a
sequence of moment sequences $\mu^{n}$, $n\geq 0$, iteratively from a
polynomial sequence that satisfies the three-term relations and show that
$\mu^{n}\to\mu$, where $\mu$ is the moment sequence that makes the polynomials
orthogonal. The second goal is to define and construct multivariate continued
fractions based on $\mu^{n}$, and to characterize the moment sequence for
which these continued fractions exist.

###  4.1. Moments from three-term recurrence

As shown in the previous section, if ${\mathcal{L}}$ is a definite moment
functional and $\mu$ is its moments, then the orthogonal polynomials satisfy
the three-term relations (3.5), whose matrix conditions in turn satisfy (3.6),
(3.16), and (3.2).

We start from the three-term recurrence relation satisfied by some polynomial
sequence. Let ${\mathbb{F}}_{n}\in(\Pi_{n}^{d})^{r_{n}^{0}}$, $n\geq 0$, be a
sequence of polynomials that satisfies

(4.1)
$x_{i}{\mathbb{F}}_{n}(x)=A_{n,i}{\mathbb{F}}_{n+1}(x)+B_{n,i}{\mathbb{F}}_{n}(x)+C_{n,i}{\mathbb{F}}_{n-1}(x),\qquad
i=1,\dots,d,$

and whose coefficient matrices satisfy (3.6), (3.16) and the rank condition
(3.2). By Theorem 3.2, there exists a quasi-determinate linear functional
${\mathcal{L}}$ such that the ${\mathbb{F}}_{n}$ form an orthogonal basis with
respect to ${\mathcal{L}}$, i.e.,

(4.2) ${\mathcal{L}}(Q{\mathbb{F}}_{n})=0,\quad
Q\in\Pi_{n-1}^{d},\qquad\det\left[{\mathcal{L}}\left({\mathbb{F}}_{n}{\mathbb{F}}_{n}^{\mathsf{T}}\right)\right]\neq
0,$

holds for any $n\geq 0$. We normalize ${\mathcal{L}}$ such that
${\mathcal{L}}(1)=1$ and denote by
$\mu_{\alpha}={\mathcal{L}}\left((\cdot)^{\alpha}\right)$,
$\alpha\in{\mathbb{N}}_{0}^{d}$, the moments of this functional. In this
setting, the moment sequence $\mu$ is well-defined but not explicitly given.

We construct a sequence $\mu^{n}$, $n\in{\mathbb{N}}_{0}$, of moment sequences
recursively in $n$ and based entirely on the recurrence (4.1). This will be
done in such a way that $\mu^{n}$ converges to $\mu$. In other words, we
provide a constructive method for the moment sequence $\mu$ from
${\mathbb{F}}_{n}$ defined by (4.1).

By Corollary 3.8, all ${\mathbb{F}}_{n}$, $n\geq 0$, are H-bases and
$\Pi^{d}/\left\langle{\mathbb{F}}_{n}\right\rangle=\Pi_{n-1}^{d}$ is the
associated interpolation space modulo
$\left\langle{\mathbb{F}}_{n}\right\rangle$. This implies that
${\mathbb{F}}_{n}$ has $r_{n-1}$ common zeros, counting multiplicities. Recall
that the multiplicity of common zeros of an ideal is not a simple number
anymore, but a structural quantity, namely a finite dimensional subspace of
$\Pi$ which is closed under differentiation or _$D$ -invariant_ for short.
More precisely, ${\mathcal{Q}}$ is $D$-invariant if $Q\in{\mathcal{Q}}$
implies that $\frac{\partial}{\partial x_{j}}Q\in{\mathcal{Q}}$ as well; in
particular, this implies that $1\in{\mathcal{Q}}$. Therefore, that
${\mathbb{F}}_{n}$ has $r_{n-1}$ common zeros is equivalent to that there
exist a set $Z_{n}$ of finitely many points in ${\mathbb{C}}^{d}$ and
$D$-invariant subspaces ${\mathcal{Q}}_{n,\zeta}$, $\zeta\in Z_{n}$,
describing the multiplicity, such that
$P\in\left\langle{\mathbb{F}}_{n}\right\rangle$ if and only if

(4.3) $\left(Q(D)P\right)(\zeta)=0,\qquad
Q\in{\mathcal{Q}}_{n,\zeta},\quad\zeta\in Z_{n},$

and $\sum_{\zeta\in Z_{n}}\dim{\mathcal{Q}}_{\zeta}=r_{n-1}$ (cf. [7]). Let
$Q_{n,\zeta}$ be a graded basis of ${\mathcal{Q}}_{n,\zeta}$ that thus
contains the constant polynomial. Then the interpolation problem with respect
to $\left(Q(D)P\right)(\zeta)$, $Q\in{\mathcal{Q}}_{n,\zeta}$, $\zeta\in
Z_{n}$, is _poised_ , i.e., uniquely solvable, in $\Pi_{n-1}$. Equivalently,
the _Vandermonde matrix_

(4.4)
$V_{n}:=\left[\left(Q(D)(\cdot)^{\alpha}\right)(\zeta):\begin{array}[]{c}Q\in{\mathcal{Q}}_{n,\zeta},\,\zeta\in
Z_{n}\\\ |\alpha|\leq n-1\end{array}\right]$

is nonsingular. Note that the case of simple zeros is described in this way as
$\dim{\mathcal{Q}}_{n,\zeta}=1$ and ${\mathcal{Q}}_{n,\zeta}={\mathbb{C}}$,
and the above Vandermonde matrix is the classical one for the _Lagrange
interpolation_ problem, in general (4.4) corresponds to a _Hermite
interpolation problem_ , however.

To reconstruct the moment sequence, the _$\Theta$ -operator_ is useful. It is
defined as

$\theta_{j}:=(\cdot)_{j}\frac{\partial}{\partial x_{j}},\quad
j=1,\dots,d,\qquad\Theta^{\alpha}:=\theta_{1}^{\alpha_{1}}\cdots\theta_{d}^{\alpha_{d}},$

and $Q(\Theta)$, $Q\in\Pi^{d}$, is defined accordingly. It has been shown, for
example in [18], that $V_{n}$ is nonsingular if and only if the _$\Theta$
-Vandermonde matrix_

$V_{\Theta,n}:=\left[\left(Q(\Theta)(\cdot)^{\alpha}\right)(\zeta):\begin{array}[]{c}Q\in{\mathcal{Q}}_{n,\zeta},\,\zeta\in
Z_{n}\\\ |\alpha|\leq n-1\end{array}\right]$

is nonsingular. Moreover, by construction,

$\left[\left(Q(\Theta)(\cdot)^{\alpha}\right)(\zeta):\begin{array}[]{c}Q\in{\mathcal{Q}}_{n,\zeta},\,\zeta\in
Z_{n}\\\ |\alpha|\leq n\end{array}\right]\widehat{\boldsymbol{F}}_{n}=0$

since ${\mathbb{F}}_{n}$ is an H-basis for the ideal defined by the dual
conditions in (4.3).

We now use the _Prony’s method_ to define a sequence of moment sequences
$\mu^{n}$, $n\geq 0$, by iteration, which is defined by using the recurrence
(4.1). To that end, we note that ${\mathbb{F}}_{0}=1$ and set, accordingly
$\mu^{0}_{\alpha}=1$, $\alpha\in{\mathbb{N}}_{0}^{d}$. Next, suppose that we
have already computed a moment sequence $\mu^{n-1}$. To define $\mu^{n}$, we
first find an exponential polynomial $g_{n}$ of the form

(4.5) $g_{n}(x):=\sum_{\zeta\in Z_{n}}G_{n,\zeta}(x)\,\zeta^{x},\qquad
G_{n,\zeta}\in{\mathcal{Q}}_{n,\zeta},$

such that

(4.6) $g_{n}(\alpha)=\mu^{n-1}_{\alpha},\qquad|\alpha|\leq n-1.$

Expanding $G_{n,\zeta}$ with respect to the basis $Q_{n,\zeta}$ of
${\mathcal{Q}}_{n,\zeta}$, we note that

$g_{n}(\alpha)=\sum_{\zeta\in
Z_{n}}\Big{(}\sum_{Q\in{\mathcal{Q}}_{n,\zeta}}c_{Q}Q(\alpha)\Big{)}\,\zeta^{\alpha}=\sum_{\zeta\in
Z_{n}}\sum_{Q\in{\mathcal{Q}}_{n,\zeta}}c_{Q}\left(Q(\Theta)(\cdot)^{\alpha}\right)(\zeta),$

hence the coefficients $c_{Q}$, when (4.6) is satisfied, are obtained by
solving the linear system

$\left[\mu^{n-1}_{\alpha}:|\alpha|\leq
n-1\right]=V_{\Theta,n}^{T}\left[c_{Q}:Q\in Q_{n,\zeta},\zeta\in
Z_{n}\right].$

Since $V_{\Theta,n}$ is square and nonsingular, the coefficients $c_{Q}$,
therefore $g_{n}$, are uniquely defined. Thus, once $\mu^{n-1}$ is present,
the exponential polynomial $g_{n}$ that satisfies (4.6) is uniquely determined
by $\mu^{n-1}$.

We are now ready to define the moment sequence $\mu^{n}$, which is simply

$\mu^{n}_{\alpha}=g_{n}(\alpha),\qquad\alpha\in{\mathbb{N}}_{0}^{d}.$

The polynomials ${\mathbb{F}}_{n}$ generate the so-called _Prony ideal_ for
the sequence $\mu^{n}$ and the associated Hankel matrices

$M^{n}_{k}:=\left[\mu^{n}_{\alpha+\beta}:\begin{array}[]{c}|\alpha|\leq k\\\
|\beta|\leq k\\\ \end{array}\right]$

can be factored into

$M^{n}_{k}=V_{\Theta,n,k}^{T}T_{n}V_{\Theta,n,k},\qquad
V_{\Theta,n,k}:=\left[\left(Q(\Theta)(\cdot)^{\alpha}\right)(\zeta):\begin{array}[]{c}Q\in
Q_{n,\zeta},\,\zeta\in Z_{n}\\\ |\alpha|\leq k\end{array}\right],$

where $T_{n}\in{\mathbb{R}}^{r_{n-1}\times r_{n-1}}$ is a nonsingular block
diagonal matrix with blocks of the size $\dim{\mathcal{Q}}_{\zeta}$, $\zeta\in
Z_{n}$, see [18]. In particular, $\det M^{n}_{n-1}\neq 0$, i.e., $\mu^{n}$ is
definite up to degree $n-1$ and $M^{n}_{n}\widehat{\boldsymbol{F}}_{n}=0$. In
the same way, the infinite Hankel matrix $M^{n}$ factorizes into

$M^{n}=V_{\Theta,n,\infty}^{T}T_{n}V_{\Theta,n,\infty},\qquad
V_{\Theta,n,\infty}:=\left[\left(Q(\Theta)(\cdot)^{\alpha}\right)(\zeta):\begin{array}[]{c}Q\in
Q_{n,\zeta},\,\zeta\in Z_{n}\\\
|\alpha|\in{\mathbb{N}}_{0}^{d}\end{array}\right],$

hence the associated Hankel operator has rank $r_{n-1}$.

######  Theorem 4.1.

For $n\geq 1$, the moments $\mu^{n}$ and $\mu$ from above satisfy

(4.7) $\mu_{\alpha}=\mu^{n}_{\alpha},\qquad|\alpha|\leq 2n-1.$

###### Proof.

For $n\geq 1$ we have that $\mu^{n}_{\alpha}=\mu^{n-1}_{\alpha}=\mu_{\alpha}$
for $|\alpha|\leq n-1$; for $n=1$ this is the fact that
$\mu^{0}_{0}=\mu_{0}=1$ and for $n\geq 2$ it immediately follows by induction.
We now set $\omega_{\alpha}:=\mu_{\alpha}-\mu^{\alpha}$,
$\alpha\in{\mathbb{N}}_{0}^{d}$, hence $\omega_{\alpha}=0$, $|\alpha|\leq
n-1$. We partition the coefficients of ${\mathbb{F}}_{n}$ as

$\widehat{\boldsymbol{F}}_{n}=\left[\begin{array}[]{c}\widehat{F}_{n,n-1}\\\
\widehat{F}_{n,n}\end{array}\right],\qquad\widehat{F}_{n,n-1}\in{\mathbb{R}}^{r_{n}^{0}\times
r_{n-1}},\,F_{n,n}\in{\mathbb{R}}^{r_{n}^{0}\times r_{n}^{0}},$

with $\det{\boldsymbol{F}}_{n,n}\neq 0$ and note that due to the orthogonality
of ${\mathbb{F}}_{n}$ we have $M_{n}\widehat{\boldsymbol{F}}_{n}=0$, hence

$0=\left(M_{n}-M_{n}^{n}\right)\widehat{\boldsymbol{F}}_{n}=W_{n}\widehat{\boldsymbol{F}}_{n}=\left[\begin{array}[]{cc}W_{n-1}^{0}&W_{n}^{0}\\\
\vdots&\vdots\\\
W_{n-1}^{n}&W_{n}^{n}\end{array}\right]\left[\begin{array}[]{c}\widehat{\boldsymbol{F}}_{n,n-1}\\\
\widehat{\boldsymbol{F}}_{n,n}\end{array}\right],$

where

$W_{n-1}^{k}=\left[\omega^{n}_{\alpha+\beta}:\begin{array}[]{c}|\alpha|=k\\\
|\beta|\leq n-1\\\ \end{array}\right],\quad
W_{n}^{k}=\left[\omega^{n}_{\alpha+\beta}:\begin{array}[]{c}|\alpha|=k\\\
|\beta|=n\\\ \end{array}\right],$

which yields that

$W_{n}^{k}=-W_{n-1}^{k}\widehat{\boldsymbol{F}}_{n,n-1}\widehat{\boldsymbol{F}}_{n,n}^{-1}.$

Then $W_{n-1}^{0}=0$ implies that $W_{n}^{1}=0$ and since
$W^{k+1}_{n-1}=\left[W^{k}_{n-1}\,W_{n}^{k}\right]$, a simple induction on $k$
yields that

$0=W_{n-1}^{n}=\left[\omega^{n}_{\alpha+\beta}:\begin{array}[]{c}|\alpha|=n\\\
|\beta|\leq n-1\\\ \end{array}\right].$

Since any $\alpha$ with $|\alpha|\leq 2n-1$ can be written as $\beta+\gamma$,
$|\beta|\leq n$, $|\gamma|\leq n-1$, the claim follows. ∎

At this point, it is worthwhile to remark that the above construction
inductively computes the initial segments $\mu^{n}$ of $\mu$ without knowing
$\mu$ at any point of this process. This holds in accordance with the nature
of continued fractions.

In the case that the Gaussian cubature rules exist, the moment sequence
$\mu^{n}$ that satisfies (4.7) has a simple formulation. Indeed, the existence
of the Gaussian cubature rule means that all zeros in $Z_{n}$ are simple, so
that $\dim{\mathcal{Q}}_{n,\zeta}=1$ and, consequently, $G_{n,\zeta}$ are
scalars. Moreover, they are exactly equal to the coefficients in the Gaussian
cubature rules, which are $K_{n}(\zeta,\zeta)^{-1}$, with $K_{n}$ being the
reproducing kernel defined in (3.9). Indeed, let $g_{n}(x)=\sum_{\zeta\in
Z_{n}}\frac{1}{K_{n}(\zeta,\zeta)}\zeta^{x}$. Then, the Gaussian cubature rule
(3.18) implies immediately that (4.6) holds for every positive integer $n$. We
state the result as a corollary.

######  Corollary 4.2.

Let ${\mathcal{L}}(f)=\int_{\Omega}f(x)w(x)\mathrm{d}x$ for a weight function
defined on $\Omega\subset{\mathbb{R}}^{d}$, for which the Gaussian cubature
rule (3.18) holds. Then the moment $\mu^{n}$ in Theorem 4.1 is precisely given
by

$\mu_{\alpha}^{n}=\sum_{\zeta\in
Z_{n}}\frac{1}{K_{n}(\zeta,\zeta)}\zeta^{\alpha},\qquad{\alpha}\in{\mathbb{N}}_{0}^{d}.$

We need the Theorem 4.1 in its full generality for the characterization of
rational approximation for the Laurent series in the last subsection.

###  4.2. Rational approximation to Laurent series

Continued fractions of univariate polynomials have rational functions as
convergents. The same holds for our three-term expansions here, provided the
recurrence matrices satisfy (3.4) and (3.14) or (3.6) and (3.16),
respectively. To that end, we use the exponential polynomial $g_{n}$ (4.5)
from which $\mu^{n}$ is generated as $\mu^{n}_{\alpha}=g_{n}(\alpha)$, and
consider the formal Laurent series

(4.8)
$\mu^{n}(z):=\sum_{\alpha\in{\mathbb{N}}_{0}^{d}}\mu_{\alpha}^{n}\,z^{-\alpha},\qquad
z\in({\mathbb{C}}\setminus 0)^{d}.$

It has already been shown by Power [15] that the Laurent series associated
with any finite rank Hankel operator induces a rational Laurent series. In
addition, one can even derive an explicit representation of this rational
function from (4.5). Recall that $G_{n,\zeta}$ are coefficients of $g_{n}(x)$
and $\Theta$ is the theta operator.

######  Proposition 4.3.

For the sequences $\mu^{n}$ given as $\mu^{n}_{\alpha}=g_{n}(\alpha)$, we have

(4.9) $\mu^{n}(z)=\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\frac{(\cdot)^{\epsilon}}{\left((\cdot)-\zeta\right)^{\epsilon}}\right)(z),\qquad\epsilon:=(1,\dots,1).$

###### Proof.

By (4.5), we have that

$\displaystyle\mu^{n}(z)$
$\displaystyle=\sum_{\alpha\in{\mathbb{N}}_{0}^{d}}\sum_{\zeta\in
Z_{n}}G_{n,\zeta}(\alpha)z^{-\alpha}\,\zeta^{\alpha}=\sum_{\alpha\in{\mathbb{N}}_{0}^{d}}\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\left(\zeta\cdot(\cdot)^{-\epsilon}\right)^{\alpha}\right)(z)$
$\displaystyle=\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\sum_{\alpha\in{\mathbb{N}}_{0}^{d}}\left(\zeta\cdot(\cdot)^{-\epsilon}\right)^{\alpha}\right)(z)=\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\frac{1}{1-\zeta\cdot(\cdot)^{-\epsilon}}\right)(z)$
$\displaystyle=\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\frac{(\cdot)^{\epsilon}}{\left((\cdot)-\zeta\right)^{\epsilon}}\right)(z),$

which is (4.9). ∎

In the case of simple zeros, or the existence of Gaussian cubature rules, we
deduce from Corollary 4.2 the following:

######  Corollary 4.4.

Under the assumption of Corollary 4.2, the rational function $\mu_{n}(z)$ in
(4.9) is given by

(4.10) $\mu^{n}(z)=\sum_{\zeta\in
Z_{n}}\frac{1}{K_{n}(\zeta,\zeta)}\frac{z^{\epsilon}}{(z-\zeta)^{\epsilon}}=\sum_{\zeta\in
Z_{n}}\frac{1}{K_{n}(\zeta,\zeta)}\frac{z_{1}\cdots
z_{d}}{(z_{1}-\zeta_{1})\cdots(z_{d}-\zeta_{d})}.$

This shows that the rational functions coming from moment sequences or finite
rank Hankel operators have some special properties: they have finitely many
poles of order $d$ and their denominator polynomials can be factored into
simple linear factors. And obviously the denominators are related to the ideal
$\left\langle{\mathbb{F}}_{n}\right\rangle$ by the fact that the denominator
is the product of the common zeros of this ideal.

Let us consider the Laurent series and its rational approximants for the
examples in Subsection 3.4 which admit the Gaussian cubature rules. Instead of
dealing with integrals over the domain $\Omega$ there, we work with integrals
over ${\mathbb{R}}^{d}$ using the mapping (3.21).

Let $w$ be a weight function defined on $[-1,1]$ and assume that it admits a
Gaussian quadrature formula. Let

${\mathcal{L}}_{-\frac{1}{2}}(f)=\frac{1}{d!}\int_{{\mathbb{R}}^{d}}f(e_{1}(x),\ldots,e_{d}(x))\prod_{i=1}^{d}w_{i}(x)\mathrm{d}x$

and

${\mathcal{L}}_{\frac{1}{2}}(f)=\frac{1}{d!}\int_{{\mathbb{R}}^{d}}f(e_{1}(x),\ldots,e_{d}(x))\prod_{1\leq
i<j\leq d}(x_{i}-x_{j})^{2}\prod_{i=1}^{d}w_{i}(x)\mathrm{d}x,$

where $e_{k}(x_{1},\ldots,x_{d})$ is the $k$-th elementary symmetric function
defined in (3.21). For the moment sequence
$\\{{\mathcal{L}}_{\pm\frac{1}{2}}(x^{\alpha})\\}$, the Laurent series (4.8)
can be easily seen to satisfy

$\mu_{-\frac{1}{2}}(z)=\frac{1}{d!}\int_{{\mathbb{R}}^{d}}\frac{z_{1}\cdots
z_{d}}{(z_{1}-e_{1}(x))\cdots(z_{d}-e_{d}(x))}\prod_{i=1}^{d}w_{i}(x)\mathrm{d}x$

and

$\mu_{\frac{1}{2}}(z)=\frac{1}{d!}\int_{{\mathbb{R}}^{d}}\frac{z_{1}\cdots
z_{d}}{(z_{1}-e_{1}(x))\cdots(z_{d}-e_{d}(x))}\prod_{1\leq i<j\leq
d}(x_{i}-x_{j})^{2}\prod_{i=1}^{d}w_{i}(x)\mathrm{d}x.$

By (3.22) and (3.23), the corresponding rational functions that approximate
$\mu_{\pm\frac{1}{2}}(z)$, as in (4.10), are given by

(4.11)
$\mu_{-\frac{1}{2}}^{n}(z)=\sum_{{\gamma}_{1}=1}^{n}\sum_{{\gamma}_{2}=1}^{{\gamma}_{1}}\cdots\sum_{{\gamma}_{d}=1}^{{\gamma}_{d-1}}\Lambda_{{\gamma},n}^{(-\frac{1}{2})}\frac{z_{1}\cdots
z_{d}}{(z_{1}-u_{{\gamma}_{1},n})\cdots(z_{d}-u_{{\gamma}_{d},n})}$

and

(4.12)
$\mu_{\frac{1}{2}}^{n}(z)=\sum_{{\gamma}_{1}=1}^{n+d-1}\sum_{{\gamma}_{2}=1}^{{\gamma}_{1}-1}\cdots\sum_{{\gamma}_{d}=1}^{{\gamma}_{d-1}-1}\Lambda_{{\gamma},n}^{(\frac{1}{2})}\frac{z_{1}\cdots
z_{d}}{(z_{1}-u_{{\gamma}_{1},n+d-1})\cdots(z_{d}-u_{{\gamma}_{d},n+d-1})},$

where $u_{{\gamma},n}$ and $u_{{\gamma},n+d-1}$ are nodes of the corresponding
Gaussian cubature rules and $\Lambda_{{\gamma},n}^{(\pm\frac{1}{2})}$ are
given by (3.24). Evidently, $\mu_{\pm\frac{1}{2}}^{n}(z)$ are rational
functions and, by Theorem 4.1, they satisfy

$\mu_{\pm\frac{1}{2}}(z)-\mu_{\pm\frac{1}{2}}^{n}(z)=O(|z|^{2n}).$

###  4.3. Continued fraction expansions of moment sequences

We can now assemble the theory of the preceding sections to give a
characterization of moment sequences that admit an associated multivariate
continued fraction expansion.

To that end, we start with a sequence
$\mu=\left(\mu_{\alpha}:\alpha\in{\mathbb{N}}_{0}\right)$ such that $\det
M_{n}\neq 0$, $n\geq 0$. Define

$\mu(z)=\sum_{{\alpha}\in{\mathbb{N}}_{0}^{d}}\mu_{\alpha}z^{-\alpha},\qquad
z\in z\in({\mathbb{C}}\setminus 0)^{d}.$

As shown in Proposition 3.1, there exists a sign-orthonormal basis
${\mathbb{P}}_{n}$, $n\geq 0$, of $\Pi^{d}$ that satisfies the recurrence
relation (3.1) and the rank condition (3.2). Moreover, the commutativity
conditions (3.4) hold for the recurrence matrices. If, in addition, the
commuting conditions (3.14) are satisfied for $n\geq 0$, then the polynomials
${\mathbb{P}}_{n}$ are H-bases and there exist sequences
$\mu^{n}=\left(\mu^{n}_{\alpha}:\alpha\in{\mathbb{N}}_{0}\right)$ such that
$\mu^{n}_{\alpha}=\mu_{\alpha}$, $|\alpha|\leq 2n-1$, and therefore

$\mu(z)-\mu^{n}(z)=O\left(z^{\alpha}:|\alpha|=2n\right),$

where $\mu^{n}(z)$ is a rational function by (4.9). This is the perfect
analogy for the continued fraction expansion of a given Laurent series.

The above procedure can be performed for _any_ vector polynomial sequence
${\mathbb{F}}_{n}\in\Pi_{n}^{r_{n}^{0}}$ given by a three-term recurrence

${\mathbb{F}}_{n+1}=A_{n}(x){\mathbb{F}}_{n}+B_{n}{\mathbb{F}}_{n}+C_{n}{\mathbb{F}}_{n-1},\qquad
A_{n}(x)=\sum_{i=1}^{d}A_{n,i}\,x_{i}$

with matrix coefficients as in the recurrence relations of orthogonal
polynomials. Moreover, we let ${\mathcal{A}}_{n}\subset{\mathbb{N}}_{0}^{d}$
be a lower set, i.e., $\alpha\in{\mathcal{A}}$ implies $\beta$ in
${\mathcal{A}}$ whenever $\beta\leq\alpha$ componentwise that has the property
that $\\{x^{\alpha}:\alpha\in{\mathcal{A}}_{n}\\}$ spans the quotient space
$\Pi^{d}/\left\langle{\mathbb{F}}_{n-1}\right\rangle$. With $Z_{n}$ denoting
the set of common zeros of ${\mathbb{F}}_{n+1}$, we then define $g_{n}$ of the
form (4.5) by requiring that

(4.13) $g_{n}(\alpha)=\mu^{n-1}_{\alpha},\qquad\alpha\in{\mathcal{A}}_{n}.$

This interpolation problem is uniquely solvable (cf. [18]) and again yields
$\mu^{n}$ by sampling $g_{n}$: $\mu^{n}_{\alpha}=g_{n}(\alpha)$,
$\alpha\in{\mathbb{N}}_{0}^{d}$. Finally, obtain by this procedure a sequence
of rational functions $\mu^{n}(z)$ whose explicit expression is given in
Proposition 4.3. Now we can give a formal definition of continued fractions
that coincides with the one from [14, 16] for $d=1$.

######  Definition 4.5.

For matrices $A_{n,i}:r_{n+1}^{0}\times r_{n}^{0}$, $i=1,\dots,d$,
$B_{n}:r_{n+1}^{0}\times r_{n}^{0}$, $C_{n}:r_{n+1}^{0}\times r_{n-1}$,
$n\in{\mathbb{N}}$, we define the _continued fraction_

(4.14)
$\frac{\left.C_{1}\right|}{\left.A_{1}(z)+B_{1}\right|}+\cdots+\frac{\left.C_{n}\right|}{\left.A_{n}(z)+B_{n}\right|}:=\sum_{\zeta\in
Z_{n}}\left(G_{n,\zeta}(-\Theta)\frac{(\cdot)^{\epsilon}}{\left((\cdot)-\zeta\right)^{\epsilon}}\right)(z),$

where $Z_{n}$ are the common zeros of ${\mathbb{F}}_{n+1}$ and the
$G_{n,\zeta}$ are determined iteratively by (4.13). The ideal
$\left\langle{\mathbb{F}}_{n+1}\right\rangle$ is called the _denominator
ideal_ of the continued fraction.

The left-hand side of (4.14) is a symbolic notation for the continued fraction
in the spirit of the classical notation of one variable. In words, with each
matrix sequence, the continued fraction expansion associates a sequence of
rational functions, which is precisely the continued fraction expansion of
rational functions, exactly like the case for $d=1$. Note, however, that the
situation is significantly more intricate in several variables since
${\mathbb{F}}_{n+1}$ may _not_ be an H-basis of the denominator ideal, and
even $Z_{n}=\emptyset$ is possible. Hence, in contrast to the univariate case,
we cannot expect these continued fractions to exist for all moment
functionals. Finally, we find it worthwhile to mention that
$\left\langle{\mathbb{F}}_{n+1}\right\rangle$ is the so-called _Prony ideal_
for the associated exponential polynomial $g_{n}$.

We are now ready to give a fundamental definition.

######  Definition 4.6.

A Laurent series $\mu(z)$ has an _associated continued fraction_ with
coefficients $A_{n,i},B_{n},C_{n}$, $n\in{\mathbb{N}}$ if

$\mu(z)-\left(\frac{\left.C_{1}\right|}{\left.A_{1}(z)+B_{1}\right|}+\cdots+\frac{\left.C_{n}\right|}{\left.A_{n}(z)+B_{n}\right|}\right)=O\left(z^{\alpha}:|\alpha|=2n\right),\qquad
n\in{\mathbb{N}}.$

In our construction of signed orthogonal polynomials above, the Laurent
ideals, hence $Z_{n}$ and ${\mathcal{Q}}_{\zeta}$, $\zeta\in Z_{n}$, depend
only on the _ideal_ $\left\langle{\mathbb{P}}_{n}\right\rangle$, hence we can
also start with a sequence of orthogonal polynomials defined by the implicit
recursion (3.1), and then refer to the explicit formula (3.8), use the
notation

$G_{n}(z)=\sum_{i=1}^{d}G_{n,i}z_{i}=\sum_{i=1}^{d}D_{n,i}^{\mathsf{T}}S_{n}z_{i},$

and consider the continued fractions

(4.15)
$\frac{-F_{1}|}{|G_{1}(z)-E_{1}}+\frac{-F_{2}|}{|G_{2}(z)-E_{2}}+\cdots\frac{-F_{n}|}{|G_{n}(z)-E_{2}},\qquad
n\in{\mathbb{N}}.$

For the continued fraction expansion defined in Definition 4.6, we can now
summarize our findings in the following way.

######  Theorem 4.7.

A Laurent series $\mu(z)=\sum\mu_{\alpha}z^{-\alpha}$ has an associated
continued fraction expansion if and only if

  1. (i)

the Hankel matrices $M_{n}$ are definite, i.e., $\det M_{n}\neq 0$, $n\geq 0$,
and

  2. (ii)

the coefficient matrices of the associated sign-orthonormal polynomials
satisfy (3.4), (3.2) and (3.14).

As an example, the rational functions $\mu_{\pm\frac{1}{2}}^{n}$ given in
(4.11) and (4.12) are the associated continued fraction expansions of the
Laurent series $\mu_{\pm\frac{1}{2}}(z)$.

Once more we want to emphasize that the critical condition is (3.14). Whenever
$\det M_{n}\neq 0$, then the rank condition (3.2) and the commuting conditions
(3.4) follow automatically, they are needed, however, to ensure that a three
term recurrence yields a well-defined multiplication.

To further understand the meaning of (3.14), we consider the _monic_
orthogonal basis $\widehat{\mathbb{P}}_{n}$ associated to the definite
sequence $\mu$, where
$\widehat{\boldsymbol{P}}_{n}=\left[\begin{array}[]{c}*\\\ I_{r_{n}^{0}\times
r_{n}^{0}}\end{array}\right]$. Here, the matrices $A_{n,i}$, $B_{n,i}$ and
$C_{n,i}$ of the recurrence relation (3.5) can be given explicitly from the
moment matrix, cf. [2], in particular

$A_{n,i}=L_{n,i}=:\sum_{|\alpha|=n}e_{\alpha}e_{\alpha+\epsilon_{i}}^{T},\qquad\text{and}\qquad
C_{n,i}=H_{n}L_{n-1,i}^{T}H_{n-1}^{-1},$

where

$H_{n}=M_{n,n}-M_{n,n-1}^{\mathsf{T}}M_{n-1}^{-1}M_{n,n-1},\qquad
M_{n}=\left[\begin{array}[]{cc}M_{n-1}&M_{n,n-1}\\\
M_{n,n-1}&M_{n,n}\end{array}\right],$

denotes the Schur complement of $M_{n-1}$ in $M_{n}$, which is nonsingular
since $M_{n}$ and $M_{n-1}$ are nonsingular. For these matrices, (3.16) (with
$n$ replaced by $n+1$ for convenience) becomes

$L_{n-1,i}H_{n}L_{n-1,j}^{T}H_{n-1}^{-1}=L_{n-1,j}H_{n}L_{n-1,i}^{T}H_{n-1}^{-1},$

or, equivalently,

$L_{n-1,i}H_{n}L_{n-1,j}^{T}=L_{n-1,j}H_{n}L_{n-1,i}^{T}.$

By [22, Lemma 4.2], this identity implies that $H_{n}$ is a Hankel matrix.
that is, $(H_{n})_{\alpha,\beta}=h_{\alpha+\beta}$. Since $M_{n,n}$ is a
Hankel matrix, we obtain the following characterization:

######  Theorem 4.8.

A sequence $\mu$ has an associated continued fraction expansion if and only if
it is definite and the matrix

$M_{n,n-1}^{\mathsf{T}}M_{n-1}^{-1}M_{n,n-1}$

is a Hankel matrix for all $n\geq 1$.

## References

  * [1] M. J. Atia, M.Benabdalla and R. S. Costas-Santos, Zeros of polynomials orthogonal with respect to a signed weight. Indag. Math. 23 (2012), 26–31.
  * [2] R. Barrio, J. M. Peña, and T. Sauer.  Three-term recurrence for the evaluation of multivariate orthogonal polynomials.  J. Approx. Theor., 162 (2010), 407–420.
  * [3] H. Berens, H. Schmid and Y. Xu, Multivariate Gaussian cubature formula. Arch. Math. 64 (1995), 26–32.
  * [4] C. Dunkl and Y. Xu, Orthogonal polynomials of several variables. Encyclopedia of Mathematics and its Applications, 155. Cambridge Univ. Press, Cambridge, 2014.
  * [5] C.F. Gauss, Methodus nova integralium valores per approximationem inveniendi.  Commentationes societate regiae scientiarum Gottingensis recentiores III (1816)
  * [6] W. Gautschi, A survey of Gauss-Christoffel quadrature formulae, in E. B. Christoffel (Aachen/Monschau, 1979), pp. 72–147, Birkhäuser, Basel-Boston, Mass., 1981.
  * [7] W. Gröbner, Über das Macaulaysche inverse System und dessen Bedeutung für die Theorie der linearen Differentialgleichungen mit konstanten Koeffizienten. Abh. Math. Sem. Hamburg 12 (1937), 127–132.
  * [8] W. Groebner, Algebraische Geometrie II. Bibliographisches Institut, 1970.
  * [9] O. Karpenkov, Geomery of Continued Fractions. Springer, 2013.
  * [10] H. Li and Y. Xu, Discrete Fourier analysis on fundamental domain of $A_{d}$ lattice and on simplex in d-variables. J. Fourier Anal. Appl. 16 (2010), 383–433.
  * [11] G. Monegato, On polynomials orthogonal with respect to particular variable-signed weight functions. Z. Angew. Math. Phys. 31 (1980), 549–555.
  * [12] B. Mourrain, A new criterion for normal form algorithms. In: Fossorier, M., Imai, H., Lin, S., Poli, A. (eds) Applied Algebra, Algebraic Algorithms and Error-Correcting Codes. AAECC 1999\. Lecture Notes in Computer Science, vol 1719, p. 430-442, Springer, 2003.
  * [13] S. Notaris, Gauss-Kronrod quadrature formulae–a survey of fifty years of research. Electron. Trans. Numer. Anal. 45 (2016), 371–404.
  * [14] O. Perron, Die Lehre von den Kettenbrüchen I/II, 3rd ed. B. G. Teubner, 1954.
  * [15] S. C. Power.  Finite rank multivariable Hankel forms.  Lin. Alg. Appl., 48 (1982), 237–244.
  * [16] T. Sauer Continued Fractions and Signal Processing. Springer, 2021.
  * [17] T. Sauer.  Gröbner bases, H–bases and interpolation.  Trans. Amer. Math. Soc., 353 (2001), 2293–2308.
  * [18] T. Sauer.  Reconstructing sparse exponential polynomials from samples: Difference operators, Stirling numbers and Hermite interpolation.  In M. Floater, T. Lyche, M.-L. Mazure, K. Moerken, and L.-L. Schumaker (eds.), Mathematical Methods for Curves and Surfaces. 9th International Conference, MMCS 2016, Tønsberg, Norway. Revised Selected Papers, Lecture Notes in Computer Science, 233–251. Springer, 2017.
  * [19] K. Schmüdgen The Moment Problem. Springer, 2017.
  * [20] F. Schweiger, Multidimensional Continued Fractions. Oxford University Press, 2000.
  * [21] G. Struble, Orthogonal polynomials: Variable-signed weight functions. Numer. Math. 5 (1963), 88–94.
  * [22] Y. Xu, On zeros of multivariate quasi-orthogonal polynomials and Gaussian cubature formulae. SIAM J. Math. Anal. 25 (1994), 991–1001.