# Connecting Exceptional Orthogonal Polynomials of Different Kind

Christiane QUESNE
Physique Nucléaire Théorique et Physique Mathématique, Université Libre de
Bruxelles,
Campus de la Plaine CP229, Boulevard du Triomphe, B-1050 Brussels, Belgium
<EMAIL_ADDRESS>
###### Abstract

The known asymptotic relations interconnecting Jacobi, Laguerre, and Hermite
classical orthogonal polynomials are generalized to the corresponding
exceptional orthogonal polynomials of codimension $m$. It is proved that
$X_{m}$-Laguerre exceptional orthogonal polynomials of type I, II, or III can
be obtained as limits of $X_{m}$-Jacobi exceptional orthogonal polynomials of
the same type. Similarly, $X_{m}$-Hermite exceptional orthogonal polynomials
of type III can be derived from $X_{m}$-Jacobi or $X_{m}$-Laguerre ones. The
quadratic transformations expressing Hermite classical orthogonal polynomials
in terms of Laguerre ones is also extended to even $X_{2m}$-Hermite
exceptional orthogonal polynomials.

Key words: exceptional orthogonal polynomials; asymptotic relations; quadratic
transformations

2020 Mathematics Subject Classification: 33C47; 34A05

##  1 Introduction

Exceptional orthogonal polynomials (EOPs) are complete families of orthogonal
polynomials that arise as eigenfunctions of a Sturm-Liouville eigenvalue
problem [6]. In contrast with the classical orthogonal polynomials (COPs) of
Jacobi, Laguerre, and Hermite types [10], there are some gaps in the sequence
of their degrees, the total number of missing “exceptional” degrees being
known as the codimension. In addition, the corresponding differential equation
contains rational coefficients instead of polynomial ones. Due to these facts,
EOPs circumvent the strong limitations of Bochner’s classification theorem,
which characterizes Sturm-Liouville COP systems.

During the last few years, a lot of research activity has been devoted to the
study of EOPs both from a mathematical viewpoint and for their applications in
mathematical physics (for a recent account see [3] and references quoted
therein).

In quantum mechanics, for instance, they have been shown to lead to some
exactly solvable rational extensions of well-known quantum potentials. The
$X_{1}$ EOPs of codimension one of Ref. [6] were indeed shown to be related to
the Darboux transformation in the context of shape invariant potentials in
supersymmetric quantum mechanics [20]. Soon after the introduction of the
first $X_{2}$ EOPs of codimension two [21], infinite families of shape
invariant potentials were introduced in relation of $X_{m}$ EOPs of arbitrary
codimension $m$ [18]. Multi-indexed families of $X_{m_{1}m_{2}\ldots m_{k}}$
EOPs, connected with multi-step Darboux transformations, were then considered
[7, 19]. Very recently, confluent Darboux transformations have been used to
generate EOP families with an arbitrary number of real parameters [4].

From a mathematical viewpoint, three main questions have been the subject of
research activity in relation to EOPs. The first one is the study of the
interlacing and asymptotic behaviour of their zeros (see, e.g., [8, 11]). The
second one has to do with the recurrence relations they satisfy and which may
be of two different types: either recurrence relations of order $2N+1$, where
$N$ is the number of Darboux steps, with coefficients that are functions of
$x$ and $n$ [5, 16], or recurrence relations of order $2m+3$, where $m$ is the
codimension, with coefficients only dependent on $n$ [2, 15, 17]. Finally, the
quest for a complete classification of EOPs is a fundamental problem [3],
which, as shown in [4], remains rather open.

The purpose of the present paper is to look into the counterpart of a well-
known property of Jacobi, Laguerre, and Hermite COPs, namely that they are
interconnected through limit relations [10]. As a first step in this inquiry,
we plan to consider here $X_{m}$-Jacobi, $X_{m}$-Laguerre, and $X_{m}$-Hermite
EOPs, the first two existing in three different types I, II, and III, while
the last ones only belong to type III.

In Section 2, we list the sets of EOPs that we are going to consider, we
provide their explicit expressions in terms of COPs, and we explain how the
latter are related to some other definitions found in the literature. In
Section 3, we present and prove limit relations connecting these EOPs. In
Section 4, we briefly comment on another type of relations generalizing the
quadratic transformations relating Laguerre and Hermite COPs. Finally, Section
5 contains the conclusion.

##  2 Sets of $X_{m}$ EOPs and their expressions in terms of COPs

The sets of $X_{m}$ EOPs that we are going to consider here mainly come from
Ref. [15] (except for type II $X_{m}$-Laguerre EOPs, which are not considered
there), where they are given in monic form. To establish relations with other
definitions where standard EOPs and COPs are used, let us remind that one goes
from standard Hermite, Laguerre, and Jacobi COPs to monic ones by the
following transformations

$\displaystyle H_{n}(x)$ $\displaystyle\to$ $\displaystyle 2^{n}H_{n}(x),$
$\displaystyle L_{n}^{(\alpha)}(x)$ $\displaystyle\to$
$\displaystyle\frac{(-1)^{n}}{n!}L_{n}^{(\alpha)}(x),$ (2.1) $\displaystyle
P_{n}^{(\alpha,\beta)}(x)$ $\displaystyle\to$
$\displaystyle\frac{\Gamma(2n+\alpha+\beta+1)}{2^{n}n!\Gamma(n+\alpha+\beta+1)}P_{n}^{(\alpha,\beta)}(x).$

For future use, we list some relations satisfied by monic COPs in Appendix A.

###  2.1 $X_{m}$-Hermite polynomials

It is well known [5, 14] that only type III $X_{m}$-Hermite EOPs exist and
that they are restricted to even $m$ values. They may be defined as

$\displaystyle\hat{H}^{({\rm III},m)}_{n}(x)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\ {\rm
i}^{m}\left[H_{m}({\rm i}x)H_{n-m}(x)+{\rm i}\frac{m}{2}H_{m-1}({\rm
i}x)H_{n-m-1}(x)\right]&\text{if $n=m+1,m+2.\ldots$},\end{cases}$ (2.2)

with $m$ missing degrees $n=1,2,\ldots,m$. This relation comes from [15],
where

$\displaystyle\hat{H}^{({\rm III},m)}_{n}(x)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\ -\frac{1}{2}{\rm
i}^{m}\left\\{H_{m}({\rm
i}x)(\partial_{x}-2x)H_{n-m-1}(x)-[\partial_{x}H_{m}({\rm
i}x)]H_{n-m-1}(x)\right\\}&\text{if $n=m+1,m+2,\ldots$},\end{cases}$

after applying (A.1) and (A.2). It is connected to the previously used
$y^{(m)}_{n}(x)$ of Ref. [14] by the relation

$y^{(m)}_{n}(x)=\begin{cases}\hat{H}^{({\rm III},m)}_{0}(x)&\text{if
$n=0$},\\\ -2^{n}\hat{H}^{({\rm III},m)}_{n}(x)&\text{if
$n=m+1,m+2,\ldots$}.\end{cases}$

###  2.2 $X_{m}$-Laguerre polynomials

$X_{m}$-Laguerre EOPS exist in three different types, which may be defined as

$\hat{L}^{({\rm
I},m)}_{n}(x;\alpha)=(-1)^{m}\left(L_{m}^{(\alpha+1)}(-x)L^{(\alpha)}_{n}(x)-nL_{m}^{(\alpha)}(-x)L^{(\alpha+1)}_{n-1}(x)\right),\qquad
n=0,1,2,\ldots,$ (2.3) $\displaystyle\hat{L}^{({\rm II},m)}_{n}(x;\alpha)$
$\displaystyle=\frac{1}{\alpha+n-m}\left(nxL_{m}^{(-\alpha)}(x)L_{n-1}^{(\alpha+1)}(x)-(m-\alpha)L_{m}^{(-\alpha-1)}(x)L_{n}^{(\alpha)}(x)\right),$
$\displaystyle\quad n=0,1,2,\ldots,$ (2.4)

and

$\displaystyle\hat{L}^{({\rm III},m)}_{n}(x;\alpha)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\
(-1)^{m}\left(L_{m}^{(-\alpha)}(-x)L_{n-m}^{(\alpha-1)}(x)-mxL_{m-1}^{(-\alpha+1)}(-x)L_{n-m-1}^{(\alpha)}(x)\right)&\text{if
$n=m+1,$}\\\ &\quad\text{$m+2,\ldots$}.\end{cases}$ (2.5)

For the first two types, the polynomial degrees are $n+m=m,m+1,m+2,\ldots$
(with degrees 0, 1, …, $m-1$ missing), while for the third type, they are
$n=0,m+1,m+2,\ldots$ (with degrees 1, 2, …, $m$ missing).

Equation (2.3) is derived from the corresponding equation of Ref. [15], which
reads

$\displaystyle\hat{L}^{({\rm I},m)}_{n}(x;\alpha)$
$\displaystyle=(-1)^{m+1}\left\\{L_{m}^{(\alpha)}(-x)(\partial_{x}-1)L_{n}^{(\alpha)}(x)-\left[\partial_{x}L_{m}^{(\alpha)}(-x)\right]L_{n}^{(\alpha)}(x)\right\\},$
$\displaystyle\quad n=0,1,2,\ldots,$

after using (A.3) and (A.4). It is directly related to the definition of type
I Laguerre EOPs $L^{{\rm I},\alpha}_{m,n}(x)$, given in Eq. (3.2) of [12] in
terms of standard Laguerre polynomials (see also [8]). After transforming the
latter into monic ones, one indeed gets

$L^{{\rm I},\alpha}_{m,n}(x)=\frac{(-1)^{n-m}}{m!(n-m)!}\hat{L}^{({\rm
I},m)}_{n-m}(x;\alpha-1),\qquad n=m,m+1,m+2,\ldots.$

To obtain Eq. (2.4), we started from the definition of type II Laguerre EOPs
$L^{{\rm II},\alpha}_{m,n}(x)$ given in [12] (see also [8]), converted
standard Laguerre COPs into monic ones and multiplied the result by an
appropriate factor to get a monic EOP. As a result, we got

$L^{{\rm
II},\alpha}_{m,n}(x)=(-1)^{n-1}\frac{\alpha+1+n-2m}{m!(n-m)!}\hat{L}^{({\rm
II},m}_{n-m}(x;\alpha+1),\qquad n=m,m+1,m+2,\ldots.$

Finally, Eq. (2.5) is derived from the corresponding equation of Ref. [15],
namely

$\displaystyle\hat{L}^{({\rm III},m)}_{n}(x;\alpha)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\
(-1)^{m+1}\Bigl{\\{}L_{m}^{(-\alpha)}(-x)(x\partial_{x}+\alpha-x)L_{n-m-1}^{(\alpha)}(x)&\\\
\quad-x\left[\partial_{x}L_{m}^{(-\alpha)}(-x)\right]L_{n-m-1}^{(\alpha)}(x)\Bigr{\\}}&\text{if
$n=m+1,m+2,\ldots$},\end{cases}$ (2.6)

by using (A.3) and (A.5). After some transformations, it can be related to the
definition of type III Laguerre EOPs $L^{{\rm III},\alpha}_{m,n}(x)$ given in
Eq. (5.12) of [12]. To show this, let us start from (2.6) and use (A.5),
thereby getting

$\displaystyle\hat{L}^{({\rm III},m)}_{n}(x;\alpha)$
$\displaystyle=(-1)^{m}\Bigl{\\{}L_{m}^{(-\alpha)}(-x)\Bigl{[}L_{n-m}^{(\alpha-1)}(x)+(\alpha-x)L_{n-m-1}^{(\alpha)}(x)\Bigr{]}$
$\displaystyle\quad-
L_{m+1}^{(-\alpha-1)}(-x)L_{n-m-1}^{(\alpha)}(x)\Bigr{]}\Bigr{\\}},\qquad
n=m+1,m+2,\ldots.$

On the right-hand side of this equation, the factor between square brackets
becomes

$\displaystyle L_{n-m}^{(\alpha-1)}(x)+(\alpha-x)L_{n-m-1}^{(\alpha)}(x)$
$\displaystyle\quad=L_{n-m}^{(\alpha)}(x)+(n-m+\alpha-x)L_{n-m-1}^{(\alpha)}(x)$
$\displaystyle\quad=-(n-m-1)\Bigl{[}L_{n-m-1}^{(\alpha)}(x)+(n-m+\alpha-1)L_{n-m-2}^{(\alpha)}(x)\bigr{]}$
$\displaystyle\quad=-(n-m-1)xL_{n-m-2}^{(\alpha+1)}(x),$

after successively using (A.4), (A.6), and (A.7). The result reads

$\displaystyle\hat{L}^{({\rm III},m)}_{n}(x;\alpha)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\
(-1)^{m+1}\Bigl{[}(n-m-1)xL_{m}^{(-\alpha)}(-x)L_{n-m-2}^{(\alpha+1)}(x)&\\\
\quad+L_{m+1}^{(-\alpha-1)}(-x)L_{n-m-1}^{(\alpha)}(x)\Bigr{]}&\text{if
$n=m+1,m+2,\ldots$},\end{cases}$ (2.7)

which may then be compared with Eq. (5.12) of [12] after rewriting there
standard Laguerre polynomials in terms of monic ones. Hence,

$L^{{\rm III},\alpha}_{m,n}(x)=\frac{(-1)^{n-m-1}}{m!(n-m-1)!}\hat{L}^{({\rm
III},m)}_{n}(x,\alpha+1),\qquad n=m+1,m+2,\ldots,$

and of course $L^{{\rm III},\alpha}_{m,0}(x)=\hat{L}^{({\rm
III},m)}_{0}(x;\alpha+1)=1$.

It is worth observing that in [12], for studying properties of EOPs, the range
of the $\alpha$ parameter has been restricted to some interval ($\alpha>0$ for
type I, $\alpha>m-1$ for type II, and $0<\alpha<1$ for type III). Here, as in
[15], we consider their formal definitions without restriction on $\alpha$.

###  2.3 $X_{m}$-Jacobi polynomials

$X_{m}$-Jacobi EOPs exist in three different types, which may be defined as

$\displaystyle\hat{P}^{({\rm I},m)}_{n}(x;\alpha,\beta)$
$\displaystyle=\frac{1}{\beta+n-m}\Bigl{\\{}P_{m}^{(\alpha,-\beta)}(x)\Bigl{[}n(1+x)P_{n-1}^{(\alpha+1,\beta+1)}(x)+\beta
P_{n}^{(\alpha,\beta)}(x)\Bigr{]}$
$\displaystyle-m(1+x)P_{m-1}^{(\alpha+1,-\beta+1)}(x)P_{n}^{(\alpha,\beta)}(x)\Bigr{\\}},\qquad
n=0,1,2,\ldots,$ (2.8) $\displaystyle\hat{P}^{({\rm
II},m)}_{n}(x;\alpha,\beta)$
$\displaystyle=\frac{1}{m-n-\alpha}\Bigl{\\{}P_{m}^{(-\alpha,\beta)}(x)\Bigl{[}n(1-x)P_{n-1}^{(\alpha+1,\beta+1)}(x)-\alpha
P_{n}^{(\alpha,\beta)}(x)\Bigr{]}$
$\displaystyle-m(1-x)P_{m-1}^{(-\alpha+1,\beta+1)}(x)P_{n}^{(\alpha,\beta)}(x)\Bigr{\\}},\qquad
n=0,1,2,\ldots,$ (2.9)

and

$\displaystyle\hat{P}^{({\rm III},m)}_{n}(x;\alpha,\beta)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\
\frac{1}{\alpha+\beta+n-2m-1}\Bigl{[}(\alpha+\beta+n-m-1)P_{m}^{(-\alpha,-\beta)}(x)P_{n-m}^{(\alpha-1,\beta-1)}(x)&\\\
\quad+m(1-x^{2})P_{m-1}^{(-\alpha+1,-\beta+1)}(x)P_{n-m-1}^{(\alpha,\beta)}(x)\Bigr{]},&\text{if
$n=m+1$},\\\ &\quad\text{$m+2,\ldots$}.\\\ \end{cases}$ (2.10)

As in the case of $X_{m}$-Laguerre EOPs, $X_{m}$-Jacobi EOPs are of degree
$n+m=m,m+1,m+2,\ldots$ (with degrees 0, 1, …, $m-1$ missing) in the case of
type I or II, while they are of degree $n=0,m+1,m+2,\ldots$ (with degrees 1,
2, …, $m$ missing) in that of type III.

These relations are derived from the corresponding definitions of [15],
namely111Note that in [15], Jacobi polynomials are denoted by
$J_{n}^{(\alpha,\beta)}(x)$ instead of $P_{n}^{(\alpha,\beta)}(x)$.

$\displaystyle\hat{P}^{({\rm I},m)}_{n}(x;\alpha,\beta)$
$\displaystyle=\frac{1}{\beta+n-m}\Bigl{\\{}P_{m}^{(\alpha,-\beta)}(x)[(1+x)\partial_{x}+\beta]P_{n}^{(\alpha,\beta)}(x)$
$\displaystyle-(1+x)\bigl{[}\partial_{x}P_{m}^{(\alpha,-\beta)}(x)\bigr{]}P_{n}^{(\alpha,\beta)}(x)\Bigr{\\}},\qquad
n=0,1,2,\ldots,$ $\displaystyle\hat{P}^{({\rm II},m)}_{n}(x;\alpha,\beta)$
$\displaystyle=\frac{1}{m-n-\alpha}\Bigl{\\{}P_{m}^{(-\alpha,\beta)}(x)[(1-x)\partial_{x}-\alpha]P_{n}^{(\alpha,\beta)}(x)$
$\displaystyle-(1-x)\bigl{[}\partial_{x}P_{m}^{(-\alpha,\beta)}(x)\bigr{]}P_{n}^{(\alpha,\beta)}(x)\Bigr{\\}},\qquad
n=0,1,2,\ldots,$

and

$\displaystyle\hat{P}^{({\rm III},m)}_{n}(x;\alpha,\beta)$
$\displaystyle\quad=\begin{cases}1&\text{if $n=0$},\\\
\frac{1}{\alpha+\beta+n-2m-1}\Bigl{\\{}P_{m}^{(-\alpha,-\beta)}(x)[(x^{2}-1)\partial_{x}+(\alpha+\beta)x+\alpha-\beta]P_{n-m-1}^{(\alpha,\beta)}(x)&\\\
\quad+(1-x^{2})\bigl{[}\partial_{x}P_{m}^{(-\alpha,-\beta)}(x)\bigr{]}P_{n-m-1}^{(\alpha,\beta)}(x)\Bigr{\\}},&\text{if
$n=m+1$},\\\ &\quad\text{$m+2,\ldots$},\\\ \end{cases}$

after using (A.8) and (A.9).

It is worth observing that types I and II $X_{m}$-Jacobi EOPs are not
independent since

$\hat{P}^{({\rm II},m)}_{n}(x;\alpha,\beta)=(-1)^{n+m}\hat{P}^{({\rm
I},m)}_{n}(-x;\beta,\alpha)$

as a direct consequence of (2.8), (2.9), and the well-known symmetry relation
$P^{(\alpha,\beta)}_{n}(x)=(-1)^{n}P^{(\beta,\alpha)}_{n}(-x)$ of Jacobi
polynomials.

As we now plan to show, it is actually type II $X_{m}$-Jacobi EOPs (2.9) that
are considered in some mathematical studies [1, 8, 13]. Specifically, we are
going to prove that $\hat{P}^{(\alpha,\beta)}_{m,n}(x)$, $n=m,m+1,m+2,\ldots$,
defined in Eq. (64) of [8] is such that

$\displaystyle\hat{P}^{(\alpha,\beta)}_{m,n}(x)$
$\displaystyle=(-1)^{m}\frac{\alpha+1+n-2m}{\alpha+1+n-m}\frac{\Gamma(2m-\alpha+\beta-1)\Gamma(2n-2m+\alpha+\beta+1)}{2^{n}m!(n-m)!\Gamma(m-\alpha+\beta-1)\Gamma(n-m+\alpha+\beta+1)}$
$\displaystyle\quad\times\hat{P}^{({\rm
II},m)}_{n-m}(x;\alpha+1,\beta-1),\qquad n=m,m+1,m+2,\ldots.$ (2.11)

Equation (64) of [8], rewritten in terms of monic Jacobi polynomials, is given
by

$\displaystyle\hat{P}^{(\alpha,\beta)}_{m,n}(x)=\frac{(-1)^{m}}{\alpha+1+n-m}\frac{\Gamma(2m-\alpha+\beta-1)\Gamma(2n-2m+\alpha+\beta+1)}{2^{n}m!(n-m)!\Gamma(m-\alpha+\beta-1)\Gamma(n-m+\alpha+\beta+1)}$
$\displaystyle\quad\times\Bigl{[}(n-m)(x-1)P^{(-\alpha-1,\beta-1)}_{m}(x)P^{(\alpha+2,\beta)}_{n-m-1}(x)+(\alpha+1-m)P^{(-\alpha-2,\beta)}_{m}(x)P^{(\alpha+1,\beta-1)}_{n-m}(x)\Bigr{]},$
$\displaystyle\qquad n=m,m+1,m+2,\ldots,$

while from (2.9), we obtain

$\displaystyle\hat{P}^{({\rm
II},m)}_{n-m}(x;\alpha+1,\beta-1)=\frac{1}{\alpha+n+1-2m}\Bigl{\\{}(n-m)(x-1)P^{(-\alpha-1,\beta-1)}_{m}(x)P^{(\alpha+2,\beta)}_{n-m-1}(x)$
$\displaystyle\quad{}+\Bigl{[}(\alpha+1)P^{(-\alpha-1,\beta-1)}_{m}(x)+m(1-x)P^{(-\alpha,\beta)}_{m-1}(x)\Bigr{]}P^{(\alpha+1,\beta-1)}_{n-m}(x)\Bigr{\\}},$
$\displaystyle\qquad n=m,m+1,m+2,\ldots.$

Hence, Eq. (2.11) is equivalent to the following identity satisfied by Jacobi
polynomials

$(\alpha+1)P^{(-\alpha-1,\beta-1)}_{m}(x)+m(1-x)P^{(-\alpha,\beta)}_{m-1}(x)=(\alpha+1-m)P^{(-\alpha-2,\beta}_{m}(x),$
(2.12)

which is proved in Appendix B.

##  3 Limit relations satisfied by $X_{m}$ EOPs

The purpose of this Section is to extend to $X_{m}$ EOPs the well-known limit
relations satisfied by Jacobi, Laguerre, and Hermite COPs [10], which in monic
form can be written as

$\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}=(-2)^{n}L^{(\alpha)}_{n}(x),$
(3.1)
$\lim_{\alpha\to\infty}\alpha^{n/2}P^{(\alpha,\alpha)}_{n}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}=H_{n}(x),$
(3.2)

and

$\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}L^{(\alpha)}_{n}\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}=H_{n}(x).$
(3.3)

###  3.1 Going from Jacobi to Laguerre EOPs

Equation (3.1) can be generalized to $X_{m}$-Jacobi and $X_{m}$-Laguerre EOPs
as follows:

$\lim_{\beta\to\infty}\beta^{n+m}\hat{P}^{({\rm
I},m)}_{n}\left(1-\frac{2x}{\beta};\alpha,\beta\right)=(-2)^{n+m}\hat{L}^{({\rm
I},m)}_{n}(x;\alpha),\qquad n=0,1,2,\ldots,$ (3.4)
$\lim_{\beta\to\infty}\beta^{n+m}\hat{P}^{({\rm
II},m)}_{n}\left(1-\frac{2x}{\beta};\alpha,\beta\right)=(-2)^{n+m}\hat{L}^{({\rm
II},m)}_{n}(x;\alpha),\qquad n=0,1,2,\ldots,$ (3.5)

and

$\lim_{\beta\to\infty}\beta^{n}\hat{P}^{({\rm
III},m)}_{n}\left(1-\frac{2x}{\beta};\alpha,\beta\right)=(-2)^{n}\hat{L}^{({\rm
III},m)}_{n}(x;\alpha),\qquad n=m+1,m+2,m+3,\ldots.$ (3.6)

Let us first consider Eq. (3.4). By using (2.8), its left-hand side can be
rewritten as

$\displaystyle\lim_{\beta\to\infty}\beta^{n+m}\hat{P}^{({\rm
I},m)}_{n}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}$
$\displaystyle\quad=\lim_{\beta\to\infty}\frac{\beta^{n+m}}{\beta+n-m}\Bigl{\\{}P^{(\alpha,-\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigl{[}2n\Bigl{(}1-\frac{x}{\beta}\Bigr{)}P^{(\alpha+1,\beta+1)}_{n-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}$
$\displaystyle\qquad{}+\beta
P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}-2m\Bigl{(}1-\frac{x}{\beta}\Bigr{)}P^{(\alpha+1,-\beta+1)}_{m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{\\}}$
$\displaystyle\quad=\lim_{\beta\to\infty}\beta^{n+m-1}\Bigl{\\{}P^{(\alpha,-\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigl{[}2nP^{(\alpha+1,\beta+1)}_{n-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}+\beta
P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad{}-2mP^{(\alpha+1,-\beta+1)}_{m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{\\}}$
$\displaystyle\quad=\Bigl{[}\lim_{\beta\to\infty}\beta^{m}P^{(\alpha,-\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}\Bigl{[}2n\lim_{\beta\to\infty}\beta^{n-1}P^{(\alpha+1,\beta+1)}_{n-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}$
$\displaystyle\qquad{}+\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}-2m\Bigl{[}\lim_{\beta\to\infty}\beta^{m-1}P^{(\alpha+1,-\beta+1)}_{m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad\times\Bigl{[}\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}.$

On employing now (3.1), as well as its corollary

$\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,-\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}=2^{n}L^{(\alpha)}_{n}(-x),$
(3.7)

we get

$\displaystyle\lim_{\beta\to\infty}\beta^{n+m}\hat{P}^{({\rm
I},m)}_{n}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}$
$\displaystyle\quad=2^{m}L^{(\alpha)}_{m}(-x)\Bigl{[}2n(-2)^{n-1}L^{(\alpha+1)}_{n-1}(x)+(-2)^{n}L^{(\alpha)}_{n}(x)\Bigr{]}$
$\displaystyle\qquad{}-2m2^{m-1}L^{(\alpha+1)}_{m-1}(-x)(-2)^{n}L^{(\alpha)}_{n}(x)$
$\displaystyle\quad=(-2)^{n+m}(-1)^{m}\Bigl{\\{}\Bigl{[}L^{(\alpha)}_{m}(-x)-mL^{(\alpha+1)}_{m-1}(-x)\Bigr{]}L^{(\alpha)}_{n}(x)-nL^{(\alpha)}_{m}(-x)L^{(\alpha+1)}_{n-1}(x)\Bigr{\\}}$
$\displaystyle\quad=(-2)^{n+m}(-1)^{m}\Bigl{[}L^{(\alpha+1)}_{m}(-x)L^{(\alpha)}_{n}(x)-nL^{(\alpha)}_{m}(-x)L^{(\alpha+1)}_{n-1}(x)\Bigr{]},$

where in the last step use is made of (A.4). Comparison with (2.3) completes
the proof of (3.4).

Next, on inserting (2.9) in the left-hand side of (3.5) and proceeding as in
the previous case, we arrive at

$\displaystyle\lim_{\beta\to\infty}\beta^{n+m}\hat{P}^{({\rm
II},m)}_{n}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}$
$\displaystyle\quad=\frac{1}{m-n-\alpha}\Bigl{\\{}2nx\Bigl{[}\lim_{\beta\to\infty}\beta^{m}P^{(-\alpha,\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\bigr{)}\Bigl{]}\Bigl{[}\lim_{\beta\to\infty}\beta^{n-1}P^{(\alpha+1,\beta+1)}_{n-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad{}-\alpha\Bigl{[}\lim_{\beta\to\infty}\beta^{m}P^{(-\alpha,\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\bigr{)}\Bigr{]}\Bigl{[}\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad{}-2mx\Bigl{[}\lim_{\beta\to\infty}\beta^{m-1}P^{(-\alpha+1,\beta+1)}_{m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}\Bigl{[}\lim_{\beta\to\infty}\beta^{n}P^{(\alpha,\beta)}_{n}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}\Bigr{\\}}$
$\displaystyle\quad=\frac{(-2)^{n+m}}{m-n-\alpha}\Bigl{[}-nxL^{(-\alpha)}_{m}(x)L^{(\alpha+1)}_{n-1}(x)-\alpha
L^{(-\alpha)}_{m}(x)L^{(\alpha)}_{n}(x)+mxL^{(-\alpha+1)}_{m-1}(x)L^{(\alpha)}_{n}(x)\Bigr{]}$
$\displaystyle\quad=\frac{(-2)^{n+m}}{\alpha+n-m}\Bigl{\\{}nxL^{(-\alpha)}_{m}(x)L^{(\alpha+1)}_{n-1}(x)+\Bigl{[}\alpha
L^{(-\alpha)}_{m}(x)-mxL^{(-\alpha+1)}_{m-1}(x)\Bigr{]}L^{(\alpha)}_{n}(x)\Bigr{\\}}$

on using (3.1). On comparing with (2.4), to prove (3.5) it remains to check
that

$mxL^{(-\alpha+1)}_{m-1}(x)-\alpha
L^{(-\alpha)}_{m}(x)=(m-\alpha)L^{(-\alpha-1)}_{m}(x).$

From (A.4), it results that such a relation is equivalent to

$mxL^{(-\alpha+1)}_{m-1}(x)-\alpha
L^{(-\alpha)}_{m}(x)=(m-\alpha)\Bigl{[}L^{(-\alpha)}_{m}(x)+mL^{(-\alpha)}_{m-1}(x)\Bigr{]}$

or

$xL^{(-\alpha+1)}_{m-1}(x)=(m-\alpha)L^{(-\alpha)}_{m-1}(x)+L^{(-\alpha)}_{m}(x),$

which is indeed satisfied owing to (A.7). This completes the proof of (3.5).

Finally, on inserting (2.10) in the left-hand side of (3.6), we obtain that
for $n=m+1,m+2,\ldots$

$\displaystyle\lim_{\beta\to\infty}\beta^{n}\hat{P}^{({\rm
III},m)}_{n}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}$
$\displaystyle\quad=\Bigl{[}\lim_{\beta\to\infty}\beta^{m}P^{(-\alpha,-\beta)}_{m}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}\Bigl{[}\lim_{\beta\to\infty}\beta^{n-m}P^{(\alpha-1,\beta-1)}_{n-m}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad{}+4mx\Bigl{[}\lim_{\beta\to\infty}\beta^{m-1}P^{(-\alpha+1,-\beta+1)}_{m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}\Bigl{[}\lim_{\beta\to\infty}\beta^{n-m-1}P^{(\alpha,\beta)}_{n-m-1}\Bigl{(}1-\frac{2x}{\beta}\Bigr{)}\Bigr{]}.$

On using (3.1) and (3.7), we obtain

$\displaystyle\lim_{\beta\to\infty}\beta^{n}\hat{P}^{({\rm
III},m)}_{n}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}$
$\displaystyle\quad=(-1)^{n-m}2^{n}\Bigl{[}L^{(-\alpha)}_{m}(-x)L^{(\alpha-1)}_{n-m}(x)-mxL^{(-\alpha+1)}_{m-1}(-x)L^{(\alpha)}_{n-m-1}(x)\Bigr{]}.$

A comparison with (2.5) then shows that Eq. (3.6) is satisfied, which
completes the proof.

###  3.2 Going from Jacobi or Laguerre to Hermite EOPs

For type III EOPs and even $m$ values, Eqs. (3.2) and (3.3) can be generalized
as follows:

$\lim_{\alpha\to\infty}\alpha^{n/2}\hat{P}^{({\rm
III},m)}_{n}\Bigl{(}\frac{x}{\sqrt{\alpha}};\alpha,\alpha\Bigr{)}=\hat{H}^{({\rm
III},m)}_{n}(x),\qquad n=m+1,m+2,\ldots,$ (3.8)

and

$\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}\hat{L}^{({\rm
III},m)}_{n}\bigl{(}\sqrt{2\alpha}x+\alpha;\alpha\bigr{)}=\hat{H}^{({\rm
III},m)}_{n}(x),\qquad n=m+1,m+2,\ldots.$ (3.9)

Let us start with the proof of (3.8). From (2.10), we obtain that for
$n=m+1,m+2,\ldots$,

$\displaystyle\lim_{\alpha\to\infty}\alpha^{n/2}\hat{P}^{({\rm
III},m)}_{n}\Bigl{(}\frac{x}{\sqrt{\alpha}};\alpha,\alpha\Bigr{)}$
$\displaystyle\quad=\lim_{\alpha\to\infty}\frac{\alpha^{n/2}}{2\alpha+n-2m-1}\Bigl{[}(2\alpha+n-m-1)P^{(-\alpha,-\alpha)}_{m}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}P^{(\alpha-1,\alpha-1)}_{n-m}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}$
$\displaystyle\qquad{}+m\Bigl{(}1-\frac{x^{2}}{\alpha^{2}}\Bigr{)}P^{(-\alpha+1,-\alpha+1)}_{m-1}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}P^{(\alpha,\alpha)}_{n-m-1}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad=\Bigl{[}\lim_{\alpha\to\infty}\alpha^{m/2}P^{(-\alpha,-\alpha)}_{m}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}\Bigr{]}\Bigl{[}\lim_{\alpha\to\infty}\alpha^{(n-m)/2}P^{(\alpha-1,\alpha-1)}_{n-m}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}\Bigr{]}$
$\displaystyle\qquad{}+\frac{m}{2}\Bigl{[}\lim_{\alpha\to\infty}\alpha^{(m-1)/2}P^{(-\alpha+1,-\alpha+1)}_{m-1}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}\Bigr{]}\Bigl{[}\lim_{\alpha\to\infty}\alpha^{(n-m-1)/2}P^{(\alpha,\alpha)}_{n-m-1}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}\Bigr{]}$
$\displaystyle\quad=(-{\rm i})^{m}H_{m}({\rm i}x)H_{n-m}(x)+\frac{m}{2}(-{\rm
i})^{m-1}H_{m-1}({\rm i}x)H_{n-m-1}(x)$ $\displaystyle\quad=\hat{H}^{({\rm
III},m)}_{n}(x),$

where Eq. (3.2) and its corollary

$\lim_{\alpha\to\infty}(-\alpha)^{n/2}P^{(-\alpha,-\alpha)}_{n}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}=H_{n}({\rm
i}x)$

or

$\lim_{\alpha\to\infty}\alpha^{n/2}P^{(-\alpha,-\alpha)}_{n}\Bigl{(}\frac{x}{\sqrt{\alpha}}\Bigr{)}=(-{\rm
i})^{n}H_{n}({\rm i}x)$

are used, as well as definition (2.2).

Let us now consider the proof of (3.9). This time, we start from Eq. (2.5) and
obtain for $n=m+1,m+2,\ldots$,

$\displaystyle\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}\hat{L}^{({\rm
III},m)}_{n}\bigl{(}\sqrt{2\alpha}x+\alpha;\alpha\bigr{)}$
$\displaystyle\quad=\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}\Bigl{[}L^{(-\alpha)}_{m}\bigl{(}-\sqrt{2\alpha}x-\alpha\bigr{)}L^{(\alpha-1)}_{n-m}\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}$
$\displaystyle\qquad{}-m\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}L^{(-\alpha+1)}_{m-1}\bigl{(}-\sqrt{2\alpha}x-\alpha\bigr{)}L^{(\alpha)}_{n-m-1}\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}\Bigr{]}$
$\displaystyle\quad=\Bigl{[}\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{m/2}}L^{(-\alpha)}_{m}\bigl{(}-\sqrt{2\alpha}x-\alpha\bigr{)}\Bigr{]}\Bigl{[}\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{(n-m)/2}}L^{(\alpha-1)}_{n-m}\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}\Bigr{]}$
$\displaystyle\qquad{}-\frac{m}{2}\Bigl{[}\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{(m-1)/2}}L^{(-\alpha+1)}_{m-1}\bigl{(}-\sqrt{2\alpha}x-\alpha\bigr{)}\Bigr{]}\Bigl{[}\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{(n-m-1)/2}}L^{(\alpha)}_{n-m-1}\bigl{(}\sqrt{2\alpha}x+\alpha\bigr{)}\Bigr{]}.$

Here, we employ (3.3), as well as its corollary

$\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}L^{(-\alpha)}_{n}\bigl{(}-\sqrt{2\alpha}x-\alpha\bigr{)}={\rm
i}^{n}H_{n}({\rm i}x),$

thus leading to

$\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{n/2}}\hat{L}^{({\rm
III},m)}_{n}\bigl{(}\sqrt{2\alpha}x+\alpha;\alpha\bigr{)}={\rm
i}^{m}H_{m}({\rm i}x)H_{n-m}(x)-\frac{m}{2}{\rm i}^{m-1}H_{m-1}({\rm
i}x)H_{n-m-1}(x),$

which amounts to $\hat{H}^{({\rm III},m)}_{n}(x)$, as given in (2.2).

In Appendix C, we provide some examples of the limit relations proved in this
Section.

##  4 Quadratic transformations relating Hermite to Laguerre EOPs

The purpose of this Section is to examine whether the quadratic
transformations relating Hermite to Laguerre COPs [10], which in monic form
can be written as

$H_{2n}(x)=L^{(-1/2)}_{n}(x^{2}),\qquad n=0,1,2,\ldots,$ (4.1)

and

$H_{2n+1}(x)=xL^{(1/2)}_{n}(x^{2}),\qquad n=0,1,2,\ldots,$ (4.2)

can be generalized to type III EOPs.

We wish to prove that the answer to this question is positive for even
$X_{2m}$-Hermite EOPs and that in such a case the relation reads

$\hat{H}^{({\rm III},2m)}_{2n}(x)=\hat{L}^{({\rm
III},m)}_{n}\Bigl{(}x^{2};\frac{1}{2}\Bigr{)},\qquad n=0,m+1,m+2,\ldots.$
(4.3)

To start with, let us note that if we formally set $m=0$ and take Eqs. (C.1)
and (C.2) into account, Eq. (4.3) reduces to the known relation (4.1). It is
also obvious that Eq. (4.3) is fulfilled for $m=1,2,3,\ldots$ and $n=0$.

Let us therefore consider $m=1,2,3,\ldots$ and $n=m+1,m+2,\ldots$. Then from
(2.2) it follows that

$\hat{H}^{({\rm III},2m)}_{2n}(x)=(-1)^{m}[H_{2m}({\rm i}x)H_{2n-2m}(x)+{\rm
i}mH_{2m-1}({\rm i}x)H_{2n-2m-1}(x)].$

Next, Eqs. (4.1) and (4.2) imply that

$\hat{H}^{({\rm
III},2m)}_{2n}(x)=(-1)^{m}\Bigl{[}L^{(-1/2)}_{m}(-x^{2})L^{(-1/2)}_{n-m}(x^{2})-mx^{2}L^{(1/2)}_{m-1}(-x^{2})L^{(1/2)}_{n-m-1}(x^{2})\Bigr{]}.$

On comparing with

$\hat{L}^{({\rm
III},m)}_{n}(x^{2};\alpha)=(-1)^{m}\Bigl{[}L^{(-\alpha)}_{m}(-x^{2})L^{(\alpha-1)}_{n-m}(x^{2})-mx^{2}L^{(-\alpha+1)}_{m-1}(-x^{2})L^{(\alpha)}_{n-m-1}(x^{2})\Bigr{]}$
(4.4)

resulting from (2.5), it follows that the right-hand sides of these relations
coincide provided we set $\alpha=1/2$ in the second one. This completes the
proof of (4.3).

Let us now turn ourselves to odd $X_{2m}$-Hermite EOPs, which from (2.2) can
be written as

$\hat{H}^{({\rm III},2m)}_{2n+1}(x)=(-1)^{m}[H_{2m}({\rm
i}x)H_{2n-2m+1}(x)+{\rm i}mH_{2m-1}({\rm i}x)H_{2n-2m}(x)]$

for $n=m,m+1,m+2,\ldots$. On taking (4.1) and (4.2) into account, this
equation reduces to

$\hat{H}^{({\rm
III},2m)}_{2n+1}(x)=(-1)^{m}x\Bigl{[}L^{(-1/2)}_{m}(-x^{2})L^{(1/2)}_{n-m}(x^{2})-mL^{(1/2)}_{m-1}(-x^{2})L^{(-1/2)}_{n-m}(x^{2})\Bigr{]}.$

Here the square bracket on the right-hand side does not look to that in (4.4)
for any value of $\alpha$.

In the special case where $n=m$, however, we directly get that

$\hat{H}^{({\rm III},2m)}_{2m+1}(x)=(-1)^{m}[xH_{2m}({\rm i}x)+{\rm
i}mH_{2m-1}({\rm i}x)],$

which by a straightforward application of (A.1) and (A.2) reduces to

$\hat{H}^{({\rm III},2m)}_{2m+1}(x)=(-1)^{m+1}{\rm i}H_{2m+1}({\rm i}x),$

or, from (4.2),

$\hat{H}^{({\rm III},2m)}_{2m+1}(x)=(-1)^{m}xL^{(1/2)}_{m}(-x^{2}).$

This result may be compared with

$\hat{L}^{({\rm
III},m-1)}_{m}(x^{2};\alpha)=(-1)^{m}L^{(-\alpha-1)}_{m}(-x^{2}),$

which is a direct consequence of (2.7). Hence one may write

$\hat{H}^{({\rm III},2m)}_{2m+1}(x)=x\hat{L}^{({\rm
III},m-1)}_{m}\Bigl{(}x^{2};-\frac{3}{2}\Bigr{)},\qquad m=1,2,3,\ldots.$

Apart from this result, the generalization of the Hermite-Laguerre quadratic
transformation for odd $n$ values remains an open problem.

##  5 Conclusion

In the present paper, we have shown that the known asymptotic relations
interconnecting Jacobi, Laguerre, and Hermite COPs can be generalized to the
corresponding EOPs of codimension $m$.

For such a purpose, we have first listed the sets of EOPs to be considered,
together with their explicit expressions in terms of COPs. We have also
provided their link with some other definitions found in the literature.

We have then stated and proved limit relations allowing to get
$X_{m}$-Laguerre EOPs of type I, II, or III from $X_{m}$-Jacobi EOPs of the
same type, as well as $X_{m}$-Hermite EOPs of type III from $X_{m}$-Jacobi or
$X_{m}$-Laguerre EOPs of type III.

Finally, we have established that a quadratic transformation applied to
$X_{m}$-Laguerre EOPs of type III enables one to obtain even $X_{2m}$-Hermite
EOPs. Whether a similar transformation would lead to odd $X_{2m}$-Hermite EOPs
remains an open problem for future investigation.

Studying an extension of the limit relations presented here to more involved
kinds of EOPs is also an open question for future work.

## Appendix A: Relations satisfied by COPs

In this Appendix, we list some relations satisfied by COPs [9, 10], rewritten
in monic form:

$\partial_{x}H_{n}(x)=nH_{n-1}(x),$ (A.1)
$(\partial_{x}-2x)H_{n}(x)=-2H_{n+1}(x),$ (A.2)
$\partial_{x}L^{(\alpha)}_{n}(x)=nL^{(\alpha+1)}_{n-1}(x),$ (A.3)
$L^{(\alpha-1)}_{n}(x)=L^{(\alpha)}_{n}(x)+nL^{(\alpha)}_{n-1}(x),$ (A.4)
$(x\partial_{x}+\alpha-x)L^{(\alpha)}_{n}(x)=-L^{(\alpha-1)}_{n+1}(x),$ (A.5)
$L^{(\alpha)}_{n+1}(x)+(2n+\alpha+1-x)L^{(\alpha)}_{n}(x)+n(n+\alpha)L^{(\alpha)}_{n-1}(x)=0,$
(A.6)
$xL^{(\alpha+1)}_{n}(x)=(n+\alpha+1)L^{(\alpha)}_{n}(x)+L^{(\alpha)}_{n+1}(x),$
(A.7)
$\partial_{x}P^{(\alpha,\beta)}_{n}(x)=nP^{(\alpha+1,\beta+1)}_{n-1}(x),$
(A.8)
$(1-x^{2})\partial_{x}P^{(\alpha,\beta)}_{n}(x)+[\beta-\alpha-(\alpha+\beta)x]P^{(\alpha,\beta)}_{n}(x)=-(n+\alpha+\beta)P^{(\alpha-1,\beta-1)}_{n+1}(x),$
(A.9)
$\displaystyle(2n+\alpha+\beta+1)(2n+\alpha+\beta+2)\Bigl{[}(1-x)P^{(\alpha+1,\beta)}_{n}(x)+P^{(\alpha,\beta)}_{n+1}(x)\Bigr{]}$
$\displaystyle\quad=2(n+\alpha+1)(n+\alpha+\beta+1)P^{(\alpha,\beta)}_{n}(x),$
(A.10)
$(2n+\alpha+\beta-1)(2n+\alpha+\beta)\Bigl{[}P^{(\alpha,\beta)}_{n}(x)-P^{(\alpha-1,\beta)}_{n}(x)\Bigr{]}=2n(n+\beta)P^{(\alpha,\beta)}_{n-1}(x),$
(A.11)
$(2n+\alpha+\beta-1)\Bigl{[}P^{(\alpha,\beta-1)}_{n}(x)-P^{(\alpha-1,\beta)}_{n}(x)\Bigr{]}=2nP^{(\alpha,\beta)}_{n-1}(x).$
(A.12)

## Appendix B: Proof of an identity satisfied by Jacobi polynomials

The purpose of this Appendix is to prove Eq. (2.12).

On using Eq. (A.10) for $n\to m-1$, $\alpha\to-\alpha-1$, $\beta\to\beta$, we
obtain

$(1-x)P^{(-\alpha,\beta)}_{m-1}(x)=\frac{2(m-\alpha-1)(m-\alpha+\beta-1)}{(2m-\alpha+\beta-2)(2m-\alpha+\beta-1)}P^{(-\alpha-1,\beta)}_{m-1}(x)-P^{(-\alpha-1,\beta)}_{m}(x),$

so the left-hand side of (2.12) becomes

$\displaystyle(\alpha+1)P^{(-\alpha-1,\beta-1)}_{m}(x)+m(1-x)P^{(-\alpha,\beta)}_{m-1}(x)=(\alpha+1)P^{(-\alpha-1,\beta-1)}_{m}(x)$
$\displaystyle\quad{}+\frac{2m(m-\alpha-1)(m-\alpha+\beta-1)}{(2m-\alpha+\beta-2)(2m-\alpha+\beta-1)}P^{(-\alpha-1,\beta)}_{m-1}(x)-mP^{(-\alpha-1,\beta)}_{m}(x).$

Equation (A.12) for $n\to m$, $\alpha\to-\alpha-1$, $\beta\to\beta$, namely

$P^{(-\alpha-1,\beta-1)}_{m}(x)=P^{(-\alpha-2,\beta)}_{m}(x)+\frac{2m}{2m-\alpha+\beta-2}P^{(-\alpha-1,\beta)}_{m-1}(x),$

enables one to transform it into

$\displaystyle(\alpha+1)P^{(-\alpha-1,\beta-1)}_{m}(x)+m(1-x)P^{(-\alpha,\beta)}_{m-1}(x)$
$\displaystyle\quad=(\alpha+1)P^{(-\alpha-2,\beta)}_{m}(x)+\frac{2m^{2}(m+\beta)}{(2m-\alpha+\beta-2)(2m-\alpha+\beta-1)}P^{(-\alpha-1,\beta)}_{m-1}(x)-mP^{(-\alpha-1,\beta)}_{m}(x).$

Finally, Eq. (A.11) for $n\to m$, $\alpha\to-\alpha-1$, $\beta\to\beta$ or

$\displaystyle(2m-\alpha+\beta-2)(2m-\alpha+\beta-1)P^{(-\alpha-1,\beta)}_{m}(x)-2m(m+\beta)P^{(-\alpha-1,\beta)}_{m-1}(x)$
$\displaystyle\quad=(2m-\alpha+\beta-2)(2m-\alpha+\beta-1)P^{(-\alpha-2,\beta)}_{m}(x)$

leads to the desired result (2.12).

## Appendix C: Examples of limit relations

To start with, let us note that if we formally set $m=0$ in the EOPs
definitions, we obtain that for $n=0,1,2,\ldots$,

$\hat{H}^{({\rm III},0)}_{n}(x)=H_{n}(x),$ (C.1) $\hat{L}^{({\rm
I},0)}_{n}(x;\alpha)=L^{(\alpha+1)}_{n}(x),\qquad\hat{L}^{({\rm
II},0)}_{n}(x;\alpha)=L^{(\alpha-1)}_{n}(x),\qquad\hat{L}^{({\rm
III},0)}_{n}(x;\alpha)=L^{(\alpha-1)}_{n}(x),$ (C.2)
$\displaystyle\hat{P}^{({\rm
I},0)}_{n}(x;\alpha,\beta)=P^{(\alpha+1,\beta-1)}_{n}(x),\qquad\hat{P}^{({\rm
II},0)}_{n}(x;\alpha,\beta)=P^{(\alpha-1,\beta+1)}_{n}(x),$
$\displaystyle\hat{P}^{({\rm
III},0)}_{n}(x;\alpha,\beta)=P^{(\alpha-1,\beta+1)}_{n}(x),$ (C.3)

after some straightforward application of identities satisfied by COPs. This
shows that for $m=0$, the new limit relations of Section 3 connecting EOPs
agree with the well-known ones satisfied by COPs.

For positive values of $m$, one gets for instance

$\hat{H}^{({\rm III},2)}_{3}(x)=x^{3}+\frac{3}{2}x,$ $\hat{L}^{({\rm
I},1)}_{2}(x;\alpha)=x^{3}-(\alpha+4)x^{2}-(\alpha+1)(\alpha+4)x+(\alpha+1)(\alpha+2)(\alpha+4),$
$\hat{L}^{({\rm
II},1)}_{2}(x;\alpha)=x^{3}-(\alpha+2)x^{2}-(\alpha-1)(\alpha+2)x+(\alpha-1)\alpha(\alpha+2),$
$\hat{L}^{({\rm
III},2)}_{3}(x;\alpha)=x^{3}-3(\alpha-2)x^{2}+3(\alpha-1)(\alpha-2)x-\alpha(\alpha-1)(\alpha-2),$
$\displaystyle\hat{P}^{({\rm
I},1)}_{2}(x;\alpha,\beta)=x^{3}+\frac{2(\alpha-\beta+2)^{2}+(\alpha+\beta)(\alpha+\beta+6)}{(\alpha+\beta+4)(\alpha-\beta+2)}x^{2}$
$\displaystyle\quad{}+\frac{(\alpha-\beta+2)^{2}+(\alpha+\beta)(2\alpha+2\beta+9)}{(\alpha+\beta+3)(\alpha+\beta+4)}x+\frac{(\alpha-\beta+2)^{2}(\alpha+\beta+2)-(\alpha+\beta)(\alpha+\beta+6)}{(\alpha-\beta+2)(\alpha+\beta+3)(\alpha+\beta+4)},$
$\displaystyle\hat{P}^{({\rm
II},1)}_{2}(x;\alpha,\beta)=x^{3}+\frac{2(\alpha-\beta-2)^{2}+(\alpha+\beta)(\alpha+\beta+6)}{(\alpha+\beta+4)(\alpha-\beta-2)}x^{2}$
$\displaystyle\quad{}+\frac{(\alpha-\beta-2)^{2}+(\alpha+\beta)(2\alpha+2\beta+9)}{(\alpha+\beta+3)(\alpha+\beta+4)}x+\frac{(\alpha-\beta-2)^{2}(\alpha+\beta+2)-(\alpha+\beta)(\alpha+\beta+6)}{(\alpha-\beta-2)(\alpha+\beta+3)(\alpha+\beta+4)},$
$\displaystyle\hat{P}^{({\rm
III},2)}_{3}(x;\alpha,\beta)=x^{3}+\frac{3(\alpha-\beta)}{\alpha+\beta-4}x^{2}+3\frac{(\alpha-\beta)^{2}+\alpha+\beta-4}{(\alpha+\beta-3)(\alpha+\beta-4)}x$
$\displaystyle\quad{}+\frac{(\alpha-\beta)[(\alpha-\beta)^{2}+3\alpha+3\beta-10]}{(\alpha+\beta-2)(\alpha+\beta-3)(\alpha+\beta-4)},$

on which one can directly check the relations

$\lim_{\beta\to\infty}\beta^{3}\hat{P}^{({\rm
I},1)}_{2}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}=-8\hat{L}^{({\rm
I},1)}_{2}(x;\alpha),$ $\lim_{\beta\to\infty}\beta^{3}\hat{P}^{({\rm
II},1)}_{2}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}=-8\hat{L}^{({\rm
II},1)}_{2}(x;\alpha),$ $\lim_{\beta\to\infty}\beta^{3}\hat{P}^{({\rm
III},2)}_{3}\Bigl{(}1-\frac{2x}{\beta};\alpha,\beta\Bigr{)}=-8\hat{L}^{({\rm
III},2)}_{3}(x;\alpha),$ $\lim_{\alpha\to\infty}\alpha^{3/2}\hat{P}^{({\rm
III},2)}_{3}\Bigl{(}\frac{x}{\sqrt{\alpha}};\alpha,\alpha\Bigr{)}=\hat{H}^{({\rm
III},2)}_{3}(x),$
$\lim_{\alpha\to\infty}\frac{1}{(2\alpha)^{3/2}}\hat{L}^{({\rm
III},2)}_{3}\bigl{(}\sqrt{2\alpha}x+\alpha;\alpha\bigr{)}=\hat{H}^{({\rm
III},2)}_{3}(x).$

## Acknowledgments

The author was supported by the Fonds de la Recherche Scientifique-FNRS under
Grant No. 4.45.10.08.

## References

  * [1] Bonneux N., Exceptional Jacobi polynomials, J. Approx. Theory 239 (2019), 72–112, arXiv:1804.01323.
  * [2] Durán A.J., Higher order recurrence relation for exceptional Charlier, Meixner, Hermite and Laguerre orthogonal polynomials, Integral Transforms Spec. Funct. 26 (2015), 357–376, arXiv:1409.4697.
  * [3] García-Ferrero M.A., Gómez-Ullate D., Milson R., A Bochner type characterization theorem for exceptional orthogonal polynomials, J. Math. Anal. Appl. 472 (2019), 584–626, arXiv:1603.04358.
  * [4] García-Ferrero M.A., Gómez-Ullate D., Milson R., Exceptional Legendre polynomials and confluent Darboux transformations, SIGMA 17 (2021), 016, 19 pages, arXiv:2008.02822.
  * [5] Gómez-Ullate D., Grandati Y., Milson R., Rational extensions of the quantum harmonic oscillator and exceptional Hermite polynomials, J. Phys. A: Math. Theor. 47 (2014), 015203, 27 pages, arXiv:1306.5143..
  * [6] Gómez-Ullate D., Kamran N., Milson R., An extended class of orthogonal polynomials defined by a Sturm-Liouville problem, J. Math. Anal. Appl. 359 (2009), 352–367, arXiv:0807.3939.
  * [7] Gómez-Ullate D., Kamran N., Milson R., Two-step Darboux transformations and exceptional Laguerre polynomials, J. Math. Anal. Appl. 387 (2012), 410–418, arXiv:1103.5724.
  * [8] Gómez-Ullate D., Marcellán F., Milson R., Asymptotic and interlacing properties of zeros of exceptional Jacobi and Laguerre polynomials, J. Math. Anal. Appl. 399 (2013), 480–495, arXiv:1204.2282.
  * [9] Gradshteyn I.S., Ryzhik I.M., Table of integrals, series, and products, Academic Press, New York, 1980.
  * [10] Koekoek R., Lesky P.A., Swarttouw R.F., Hypergeometric orthogonal polynomials and their $q$-analogues, Springer-Verlag, Berlin Heidelberg, 2010.
  * [11] Kuijlaars A.B.J., Milson R., Zeros of exceptional Hermite polynomials, J. Approx. Theory 200 (2015), 28–39, arXiv:1412.6364.
  * [12] Liaw C., Littlejohn L.L., Milson R., Stewart J., The spectral analysis of three families of exceptional Laguerre polynomials, J. Approx. Theory 202 (2016), 5–41.
  * [13] Liaw C., Littlejohn L., Stewart Kelly J., Spectral analysis for the exceptional $X_{m}$-Jacobi equation, Electron. J. Differential Equations 2015 (2015) No. 194, pp. 1–10, arXiv:1501.04698.
  * [14] Marquette I., Quesne C., Two-step rational extensions of the harmonic oscillator: exceptional orthogonal polynomials and ladder operators, J. Phys. A: Math. Theor. 46 (2013), 155201, 14 pages, arXiv:1212.3474.
  * [15] Miki H., Tsujimoto S., A new recurrence formula for generic exceptional orthogonal polynomials, J. Math. Phys. 56 (2015), 033502, 13 pages, arXiv:1410.0183.
  * [16] Odake S., Recurrence relations of the multi-indexed orthogonal polynomials, J. Math. Phys. 54 (2013), 083506, 18 pages, arXiv:1303.5820.
  * [17] Odake S., Recurrence relations of the multi-indexed orthogonal polynomials : III, J. Math. Phys. 57 (2016), 023514, 24 pages, arXiv:1509.08213.
  * [18] Odake S., Sasaki R., Infinitely many shape invariant potentials and new orthogonal polynomials, Phys. Lett. B 679 (2009), 414–417, arXiv:0906.0142.
  * [19] Odake S., Sasaki R., Exactly solvable quantum mechanics and infinite families of multi-indexed orthogonal polynomials, Phys. Lett. B 702 (2011), 164–170, arXiv:1105.0508.
  * [20] Quesne C., Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry, J. Phys. A: Math. Theor. 41 (2008), 392001, 6 pages, arXiv:0807.4087.
  * [21] Quesne C., Solvable rational potentials and exceptional orthogonal polynomials in supersymmetric quantum mechanics, SIGMA 5 (2009), 084, 24 pages, arXiv:0906.2331.


# Learning the Proximity Operator in Unfolded ADMM
for Phase Retrieval††thanks: This work is supported by the European Research
Council (ERC FACTORY-CoG-6681839) and by ANITI under grant agreement
ANR-19-PI3A-0004.

Pierre-Hugo Vial , Paul Magron , Thomas Oberlin , Cédric
Févotte22footnotemark: 2 IRIT, Université de Toulouse, CNRS, Toulouse, France
(e-mail: firstname.lastname@irit.fr).Université de Lorraine, CNRS, Inria,
LORIA, F-54000 Nancy, France (e-mail: firstname.lastname@inria.fr).ISAE-
SUPAERO, Université de Toulouse, France (e-mail: firstname.lastname@isae-
supaero.fr).

###### Abstract

This paper considers the phase retrieval (PR) problem, which aims to
reconstruct a signal from phaseless measurements such as magnitude or power
spectrograms. PR is generally handled as a minimization problem involving a
quadratic loss. Recent works have considered alternative discrepancy measures,
such as the Bregman divergences, but it is still challenging to tailor the
optimal loss for a given setting. In this paper we propose a novel strategy to
automatically learn the optimal metric for PR. We unfold a recently introduced
ADMM algorithm into a neural network, and we emphasize that the information
about the loss used to formulate the PR problem is conveyed by the proximity
operator involved in the ADMM updates. Therefore, we replace this proximity
operator with trainable activation functions: learning these in a supervised
setting is then equivalent to learning an optimal metric for PR. Experiments
conducted with speech signals show that our approach outperforms the baseline
ADMM, using a light and interpretable neural architecture.

Keywords— Phase retrieval, audio, proximity operator learning, deep unfolding,
ADMM.

##  1 Introduction

Phase retrieval (PR) consists in reconstructing data from phaseless
nonnegative measurements. This problem finds practical applications in various
areas such as optical imaging [1, 2] and audio signal processing [3, 4], which
is the domain of interest of this paper. Traditionally, PR is formulated as
the following optimization problem:

$\underset{\mathbf{x}\in\mathbb{R}^{L}}{\text{min}}\||\mathbf{A}\mathbf{x}|^{d}-\mathbf{r}\|^{2},$
(1)

where $\mathbf{A}\in\mathbb{C}^{K\times L}$ is the measurement operator,
$\mathbf{r}\in\mathbb{R}_{+}^{K}$ are the phaseless measurements, and $||.||$
denotes the Euclidean norm. In audio signal processing, $\mathbf{A}$ is often
the short-time Fourier transform (STFT), and the measurements are either
magnitude ($d=1$) or power ($d=2$) spectrograms. In the seminal work [5], the
authors consider STFT magnitude measurements and propose an iterative
procedure, known as the Griffin-Lim algorithm (GLA) which is proved to
converge to a critical point of (1). Other optimization algorithms were
proposed to tackle Problem (1), such as majorization-minimization [6],
gradient descent [7] or alternating direction method of multipliers (ADMM)
[8]. PR has also been addressed by replacing the quadratic loss in (1) with
alternative discrepancy measures such as Bregman divergences [9], which have
been shown to lead to superior results in many audio applications [10, 11].
However, prescribing a loss that is optimal for all signal processing problems
and classes of audio signals remains challenging [9].

On the other hand, recent PR approaches have leveraged deep neural networks
(DNNs) [12, 13, 14, 15]. Despite their successful performances in a large
number of tasks, the enthusiasm for DNNs can be tempered by a general lack of
explainability due to their black box structure, and by their limited ability
to generalize to unseen data or experimental conditions. Deep unfolding (or
unrolling) [16, 17] is a promising attempt to alleviate these limitations with
model-based architectures derived from iterative algorithms. This strategy
consists in considering each iteration of an optimization algorithm as a
(trainable) layer of a DNN. It has proved very promising when applied to many
algorithms in signal processing [17, 18, 19], including audio PR [20, 21, 22].
However, these approaches treat the linear layers as trainable parameters
while activation functions remain fixed, although recent studies show that
trainable activation functions can upgrade performance [23, 24, 25]. Besides,
they do not address the problem of learning the metric involved in the PR
formulation problem.

In this paper, we propose to unfold the ADMM algorithm for PR proposed in [9].
Our method builds upon observing that the choice of the discrepancy measure
only affects the computation of a proximity operator in the ADMM updates.
Therefore, we can recast the problem of metric learning as a problem of
proximity operator learning in the unfolded ADMM. To that end, we replace this
proximity operator with a trainable activation function. We show that the
proposed parametrization of the network is connected to the metric involved in
the original optimization problem, which yields an interpretable architecture.
Experiments performed on speech signals demonstrate the efficiency of our
method, which outperforms a baseline ADMM [8] with a number of iterations
equal to the number of layers in the unfoldeed ADMM.

The rest of this paper is structured as follows. Section 2 presents the
related work. Then, the proposed method is introduced in Section 3 and tested
experimentally in Section 4. Finally, Section 5 draws some concluding remarks.

Mathematical notation:

  * •

$\mathbf{A}$ (capital, bold font): matrix.

  * •

$\mathbf{x}$ (lower case, bold font): vector.

  * •

$z$ (regular): scalar.

  * •

$|.|$, $\angle(.)$: magnitude and complex angle, respectively.

  * •

$\mathbf{x}^{\mathsf{H}}$: Hermitian transpose.

  * •

$\odot$, $(.)^{d}$, fraction bar: element-wise matrix or vector
multiplication, power, and division, respectively.

##  2 Related work

###  2.1 PR with Bregman divergences

In [9], we proposed to reformulate the PR problem by substituting the
quadratic loss in (1) with a Bregman divergence. Bregman divergences encompass
beta-divergences [26], which include the Kullback-Leibler and Itakura-Saito
divergences, as well as the quadratic loss. They are defined as follows:

$\mathcal{D}_{\psi}(\mathbf{p}\,\bm{|}\,\mathbf{q})=\sum_{k=1}^{K}\left[\psi(p_{k})-\psi(q_{k})-\psi^{\prime}(q_{k})(p_{k}-q_{k})\right],$
(2)

where $\psi$ is a strictly-convex, continuously-differentiable generating
function (with derivative $\psi^{\prime}$). Since Bregman divergences are not
symmetric in general, we considered the two following problems, respectively
termed “right” and “left”:

$\underset{\mathbf{x}\in\mathbb{R}^{L}}{\text{min}}\,\mathcal{D}_{\psi}(\mathbf{r}\,\bm{|}\,|\mathbf{A}\mathbf{x}|^{d})\text{
and
}\underset{\mathbf{x}\in\mathbb{R}^{L}}{\text{min}}\,\mathcal{D}_{\psi}(|\mathbf{A}\mathbf{x}|^{d}\,\bm{|}\,\mathbf{r}).$
(3)

Two algorithms were derived for solving these problems, based on gradient
descent and ADMM [9].

###  2.2 PR via ADMM

In this work we focus on the “left” PR problem using ADMM, since it produced
the best reconstruction results in [9]. Drawing on [8], this approach consists
in reformulating (3) into the following constrained optimization problem:

$\underset{\mathbf{x}\in\mathbb{R}^{L},\mathbf{u}\in\mathbb{R}^{K}_{+},\theta\in\left[0;2\pi\right[^{K}}{\text{min}}\,\mathcal{D}_{\psi}(\mathbf{u}\,\bm{|}\,\mathbf{r})\quad\mbox{s.t.}\quad(\mathbf{A}\mathbf{x})^{d}=\mathbf{u}\odot
e^{i\bm{\theta}},$ (4)

where $\mathbf{u}$ and $\bm{\theta}$ are auxiliary variables for the magnitude
and phase of $\mathbf{A}\mathbf{x}$. From (4) we obtain the augmented
Lagrangian:

$\mathcal{L}(\mathbf{x},\mathbf{u},\bm{\theta},\bm{\lambda})=\mathcal{D}_{\psi}(\mathbf{u}\,\bm{|}\,\mathbf{r})+\mathfrak{R}\left(\bm{\lambda}^{\mathsf{H}}((\mathbf{A}\mathbf{x})^{d}-\mathbf{u}\odot
e^{i\bm{\theta}})\right)+\frac{\rho}{2}\left\|(\mathbf{A}\mathbf{x})^{d}-\mathbf{u}\odot
e^{i\bm{\theta}}\right\|^{2},$

where $\mathfrak{R}$ is the real part function,
$\bm{\lambda}\in\mathbb{C}^{K}$ are the Lagrange multipliers and $\rho>0$ is
the augmentation parameter. Alternate minimization of $\mathcal{L}$ leads to
the following update rules [9]:

$\displaystyle\mathbf{h}_{t+1}$
$\displaystyle=(\mathbf{A}\mathbf{x}_{t})^{d}+\frac{\bm{\lambda}_{t}}{\rho}$
(5) $\displaystyle\mathbf{u}_{t+1}$
$\displaystyle=\text{prox}_{\rho^{-1}\mathcal{D}_{\psi}(\cdot\,\bm{|}\,\mathbf{r})}(|\mathbf{h}_{t+1}|)$
(6) $\displaystyle\bm{\theta}_{t+1}$ $\displaystyle=\angle\mathbf{h}_{t+1}$
(7) $\displaystyle\mathbf{x}_{t+1}$
$\displaystyle=\mathbf{A}^{\mathsf{H}}\big{(}\mathbf{u}_{t+1}\odot
e^{i\bm{\theta}_{t+1}}-\frac{\bm{\lambda}_{t}}{\rho}\big{)}^{1/d}$ (8)
$\displaystyle\bm{\lambda}_{t+1}$
$\displaystyle=\bm{\lambda}_{t}+\rho(\mathbf{A}\mathbf{x}_{t+1}-\mathbf{u}_{t+1}\odot
e^{i\bm{\theta}_{t+1}}),$ (9)

where $\text{prox}_{f}$ denotes the proximity operator of a function $f$,
defined as the mapping of a vector $\mathbf{y}$ to the solution of the
minimization problem:

$\text{prox}_{\rho^{-1}f}(\mathbf{y}):=\text{arg}\underset{\mathbf{z}\in\mathbb{R}^{K}}{\text{min}}f(\mathbf{z})+\frac{\rho}{2}\|\mathbf{z}-\mathbf{y}\|^{2}.$
(10)

However, a closed-form expression of this operator is not available for every
Bregman divergence. To alleviate this issue, we propose in this paper to
replace this proximity operator with a learnable neural activation function,
as detailed in the next section.

##  3 Proposed method

###  3.1 General architecture

Figure 1: One layer of the proposed unfolded architecture.

The ADMM updates detailed in 2.2 consist in successive linear and nonlinear
computations. As such, this algorithm can be viewed as a neural network
$\bm{\mathsf{U}}$ via unfolding:

$(\mathbf{x}_{T},\,\bm{\lambda}_{T})=\bm{\mathsf{U}}(\mathbf{x}_{0},\,\bm{\lambda}_{0})=\mathsf{U}_{1}\circ\dots\circ\mathsf{U}_{T}(\mathbf{x}_{0},\,\bm{\lambda}_{0}),$
(11)

where $\mathsf{U}_{t}$ denotes the $t$-th layer of the network, mimicking the
$t$-th iteration of the ADMM algorithm, as illustrated in Fig. 1. The layer
$\mathsf{U}_{t}$ can be decomposed into two linear parts denoted by
$\mathsf{L}_{t}^{(1)}$ and $\mathsf{L}_{t}^{(2)}$, and a nonlinear part
$\mathsf{NL}_{t}$ as follows:

$\displaystyle\mathsf{L}_{t}^{(1)}$
$\displaystyle:(\mathbf{x}_{t-1},\,\bm{\lambda}_{t-1})\mapsto\mathbf{h}_{t}$
(12) $\displaystyle\mathsf{NL}_{t}$
$\displaystyle:\mathbf{h}_{t}\mapsto(\mathbf{u}_{t},\,\bm{\theta}_{t})=\left(\mathsf{F}_{t}(|\mathbf{h}_{t}|,\,\mathbf{r}),\,\angle\mathbf{h}_{t}\right)$
(13) $\displaystyle\mathsf{L}_{t}^{(2)}$
$\displaystyle:(\mathbf{x}_{t-1},\,\bm{\lambda}_{t-1},\mathbf{u}_{t},\,\bm{\theta}_{t})\mapsto(\mathbf{x}_{t},\,\bm{\lambda}_{t}),$
(14)

with $\mathbf{h}_{t},\mathbf{x}_{t},\,\bm{\lambda}_{t}$ respectively defined
as in (5), (8) and (9). $\mathsf{F}_{t}$ denotes a parameterized sublayer
modeling the proximity operator of equation (6). Since the choice of the
discrepancy measure $\mathcal{D}_{\psi}$ only affects the proximity operator
(6) in the updates, we can recast the problem of metric learning as the
problem of proximity operator learning. We propose to leverage a trainable
activation function in order to model this layer and learn the proximity
operator.

###  3.2 Proposed parameterization

To build the non-linear sublayers $\mathsf{F}_{t}$ that model
$\text{prox}_{\rho^{-1}\mathcal{D}_{\psi(\cdot\,|\,\mathbf{r})}}$, we first
reformulate this operator as follows. Let $\mathbf{v}\in\mathbb{R}^{K}$ and
$f(\mathbf{z})=\sum_{k=1}^{K}[\psi(z_{k})+v_{k}z_{k}]$. We have [27]:

$\text{prox}_{\rho^{-1}f}(\mathbf{y})=\text{prox}_{\rho^{-1}\tilde{\psi}}\left(\mathbf{y}-\rho^{-1}\mathbf{v}\right),$
(15)

where $\tilde{\psi}(\mathbf{z})=\sum_{k}\psi(z_{k})$. Setting
$\mathbf{v}=-\psi^{\prime}(\mathbf{r})$ in (15), with $\psi^{\prime}$ applied
entrywise, it is straightforward to see that:

$\text{prox}_{\rho^{-1}\mathcal{D}_{\psi}(\cdot\,\bm{|}\,\mathbf{r})}(\mathbf{y})=\text{prox}_{\rho^{-1}\tilde{\psi}}(\mathbf{y}+\rho^{-1}\psi^{\prime}(\mathbf{r})).$
(16)

This formulation of the proximity operator is more convenient than (6) since
the measurements $\mathbf{r}$ no longer appear in the input function of the
proximity operator, but instead in the argument of the latter (with
$\mathbf{y}$). This leads to a more natural parametrization for unrolling.

Let us first derive the proximity operator (16) in a simple scenario, namely
the quadratic loss ($\tilde{\psi}=\frac{1}{2}\|\cdot\|^{2}$,
$\psi^{\prime}(\mathbf{r})=\mathbf{r}$). In this case we have [8]:

$\text{prox}_{\rho^{-1}\frac{1}{2}\|\cdot-\mathbf{r}\|^{2}}(\mathbf{y})=\frac{\mathbf{y}+\rho^{-1}\mathbf{r}}{1+\rho^{-1}}.$
(17)

As a result, a first simple approach for proximity operator learning would
consist in treating $\rho$ as a learnable parameter. However, early
experiments have shown poor performance with this approach, which is due to
the very low expressive power of such a model (only one scalar value). More
generally, one can consider a beta-divergence with shape parameter $\beta$,
for which $\psi^{\prime}(\mathbf{r})=\frac{\mathbf{r}^{\beta-1}}{\beta-1}$
[26]. However, the proximity operator of $\tilde{\psi}$ is not available for
every beta-divergence.

To alleviate this issue, we model this unavailable proximity operator using
Adaptive Piecewise Linear (APL) activations [28]. They are defined by:

$\textsf{APL}(\mathbf{y}):=\max(\mathbf{y},0)+\sum_{c=1}^{C}w_{c}\max(-\mathbf{y}+b_{c},0),$
(18)

where $w_{c}$ and $b_{c}$ are learnable parameters controlling the slopes and
biases of the linear segments, and the $\max$ is applied entry-wise. Then, we
propose the following parametrization of the nonlinear layer $\mathsf{F}_{t}$:

$\mathsf{F}_{t}(\mathbf{y},\,\mathbf{r})=\textsf{APL}_{t}\left(\gamma^{(1)}_{t}\mathbf{y}+\gamma^{(2)}_{t}\frac{\mathbf{r}^{\beta_{t}-1}}{\beta_{t}-1}\right),$
(19)

with learnable parameters $w_{c,t},b_{c,t}$, $\gamma^{(1)}_{t}$,
$\gamma^{(2)}_{t}$, and $\beta_{t}$. Even though ad hoc, this parametrization
is motivated by the following considerations:

  * •

APL can represent any continuous piecewise linear function over a subset of
real numbers. As such, it generalizes the proximity operator obtained in the
quadratic case [8].

  * •

The term in the form of $\frac{\mathbf{r}^{\beta-1}}{\beta-1}$ in (19) is
reminiscent of $\psi^{\prime}$ for beta-divergences, as mentioned above.

  * •

Introducing learnable weights $\gamma^{(1)}_{t}$ and $\gamma^{(2)}_{t}$ allows
to increase the model capacity, as it was shown beneficial in our preliminary
experiments.

Note that when
$w_{c}=0,\,\gamma_{t}^{(1)}=\frac{1}{1+\rho^{-1}},\,\gamma_{t}^{(2)}=\frac{\rho^{-1}}{1+\rho^{-1}}$
and $\beta_{t}=2$, $\mathsf{F}_{t}$ is equal to the proximity operator for the
quadratic loss (17). Overall, our parametrization (19) is a good trade-off
between tractability, interpretability, and expressiveness.

Two variants of the proposed architecture will be considered in our
experiments. In the “untied” variant, each layer uses different parameters,
and the global set of parameters is
$\left\\{\\{w_{c,t},b_{c,t}\\}_{c=1}^{C},\gamma^{(1)}_{t},\gamma^{(2)}_{t},\beta_{t}\right\\}_{t=1}^{T}$,
while in the “tied” variant, the parameters are shared among layers, i.e.,
constant with $t$.

In the end, after learning these parameters (see Section 4.1.2), the proposed
method, termed unfolded ADMM (UADMM), estimates a signal $\mathbf{x}_{T}$ via
$(\mathbf{x}_{T},\bm{\lambda}_{T})=\bm{\mathsf{U}}(\mathbf{x}_{0},\bm{\lambda}_{0})$,
where $\mathbf{x}_{0}$ is some initial estimate.

###  3.3 Discussion about interpretability

Under mild assumptions and as detailed in the supplementary material,111In
particular, the APL unit must be strictly increasing, something we ensure by
imposing the weights to be negative, using $w_{c,t}=-\tilde{w}_{c,t}^{2}$. we
can prove that there exists a close-form function
$f_{\mathbf{r},t}:\mathbb{R}^{K}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ such
that
${\mathsf{F}_{t}(\mathbf{y},\mathbf{r})=\text{prox}_{f_{\mathbf{r},t}}(\mathbf{y})}$,
(see (31) in the supplementary material). In the “tied” variant, where
$f_{\mathbf{r},t}=f_{\mathbf{r}}$, reconstructing $f_{\mathbf{r}}$ from the
learned parameters is analogous to identifying the metric
$\mathcal{D}_{\psi}(\cdot\,|\,\mathbf{r})$ involved in the PR optimization
problem. With the relaxation proposed in the “untied” case, this
interpretation is more limited as $f_{\mathbf{r},t}$ is different in each
layer of the network.

Note that when replacing (16) with (19), we have disentangled the proximity
operator of $\tilde{\psi}$ and the derivative $\psi^{\prime}$, in addition to
introducing weights $\gamma^{(1)}_{t}$ and $\gamma^{(2)}_{t}$. As a result,
the function $f_{\mathbf{r}}$ is no longer guaranteed to be a Bregman
divergence, strictly speaking. Nevertheless, we can still interpret it as a
measure of discrepancy between $\mathbf{y}$ and $\mathbf{r}$.

##  4 Experiments

In this section, we assess the potential of UADMM for PR of speech signals.
Our code is implemented using the PyTorch framework [29] and will be made
available online for reproducibility, as soon as the paper is accepted.

###  4.1 Protocol

####  4.1.1 Data

We build a set of speech signals from the TIMIT dataset [30]. The dataset is
split into training, validation, and test subsets containing $1000$, $10$, and
$50$ utterances, respectively. The signals are mono, sampled at $22,500$ Hz
and cropped to 2 seconds. The STFT is computed with a 1024 samples-long (46
ms) self-dual sine window [9] and 50$\%$ overlap. STFT magnitudes ($d=1$) are
considered as nonnegative observations $\mathbf{r}$.

####  4.1.2 Training

The network is trained with the ADAM algorithm [31] using a learning rate of
$10^{-4}$. We use a structure with $T=15$ layers and $C=3$, as these values
have shown to be a good trade-off between performance and number of parameters
in preliminary experiments. Batches of $10$ signals with a maximum of $200$
epochs are used for training. Training is stopped when the loss function on
the validation subset starts increasing. Given that we consider speech
signals, we train the network by minimizing the negative short-term objective
intelligibility measure (STOI) [32] between the estimated and ground truth
signals. Indeed, this strategy was shown to be efficient for speech
enhancement applications [33, 34]. The negative STOI metric used for training
the network is implemented in PyTorch via the pytorch-stoi library [35].

####  4.1.3 Methods

As baselines, we consider GLA [5] (run for $1500$ iterations), and ADMM using
a quadratic loss and ${\rho=10^{-3}}$, since this setup has exhibited good
performance in our previous study under similar conditions [9]. ADMM is run
for a variable number of iterations, with a maximum of $1500$ (performance
does not further improve beyond). For fairness, the linear parts of UADMM use
the same value for $\rho$, and it is initialized such that $\mathsf{F}_{t}$
replicates the quadratic proximity operator (cf. Section 3.2). All methods use
the same initial signal estimate $\mathbf{x}_{0}$ computed using the ground
truth magnitudes $\mathbf{r}$ and a random uniform phase, and
$\bm{\lambda}_{0}=0$.

####  4.1.4 Evaluation

Reconstruction performance is assessed with the STOI metric (which ranges
between $0$ and $1$, higher is better) computed on the test set with the
pystoi library [36].

###  4.2 Results

Figure 2: Training loss (negative STOI) over epochs. Note that pytorch-stoi
implementation does not exactly replicate the original metric and consequently
yields values lower than $-1$. Figure 3: Performance on the test set. Each
box-plot is made up of a central line indicating the median, box edges
indicating the $1^{\text{st}}$ and $3^{\text{rd}}$ quartiles, whiskers
indicating the extremal values, and circles representing the outliers.

First, we display the training loss over epochs in Fig. 2. Both UADMM variants
outperform the baseline ADMM with $15$ iterations on the training set. UADMM-
untied reaches a lower loss value than its tied counterpart, which was
expected since this variant contains more learnable parameters. They reach a
performance comparable to that of ADMM using $150$ and $75$ iterations,
respectively. The results on the test set presented in Fig. 3 confirm that the
proposed UADMM approach significantly outperforms the classical ADMM using the
same number of iterations, as well as the GLA baseline. A fully-converged ADMM
algorithm (using $1500$ iterations) exhibits a higher STOI than our $15$
layers-based approach. Nonetheless, a more fair comparison would involve that
both approaches use the same total number of iterations/layers.

Figure 4: Evaluation with STOI over test dataset with iterated model. The
solid lines denote the mean STOI and the the light colored areas the values
between the first and the third quartile.

To that end, we consider an _ad hoc_ variant of our method, where we apply the
whole UADMM network several times: this allows to increase the total number of
“iterations” without the need for training a larger network from scratch. The
results presented in Fig. 4 show that this method consistently and
significantly outperforms ADMM for any number of iterations. In particular,
the performance of the fully-converged ADMM (after $1500$ iterations) is
reached at only $30$ “iterations” for UADMM-untied (i.e., twice the number of
trained layers), which exhibits the computational advantage of the proposed
approach. Note that UADMM-tied with $T$ layers is equivalent to applying $T$
iterations of a standard ADMM algorithm using a learned metric
$f_{\mathbf{r}}$. We display these learned metrics ($f_{\mathbf{r}}$ in the
tied case and $f_{\mathbf{r},t}$ in the untied case) in Fig. 5. These resemble
beta-divergences with $\beta\in[1.5,2.5]$. This is consistent with previous
results from the literature [37], where this range of values has shown good
performance for audio spectral decomposition.

Figure 5: Learned metrics $f_{r,t}(y)$ with $r=1$. The quadratic loss and
Kullback-Leibler divergence $\mathcal{D}_{KL}(y\,|\,r)$ are also displayed for
the sake of comparison. In the “tied” case, $f_{r}$ in analogous to
$\mathcal{D}_{\psi}(\cdot\,|\,r)$ involved in the PR optimization problem. For
clarity, only 3 of the 15 trained layers $f_{r,t}$ are displayed for the
“untied” case.

##  5 Conclusion

In this paper, we have addressed the problem of metric learning for phase
retrieval by unfolding the recently proposed ADMM algorithm [9] into a neural
network. We proposed to replace the proximity operator involved in this
algorithm with learnable activation functions, since this operator conveys the
information about the discrepancy measure used in formulating the PR problem.
Experiments conducted on speech signals show that this approach outperforms
the ADMM algorithm while keeping a light and interpretable structure.

## References

  * [1] A. Walther,  “The question of phase retrieval in optics,”  Optica Acta: International Journal of Optics, vol. 10, no. 1, pp. 41–49, 1963.
  * [2] R. W. Harrison,  “Phase problem in crystallography,”  Journal of the Optical Society of America A, vol. 10, no. 5, pp. 1046–1055, 1993.
  * [3] T. Gerkmann, M. Krawczyk-Becker, and J. Le Roux,  “Phase processing for single-channel speech enhancement: History and recent advances,”  IEEE Signal Processing Magazine, vol. 32, no. 2, pp. 55–66, March 2015.
  * [4] P. Mowlaee, R. Saeidi, and Y. Stylianou,  “Advances in phase-aware signal processing in speech communication,”  Speech Communication, vol. 81, pp. 1–29, July 2016.
  * [5] D. Griffin and J. Lim,  “Signal estimation from modified short-time Fourier transform,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 32, no. 2, pp. 236–243, April 1984.
  * [6] T. Qiu, P. Babu, and D. P. Palomar,  “PRIME: Phase retrieval via majorization-minimization,”  IEEE Transactions on Signal Processing, vol. 64, no. 19, pp. 5174–5186, October 2016.
  * [7] E. J Candès, X. Li, and M. Soltanolkotabi,  “Phase retrieval via Wirtinger flow: Theory and algorithms,”  IEEE Transactions on Information Theory, vol. 61, no. 4, pp. 1985–2007, April 2015.
  * [8] J. Liang, P. Stoica, Y. Jing, and J. Li,  “Phase retrieval via the alternating direction method of multipliers,”  IEEE Signal Processing Letters, vol. 25, no. 1, pp. 5–9, January 2018.
  * [9] P.-H. Vial, P. Magron, T. Oberlin, and C. Févotte,  “Phase retrieval with Bregman divergences and application to audio signal recovery,”  IEEE Journal of Selected Topics in Signal Processing, vol. 15, no. 1, pp. 51–64, January 2021.
  * [10] C. Févotte, N. Bertin, and J.-L. Durrieu,  “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis,”  Neural computation, vol. 21, no. 3, pp. 793–830, March 2009.
  * [11] J. Le Roux, H. Kameoka, N. Ono, A. de Cheveigné, and S. Sagayama,  “Computational auditory induction as a missing-data model-fitting problem with Bregman divergence,”  Speech Communication, vol. 53, no. 5, pp. 658 – 676, May-June 2011\.
  * [12] S. Takamichi, Y. Saito, N. Takamune, D. Kitamura, and H. Saruwatari,  “Phase reconstruction from amplitude spectrograms based on von-Mises-distribution deep neural network,”  in Proc. International Workshop on Acoustic Signal Enhancement (IWAENC), September 2018, pp. 286 – 290.
  * [13] L. Thieling, D. Wilhelm, and P. Jax,  “Recurrent phase reconstruction using estimated phase derivatives from deep neural networks,”  in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), June 2021, pp. 7088–7092.
  * [14] N. B. Thien, Y. Wakabayashi, K. Iwai, and T. Nishiura,  “Two-stage phase reconstruction using DNN and von Mises distribution-based maximum likelihood,”  in Proc. Annual Summit and Conference of the Asia-Pacific Signal and Information Processing Association (APSIPA ASC), December 2021, pp. 995–999.
  * [15] S. Ö. Arık, H. Jun, and G. Diamos,  “Fast spectrogram inversion using multi-head convolutional neural networks,”  IEEE Signal Processing Letters, vol. 26, no. 1, pp. 94–98, January 2018.
  * [16] K. Gregor and Y. LeCun,  “Learning fast approximations of sparse coding,”  in Proc. International Conference on Machine Learning (ICML), June 2010, pp. 399–406.
  * [17] J. R. Hershey, J. Le Roux, and F. Weninger,  “Deep unfolding: Model-based inspiration of novel deep architectures,”  arXiv preprint arXiv:1409.2574, 2014.
  * [18] C. Bertocchi, E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet, and M. Prato,  “Deep unfolding of a proximal interior point method for image restoration,”  Inverse Problems, vol. 36, no. 3, pp. 034005, February 2020.
  * [19] V. Monga, Y. Li, and Y. C. Eldar,  “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,”  IEEE Signal Processing Magazine, vol. 38, no. 2, pp. 18–44, March 2021.
  * [20] Y. Masuyama, K. Yatabe, Y. Koizumi, Y. Oikawa, and N. Harada,  “Deep Griffin–Lim iteration,”  in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2019, pp. 61–65.
  * [21] Z.-Q. Wang, J. Le Roux, D. Wang, and J. R. Hershey,  “End-to-end speech separation with unfolded iterative phase reconstruction,”  in Proc. Interspeech, September 2018.
  * [22] G. Wichern and J. Le Roux,  “Phase reconstruction with learned time-frequency representations for single-channel speech separation,”  in Proc. International Workshop on Acoustic Signal Enhancement (IWAENC), September 2018, pp. 396–400.
  * [23] S. Qian, H. Liu, C. Liu, S. Wu, and H. San Wong,  “Adaptive activation functions in convolutional neural networks,”  Neurocomputing, vol. 272, pp. 204–212, January 2018.
  * [24] A. Apicella, F. Isgrò, and R. Prevete,  “A simple and efficient architecture for trainable activation functions,”  Neurocomputing, vol. 370, pp. 1–15, December 2019.
  * [25] S. Scardapane, S. Van Vaerenbergh, S. Totaro, and A. Uncini,  “Kafnets: Kernel-based non-parametric activation functions for neural networks,”  Neural Networks, vol. 110, pp. 19–32, February 2019.
  * [26] R. Hennequin, B. David, and R. Badeau,  “Beta-divergence as a subclass of Bregman divergence,”  IEEE Signal Processing Letters, vol. 18, no. 2, pp. 83–86, February 2011.
  * [27] P. L. Combettes and V. R. Wajs,  “Signal recovery by proximal forward-backward splitting,”  Multiscale Modeling & Simulation, vol. 4, no. 4, pp. 1168–1200, November 2005.
  * [28] F. Agostinelli, M. Hoffman, P. Sadowski, and P. Baldi,  “Learning activation functions to improve deep neural networks,”  in Proc. International Conference on Learning Representations (ICLR) workshop, May 2015.
  * [29] “PyTorch,” https://pytorch.org/.
  * [30] J. S. Garofolo, L. F. Lamel, W. M. Fisher, J. G. Fiscus, D. S. Pallett, N. L. Dahlgren, and V. Zue,  “TIMIT acoustic-phonetic continuous speech corpus,”  Linguistic data consortium, 1993.
  * [31] D. P. Kingma and J. L. Ba,  “Adam: a method for stochastic optimization,”  in Proc. International Conference on Learning Representations (ICLR), May 2015.
  * [32] C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen,  “An algorithm for intelligibility prediction of time–frequency weighted noisy speech,”  IEEE Transactions on Audio, Speech, and Language Processing, vol. 19, no. 7, pp. 2125–2136, February 2011.
  * [33] Y. Zhao, B. Xu, R. Giri, and T. Zhang,  “Perceptually guided speech enhancement using deep neural networks,”  in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), April 2018, pp. 5074–5078.
  * [34] G. Naithani, J. Nikunen, L. Bramslow, and T. Virtanen,  “Deep neural network based speech separation optimizing an objective estimator of intelligibility for low latency applications,”  in Proc. International Workshop on Acoustic Signal Enhancement (IWAENC), September 2018, pp. 386–390.
  * [35] M. Pariente,  “PyTorch implementation of STOI,” https://github.com/mpariente/pytorch_stoi, 2021,  [Online; accessed 13-February-2022].
  * [36] M. Pariente,  “PySTOI,” https://github.com/mpariente/pystoi, 2018.
  * [37] D. Fitzgerald, M. Cranitch, and E. Coyle,  “On the use of the beta divergence for musical source separation,”  in Proc. of IET Irish Signals and Systems Conference (ISSC), June 2008, pp. 1–6.
  * [38] J.-J. Moreau,  “Proximité et dualité dans un espace hilbertien,”  Bulletin de la Société mathématique de France, vol. 93, pp. 273–299, 1965.
  * [39] P. L. Combettes and J.-C. Pesquet,  “Proximal thresholding algorithm for minimization over orthonormal bases,”  SIAM Journal on Optimization, vol. 18, no. 4, pp. 1351–1376, 2008\.
  * [40] R. Gribonval and M. Nikolova,  “A characterization of proximity operators,”  Journal of Mathematical Imaging and Vision, vol. 62, no. 6, pp. 773–789, July 2020.

## Supplementary: Characterization of $\mathsf{F}(\mathbf{y},\mathbf{r})$ as a
proximity operator.

In this supplementary, we address the problem of identifying a function
$f_{\mathbf{r}}$ $:\mathbb{R}^{K}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ such
that
${\mathsf{F}(\mathbf{y},\mathbf{r})=\text{prox}_{f_{\mathbf{r}}}(\mathbf{y})}$,
with $\mathsf{F}$ defined in (19) of the main paper. Note that we ignore here
the layer index $t$ for simplicity.

First, we present in Section A sufficient conditions for the existence of such
a function in a general setting. Then, in Section B we consider more
particularly the case of APL units. Finally, in Section C we recover the
function $f_{\mathbf{r}}$ corresponding to our sublayer $\mathsf{F}$.

###  A Proximity operators

First, let us recall that the proximity operator of a real-valued convex
function $\sigma:\mathbb{R}^{K}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ is
defined as the mapping of a vector $\mathbf{z}$ to the solution of the
minimization problem:

$\text{prox}_{\rho^{-1}\sigma}(\mathbf{z}):=\text{arg}\underset{\mathbf{p}\in\mathbb{R}^{K}}{\text{min}}\sigma(\mathbf{p})+\frac{\rho}{2}\|\mathbf{p}-\mathbf{z}\|^{2}.$
(20)

This definition can be extended to nonconvex functions, resulting in a
possibly set-valued operator. Following [38, 39, 40], a function
$g:\mathcal{Z}\subset\mathbb{R}^{K}\rightarrow\mathbb{R}^{K}$ can be
characterized as the proximity operator of a convex function when it is
strictly increasing over each coordinate. More precisely, the authors in [40]
show that if there exists a convex, lower semi-continuous function
$\tilde{g}:\mathbb{R}^{K}\rightarrow\mathbb{R}$ such that $g(\mathbf{z})$ is a
subgradient of $\tilde{g}(\mathbf{z})$ for all $\mathbf{z}\in\mathcal{Z}$,
i.e.,

$g(\mathbf{z})\in\\{\mathbf{v}\in\mathbb{R}^{K}\,|\,\forall\mathbf{p}\in\mathbb{R}^{K},\,\langle\mathbf{p}-\mathbf{z},\,\mathbf{v}\rangle+\tilde{g}(\mathbf{z})\leq\tilde{g}(\mathbf{p})\\},$
(21)

then there exists $\sigma$ such that
$g(\mathbf{z})\in\text{prox}_{\sigma}(\mathbf{z})$ and the following relation
stands:

$\forall\mathbf{z}\in\mathcal{Z},\,\sigma(g(\mathbf{z}))=\langle\mathbf{z},\,g(\mathbf{z})\rangle-\frac{1}{2}\|g(\mathbf{z})\|^{2}-\tilde{g}(\mathbf{z}).$
(22)

Finally, when $g$ is invertible and with $\mathbf{y}=g(\mathbf{z})$, we have:

$\sigma(\mathbf{y})=\langle
g^{-1}(\mathbf{y}),\,\mathbf{y}\rangle-\frac{1}{2}\|\mathbf{y}\|^{2}-\tilde{g}(g^{-1}(\mathbf{y})).$
(23)

###  B Characterization with APL

We now address the case of (strictly increasing) APL functions as defined in
(18), with negative weights (cf. 3.3) and at least one nonnegative bias
$b_{c}$. Let us consider the following convex, lower semi-continuous function
$\overline{\textsf{APL}}$ such that $\forall z\in\mathbb{R}$:

$\overline{\textsf{APL}}(z)=\frac{z^{2}}{2}\chi_{[0;+\infty]}(z)\\\
+\sum_{c=1}^{C}w_{c}\left(\frac{-z^{2}}{2}+b_{c}z\right)\chi_{]-\infty;b_{c}]}(z),$
(24)

where $\chi_{\Pi}$ is the indicator function of set $\Pi$. Since for any
$z\in\mathbb{R}$, $\textsf{APL}(z)$ is a subgradient of
$\overline{\textsf{APL}}(z)$, and denoting
$\widetilde{\textsf{APL}}(\mathbf{z})=\sum_{k=1}^{K}\overline{\textsf{APL}}(z_{k})$,
it is straightforward to show that the relation (21) stands for
${g(\mathbf{z})=\textsf{APL}(\mathbf{z})}$ and
${\tilde{g}(\mathbf{z})=\widetilde{\textsf{APL}}(\mathbf{z})}$.

Besides, since APL is invertible we can use the relation (23) to identify
$\sigma$:

$\sigma(\mathbf{y})=\langle\textsf{APL}^{-1}(\mathbf{y}),\,\mathbf{y}\rangle-\frac{1}{2}\|\mathbf{y}\|^{2}-\widetilde{\textsf{APL}}(\textsf{APL}^{-1}(\mathbf{y})),$
(25)

with:

$\textsf{APL}^{-1}(\mathbf{y})=\frac{\mathbf{y}-\sum_{c=1}^{C}w_{c}b_{c}\chi_{]-\infty,\,\textsf{APL}(b_{c})]}(\mathbf{y})}{\chi_{[\textsf{APL}(0),\,+\infty[}(\mathbf{y})-\sum_{c=1}^{C}w_{c}\chi_{]-\infty,\,\textsf{APL}(b_{c})]}(\mathbf{y})}.$
(26)

###  C Characterization with $\mathsf{F}$

Finally, let us retrieve
$f_{\mathbf{r}}:\mathbb{R}^{K}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ such
that
${\mathsf{F}(\mathbf{y},\mathbf{r})=\text{prox}_{f_{\mathbf{r}}}(\mathbf{y})}$.
Drawing on the previous section and using the definition of $\mathsf{F}$ from
(19), we have:

$\mathsf{F}(\mathbf{y},\,\mathbf{r})=\text{prox}_{\sigma}\left(\gamma^{(1)}\mathbf{y}+\gamma^{(2)}\frac{\mathbf{r}^{\beta-1}}{\beta-1}\right).$
(27)

To fully identify $f_{\mathbf{r}}$, we first need to reformulate (27) so that
the argument of the right hand side term simply becomes $\mathbf{y}$. To that
end, we leverage a property from [27], which consists in first rewriting (27)
as follows:

$\mathsf{F}(\mathbf{y},\,\mathbf{r})=\text{prox}_{\sigma}\left(\frac{\mathbf{y}-\mathbf{q}}{2\alpha+1}\right),$
(28)

with $\alpha=\displaystyle\frac{1-\gamma^{(1)}}{2\gamma^{(1)}}$ and
$\mathbf{q}=-\displaystyle\frac{\gamma^{(2)}}{\gamma^{(1)}}\displaystyle\frac{\mathbf{r}^{\beta-1}}{\beta-1}$.
The property from [27] then states that:

$\text{prox}_{\varphi+\alpha\|\cdot\|^{2}+\langle\mathbf{q},\,\cdot\rangle}(\mathbf{y})=\text{prox}_{\varphi/(2\alpha+1)}\left(\frac{\mathbf{y}-\mathbf{q}}{2\alpha+1}\right).$
(29)

Let $\varphi=(2\alpha+1)\sigma$. Combining (28) and (29) yields:

$\mathsf{F}(\mathbf{y},\,\mathbf{r})=\text{prox}_{(2\alpha+1)\sigma+\alpha\|\cdot\|^{2}+\langle\mathbf{q},\,\cdot\rangle}(\mathbf{y}).$
(30)

As a result, from (30) we can identify $f_{\mathbf{r}}$ such that its
proximity operator is $\mathsf{F}$. If we further exploit the definition of
$\sigma$ from (25), we finally have:

$f_{\mathbf{r}}(\mathbf{y})=\frac{1}{\gamma^{(1)}}\left\langle\textsf{APL}^{-1}(\mathbf{y})-\gamma^{(2)}\frac{\mathbf{r}^{\beta-1}}{\beta-1},\,\mathbf{y}\right\rangle\\\
-\frac{1}{2}\|\mathbf{y}\|^{2}-\frac{1}{\gamma^{(1)}}\widetilde{\textsf{APL}}(\textsf{APL}^{-1}(\mathbf{y})).$
(31)

Therefore, using (31) one can recover the function associated with the learned
proximity operator, and consequently identify the metric involved in the
formulation of the PR problem.


# Deep convolutional neural networks for cyclic sensor data

1st Payman Goodarzi  Lab for Measurement Technology
Saarland University
Saarbruecken, Germany
<EMAIL_ADDRESS>2rd Yannick Robin  Lab for Measurement
Technology
Saarland University
Saarbruecken, Germany
<EMAIL_ADDRESS>3th Andreas Schütze  Lab for Measurement
Technology
Saarland University
Saarbruecken, Germany
<EMAIL_ADDRESS>4th Tizian Schneider  Lab for Measurement
Technology
Saarland University
Saarbruecken, Germany
<EMAIL_ADDRESS>
###### Abstract

Predictive maintenance plays a critical role in ensuring the uninterrupted
operation of industrial systems and mitigating the potential risks associated
with system failures. This study focuses on sensor-based condition monitoring
and explores the application of deep learning techniques using a hydraulic
system testbed dataset. Our investigation involves comparing the performance
of three models: a baseline model employing conventional methods, a single CNN
model with early sensor fusion, and a two-lane CNN model (2L-CNN) with late
sensor fusion. The baseline model achieves an impressive test error rate of 1%
by employing late sensor fusion, where feature extraction is performed
individually for each sensor. However, the CNN model encounters challenges due
to the diverse sensor characteristics, resulting in an error rate of 20.5%. To
further investigate this issue, we conduct separate training for each sensor
and observe variations in accuracy. Additionally, we evaluate the performance
of the 2L-CNN model, which demonstrates significant improvement by reducing
the error rate by 33% when considering the combination of the least and most
optimal sensors. This study underscores the importance of effectively
addressing the complexities posed by multi-sensor systems in sensor-based
condition monitoring.

###### Index Terms:

Predictive maintenance, hydraulic system, deep learning, convolutional neural
network

##  I Introduction

Industrial systems and factories operate continuously, necessitating
uninterrupted performance to avoid process downtime, significant financial
losses, and potential safety hazards. To mitigate these risks, companies
employ various maintenance approaches, including corrective maintenance,
preventive maintenance, and predictive maintenance. Predictive maintenance
(PdM) heavily relies on monitored signals from diverse sensors, with machine
learning methods playing a pivotal role in data-driven PdM. These methods can
be categorized into two groups: conventional approaches and deep learning
techniques. Conventional methods involve preprocessing, feature extraction
(FE), feature selection (FS), and the subsequent application of classification
or regression algorithms [1], commonly referred to as FESC/FESR in this study.
In contrast, modern deep neural networks have demonstrated exceptional
performance across various applications, including PdM [2, 3, 4].

In line with these advancements, gas mixture measurement has emerged as a
promising application for deep neural networks in recent research. Notably,
Robin et al. [5] introduced a convolutional neural network (CNN) specifically
designed for indoor air quality monitoring (TCOCNN), accurately predicting
volatile organic compounds using temperature-cycled operation sensors. The
proposed method surpassed existing data evaluation techniques, underscoring
the effectiveness of CNNs in this domain. It is worth noting that the signals
utilized in their study bear resemblance to the typical data encountered in
condition monitoring and predictive maintenance applications, i.e. multiple
sensors with periodic or cyclic data.

Motivated by the aforementioned findings, the objective of our study is to
compare our previously published method, TCOCNN, with a benchmark method in a
different application context. To achieve this, we utilize a publicly
available dataset from a hydraulic system testbed [6]. Recent research has
applied various deep learning techniques to the dataset under investigation
[7, 8, 9, 10, 11]. Prakash et al. [7] employed a 1D CNN model to analyze the
pressure difference between two pressure sensors. Huang et al. [12] took a
parallel approach by utilizing multiple independent convolutional neural
networks to extract features from individual sensors. Furthermore, Berghout et
al. [10] introduced a novel neural network model specifically designed to
process the extracted features. In a distinct approach, Zhang et al. [9]
demonstrated the application of a Transformer model with self-attention,
originally trained on natural language, to the task of sensor fusion.
Collectively, these studies contribute to the exploration of diverse
methodologies for analyzing sensor data and extracting meaningful insights.
The primary goal of our comprehensive evaluation is to assess the performance
of the air quality model when applied to the field of condition monitoring. In
doing so, we aim to address potential challenges and difficulties associated
with multiple sensors of different types, thereby providing valuable insights
for future research in this area. The remaining structure of this paper is
outlined as follows: Section II describes the materials and methods, including
details about the dataset and the utilized convolutional neural network.
Section III reports the results, and finally, Section IV presents the
conclusions derived from this study.

##  II Materials and methods

###  II-A Dataset

This study utilizes a dataset that captures the behavior of a hydraulic system
(HS) testbed, which has been specifically designed to simulate various common
faults encountered in such systems [6]. The ZeMA111Zentrum für Mechatronik und
Automatisierungstechnik gemeinnützige GmbH dataset includes simulated faults
such as decreased cooler performance, main valve switching degradation,
internal pump leakage, and accumulator pre-charge pressure reduction with the
control system enabling independent adjustment of each fault condition. Fig.
1a provides an illustration of the conditions of the cooler, valve, pump, and
accumulator within the dataset, which consists of recordings from 17 sensors
over a constant operating cycle lasting 60 seconds as would be typical, e.g.
for a hydraulic press operation. These sensors measure process values,
including pressure (PS1 - PS6), flow (FS1, FS2), temperature (TS1 - TS5),
electrical power (EPS1), and vibration (VS1). Additionally, the dataset
includes three virtual sensors, namely cooling efficiency (CE), cooling power
(CP), and system efficiency (SE). These virtual sensors are calculated using a
physical model that combines various measured values. The dataset comprises
sensors of various types, and their sampling rates vary based on the measured
parameter. The sampling frequencies range from 1 to 100 Hz, resulting in
observations with 60 to 6000 data samples per sensor per cycle. Fig. 1b
effectively showcases the distinct characteristics exhibited by two sensors
through three cycles, highlighting the multimodality of the sensor data [13].
In the context of the present study, the target is to predict the condition of
the accumulator. Specifically, the model aims to classify the hydraulic
accumulator pre-charge pressure into categories, i.e., ”optimal pressure,”
”lightly reduced pressure,” ”severely reduced pressure,” and ”close to total
failure.”

(a)

(b)

Figure 1: ZeMA dataset control variables (a), and three cycles of two selected
sensors (b).

###  II-B Algorithms

Conventional ML: The baseline model for this study is selected from [14],
where an AutoML toolbox is employed to analyze the dataset. The AutoML toolbox
explores various combinations of FESC methods to identify the most effective
approach. In this study, we utilize the method identified by the toolbox,
which achieved the highest cross-validation accuracy. The selected method
involves extracting statistical moments (mean, standard deviation, skewness,
and kurtosis) from the raw data as the features. Pearson correlation is used
as the feature selector, and linear discriminant analysis with Mahalanobis
distance serves as the classifier.

Deep learning methods: Deep learning techniques were employed as the second
approach to construct the models. Specifically, two CNN models were utilized:
TCOCNN and 2L-CNN, as depicted in Figure 2. TCOCNN is a deep network
comprising 10 convolutional layers, and we also evaluated the late fusion
version of TCOCNN, which has demonstrated its effectiveness in various
multimodal systems [15, 12]. In the 2L-CNN model, each sensor was assigned its
own set of convolutional filters. The resulting feature maps from each
convolutional lane were concatenated, and a fully connected layer was then
applied.

Figure 2: Architecture of the TCOCNN and 2L-CNN.

In a previous study [5], the employed model was utilized to predict multiple
gas concentrations, such as acetone, ethanol, formaldehyde, toluene, the total
concentration of all volatile organic compounds (VOCsum), as well as the
inorganic gases carbon monoxide and hydrogen. However, in the present study, a
different approach is adopted as it focuses on a classification task.
Consequently, the last fully connected layer and output function differ from
those used in [5]. To determine the network’s hyperparameters, a
hyperparameter (HP) tuning process is employed. Due to computational resource
limitations, the HP tuning involves conducting 50 iterations within a
predefined search space for the HPs, as detailed in Table I.

TABLE I: The hyperparameter ranges for the CNN. Initial Learning Rate (Log Scale) |  Number of Filters (First Two Layers) |  Kernel Size (First Two Layers) |  Stride Size (First Layer) |  Dropout Number of Neurons (FC) |  Number of Neurons (FC)
---|---|---|---|---|---
$1\times 10^{-7}-1\times 10^{-4}$ |  10-100 |  100-300 |  100-175 |  30-50% |  500-2500

To address the issue of multiple sensors having different sampling rates and,
thus, lengths, a preprocessing step is employed where all sensors are
upsampled to 6000 data points, ensuring the same number of raw data for each
sensor. These upsampled sensor readings are then combined to form a 2D matrix,
enabling the application of a 2D convolutional network. Following the
preprocessing step, a random split is performed on the dataset. The data is
divided into three subsets: 70% of the data is allocated for training, 10% for
validation, and 20% for testing. This partitioning ensures that the model is
trained on a substantial portion of the data while having separate datasets
for validation and testing. The HP tuning is then exclusively performed on the
training data, enabling the identification of optimal HPs based on minimizing
the validation loss. This ensures the model is fine-tuned and optimized for
performance on unseen data.

##  III Results

The results obtained from the baseline model, which utilizes conventional
methods, showcased excellent performance for the defined task. When all 17
sensors were used as input, the model achieved a test error rate of 1%. The
conventional method applies feature extraction to each sensor individually,
thereby avoiding any challenges associated with multiple variant sensors.

In contrast, the CNN model incorporating all 17 sensors and the mentioned
preprocessing steps exhibited lower performance, achieving an error rate of
20.5%. This outcome highlighted the challenges arising from the multimodal
characteristics of the sensors. To highlight this issue, we trained the
network separately for each sensor and observed varying error rates. Fig. 3a
displays the error results, with PS1 achieving the best error rate of 1.7% and
PS2 obtaining the highest error of 73.4%.

(a)

(b)

Figure 3: Test error when only single sensors are used (a), and when
combinations of sensors are used (b). Each HP tuning is repeated 5 times, and
the error bars represent the range for the standard deviation of the results.

To further investigate the impact of dissimilar sensors on the deep neural
network model, we conducted a series of tests using different combinations of
sensors. In order to keep the experiments manageable, we focused on two
sensors and trained the network using all possible combinations involving the
best sensor (PS1) and the remaining sensors. The results are illustrated in
Figure 1, where the most notable findings are highlighted. The worst
combination, involving PS1 and PS5, exhibited a significantly higher error
rate of 39.8% (B&W). Moreover, when comparing the performance of using the
best sensor alone to using it twice (B&B), we observed that increasing the
input size (2x6000) did not have a significant impact on the results. However,
introducing a second sensor with uniform noise (B&N) as an input resulted in a
substantial increase in the error rate, reaching 30.3%.

Fig. 3b provides an illustration of the result achieved by this new model,
denoted as ”2L-CNN,” using the best and worst sensors as input. The new model
demonstrated improved performance compared to the normal single-lane network,
achieving error rates of 39.8% and 6.5% respectively. This finding provides
additional support for the importance of effectively addressing the challenges
associated with dissimilar sensors and multimodal learning in order to
maintain and improve performance.

##  IV Conclusion

In conclusion, this study aimed to assess the performance of a network
originally designed for temperature-cycled gas sensors in a different
application, namely the fault classification of a hydraulic system, and
compare it with conventional methods. The conventional baseline model
demonstrated impressive performance in handling the 17 diverse sensors,
achieving an outstanding error rate of 1%. This method utilizes independent
late data fusion after extracting features individually from each sensor.

In contrast, the CNN model that incorporated all 17 sensors and utilized
preprocessing to achieve the same raw data length demonstrated significantly
lower performance, resulting in an error rate of 20.5%. This outcome
highlighted the challenges arising from the dissimilar characteristics of the
sensors, which hindered the network’s ability to effectively handle them.
Notably, training the network separately for each sensor revealed substantial
variations in error rates, ranging from 1.7% for PS1 to 73.4% for PS2.

To emphasize the importance of pertinent input data, we performed experiments
using various sensor combinations. These tests provided clear evidence that
incorporating irrelevant sensors in the input data significantly compromised
the results. Additionally, we evaluated the performance of a 2L-CNN model that
utilizes late-sensor fusion, which proved to be a viable strategy for
addressing the issues caused by irrelevant sensors. This discovery underscores
the crucial significance of meticulously choosing and prioritizing relevant
sensors to enhance the model’s performance.

Overall, this study underscores the challenges posed by the multimodal nature
of sensor data. It emphasizes the significance of effectively addressing these
challenges to unlock the full potential of sensor-based applications and
enhance their overall performance. Future research endeavors can focus on
exploring advanced techniques to overcome these obstacles and further enhance
the accuracy of fault detection models.

## References

  * [1] T. Schneider, S. Klein, and A. Schütze, “Machine learning in industrial measurement technology for detection of known and unknown faults of equipment and sensors,” _Technisches Messen_ , vol. 86, pp. 706–718, 11 2019.
  * [2] S. Namuduri, B. N. Narayanan, V. S. P. Davuluru, L. Burton, and S. Bhansali, “Review—deep learning methods for sensor based predictive maintenance and future perspectives for electrochemical sensors,” _Journal of The Electrochemical Society_ , vol. 167, p. 037552, 2 2020.
  * [3] R. Zhao, R. Yan, Z. Chen, K. Mao, P. Wang, and R. X. Gao, “Deep learning and its applications to machine health monitoring,” _Mechanical Systems and Signal Processing_ , vol. 115, pp. 213–237, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0888327018303108
  * [4] R. Magar, L. Ghule, J. Li, Y. Zhao, and A. B. Farimani, “Faultnet: A deep convolutional neural network for bearing fault classification,” _IEEE Access_ , vol. 9, pp. 25 189–25 199, 2021.
  * [5] Y. Robin, J. Amann, T. Baur, P. Goodarzi, C. Schultealbert, T. Schneider, and A. Schütze, “High-performance voc quantification for iaq monitoring using advanced sensor systems and deep learning,” _Atmosphere_ , vol. 12, p. 1487, 11 2021.
  * [6] N. Helwig, E. Pignanelli, and A. Schutze, “Condition monitoring of a complex hydraulic system using multivariate statistics.” IEEE International Instrumentation and Measurement Technology Conference (I2MTC) Proceedings, 5 2015, pp. 210–215.
  * [7] J. Prakash, S. Singh, A. Miglani, and P. K. Kankar, “Pressure signal-based analysis of anomalies in switching behavior of a two-way directional control valve,” _ASME Open Journal of Engineering_ , vol. 2, 1 2023.
  * [8] S. Pillai and P. Vadakkepat, “Deep learning for machine health prognostics using kernel-based feature transformation,” _Journal of Intelligent Manufacturing_ , vol. 33, pp. 1665–1680, 8 2022.
  * [9] Z. Zhang, M. Farnsworth, B. Song, D. Tiwari, and A. Tiwari, “Deep transfer learning with self-attention for industry sensor fusion tasks,” _IEEE Sensors Journal_ , vol. 22, pp. 15 235–15 247, 8 2022.
  * [10] T. Berghout, M. Benbouzid, S. M. Muyeen, T. Bentrcia, and L.-H. Mouss, “Auto-nahl: A neural network approach for condition-based maintenance of complex industrial systems,” _IEEE Access_ , vol. 9, pp. 152 829–152 840, 2021.
  * [11] X. Ma, P. Wang, B. Zhang, and M. Sun, “A multirate sensor information fusion strategy for multitask fault diagnosis based on convolutional neural network,” _Journal of Sensors_ , vol. 2021, pp. 1–17, 6 2021.
  * [12] K. Huang, S. Wu, F. Li, C. Yang, and W. Gui, “Fault diagnosis of hydraulic systems based on deep learning model with multirate data samples,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 33, pp. 6789–6801, 11 2022.
  * [13] T. Baltrusaitis, C. Ahuja, and L.-P. Morency, “Multimodal machine learning: A survey and taxonomy,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 41, pp. 423–443, 2 2019.
  * [14] T. Schneider, N. Helwig, and A. Schütze, “Industrial condition monitoring with smart sensors using automated feature extraction and selection,” _Measurement Science and Technology_ , vol. 29, p. 94002, 8 2018. [Online]. Available: https://doi.org/10.1088/1361-6501/aad1d4
  * [15] H. R. V. Joze, A. Shaban, M. L. Iuzzolino, and K. Koishida, “Mmtm: Multimodal transfer module for cnn fusion.” Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 13 289–13 299.