Out-of-distribution detection using normalizing flows on the data manifold

Seyedeh Fatemeh Razavia (razavi_f@ut.ac.ir), Mohammad Mahdi Mehmanchia
(mahdi.mehmanchi@ut.ac.ir), Reshad Hosseinia (reshad.hosseini@ut.ac.ir),
Mostafa Tavassolipoura<EMAIL_ADDRESS>
a School of ECE, College of Engineering, University of Tehran, Tehran, Iran

Corresponding Author:
Reshad Hosseini
School of ECE, College of Engineering, University of Tehran, Tehran, Iran
Tel: (+98) 21-82085099
Email<EMAIL_ADDRESS>
###### Abstract

One possible solution to detect out-of-distribution data is by fitting a
distribution and using the exact likelihood value as a measurement.
Normalizing flows are likelihood-based generative models providing a tractable
density estimation via dimension-preserving invertible transformations.
Conventional normalizing flows are prone to fail in out-of-distribution
detection, because of the well-known curse of dimensionality problem of the
likelihood-based models. According to the manifold hypothesis, real-world data
often lie on a low-dimensional manifold. Using normalizing flows that jointly
estimate the density and the manifold of data, this paper proposes a new
criterion for out-of-distribution detection by combining the log-likelihood of
normalizing flows for given data and the distance from the manifold. However,
despite the successful combination of these two criteria, each one on its own
is ineffective for out-of-distribution detection. The extensive experimental
results show that the manifold learning improves out-of-distribution detection
without modifying the model structure or using auxiliary out-of-distribution
data during training. The extensive experimental results show that manifold
learning improves the out-of-distribution detection ability of a class of
likelihood-based models known as normalizing flows. This improvement is
achieved without modifying the model structure or using auxiliary out-of-
distribution data during training.

###### keywords:

Generative models , Manifold learning , Normalizing flows , Out-of-
distribution

††journal: Expert Systems with Applications

Abbreviations:  AR, Auto-Regressive; AUROC, Area Under the Receiver Operating
Characteristic Curve; BPD, Bit-Per-Dim; DNF, Denoising Normalizing Flow; MSE,
Mean Squared Error; NF, Normalizing Flow; NLL, Negative Log-Likelihood; OOD,
Out-Of-Distribution; IC, Input Complexity; ID, In-Distribution; VAE,
Variational Auto-Encoder.

##  1 Introduction

Out-Of-Distribution (OOD) detection classifies test data into in-/out- of
distribution data [yang2021generalized]. Generative models that rely on the
likelihood estimation are a prominent candidate for detecting OOD data. It is
expected they can assign low likelihood value to OOD data. However, high-
dimensional likelihood-based generative models are susceptible to failure in
OOD detection. However, high-dimensional likelihood-based generative models
are susceptible to failure in OOD detection due to the tendency of the
likelihood to be extremely small across significant parts of the space.

Although Normalizing Flows (NFs) as likelihood-based generative models with
tractable likelihood appear to be well-suited for addressing the OOD detection
problem, their likelihood does not provide a significant distinction for OOD
data. In fact, NFs tend to assign high likelihood to OOD data
[nalisnick2018deep, kirichenko2020normalizing]. Therefore, it has been stated
that they learn local transformations instead of semantics. This problem is
fundamental because the cost function of NFs (Negative Log-Likelihood (NLL))
for high dimensional data tends to cover the entire space regardless of its
low-/high- density regions due to the assignment of equal weights to samples
theis2015note. Minimizing the NLL (or equivalently minimizing the Kullback-
Leibler divergence) over the entire high-dimensional space is challenging due
to the assignment of equal weights to samples theis2015note. Using alternative
measurement instead of pure likelihood is a candidate solution followed by
researchers such as likelihood ratio [ren2019likelihood], likelihood regret
[xiao2020likelihood], and Input Complexity (IC) [Serrà2020Input]. There are
some other alternative measurements instead of pure likelihood such as
likelihood ratio [ren2019likelihood], likelihood regret [xiao2020likelihood],
and Input Complexity (IC) [Serrà2020Input]. To the best of our knowledge, the
role of manifold learning to tackle the aforementioned problem has not yet
been studied for NFs.

The important structural limitation of common NFs is that they cannot learn
the embedded manifolds of the data. Real data are usually embedded in low-
dimension manifolds, and powerful generating models use this intuition
[kingma2013auto, goodfellow2014generative]. NFs use bijective transformation,
and therefore they preserve dimensionality between the input and transformed
spaces. One can think that this can be easily solved by using non-bijective
transformation that can transform the data space to a lower dimensional one,
and maximizing the likelihood of the data on the low-dimensional manifold is
obtained by such a transformation. But unfortunately, this optimization
problem cannot be solved exactly. Some other well-known likelihood-based
methods on the manifold also approximate a lower-bound of the likelihood
[kingma2013auto].

Recently, several researchers proposed solutions to the problem of maximizing
the likelihood of the data on a manifold using NFs. Some use an injective
transformation in NFs [brehmer2020flows, caterini2021rectangular,
cunningham2020normalizing]. By using injective transformations, non-square
Jacobian matrices appear in the likelihood term making the optimization
computationally expensive. Using injective transformations makes the
optimization computationally expensive. brehmer2020flows used two-step
training procedure, first the manifold is learned using an NF, then the
density is estimated using another NF. This procedure simplifies the training,
but it can lead to poor density estimation. Learning only the manifold and not
the Jacobian of the initial NF during the first training step can lead to
unexpected alterations in the space deformation by NF, potentially making
density estimation challenging during the second phase. horvat2021denoising
introduced a method that learns a low-dimensional manifold from NFs in a
single-step training procedure without using injective transformations. The
transformed space of the NF is split into two parts, manifold and noise,
modeled by another NF and a low-variance Gaussian distribution, respectively.
They also add this low-variance noise to the input.

Inspired by the previously mentioned methods, we provide an approach to
overcome the aforementioned inherent shortcoming of NFs as a likelihood-based
generative model to estimate a manifold, particularly in high dimensions.
Then, we employ the model for OOD detection. Fig. LABEL:fig:idea is a toy
example that illustrates the intuition behind our proposed method. Suppose
that we have a semicircle dataset with unit radius where the data points are
generated according to a truncated Gaussian distribution
$\mathcal{N}(\mu=\frac{\pi}{2},\sigma=0.85,a=\frac{-\pi}{1.7},b=\frac{\pi}{1.7})$.
The data in this two-dimensional representation clearly resides on a one-
dimensional manifold (i.e., a straight finite line). To find this line, for
every $x\in\mathbb{R}$ and $y\in\mathbb{R}^{+}$, we can project point $(x,y)$
to $\theta$, where $0\leq\theta\leq\pi$. By using the equations
$x^{\prime}=\cos(\theta)$ and $y^{\prime}=\sin(\theta)$, we can then transform
this line back to the semicircle and obtain the reconstructed point
$(x^{\prime},y^{\prime})$. If the initial point $(x,y)$ lies on the arc, then
$(x^{\prime},y^{\prime})=(x,y)$. Otherwise, $(x,y)$ will be mapped to the
nearest point $(x^{\prime},y^{\prime})$ on the arc. So, we can define a
reconstruction loss between each point and its reconstructed point based on a
distance function. Fig. LABEL:fig:rec shows the reconstruction loss for points
on a grid, where we used the Mean Squared Error (MSE) loss (square Euclidean
distance). The NLL on the manifold for points on the grid is shown in Fig.
LABEL:fig:nll. Now consider points $A$ and $C$. Point $A$ has low likelihood
and point $C$ has high likelihood (see Fig. LABEL:fig:nll). If we only
consider the likelihood on the manifold, we may incorrectly consider $A$ an
OOD sample and $C$ an In-Distribution (ID) sample. To address this issue, we
add the reconstruction loss to the NLL term with a positive constant factor.
This way, our OOD score includes both the NLL and the reconstruction loss. We
have presented the proposed OOD score in Fig. LABEL:fig:score. Based on the
provided example, it may be tempting to use only the reconstruction loss for
OOD detection. However, we should consider points $B$ and $D$ to fully
understand the importance of the likelihood term. Although points $B$ and $D$
have the same reconstruction loss, point $D$ is near a region with a high
density, unlike point $B$. Therefore, it is reasonable to consider point $D$
as an ID sample and point $B$ as an OOD sample. In other words, we can
tolerate a bigger reconstruction loss for samples near high-density regions
compared to those near low-density regions. This intuition can be clearly
observed in Fig. LABEL:fig:score.

To jointly estimate the density and learn the manifold, the objective would be
to optimize the likelihood within the on-manifold region while imposing
penalties on the off-manifold region. To achieve this, a penalty function is
employed, which encourages the inverse mapping to reconstruct the original
data only based on the on-manifold subspace. We combine the achieved NLL and
the manifold score (reconstruction loss) to create an OOD detection score. We
find this approach promising as it enhances the performance of common NFs in
OOD detection tasks without increasing architectural complexity or using
auxiliary OOD training data. Accordingly, our main contributions are
summarized as follows:

  1. 1.

We explore the impact of incorporating manifold learning into NFs for OOD
detection and demonstrate that the use of NFs with manifold learning can
enhance OOD detection in some cases.

  2. 2.

While manifold learning with NFs algorithms typically employs the MSE penalty
function, we also investigate the impact of using the Huber penalty function
[huber] during the manifold learning phase on OOD detection.

The paper is arranged as follows. Section 2 reviews the recently proposed
methods related to manifold learning in NFs and OOD detection. The
preliminaries are introduced in Section 3. Section 4 introduces the proposed
method. Experimental results are presented in Section 5. In the end, the
conclusion is discussed in Section 6.

##  2 Related Work

Studies on manifold learning in NFs methods and OOD detection methods are
provided in Section 2.1 and Section 2.2, respectively.

###  2.1 Manifold learning in NFs

This section provides a brief overview of the most relevant research about the
manifold learning in NFs. In general, the methods are categorized into two
classes: NFs on prescribed manifolds and NFs on learnable manifolds. The first
class focuses on learning flows on prescribed manifolds, while the second one
includes learning the manifold, which aligns closely with our study. We
discuss the most related research in the second category in the next.

Recently, solutions to the problem of manifold learning and density estimation
with NFs have been proposed [kim2020softflow, kothari2021trumpets,
cunningham2020normalizing, brehmer2020flows, caterini2021rectangular,
horvat2021denoising, kalatzis2021multi, ross2021conformal]. In the following,
we review the most related works to our research and highlight their key
properties. A proposed scientific path is about separating manifold learning
and density estimation. The pioneer of this study in NF is a method named
$\mathcal{M}$-flow [brehmer2020flows]. In $\mathcal{M}$-flow, an injective NF
transformation is used to transform the data into a manifold. After that,
another NF is used to estimate the density on the manifold. Manifold learning
and density estimation are done separately to avoid the computation of
Jacobian term induced by the injective transformation (two-phase training). An
extension of this study called multi-chart flow [kalatzis2021multi], which
employs multiple mapping instead of one to find a manifold with multiple
charts. It also suffers from two-phase training. Another followed research of
$\mathcal{M}$-flow, named Rectangular flow [caterini2021rectangular],
overcomes the calculation of injective transformation by relying on automatic
differentiation and linear algebra tricks. A recent study, named Denoising NF
(DNF) [horvat2021denoising], overcomes the limitations of using injective
transformations and separating model training by splitting the transformed
space of an NF into two parts, noise-insensitive and noise-sensitive. These
parts are modeled by another NF and a low-variance Gaussian distribution,
respectively. The noise is also added to the input training data. Using this
structure, two-phase training is no longer needed and the two NFs are trained
simultaneously.

###  2.2 Out-of-distribution detection

A brief overview of recent and prominent developments in OOD detection is
provided in the following. It should be noted that while there is a wide range
of topics for this discussion [yang2021generalized], the main focus of the
current research is OOD detection through density estimation.

##### Non-density-based methods

Contrary to the main focus of this paper, OOD detection methods are not
limited to density-based models. Using data labels can be a guide for this
problem. Applying an appropriate threshold on a pre-trained classifier,
decreasing the model’s overconfidence, changing the training approach, and
using an auxiliary OOD training data are common approaches [lee2017training,
hendrycks2018deep, devries2018learning, liang2018enhancing].

lee2017training proposes a pre-training algorithm by using generated auxiliary
OOD training data. Since classifiers are not trained for the purpose of OOD
detection, but instead to increase accuracy in inference time, a targeted pre-
training approach can help the model enrich its features and decrease the
model’s confidence in OOD data. Same to the mentioned research, one of the
related researches (called outlier exposure) defines a specific loss function
to discriminate between ID data and an auxiliary OOD data with a ranking loss
to enhance the pre-training approach [hendrycks2018deep]. Confidence
estimation along with the main objective function [devries2018learning] or
using a temperature-based Softmax function [liang2018enhancing] are other
solutions pursued in the literature.

##### Density-based methods

It seems likelihood-based generative models can be good candidates for OOD
detection by assigning less probability to OOD data. However, a frequent point
about likelihood-based models such as Auto-Regressive (AR) models
[ren2019likelihood], Variational Auto-Encoders (VAEs) [xiao2020likelihood],
and NFs [nalisnick2018deep, kirichenko2020normalizing] is that the likelihood
is not a distinguishing score for OOD data [UnderstandingBackground,
UnderstandingFailure]. The common weakness of all the mentioned models is that
they also learn the background and irrelevant information. In other words,
they equally attend to whole space.

Using the likelihood ratio instead of pure likelihood in AR models is a well-
known solution that was proposed for the first time to solve the sequential
genomics OOD detection [ren2019likelihood]. It considers a background/semantic
decomposition of data. The ratio score is obtained by dividing the likelihood
value of the model by the likelihood value of a background model. The
background model is trained by perturbing data to learn the background
information. The defined ratio implicitly ignores irrelevant information.
Despite the success of this model in genomic OOD detection, it has not
achieved state-of-the-art results on image datasets. A remarkable solution for
VAE is named likelihood regret [xiao2020likelihood]. This method is based on
the principle that OOD data can drastically shift the likelihood after several
fine-tuning steps. The authors employ a VAE to avoid overfitting in fine-
tuning steps with freezing the decoder block.

According to the literature, NFs only learn local transformations rather than
semantics [nalisnick2018deep, kirichenko2020normalizing]. Therefore, despite
having tractable likelihood, they suffer from shortcomings in OOD detection.
Previously proposed methods try to integrate high-level information in any way
to avoid the OOD fault. At first, the results of training NFs on the embedded
data instead of raw one confirm this claim [kirichenko2020normalizing].
Moreover, improving OOD detection can be achieved by redesigning the model’s
architecture. Using attention units [kumar2021inflow] or deep residual flows
[zisselman2020deep] can be mentioned as pioneer works. Same to non-density-
based methods, changing training strategies (such as using information theory
[ardizzone2020training]) can be effective for OOD detection. Employing a test-
time OOD score based on input complexity is a valuable measurement
[Serrà2020Input]. It is considered as a likelihood ratio between the learned
density and a universal lossless compressor. Simple data are coded in fewer
bits. So, the ratio between the trained model and the compressor achieves a
likelihood ratio to reach the actual number of bits. A line of research
overcomes this limitation by ensembling models [choi2018waic]. Surprisingly,
this approach was successful despite its over-parametrization. Although every
density estimation model is not able to make an accurate diagnosis, but the
mentioned research shows that the combination of such models in aggregate form
is successful contrary to expectations.

Best of our knowledge, investigating the effect of manifold learning in NFs on
OOD detection has not been covered in the literature.

##  3 Preliminaries

This section introduces our notations and related preliminary works to make it
easier for the reader to follow the subject. The rest of this section is
arranged as follows: At first, the standard normalizing flow is discussed in
Section 3.1. After that, a robust reconstruction loss function is presented in
Section 3.2.

###  3.1 Normalizing Flow

There are a variety of well-known deep likelihood-based generative methods,
like VAEs [kingma2013auto], NFs [rezende2015variational], AR models
[murphy2022probabilistic], energy base models [murphy2022probabilistic], and
diffusion models [pmlr-v37-sohl-dickstein15]. But among the mentioned models,
only ARs and NFs can exactly compute the likelihood. VAEs and diffusion models
find a lower bound on the likelihood, while energy base models approximate it.
Sampling in common AR models is computationally expensive, due to sequential
nature of these models. Sampling in NFs is not sequential, but they have a
structural limitation (dimension-preserving) that limits their applicability
and generating power.

NF is a parametric diffeomorphism transformation,
$f_{\phi}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$. So, it is a two-side
differentiable bijective transformation. By choosing a random variable from a
pre-defined distribution in $Z\in\mathbb{R}^{D}$ and transferring back it by
$f_{\phi}^{-1}$, the distribution of data can be found by a change-of-variable
formula like

$p_{X}(X)=p_{Z}(f_{\phi}(X))|\det J_{f_{\phi}}(f_{\phi}(X))|,$ (1)

where $J_{f_{\phi}}\in\mathbb{R}^{D\times D}$ is the Jacobian matrix of the
transformation $f_{\phi}$. It is clear that the practical computability
constraint for this term is the implicit constraint imposed on NFs. Moreover,
the model parameters are estimated by the NLL criterion

$\phi^{*}=\arg\underset{\phi}{\min}\big{(}-\log p_{X}(x)\big{)},$ (2)

where $x=\\{x_{n}\\}_{n=1}^{n=N}$ is the available training data from the
distribution $p_{X}$.

###  3.2 Reconstruction loss functions

During the training of NFs on a manifold, all off-manifold data have a
likelihood of 0. In the literature, penalizing the off-manifold part through a
quadratic reconstruction function, typically MSE, is a common approach to
addressing this problem [brehmer2020flows, caterini2021rectangular,
horvat2021denoising, ross2021conformal]. A weakness of quadratic functions is
its sensitivity to off-manifold data points. However, a regularized
reconstruction function utilizes both linear and quadratic functions, with the
function type changing depending on a comparison between the error value and a
threshold $\delta$. Accordingly, this switching leads the off-manifolds are
not penalized as much as on-manifold data.

Eq. 3 defines a switching function named the Huber function [huber] with a
strong background in robust statistics.

$H_{\delta}(x,y)=\begin{cases}0.5(x-y)^{2},&\text{if }|x-y|<\delta\\\
\delta(|x-y|-0.5\delta),&\text{otherwise}\\\ \end{cases}$ (3)

##  4 Proposed method

By considering a standard NF like $f_{\phi}$, the density of data is estimated
by Eq. 1. It would be exciting if we could calculate the density of data on
the manifold and measure its impact on OOD detection. However, density
estimation on manifolds is not easy. In this paper, we take inspiration from
existing methods of manifold learning in NFs ($\mathcal{M}$-flow
[brehmer2020flows] and DNF [horvat2021denoising]) and propose a new one to
detect OOD data.

From the manifold perspective, the transformed (or latent) space ($z$) can be
disentangled into two spaces, on-manifold ($u$), and off-manifold ($v$) space.
This is because placing all real-value data on a manifold is not necessarily a
correct assumption, and some of the data is generally outside the manifold.
The distribution of the transformed data is a joint distribution of these two
parts like $p_{Z}(z)=p_{Z}(u,v)$. With assuming independence among two spaces
$u\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}v$, the distribution of the transformed space
is decomposed into the product of two sub-space’s distributions as
$p_{Z}(z)=p_{U}(u)\times p_{V}(v)$, where
$U:\mathcal{M}_{M}\subseteq\mathbb{R}^{d}$ and
$V:\mathcal{M}_{O}\subseteq\mathbb{R}^{D-d}$ are the on-manifold and the off-
manifold sub-spaces, respectively.

Our goal is to rearrange the transformed space of an NF in such a manner that
a portion corresponds to the on-manifold part, while the other corresponds to
the off-manifold part. Let $f_{\phi}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$
be a standard NF. We denote the first $d$ components of the output of
$f_{\phi}$ as $u$ and the remaining ones as $v$ where $u$ stands for an on-
manifold part of data and $v$ for an off-manifold part. Formally,
$z=f_{\phi}(x)=(z_{1},z_{2},...,z_{D}),u=(z_{1},z_{2},...,z_{d}),v=(z_{d+1},z_{d+2},...,z_{D}),$
where $x$ is the input, and $z$ is the corresponding transformed data by the
NF. We choose $p_{Z}(z)=\mathcal{N}(0,I_{D})$ and factorize $p_{Z}(z)$ such
that $p_{Z}(z)=p_{U}(u)\times p_{V}(v)$.

It should be mentioned that we can use another NF $h_{\theta}$ to fit a more
complex distribution on the on-manifold part ($u$) rather than a Gaussian
distribution. This idea was followed by DNF [horvat2021denoising] and
$\mathcal{M}$-flow [brehmer2020flows] for density estimation. Briefly, in this
case, $p_{U}(u)$ is replaced by $p_{U^{\prime}}(u^{\prime})|\det
G_{h_{\theta}}(u^{\prime})|$, where $u^{\prime}$ is the transformed space by
$h_{\theta}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, and
$p_{U^{\prime}}(u^{\prime})=\mathcal{N}(0,I_{d})$. Moreover, DNF fits a tight
Gaussian distribution $p_{V}\sim\mathcal{N}(0,\epsilon I_{D-d})$ on the off-
manifold part based on the added noise value $\epsilon$ to the data. However,
the influence of the added noise in DNF was not reported to have a significant
effect. This is likely because the model is trained on real-valued data that
have inherent quantization noise, making the addition of more noise
unnecessary. Consequently, we do not add noise to data in our framework, and
so there are some minor differences with DNF.

The mentioned changes on standard NFs are not sufficient for manifold
learning. It is necessary to penalize the off-manifold sub-space so that the
model is able to reconstruct the data from the on-manifold sub-space.
Accordingly, a constrained optimization problem on NLL appears like

$\underset{\phi}{\min}\;\textrm{NLL}(x;\phi)\;\textrm{s.t.}\;\mathcal{C}(x,\tilde{x})\leq\tau,$
(4)

where $x=(x_{1},x_{2},...,x_{D})$,
$\tilde{x}=(\tilde{x}_{1},\tilde{x}_{2},...,\tilde{x}_{D})=f_{\phi}^{-1}(\text{proj}(f_{\phi}(x)))$,
$\mathcal{C}$, and $\tau$ are the input data, corresponding reconstruction
(computed through the inverse of normalizing flow $f_{\phi}$), the penalty
function, and the penalization threshold, respectively. Moreover,
$\text{proj}(f_{\phi}(x))$ is the first $d$ components of $f_{\phi}(x)$
corresponding to the data manifold, padded with zero as
$\text{proj}(f_{\phi}(x))=(u,\vec{0}_{D-d})$.

Existing methods such as $\mathcal{M}$-flow [brehmer2020flows], DNF
[horvat2021denoising], and Rectangular flow [caterini2021rectangular] penalize
the off-manifold with a quadratic constraint. An important characteristic of a
quadratic function is its sensitivity to outliers, that can be interpreted as
off-manifold parts of the data. One aspect of our contribution involves using
a switching penalization function, containing both linear and quadratic forms,
which is determined based on the distance to the manifold. The switching is
modeled with a Huber function (Eq. (3)). In the proposed method, an element-
wise Huber function with threshold $\delta$, denoted $H_{\delta}$, is applied
to the difference between the input data ($x$) and the reconstructed one
($\tilde{x}$). The reconstructed data is a function of manifold
($\tilde{x}=f_{\phi}^{-1}(\text{proj}(f_{\phi}(x)))$) rather than joint on-
manifold and off-manifold parts. Therefore, the penalization term is computed
by averaging an element-wise function like

$\mathcal{C}(x,\tilde{x};\delta)=\frac{1}{D}\sum_{i=1}^{D}H_{\delta}(|x_{i}-\tilde{x}_{i}|,0).$
(5)

Inspired by the Lagrange multipliers method, the proposed constrained
optimization problem can be converted to its unconstrained equivalent one

$\ell=-\log P_{u}(u)-\log P_{v}(v)+\log|\det
J_{f_{\phi}}(x)|+\lambda\mathcal{C}(x,\tilde{x};\delta),$ (6)

where $\lambda$ is the penalization term hyper-parameter, and the first three
terms are the objective of a standard NF (called NLL).

The proposed approach not only facilitates the estimation of manifold density
for NFs but also addresses the limitations of relying solely on likelihood or
reconstruction for OOD detection, as extensively explained in Section 1 (see
Fig. LABEL:fig:idea). In other words, high likelihood or low reconstruction
error alone are not reliable indicators for identifying ID data. However, by
combining these two factors using the introduced cost function (Eq. 6), we not
only estimate the manifold using NFs but also introduce an appropriate
indicator for OOD detection that is validated through comprehensive
experiments.

##  5 Results

The primary goal of the experiments is to explore the significance of manifold
learning in NFs for enhancing OOD detection. Moreover, the performance of the
proposed method is also evaluated in terms of image generation because it is
fundamentally a generative model. As previously mentioned, $\mathcal{M}$-flow
[brehmer2020flows] and DNF [horvat2021denoising] are existing methods that
closely resemble our approach for manifold learning in NFs. Our proposed OOD
detection score, combines the precise likelihood value with the reconstruction
loss. However, $\mathcal{M}$-flow cannot effectively achieve high-dimensional
likelihood due to using injective transformations. Therefore, our focus is on
comparing our approach with the DNF framework in terms of image generation
performance. Subsequently, we compare our OOD detection results with some
other state-of-the-art likelihood-based OOD detection methods such as
likelihood ratio [ren2019likelihood] and IC [Serrà2020Input]. All in all, we
observed that penalizing the learned manifold in the proposed method leads to
consider the dominant manifold as the ID manifold and the remains (like less
frequent details) as OOD regions.

###  5.1 Experimental setting

The used dataset and experimental setting are described in Section 5.1.1 and
Section 5.1.2, respectively.

####  5.1.1 Dataset

The seven used datasets in the current study are introduced in Table 1. We
could increase the number of layers and parameters in proportion to the
original image’s dimensions. Still, due to resource allocation limitations, we
resize the color images to $32\times 32$ in our experiments.

Table 1: The information of the used datasets in the paper.

####  5.1.2 Architecture

For two-blocks models such as standard DNF or architecture searching
experiments, we use the Glow model [kingma2018glow] for $f_{\phi}$ and the
RealNVP model [dinh2016density] for $h_{\theta}$. Meanwhile, for one-block
models, we use either one of the models. Our Glow model has three blocks, each
with 32 steps of the flow. Each step consists of an actnorm layer, an
invertible $1\times 1$ convolution layer, and an affine coupling layer (see
kingma2018glow for the details of each layer). There is also a squeeze layer
before each step and a split layer after each step. The RealNVP model
comprises six affine coupling layers and six masking layers before each affine
coupling layer. While DNF has chosen a specific architecture in its
experiments, the employed architectures in this paper are based on two well-
known NFs, Glow and RealNVP. Our main goal is to show the importance of
manifold learning in NFs for OOD detection, regardless of the specific
architecture.

####  5.1.3 Training

All implementations are done in the PyTorch111https://pytorch.org/(version
1.10.0+cu102) framework. We use Adam [adam] as our optimization algorithm with
a learning rate of 1e-5 and a batch size of 64. Also, we set the value of the
penalization term coefficient ($\lambda$) to 100. All training experiments
were done on a single GPU (GTX 1080 Ti), with Cuda version 11.4 and 11019GB
memory. The implementation codes of the proposed methods are attached to this
file.

###  5.2 Architecture searching

In this section, we consider four structures: DNF, MSE-based off-manifold
penalization with two blocks (called Proposed-M), Huber-based off-manifold
penalization with two blocks (called Proposed-H), and Huber-based off-manifold
penalization with one block (called Proposed-D). In the first three, two NFs
(named $f_{\phi}$ and $h_{\theta}$, Glow followed by RealNVP) are used, and in
the last structure, a single NF (named $f_{\phi}$) is used. Briefly, we employ
two blocks together as an NF to be more comparable in structure with DNF.
Therefore, the two-blocks method means using a dimension preserving NF
(RealNVP in our setting) to estimate the density of the manifold part (same as
DNF) instead of the Gaussian distribution. We choose all DNF backbone models
the same as ours in all experiments for a fair comparison.

We only present the results of CelebA here and the rest are available in the
supplementary files. Based on the reported results for CelebA dataset by
horvat2021denoising (DNF), we consider $d=|\mathcal{M}_{M}|\sim 500$ for
architecture searching. It is worth noting that, we assume prior knowledge of
the manifold’s intrinsic dimension from existing works in the field. Our goal
is not to determine the intrinsic dimension of the manifold. Readers seeking
dimensionality determination can refer to NEURIPS2022_4f918fa3 and
DBLP:conf/iclr/PopeZAGG21

The results of experiment on CelebA dataset are shown in Table 2 and Fig.
LABEL:fig:figure_face_manifold_CelebA. The main goal of the designed
experiment is to find an appropriate architecture for generating data from the
manifold as long as it does not decrease the likelihood of data. The reported
results in Table 2 confirm that our proposed methods (Proposed-H, Proposed-M,
and Proposed-D) outperform DNF in terms of Bit-Per-Dim (BPD) (calculated by
$\frac{\text{NLL}}{D\times\log 2}$), whereas their generated images look
visually similar, according to Fig. LABEL:fig:figure_face_manifold_CelebA. An
important aspect of DNF is its reliance on incorporating noise during
training. It is noteworthy that the reported inference BPD for the DNF is
computed using noiseless data, and making it an approximation. It is worth
mentioning that best generated images in horvat2021denoising (DNF) are based
on a well-defined StyleGAN manifold [DBLP:conf/cvpr/KarrasLAHLA20], not on raw
data. Additionally, their reported generated images on the CelebA dataset are
not visually perfect.

Table 2: The best MSE/BPD scores for data lie in
$\mathcal{M}\subset\mathbb{R}^{500}$. Dimension changing of the one-block
(Proposed-D) and two-block (DNF, Proposed-M, Proposed-H) methods are
$\mathbb{R}^{3072}\rightarrow\mathbb{R}^{500}$ and
$\mathbb{R}^{3072}\rightarrow\mathbb{R}^{500}\rightarrow\mathbb{R}^{500}$.
#Params means the number of trainable parameters.

Moreover, since the Huber function applies a linear function to penalize the
off-manifold pixels rather than a quadratic one, we anticipated that the MSE
would be higher in the models that were penalized with the Huber function. In
addition to the aforementioned factors, the integration of reconstruction loss
and BPD value is a crucial component in our study. Although the reconstruction
loss may slightly lag behind the DNF method in Table 2, our approach does not
just depend on reconstruction or BPD (as scaled NLL) for OOD detection. Based
on the evaluation of the results in the latter three columns, we can
confidently assert that the proposed method offers a promising advantage in
this regard. Another strength is that the Proposed-D method achieves same
images (in terms of visual quality in Fig.
LABEL:fig:figure_face_manifold_CelebA) to DNF and other two-blocks proposed
methods, with fewer parameters while maintaining the likelihood. Considering
the outcomes, the proposed-D model is chosen for further experiments in the
following. To summarize, we use the term "proposed" to denote the model
referred to as "proposed-D".

By choosing Proposed-D as the base model, more experiments are presented here
to evaluate the performance of the model in terms of generation and manifold
learning on ID/OOD data. The order of experiments is in such a way that the
ability of the models is measured for different ID/OOD data and different
manifold dimensions $(10,50,100,500,1000)$. The results for CelebA test data
are presented here. The results of other datasets are available in the
supplementary file.

Fig. LABEL:fig:figure_CelebA_generation_reconstruction contains the generated
and reconstructed data (learned manifold) for a model trained on CelebA. In
low-dimensional manifolds, images generated by the model are limited to
displaying the main manifold of the face and the background details are not
generated. However, as the manifold dimensions increase, the inclusion of less
frequent details (mostly related to the background) in the generated image is
observed. If the data manifold has simplistic visual attributes (e.g., as
observed in the SVHN dataset results presented in the supplementary file), the
likelihood can be better fit to the data, and thus sharpened images can be
generated as manifold dimension increases. This trend is also observed in more
complex datasets such as CIFAR10 (as noted in the supplementary file) when the
likelihood is fitted using various dimensions of the manifold.

The second row of Fig. LABEL:fig:figure_CelebA_generation_reconstruction
contains the reconstruction of CelebA test data on the model trained with
CelebA for the corresponding manifold dimension. The reconstruction ability is
improved by increasing the dimension. As it is clear, the model’s
invertibility power will appear by increasing manifold dimension, the same as
standard NF. The model trained on CelebA not only excels at generating and
reconstructing ID data but also demonstrates great performance in
discriminating OOD data. In case of reconstructiong OOD data in third and
forth rows of Fig. LABEL:fig:figure_CelebA_generation_reconstruction,
especially from low dimensions, it cannot reconstruct OOD data well. In other
words, learning low dimension manifold for ID data leads to focusing on ID
specific features instead of global features. It is worth noting that our
point is about manifold learning while in case of dimension preserving (like
standard NFs), the reconstruction is done without error.

###  5.3 Image OOD detection

In this section, we present the results of OOD detection for gray-scale
images, and we also expand our analysis to include high-dimensional color
images and hyper-parameter ($d$ and $\delta$) search.

At first, the pre-defined RealNVP is employed as a backbone model for manifold
learning and OOD detection of gray-scale images. The purpose of the designed
experiment is to evaluate the performance of various penalty functions in
cases where distinguishing between ID and OOD data becomes challenging due to
visual similarity. The model is trained on the MNIST dataset as ID data with
an embedded manifold in $\mathbb{R}^{3}$. Then, it is evaluated with $28\times
28$ gray-scale SVHN as OOD data. A hard OOD detection threshold is determined
as the maximum score obtained from the MNIST test data. Any sample with higher
score than this threshold is classified as OOD data. The reported results in
Fig. LABEL:fig:grayscale_OOD confirm the considerable superiority of using the
Huber function compared to MSE in successfully identifying OOD data with high
overlap.

Based on Fig. LABEL:fig:grayscale_OOD, the Huber function exhibits
significantly fewer mistakes on the SVHN test data compared to MSE.
Furthermore, the mistakes made by the Huber function are more meaningful,
particularly in situations where the data visually resembles the ID data
(MNIST). To provide further clarification, the Huber function nearly almost
tends to make errors on images resembling data within a distribution that
consists of a single digit on a dark background. It is important to highlight
that MNIST is structured in such a manner that each image contains a single
digit placed on a black background. Nevertheless, when the data deviates from
this pattern, as seen in cases like SVHN (including house license plate
numbers) where multiple digits may appear in a single image, the Huber
function may exhibit fewer errors. On the other hand, the MSE function
exhibits a higher number of errors and even mistakenly recognizes images
containing multiple digits as ID data.

The effect of manifold learning on color image OOD detection for NFs in
comparison with other methods [kingma2018glow, ren2019likelihood,
Serrà2020Input] is evaluated in terms of AUROC (Area Under the Receiver
Operating Characteristic Curve) criterion in Table 3 (trained on SVHN), Table
4 (trained on CIFAR10), and Table 5 (trained on CelebA). AUROC computes the
area under the ROC curve by using un-thresholded normalized predictions to the
range of $[0,1]$ on the data along with binary true labels. In our
configuration, we assign a true label of 0 for ID data and a true label of 1
for OOD data. A higher AUROC value indicates better performance.

The most important reason behind selecting the methods for comparison in Table
3, Table 4, and Table 5 is their consistency in the assumptions with our
proposed method, meaning that auxiliary OOD training data is not employed, and
the model structure is not tuned for this task, despite still achieving state-
of-the-art results. Each reported table consists of two separate parts (first
for the Huber and the next for MSE penalization). The first column includes
the name of the OOD data used during testing. The next columns are associated
with different manifold dimensions. In the case of dimension-preserving
($d=D$), the baseline method [kingma2018glow] (reported architecture in
Section 5.1.2), the likelihood ratio method [ren2019likelihood] (with the same
architecture as reported Glow network in Section 5.1.2), and the input
complexity method [Serrà2020Input] (with the same architecture as reported
Glow network in Section 5.1.2) are evaluated. The next columns for different
manifold dimensions have the same partitioning: the proposed method (named P),
the proposed method with considering likelihood ratio (named
$\text{P}+\text{ratio}$), and the proposed method with considering input
complexity (named $\text{P}+\text{IC}$). To summarize the approach of
combining the proposed method with the two mentioned methods (likelihood ratio
and input complexity), it can be stated that, based on the details provided in
Section 2, as both methods are test time scores, we incorporate only the
compressor loss from the input complexity method and the background loss from
the likelihood ratio method into our score.

Table 3 shows the superiority of the proposed method for different manifold
dimensions in complicated datasets. However, for simple-manifold family
datasets, the combination of the proposed method and IC is more effective.
Additionally, the performance of the proposed method appears to be nearly
indistinguishable when using both types of penalty functions. The Huber
function only demonstrates a slight advantage over MSE in very low dimensions.
In short, regardless of the specific penalty function used, the concept of
proposed method demonstrates some advantageous in OOD detection. The presented
results in Table 4 and Table 5 demonstrate similar performance to those in
Table 3. The key difference is that the proposed method alone or its
combination with the two mentioned methods (ratio, IC) outperform the
likelihood ratio and IC methods in terms of AUROC.

Table 3: AUROC score for Glow [kingma2018glow], likelihood ratio method (named
ratio) [ren2019likelihood], input complexity method (named IC)
[Serrà2020Input] with a lossless PNG compressor, and all cases of the proposed
method, namely the pure case (named P), including ratio (P+ratio) and
including IC (P+IC). Models are trained on SVHN. Table 4: AUROC score for Glow
[kingma2018glow], likelihood ratio method (named ratio) [ren2019likelihood],
input complexity method (named IC) [Serrà2020Input] with a lossless PNG
compressor, and all cases of the proposed method, namely the pure case (named
P), including ratio (P+ratio) and including IC (P+IC). Models are trained on
CIFAR10. Table 5: AUROC score for Glow [kingma2018glow], likelihood ratio
method (named ratio) [ren2019likelihood], input complexity method (named IC)
[Serrà2020Input] with a lossless PNG compressor, and all cases of the proposed
method, namely the pure case (named P), including ratio (P+ratio) and
including IC (P+IC). Models are trained on CelebA.

At end, one notable observation derived from analyzing the outcomes of the
trained model on the SVHN dataset, which can be generalized to other datasets,
is the effect of adjusting the parameters $d$ and $\delta$ on the AUROC value.
These results emphasize the significance of $\delta$ and $d$ in determining
the model’s performance, as depicted in Fig. LABEL:fig:search_parameter.

The left figure in Fig. LABEL:fig:search_parameter illustrates the changes in
the model’s performance as $\delta$ is increased. The mean and standard
deviation reported for each $\delta$ are based on models trained with that
specific threshold but with varying values of $d$. It is obvious that
increasing $\delta$ leads to higher AUROC values for datasets with simpler
data manifolds (like MNIST, FMNIST, and constant). However, the impact of
$\delta$ on the AUROC value appears to be unnoticeable for datasets with
complex manifolds (like CelebA, CIFAR10, and LSUN). It is important to note
that excessively high values of $\delta$ may deviate from the primary
objective of the problem. The mean and standard deviation reported for each
$\delta$ are based on models trained with that specific threshold but with
varying values of $d$. The right plot in Fig. LABEL:fig:search_parameter
demonstrates the effect of $d$ on the AUROC value. Each star point on the
diagram represents the mean and standard deviation of the output from trained
models with different $\delta$ for a fixed $d$. We observe that as the
dimensionality increases, the model’s ability for OOD detection decreases,
particularly in datasets with simpler manifolds. However, it is clear that
excessively reducing the dimensionality by choosing a value of $d$ smaller
than the actual data manifold is not suitable.

##  6 Conclusion

NFs cannot directly provide an estimation of the likelihood on the data
manifold. To address this, we adopted an algorithmic approach to jointly
estimate the likelihood and data manifold, and used it for OOD detection. Our
proposed method is based on the intuition that relying only on likelihood
value on the manifold or reconstruction loss is insufficient for OOD
detection. Accordingly, we found that usings these two indicators, namely
likelihood and reconstruction loss, jointly significantly improves OOD
detection. Through extensive experiments, we showed the strength of the
proposed method in finding data manifold and estimating its likelihood,
besides showing excellent performance for image OOD detection. The proposed
method can also generate samples similar to the latest state-of-the-art
manifold learning using normalizing flows models in terms of visual quality. A
per-pixel (or element-wise) Huber loss function was used in this paper for
manifold learning. The Huber function was employed here due to its transition
mode from linear to quadratic. It penalizes off-manifold parts linearly, and
the model focuses on learning the on-manifold sub-space well. Per-object (a
group of pixels) and per-sample (all pixels) loss can be an interesting work
to follow in a future work. It should be noted that a primary assumption of
the current study is that data lie on a manifold with one chart. In future
work, we aim to extend a multi-chart version of the current method with
possibly improved ability of OOD detection.


# Evaluations of certain Catalan-Hankel Pfaffians via classical skew
orthogonal polynomials

Bo-Jian Shen  School of Mathematical Sciences, Shanghai Jiaotong University,
People’s Republic of China<EMAIL_ADDRESS>,  Shi-Hao Li
Department of Mathematics, Sichuan University, Chengdu, 610064, China
<EMAIL_ADDRESS>and  Guo-Fu Yu  School of Mathematical Sciences,
Shanghai Jiaotong University, People’s Republic of China<EMAIL_ADDRESS>
###### Abstract.

This paper is to evaluate certain Catalan-Hankel Pfaffians by the theory of
skew orthogonal polynomials. Due to different kinds of hypergeometric
orthogonal polynomials underlying the Askey scheme, we explicitly construct
the classical skew orthogonal polynomials and then give different examples of
Catalan-Hankel Pfaffians with continuous and $q$-moment sequences.

###### Key words and phrases:

Catalan-Hankel Pfaffian; Classical orthogonal polynomials; Skew orthogonal
polynomials

###### 2010 Mathematics Subject Classification:

15A15, 33E20

##  1\. Introduction

Hankel determinants as a specific determinant has attracted much attention
from researchers in many different subjects (see, for example [5, 12, 29] and
references therein). It is well known that if $\\{\mu_{m}\\}_{m\geq 0}$ is the
moment sequence taking the form

$\displaystyle\mu_{n}=\int_{\mathbb{R}}x^{n}\omega(x)dx,$

then the Hankel determinant $\det(\mu_{i+j})_{i,j=0}^{n-1}$ has a nice
integral formula due to the Andreiéf formula [14]

$\displaystyle\det(\mu_{i+j})_{i,j=0}^{n-1}=\frac{1}{n!}\int_{\mathbb{R}^{n}}\prod_{1\leq
i<j\leq n}|x_{i}-x_{j}|^{2}\prod_{i=1}^{n}\omega(x_{i})dx_{i}.$

Moreover, if $\omega(x)$ is a classical weight, e.g. the Gaussian weight,
Laguerre weight, Jacobi weight and circular Jacobi weight, then Hankel
determinant with classical moments can be explicitly evaluated due to the
Selberg integral [13, Chapter 4]. Moreover, $q$-versions of Selberg integrals
which relate to the moment sequence with discrete measure play important roles
in modern mathematical physics as well [17].

In recent years, the Pfaffian version of Hankel determinant, called the
Catalan-Hankel Pfaffian, was proposed due to its potential applications in the
theory of combinatorics [19]. The Pfaffian version of Andreiéf formula, which
is called the de Bruijn formula, has been applied into the evaluation of
certain Catalan-Hankel Pfaffian with the help of the Selberg integral and its
$q$-version [19, 20]. Catalan-Hankel Pfaffian has the form

$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j-1}\right]_{i,j=0}^{2N-1}$ (1.1)

and its $q$-version

$\displaystyle\text{Pf}\left[([j]_{q}-[i]_{q})\mu_{i+j-1}\right]_{i,j=0}^{2N-1},\quad[j]_{q}=\frac{1-q^{j}}{1-q},$
(1.2)

where $\\{\mu_{n}\\}_{n\geq 0}$ is a moment sequence related to discrete
$q$-measure. The specific little $q$-Jacobi case were considered in [19] and
some recent developments like the Al-Salam $\&$ Carlitz I case was given in
[20].

Regarding evaluations of ($q$-)Catalan-Hankel Pfaffian, one application is to
give formula to the weighted enumeration of some specific plane partitions.
For example, in [19], the authors considered the little $q$-Jacobi case and
enumerated a special family of shifted reverse plane partitions with weights
that resemble the weight in the inner product of little $q$-Jacobi
polynomials. The importance of these evaluations lies in the random matrix
theory as well. In an earlier work of Mehta and Wang [25], they gave an
evaluation of the Hankel Pfaffian

$\displaystyle\text{Pf}[(j-i)\Gamma(i+j+b)]_{i,j=0}^{2N-1},$ (1.3)

which was related to the skew orthogonal Laguerre polynomials given in [1].
Moreover, the shifted Catalan-Hankel Pfaffian

$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j+r}\right]_{i,j=0}^{2N-1}\quad\text{or}\quad\text{Pf}\left[([j]_{q}-[i]_{q})\mu_{i+j+r}\right]_{i,j=0}^{2N-1}$

is closely related to the adjacent families of skew-orthogonal polynomials
considered in [23]. These evaluations may give more hints in exploring novel
examples in integrable systems and random matrices.

As mentioned, one way to evaluate these Pfaffians is based on the de Bruijn
formula and Selberg integral [19, 20]. In this paper, we will investigate
another way, namely by using the theory of classical ($q$-)skew orthogonal
polynomials, to make these evaluations and demonstrate the effectiveness—by
simply considering different classsical weights underlying Askey scheme, we
can get different examples of ($q$-)Catalan-Hankel Pfaffians. Since the Hankel
determinant is closely related to the theory of classical orthogonal
polynomials, then if we know normalisation constants of these classical
orthogonal polynomials, the evaluations of Hankel determinant can be made. It
is natural to ask whether we can apply the theory of classical ($q$-)skew
orthogonal polynomials into evaluations of the certain ($q$-)Catalan-Hankel
Pfaffians. The answer is affirmative. If we know the normalisation factors of
skew orthogonal polynomials, then these Pfaffians can be explicitly given by
those normalisation factors. Due to different expressions in discrete and
continuous cases (c.f. (1.1) and (1.2)), we discuss them separately.
Continuous cases including Hermite-type, Laguerre-type and Jacobi-type skew
orthogonal polynomials were firstly constructed by Adler et al [1] and we give
a brief review in section 2. Some evaluations of Catalan-Hankel Pfaffians
related to these weights including formula (1.3) will be demonstrated.
Moreover, we consider a Catalan-Hankel Pfaffian with Cauchy weight in that
part. Evaluations of $q$-Catalan-Hankel Pfaffians will be based on the theory
of classical $q$-skew orthogonal polynomials. With the help of discrete
Pearson relation, we construct different examples of classical ($q$-)skew
orthogonal polynomials, thus obtaining different kinds of skew orthogonal
polynomials and normalisation factors with respect to different classical
$q$-weights, which enlarge the results given by [16, 20].

##  2\. Continuous measure

Let’s consider the skew inner product related to (1.1)

$\displaystyle\langle\phi(x),\psi(x)\rangle_{4,\omega}=\frac{1}{2}\int_{\mathbb{R}}[\phi(x)\psi^{\prime}(x)-\phi^{\prime}(x)\psi(x)]\omega(x)dx,$
(2.1)

then the skew moments are given by

$\displaystyle m_{i,j}:=\langle
x^{i},x^{j}\rangle_{4,\omega}=\frac{1}{2}(j-i)\int_{\mathbb{R}}x^{i+j-1}\omega(x)dx.$

The skew LU decomposition (or so-called skew Borel decomposition [2]) of the
moment matrix $(m_{i,j})_{i,j\in\mathbb{N}}$ could give rise to the monic skew
orthogonal polynomials $\\{Q_{j}(x)\\}_{j\in\mathbb{N}}$ satisfying the skew
orthogonal relation

$\displaystyle\langle
Q_{2n}(x),Q_{2m+1}(x)\rangle_{4,\omega}=u_{n}\delta_{n,m},\quad\langle
Q_{2n}(x),Q_{2m}(x)\rangle_{4,\omega}=\langle
Q_{2n+1}(x),Q_{2m+1}(x)\rangle_{4,\omega}=0$ (2.2)

for certain $u_{n}>0$. Moreover, these skew orthogonal polynomials have the
following Pfaffian expressions [8]

$\displaystyle Q_{2n}(x)=\frac{1}{\tau_{2n}}\text{Pf}(0,\cdots,2n,x),\quad
Q_{2n+1}(x)=\frac{1}{\tau_{2n}}\text{Pf}(0,\cdots,2n-1,2n+1,x)$

with $\tau_{2n}=\text{Pf}(0,\cdots,2n-1)$ and the Pfaffian elements are given
by $\text{Pf}(i,j)=m_{i,j}$ and $\text{Pf}(i,x)=x^{i}$. By putting $m_{i,j}$
into the expression of $\tau_{2n}$, we have

$\displaystyle\tau_{2n}=\frac{1}{2^{n}}\text{Pf}\left[(j-i)\mu_{i+j-1}\right]_{i,j=0}^{2n-1},\quad\mu_{n}=\int_{\mathbb{R}}x^{n}\omega(x)dx.$
(2.3)

Interestingly, $u_{n}$ in the skew orthogonal condition (2.2) is the ratio
$\tau_{2n+2}/\tau_{2n}$, and therefore, if we can explicitly compute
$\\{u_{n}\\}_{n\in\mathbb{N}}$ in the skew orthogonal condition, then we
directly have $\tau_{2N}=\prod_{i=0}^{N-1}u_{i}.$ In fact, normalisation
factors $\\{u_{n}\\}_{n\in\mathbb{N}}$ can be explicitly computed when
$\\{Q_{n}(x)\\}_{n\in\mathbb{N}}$ are classical skew orthogonal polynomials.
For details, please refer to [1] and we give a brief review here.

Starting from the inner product

$\displaystyle\langle\phi(x),\psi(x)\rangle_{2,\rho}=\int_{\mathbb{R}}\phi(x)\psi(x)\rho(x)dx,$

one can construct a family of monic orthogonal polynomials
$\\{p_{j}(x)\\}_{j\in\mathbb{N}}$ satisfying the orthogonal relation

$\displaystyle\langle p_{j}(x),p_{k}(x)\rangle_{2,\rho}=h_{j}\delta_{j,k}.$

We call these orthogonal polynomials classical if the weight function
$\rho(x)$ satisfies the following Pearson equation

$\displaystyle\frac{\rho^{\prime}(x)}{\rho(x)}=-\frac{g(x)}{f(x)}\quad\text{with
deg $f(x)\leq 2$ and deg $g(x)\leq 1$}.$

From the relation, one can construct an operator

$\displaystyle\mathcal{A}=f\partial_{x}+\frac{f^{\prime}-g}{2},$ (2.4)

such that

$\displaystyle\langle\phi(x),\mathcal{A}\psi(x)\rangle_{2,\rho}=\langle\phi(x),\psi(x)\rangle_{4,\omega},\quad\omega(x)=\rho(x)f(x),$

where the skew inner product $\langle\cdot,\cdot\rangle_{4,\omega}$ was given
in (2.1). This is the key formula to establish the relation between classical
orthogonal polynomials and classical skew orthogonal polynomials. Moreover,
the normalisation factor $u_{n}$ could be computed via

$\displaystyle\mathcal{A}p_{k}(x)=-\frac{c_{k}}{h_{k+1}}p_{k+1}(x)+\frac{c_{k-1}}{h_{k-1}}p_{k-1}(x),\quad
u_{n}=c_{2n},$ (2.5)

with $h_{k}$ the normalisation factor of orthogonal polynomials. Therefore, to
obtain normalisation factors of skew orthogonal polynomials, our attention
should be paid to the computation of $c_{k}$ in the above equation. The
fastest method to compute $c_{k}$ is to compare the coefficients of $x^{k+1}$
on both sides and the following are the examples of the continuous classical
weights including Hermite, Laguerre, Jacobi and Cauchy weights.

From weight functions

$\displaystyle\rho(x)=\left\\{\begin{array}[]{ll}e^{-x^{2}},&\text{Hermite},\\\
x^{a}e^{-x},&\text{Laguerre},\\\ x^{a}(1-x)^{b},&\text{Jacobi},\\\
(1+x^{2})^{-a},&\text{Cauchy},\end{array}\right.$

one can obtain the Pearson pair

$\displaystyle(f,g)=\left\\{\begin{array}[]{ll}(1,2x),&\text{Hermite},\\\
(x,x-a),&\text{Laguerre},\\\ (x(1-x),(a+b)x-a),&\text{Jacobi},\\\
(1+x^{2},2ax),&\text{Cauchy}.\end{array}\right.$

This table was given by [13, eq. (5.58)] and it should be remarked that we use
the weight function of Jacobi as $x^{a}(1-x)^{b}$ supported in $[0,1]$ to make
the moments easily written down. There are two quantities obtained from this
table. One is the Hankel sequences given in (2.3). From
$\omega(x)=f(x)\rho(x)$, one knows that the weights of classical skew
orthogonal polynomials are

$\displaystyle\omega^{(H)}(x)=e^{-x^{2}},\quad\omega^{(L)}(x)=x^{a+1}e^{-x},\quad\omega^{(J)}(x)=x^{a+1}(1-x)^{b+1},\quad\omega^{(C)}(x)=(1+x^{2})^{-a+1}.$

Therefore, we have the following moments form

$\displaystyle\mu_{n}^{(H)}=\frac{1+(-1)^{n}}{2}\Gamma\left(\frac{n+1}{2}\right),\quad\mu_{n}^{(L)}=\Gamma(n+a+2),$
$\displaystyle\mu_{n}^{(J)}=\frac{\Gamma(n+a+2)\Gamma(b+2)}{\Gamma(n+a+b+4)},\quad\,\,\,\mu_{n}^{(C)}=\frac{1+(-1)^{n}}{2}\frac{\Gamma((n+1)/2)\Gamma(a-2-(n+1)/2)}{\Gamma(a-1)}.$

On the other hand, we can obtain

$\displaystyle
c_{k}=\left\\{\begin{array}[]{ll}h_{k+1}^{(H)},&\text{Hermite},\\\
h_{k+1}^{(L)}/2,&\text{Laguerre},\\\
(k+1+(a+b)/2)h_{k+1}^{(J)},&\text{Jacobi},\\\
(a-1-k)h_{k+1}^{(C)},&\text{Cauchy}\end{array}\right.$

where $\\{h_{k}\\}_{k\in\mathbb{N}}$ are normalisation factors with respect to
different weights [13, Chap. 5]

$\displaystyle h_{k}^{(H)}=\pi^{1/2}2^{-k}k!,\quad
h_{k}^{(L)}=\Gamma(k+1)\Gamma(a+k+1),$ $\displaystyle
h_{k}^{(J)}=2^{a+b+1+2k}\frac{\Gamma(k+1)\Gamma(a+b+1+k)\Gamma(a+1+k)\Gamma(b+1+k)}{\Gamma(a+b+2k+1)\Gamma(a+b+2k+2)},$
$\displaystyle h_{k}^{(C)}=\pi
2^{2k+2a+2}\frac{\Gamma(k+1)\Gamma(-2k-2a)\Gamma(-2a-2k-1)}{\Gamma(-2a-k)(\Gamma(-a-k))^{2}}.$

By using above results, we could state the following proposition.

######  Proposition 2.1.

We have the following evaluations of certain Catalan-Hankel Pfaffians

$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j-1}^{(H)}\right]_{i,j=0}^{2N-1}=2^{-N(N-1)}(\sqrt{\pi})^{N}\prod_{i=0}^{N-1}\Gamma(2i+2),$
$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j-1}^{(L)}\right]_{i,j=0}^{2N-1}=\prod_{i=0}^{N-1}\Gamma(2i+2)\Gamma(2i+a+2),$
$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j-1}^{(J)}\right]_{i,j=0}^{2N-1}=4^{N}\frac{\Gamma(N+(a+b+2)/4)}{\Gamma((a+b+2)/4)}\prod_{i=0}^{N-1}\frac{\Gamma(a+2i+2)\Gamma(b+2i+2)\Gamma(2i+3)}{\Gamma(a+b+4i+3)},$
$\displaystyle\text{Pf}\left[(j-i)\mu_{i+j-1}^{(C)}\right]_{i,j=0}^{2N-1}=(-1)^{N}2^{2N^{2}-N(a-3)}\frac{\Gamma(N+(1-a)/2)}{\Gamma((1-a)/2)}\prod_{i=0}^{N-1}\frac{\Gamma(2a-4i-2)\Gamma(2i+3)}{(\Gamma(a-2i-1))^{2}}.$

######  Remark 2.2.

The second formula is exactly equation (1.3) with the shift $a\to a-1$.

##  3\. Discrete measure: $q$-case

This part is devoted to the evaluations of $q$-Catalan Hankel Pfaffians given
by formula (1.2). The $q$-case corresponds to a special discrete measure
distributed on exponential lattices $x(i)=q^{i}$ with $0<q<1$ and
$i\in\mathbb{Z}$. By using the definition of Jackson’s $q$-integral 111The
definition of $q$-integral on the interval $[0,\infty)$ is different from the
one defined on $[0,1]$, see [21]. However, these two cases can be treated
similarly and we just consider the former one here.

$\displaystyle\int_{0}^{\infty}f(x)d_{q}x=(1-q)\sum_{s=-\infty}^{\infty}f(q^{s})q^{s},$

one can define the following inner product

$\displaystyle\langle\phi(x),\psi(x)\rangle_{2,\rho}$
$\displaystyle:=\int_{0}^{\infty}\phi(x)\psi(x)\rho(x)d_{q}x=(1-q)\sum_{s\in\mathbb{Z}}\phi(q^{s})\psi(q^{s})\rho(q^{s})q^{s}.$
(3.1)

With this inner product, $q$-orthogonal polynomials
$\\{p_{n}(x;q)\\}_{n\in\mathbb{N}}$ are defined by the orthogonal relation

$\displaystyle\langle
p_{i}(x;q),p_{j}(x;q)\rangle_{2,\rho}=h_{n}(q)\delta_{n,m}$ (3.2)

with respect to the weight $\rho(x;q)$. Moreover, if $\rho(x;q)$ is a
classical weight, we call the corresponding orthogonal polynomials the
classical $q$-orthogonal polynomials. Classical $q$-orthogonal polynomials
include many interesting examples underlying the $q$-Askey scheme. For
details, please refer to [3, 21, 22]. One important property of classical
$q$-orthogonal polynomials is that the weight function satisfies an analogy of
Pearson relation given by Nikiforov and Suslov [26]

$\displaystyle\frac{\rho(qx)}{\rho(x)}=\frac{f(x)-q^{-\frac{1}{2}}(1-q)xg(x)}{f(qx)}\quad\text{with
deg $f(x)\leq 2$ and deg $g(x)\leq 1$}.$ (3.3)

In the following, we will demonstrate how to connect $q$-inner product with
$q$-skew inner product.

Let’s define a $q$-analogy of the skew inner product (2.1)

$\displaystyle\langle\phi(x),\psi(x)\rangle_{4,\omega}=\int_{0}^{\infty}[\phi(x)D_{q}\psi(x)-D_{q}\phi(x)\psi(x)]\omega(x)d_{q}x$
(3.4)

with $q$-difference operator

$D_{q}f(x)=\frac{f(x)-f(qx)}{(1-q)x}.$

Then by defining an operator $\mathcal{A}_{q}$ [16]

$\displaystyle\mathcal{A}_{q}=q^{-\frac{1}{2}}g(x)T_{q}+q^{-1}f(x)D_{q^{-1}}+f(x)D_{q},\quad
T_{q}f(x)=f(qx),$ (3.5)

one can find such a connection formula

$\displaystyle\langle\phi(x;q),\mathcal{A}_{q}\psi(x;q)\rangle_{2,\rho}=\langle\phi(x;q),\psi(x;q)\rangle_{4,\omega},\quad\omega(x)=\rho(qx)f(qx).$
(3.6)

On the other hand, from the skew inner product (3.1), one has the following
($q$-)skew moments

$\displaystyle
m_{i,j}:=([j]_{q}-[i]_{q})\int_{0}^{\infty}x^{i+j-1}\omega(x)d_{q}x,\quad[j]_{q}=\frac{1-q^{j}}{1-q}.$
(3.7)

Similarly, it is known that a family of monic skew orthogonal polynomials
$\\{Q_{i}(x;q)\\}_{i\in\mathbb{N}}$ could be constructed from those moments if
even-order moment matrices $\text{Pf}(m_{i,j})_{i,j=0}^{2n-1}\neq 0$ for all
$n\in\mathbb{N}_{+}$ [16]. Furthermore, polynomials
$\\{Q_{i}(x;q)\\}_{i\in\mathbb{N}}$ admit the following Pfaffian expressions

$\displaystyle
Q_{2n}(x;q)=\frac{1}{\tau_{2n}(q)}\text{Pf}(0,\cdots,2n,x),\quad
Q_{2n+1}(x;q)=\frac{1}{\tau_{2n}(q)}\text{Pf}(0,\cdots,2n-1,2n+1,x)$
$\displaystyle\tau_{2n}(q)=\text{Pf}(0,\cdots,2n-1),\quad\text{Pf}(i,j)=m_{i,j},\quad\text{Pf}(i,x)=x^{i},$

and they satisfy the following skew orthogonal relations

$\displaystyle\langle
Q_{2n}(x;q),Q_{2m+1}(x;q)\rangle_{4,\omega}=\frac{\tau_{2n+2}(q)}{\tau_{2n}(q)}:=u_{n}(q)\delta_{n,m},$
$\displaystyle\langle Q_{2n}(x;q),Q_{2m}(x;q)\rangle_{4,\omega}=\langle
Q_{2n+1}(x;q),Q_{2m+1}(x;q)\rangle_{4,\omega}=0.$

Therefore, one knows that

$\displaystyle\text{Pf}\left[([j]_{q}-[i]_{q})\int_{0}^{\infty}x^{i+j-1}\omega(x)d_{q}x\right]_{i,j=0}^{2N-1}=\text{Pf}(0,\dots,2N-1)=\prod_{i=0}^{N-1}u_{i}(q).$

Interestingly, there is a method to evaluate the value of $u_{n}(q)$ quickly
by taking advantage of $q$-Pearson relation.

As an analogy of equation (2.5), there holds the formula

$\displaystyle\mathcal{A}_{q}p_{k}(x;q)=-\frac{c_{k}(q)}{h_{k+1}(q)}p_{k+1}(x;q)+\frac{c_{k-1}(q)}{h_{k-1}(q)}p_{k-1}(x;q)$
(3.8)

in the $q$-case due to the property of $\mathcal{A}_{q}$ [16, eq. (4.28)].
Moreover, the quantity $h_{k}(q)$ is the normalisation factor of the
orthogonal relation (3.2) and $c_{k}(q)$ is closely related to $u_{k}(q)$ via
the relation $u_{k}(q)=c_{2k}(q)$. Therefore, it is the point to compute
$c_{k}(q)$ from above equation and then the exact value of $q$-Catalan Hankel
Pfaffian could be obtained. In the following part, we discuss several
different cases including the examples in [19, 20].

###  3.1. Al-Salam $\&$ Carlitz I case

The first example considered here is the Al-Salam $\&$ Carlitz I case with the
weight function

$\displaystyle\rho(x;q)=(qx,a^{-1}qx;q)_{\infty},\quad a<0.$

It is well known that the Al-Salam $\&$ Carlitz polynomials
$\\{U_{n}^{(a)}(x;q)\\}_{n\geq 0}$ have the orthogonality

$\displaystyle\int_{a}^{1}U_{m}^{(a)}(x;q)U_{n}^{(a)}(x;q)\rho(x;q)d_{q}x=(-a)^{n}(1-q)(q;q)_{n}(q,a,a^{-1}q;q)_{\infty}q^{n\choose
2}\delta_{n,m}:={h}_{n}\delta_{n,m},$

and canonical moments were given by [27]

$\displaystyle\mu_{n}=\int_{a}^{1}x^{n}\rho(x;q)d_{q}x=(1-q)\left(q,a,a^{-1}q;q\right)_{\infty}\sum_{i=0}^{n}\left[\begin{array}[]{c}n\\\
i\end{array}\right]_{q}a^{i}$

Specifically, Al-Salam $\&$ Carlitz I polynomials have $q$-hypergeometric
function expressions

$\displaystyle U_{n}^{(a)}(x;q)=(-a)^{n}q^{n\choose
2}{{}_{2}}\phi_{1}\left(\left.\begin{array}[]{c}q^{-n},x^{-1}\\\
0\end{array}\right|q;\frac{qx}{a}\right).$

As was shown in [16, Sec. 4.4.1(1)], the Pearson pair in this case is

$\displaystyle(f,g)=\left(x^{2}-(1+a)x+a,\,\frac{q^{1/2}}{1-q}(x-(1+a))\right),$

and the coefficient $c_{n}$ in (3.8) can be explicitly computed as
$c_{n}=-q^{-n}{h}_{n+1}/(1-q)$. Thus the moment Pfaffian can be written as

$\displaystyle\text{Pf}(m_{i,j})_{i,j=0}^{2N-1}$
$\displaystyle=\prod_{k=0}^{N-1}a^{2k+1}(q;q)_{2k+1}q^{\binom{2k}{2}}\left(q,a,a^{-1}q;q\right)_{\infty}$
$\displaystyle=a^{N^{2}}q^{\frac{1}{6}N(N-1)(4N-5)}\left(q,a,a^{-1}q;q\right)_{\infty}^{N}\prod_{k=0}^{N-1}(q;q)_{2k+1},$

where $m_{i,j}=([j]_{q}-[i]_{q})\int_{a}^{1}x^{i+j-1}\omega(x;q)d_{q}x$. Since
$\omega(x;q)=f(qx;q)\rho(qx;q)=a\rho(x;q)$, we have

$\displaystyle
m_{i,j}=a([j]_{q}-[i]_{q})\mu_{i+j-1}=(q^{i}-q^{j})\left(q,a,a^{-1}q;q\right)_{\infty}\sum_{k=0}^{i+j-1}\left[\begin{array}[]{c}i+j-1\\\
k\end{array}\right]_{q}a^{k+1}.$

Dividing both sides by $\left[a(q,a,a^{-1}q;q)_{\infty}\right]^{N}$ leads to

$\displaystyle\text{Pf}\left((q^{i}-q^{j})\sum_{k=0}^{i+j-1}\left[\begin{array}[]{c}i+j-1\\\
k\end{array}\right]_{q}a^{k}\right)_{i,j=0}^{2N-1}=a^{N(N-1)}q^{\frac{1}{6}N(N-1)(4N-5)}\prod_{k=0}^{N-1}(q;q)_{2k+1}.$

######  Remark 3.1.

With $a=-1$, Al-Salam $\&$ Carlitz I polynomials reduce to the $q$-Hermite I
polynomials [22, Sec. 14.28], therefore the above mentioned method can be
applied to the $q$-Hermite I case as well. The moments of q-Hermite I
polynomials take the form [27]

$\displaystyle\mu_{n}=(1-q)(q,-1,-q;q)_{\infty}\frac{1+(-1)^{n}}{2}\left(q;q^{2}\right)_{n/2}.$

Thus we have the following evaluation of $q$-Catalan-Hankel Pfaffian

$\displaystyle\text{Pf}\left((q^{i}-q^{j})\frac{1+(-1)^{i+j-1}}{2}(q;q^{2})_{\frac{i+j-1}{2}}\right)_{i,j=0}^{2N-1}=(-1)^{N}q^{\frac{1}{6}N(N-1)(4N-5)}\prod_{k=0}^{N-1}(q;q)_{2k+1}.$

###  3.2. Stieltjes-Wigert case

The Stieltjes-Wigert polynomials are well studied in random matrix theory of
the so-called Stieltjes-Wigert ensemble, which firstly appeared in the study
of non-intersecting Brownian walkers, and subsequently in quantum many body
systems etc.. For a detailed review, please refer to [15] and references
therein. The weight of Stieltjes-Wigert polynomials was first given by
Stieltjes as an example of indeterminate moment problems [28] and further
studied by Wigert [30]. Usually there are several different expressions for
the Stieltjes-Wigert’s weight function [11]. The original one corresponding to
a continuous measure is

$\displaystyle w(x)=\frac{1}{\sqrt{\pi}}kx^{-k^{2}\log x},\quad x>0$

with the moment $\mu_{n}=\int_{0}^{\infty}x^{n}w(x)dx=e^{(n+1)^{2}/4k^{2}}$.
If we set $q=e^{-1/2k^{2}}$, then $\mu_{n}$ can be written as
$q^{-(n+1)^{2}/2}$. With this measure, Wigert found the following expression
for Stieltjes-Wigert polynomials [30]

$\displaystyle
P_{n}(x)=(-1)^{n}\frac{q^{n/2+1/4}}{\sqrt{(q;q)_{n}}}\sum_{k=0}^{n}\left[\begin{array}[]{l}n\\\
k\end{array}\right]_{q}(-1)^{k}q^{k^{2}+k/2}x^{k},$

with orthogonality

$\displaystyle\int_{0}^{\infty}P_{n}(x)P_{m}(x)w(x)dx=\delta_{mn}.$

In [9], Chihara proposed a discrete weight on the exponential lattices
admitting the form

$\displaystyle\xi(x)=\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{q}M}q^{n+n^{2}/2}&,x=q^{n}\\\
0&,x\neq q^{n}\end{array}\right.,\quad
M=(-q\sqrt{q},-q^{-1/2},q;q)_{\infty},\quad n\in\mathbb{Z}.$

Then the corresponding inner product becomes

$\displaystyle\langle
p_{n}(x),p_{m}(x)\rangle_{2,\xi}=\frac{1}{\sqrt{q}M}\sum_{k=-\infty}^{\infty}p_{n}(q^{k})p_{m}(q^{k})q^{k+k^{2}/2}$

We can prove that this discrete weight is equivalent to the continuous one by
showing that they have same moments. For completeness, we give a short proof
to this fact.

One can easily check that

$\displaystyle\xi(qx)=q^{\frac{3}{2}}x\xi(x),\quad x\in\mathbb{R}$

and therefore obtain a recurrence relation for the moments
$\\{\mu_{n}\\}_{n\in\mathbb{N}}$

$\displaystyle\mu_{n}$ $\displaystyle:=\langle
x^{n},1\rangle_{2,\xi}=\sum_{k=-\infty}^{\infty}q^{nk}\xi(q^{k})=\sum_{k=-\infty}^{\infty}q^{n(k+1)}\xi(q^{k+1})$
$\displaystyle\quad=q^{n+\frac{3}{2}}\sum_{k=-\infty}^{\infty}q^{(n+1)k}\xi(q^{k})=q^{n+\frac{3}{2}}\mu_{n+1}.$

As a result, we have $\mu_{n}=q^{-n^{2}/2-n}\mu_{0}$. By making use of the
Jacobi triple product identity [18, Thm 352]

$\displaystyle\sum^{\infty}_{n=-\infty}(-1)^{n}q^{\binom{n}{2}}x^{n}=(x,q/x,q;q)_{\infty},$

$\mu_{0}$ can be directly computed as

$\displaystyle\mu_{0}=\frac{1}{\sqrt{q}M}\sum^{\infty}_{k=-\infty}q^{k+k^{2}/2}=\frac{1}{\sqrt{q}}.$

Thus $\mu_{n}=q^{-\frac{(n+1)^{2}}{2}}$, which is the same as that of the
continuous measure.

For the present need, we would consider the discrete measure. Denote

$\displaystyle\rho(x;q)=\frac{1}{\sqrt{q}M}x^{\frac{\ln x}{2\ln q}}$

to be the weight function and define

$\displaystyle
p_{n}(x)=\frac{\sqrt{(q;q)_{n}}}{q^{n^{2}+n+\frac{1}{4}}}P_{n}(x)$

to be the monic Stieltjes-Wigert polynomials. We can rewrite the orthogonality
in terms of $\rho(x)$ and $\\{p_{n}(x)\\}_{n\geq 0}$ as

$\displaystyle\int_{0}^{\infty}p_{n}(x)p_{m}(x)\rho(x)d_{q}x=(1-q)\frac{(q;q)_{n}}{q^{2n^{2}+2n+\frac{1}{2}}}\delta_{mn}:=h_{n}\delta_{n,m}.$

By assuming $f(x)=x$ and solving the Pearson equation (3.3), we get
$g(x)={(q^{2}x-q^{1/2})}/{(q-1)}$ which leads to

$\displaystyle c_{n}(q)=\frac{q^{n+\frac{3}{2}}}{1-q}h_{n+1}.$

As a consequence, we have the following evaluation of the moment Pfaffian

$\displaystyle\text{Pf}(m_{i,j})_{i,j=0}^{2N-1}=q^{-\frac{1}{6}N(N+1)(8N-5)}\prod_{k=0}^{N-1}(q;q)_{2k+1}.$

On the other hand, we can compute the skew moments (3.7) as a scalar product
of canonical moments $\mu_{n}$

$\displaystyle m_{i,j}=$
$\displaystyle([j]_{q}-[i]_{q})\int_{0}^{\infty}x^{i+j-1}\omega(x;q)d_{q}x=([j]_{q}-[i]_{q})\int_{0}^{\infty}x^{i+j+1}q^{\frac{3}{2}}\rho(x;q)d_{q}x$
$\displaystyle=$
$\displaystyle([j]_{q}-[i]_{q})q^{{3}/{2}}\mu_{i+j+1}=([j]_{q}-[i]_{q})q^{-[(i+j+2)^{2}-3]/2}.$

By combining above results, we have the evaluation

$\displaystyle\text{Pf}(([j]_{q}-[i]_{q})q^{-(i+j+2)^{2}/2})_{i,j=0}^{2N-1}=q^{-\frac{1}{6}N(2N+1)(8N-1)}\prod_{k=0}^{N-1}(q;q)_{2k+1}.$

###  3.3. Little $q$-Jacobi case

Little $q$-Jacobi polynomials are important in many mathematical fields such
as polynomials theory [22] and quantum group [24]. These polynomials have the
following series form (c.f. [24, eq. (2.21)])

$\displaystyle p_{n}^{(\alpha,\beta)}(z;q)=\sum_{r\geq
0}\frac{(q^{-n};q)_{r}(q^{\alpha+\beta+n+1};q)_{r}}{(q;q)_{r}(q^{\alpha+1};q)_{r}}(qz)^{r}$

and they obey the following orthogonal relation [24, Prop. 3.9]

$\displaystyle\int_{0}^{1}p_{m}^{(\alpha,\beta)}(z;q)p_{n}^{(\alpha,\beta)}(z;q)z^{\alpha}(qz;q)_{\beta}d_{q}z=\delta_{n,m}q^{(\alpha+1)n}\frac{(1-q)(q;q)_{\alpha}^{2}(q;q)_{\beta+n}(q;q)_{n}}{(1-q^{\alpha+{\beta}+2n+1})(q;q)_{\alpha+n}(q;q)_{\alpha+\beta+n}}.$

In this case, the canonical moments are given by [27]

$\displaystyle\mu_{n}^{(\alpha,\beta)}=\int_{0}^{1}z^{{\alpha}+n}(qz;q)_{\beta}d_{q}z=\frac{(1-q)(q^{\alpha+\beta+1};q)_{\infty}(q;q)_{\infty}(q^{\alpha+1};q)_{n}}{(1-q^{\alpha+\beta+1})(q^{\alpha+1};q)_{\infty}(q^{\beta+1};q)_{\infty}(q^{\alpha+\beta+2};q)_{n}}.$

If we define monic little $q$-Jacobi polynomials
$\\{\tilde{p}_{n}^{({\alpha},{\beta})}(z;q)\\}_{n\geq 0}$

$\displaystyle\tilde{p}_{n}^{({\alpha},{\beta})}(z;q)=(-1)^{n}q^{\frac{n(n-1)}{2}}\frac{(q^{{\alpha}+1};q)_{n}}{(q^{{\alpha}+{\beta}+n+1};q)_{n}}p_{n}^{({\alpha},{\beta})}(z;q),$

then they satisfy the orthogonal relation

$\displaystyle\begin{aligned}
\int_{0}^{1}&\tilde{p}_{n}^{({\alpha},{\beta})}(z;q)\tilde{p}_{m}^{({\alpha},{\beta})}(z;q)z^{\alpha}(qz;q)_{\beta}d_{q}z\\\
&=q^{n(n+{\alpha})}[{\alpha}+{\beta}+2n+1]_{q}^{-1}\frac{(q;q)_{n}(q;q)_{n+{\alpha}}(q;q)_{n+{\beta}}(q;q)_{n+{\alpha}+{\beta}}}{(q;q)_{2n+{\alpha}+{\beta}}^{2}}\delta_{n,m}:=h_{n}\delta_{n,m}.\end{aligned}$
(3.9)

Regarding the weight of little $q$-Jacobi polynomials

$\displaystyle\rho(z;q)=z^{\alpha}(qz;q)_{\beta},$

by solving (3.3), it admits the Pearson pair

$\displaystyle(f,g)=\left(-x^{2}+x,-q^{\frac{1}{2}}\left([\alpha+\beta+2]_{q}x-[\alpha+1]_{q}\right)\right).$

Therefore, the coefficient $c_{n}$ has the expression

$\displaystyle c_{n}=q^{-n}[2n+2+\alpha+\beta]_{q}h_{n+1},$

where $h_{n}$ is the nomalization constant in the orthogonal relation of monic
polynomials (3.9). Then we have the following expression for the moment
Pfaffian

$\displaystyle\text{Pf}(m_{ij}^{(\alpha,\beta)})_{i,j=0}^{2N-1}$
$\displaystyle=q^{\frac{1}{3}N(4N^{2}-3N+2)+\alpha N^{2}}$
$\displaystyle\times\prod_{k=0}^{N-1}\frac{(1-q^{\alpha+\beta+4k+2})(q,q^{\alpha+1},q^{\beta+1},q^{\alpha+\beta+1};q)_{2k+1}(q;q)_{\infty}(q^{\alpha+\beta+1};q)_{\infty}}{(q^{\alpha+1},q^{\beta+1};q)_{\infty}(q^{\alpha+\beta+1};q)_{4k+2}(q^{\alpha+\beta+1};q)_{4k+3}}.$

Since the weight of the corresponding skew orthogonal little $q$-Jacobi
polynomials is

$\displaystyle\omega(x;q)=f(qx;q)\rho^{({\alpha},{\beta})}(qx;q)=q^{\alpha+1}\rho^{(\alpha+1,\beta+1)}(x;q)=q^{\alpha+1}x^{\alpha+1}(qx;q)_{\beta+1},$

the skew moments $\\{m_{i,j}^{(\alpha,\beta)}\\}_{i,j\in\mathbb{N}}$ are
related to the canonical moments by

$\displaystyle m_{i,j}^{(\alpha,\beta)}=$
$\displaystyle([j]_{q}-[i]_{q})\int_{0}^{1}x^{i+j-1}\omega(x;q)d_{q}x=([j]_{q}-[i]_{q})q^{\alpha+1}\mu_{i+j-1}^{(\alpha+1,\beta+1)}$
$\displaystyle=$
$\displaystyle([j]_{q}-[i]_{q})q^{\alpha+1}\frac{(1-q)(q^{\alpha+\beta+1};q)_{\infty}(q;q)_{\infty}(q^{\alpha+1};q)_{n}}{(1-q^{\alpha+\beta+1})(q^{\alpha+1};q)_{\infty}(q^{\beta+1};q)_{\infty}(q^{\alpha+\beta+2};q)_{n}}.$

After eliminating the constant
$\left[q^{\alpha+1}{(q^{\alpha+\beta+1};q)_{\infty}(q;q)_{\infty}}/{\left((1-q^{\alpha+\beta+1})(q^{\alpha+1},q^{\beta+1};q)_{\infty}\right)}\right]^{N}$,
we have

$\displaystyle\text{Pf}\left(([j]_{q}-[i]_{q})\frac{(1-q)(q^{{\alpha}+2};q)_{i+j-1}}{(q^{{\alpha}+{\beta}+4};q)_{i+j-1}}\right)_{i,j=0}^{2N-1}$
(3.10)
$\displaystyle\qquad=q^{\frac{1}{3}N(N-1)(4N+1)+{\alpha}N(N-1)}\prod_{k=0}^{N-1}\frac{(q;q)_{2k+1}(q^{{\alpha}+2},q^{{\beta}+2};q)_{2k}}{(q^{{\alpha}+{\beta}+4};q)_{2k-2}(q^{{\alpha}+{\beta}+2k+2};q)_{2k}(q^{{\alpha}+{\beta}+2k+2};q)_{2k+2}},$

which coincides with the result in [20, Thm. 5.1].

###  3.4. Big $q$-Jacobi case

Big $q$-Jacobi polynomials were introduced by Andrews and Askey as an
infinite-dimentional version of the $q$-Hahn polynomials [4]. In addition, big
$q$-Jacobi polynomials were also contained in the Bannai-Ito scheme of dual
systems of orthogonal polynomials as an infinite dimension analogue of the
q-Racah polynomials [6]. These polynomials take the hypergeometric function
form

$\displaystyle
P_{n}(x;a,b,c;q)={}_{3}\phi_{2}\left(\left.\begin{array}[]{c}q^{-n},abq^{n+1},x\\\
aq,cq\end{array}\right|q;q\right),$

orthogonal with respect to the weight

$\displaystyle\rho^{(a,b,c)}(x;q)=\frac{(a^{-1}x,c^{-1}x;q)_{\infty}}{(x,bc^{-1}x;q)_{\infty}}$

and have the orthogonality

$\displaystyle\int_{cq}^{aq}P_{m}(x;a,b,c;q)P_{n}(x;a,b,c;q)\rho^{(a,b,c)}(x;q)d_{q}x$
$\displaystyle\quad=aq(1-q)\frac{\left(q,a^{-1}c,ac^{-1}q,abq^{2};q\right)_{\infty}}{\left(aq,bq,cq,abc^{-1}q;q\right)_{\infty}}\frac{(1-abq)}{\left(1-abq^{2n+1}\right)}\frac{\left(q,bq,abc^{-1}q;q\right)_{n}}{(abq,aq,cq;q)_{n}}\left(-acq^{2}\right)^{n}q^{\binom{n}{2}}\delta_{mn}.$

One can easily check that the normalisation constant for the monic polynomials
is

$\displaystyle
h_{n}=\frac{aq(1-q)(-acq^{2})^{n}q^{\binom{n}{2}}(q,a^{-1}c,ac^{-1}q,abq;q)_{\infty}(q,aq,bq,cq,abc^{-1}q;q)_{n}}{(1-abq^{2n+1})(aq,bq,cq,abc^{-1}q;q)_{\infty}(abq;q)_{n}(abq^{n+1};q)_{n}^{2}}$

and the canonical moments are given by [27]

$\displaystyle\mu_{n}^{(a,b,c)}$
$\displaystyle=\int_{cq}^{aq}x^{n}\rho^{(a,b,c)}(x;q)d_{q}x$
$\displaystyle=aq\frac{\left(abq^{2},a^{-1}c,ac^{-1}q;q\right)_{\infty}}{\left(aq,bq,cq,abc^{-1}q;q\right)_{\infty}}\sum_{m=0}^{n}(-1)^{m}\left[\begin{array}[]{c}n\\\
m\end{array}\right]_{q}q^{-nm+\binom{m+1}{2}}\frac{(aq,cq;q)_{m}}{\left(abq^{2};q\right)_{m}}.$

Solving equation (3.3) gives the Pearson pair

$\displaystyle(f,g)=\left((1-\frac{x}{aq})(1-\frac{x}{cq}),\frac{q^{\frac{1}{2}}}{1-q}\left((\frac{1}{acq^{2}}-\frac{b}{c})x+\frac{b}{c}+1-\frac{1}{aq}-\frac{1}{cq}\right)\right).$

Then, by comparing the leading coefficients on both sides of (3.8) we have

$\displaystyle c_{n}=\frac{abq^{2n+2}-1}{ac(1-q)q^{n+2}}h_{n+1},$

which leads to the following evaluation of the moment Pfaffian

$\displaystyle\text{Pf}(m_{ij}^{(a,b,c)})_{i,j=0}^{2N-1}=$
$\displaystyle(-1)^{N}c^{N(N-1)}a^{N^{2}}q^{\frac{1}{6}N(N+1)(4N-1)}$
$\displaystyle\times$
$\displaystyle\prod_{k=0}^{N-1}\frac{(q,a^{-1}c,ac^{-1}q,abq;q)_{\infty}(q,aq,bq,cq,abc^{-1}q;q)_{2k+1}}{(aq,bq,cq,abc^{-1}q;q)_{\infty}(abq;q)_{2(k+N)+1}}.$

where $m_{i,j}^{(a,b,c)}$ is given by

$\displaystyle
m_{i,j}^{(a,b,c)}=([j]_{q}-[i]_{q})\int_{cq}^{aq}x^{i+j-1}\omega(x;q)^{(a,b,c)}d_{q}x,$

and $\omega(x;q)$ is expressed by

$\displaystyle\omega(x;q)^{(a,b,c)}=f(qx;q)\rho^{(a,b,c)}(qx;q)=(1-x)(1-bc^{-1}x)\rho(x;q)^{(a,b,c)}=\rho^{(qa,qb,qc)}(qx;q).$

Note that

$\displaystyle\mu_{n}^{(a,b,c)}=\int_{cq}^{aq}x^{n}\rho^{(a,b,c)}(x;q)d_{q}x=q^{n+1}\int_{c}^{a}x^{n}\rho^{(a,b,c)}(qx;q)d_{q}x,$

we have the expression

$\displaystyle m^{(a,b,c)}_{i,j}=$
$\displaystyle([j]_{q}-[i]_{q})q^{-(i+j)}\mu_{i+j-1}^{(qa,qb,qc)}$
$\displaystyle=$
$\displaystyle([j]_{q}-[i]_{q})aq^{-(i+j-1)}\frac{\left(abq^{2},a^{-1}c,ac^{-1}q;q\right)_{\infty}}{\left(aq,bq,cq,abc^{-1}q;q\right)_{\infty}}\left(\sum_{m=0}^{n}(-1)^{m}\left[\begin{array}[]{c}n\\\
m\end{array}\right]_{q}q^{-nm+\binom{m+1}{2}}\frac{(aq,cq;q)_{m}}{\left(abq^{2};q\right)_{m}}\right).$

Dividing both sides by
$\left[aq{\left(abq^{2},a^{-1}c,ac^{-1}q;q\right)_{\infty}}/{\left(aq,bq,cq,abc^{-1}q;q\right)_{\infty}}\right]^{N}$
and changing variables $(a,b,c)\rightarrow(q^{-1}a,q^{-1}b,q^{-1}c)$, we get

Pf
$\displaystyle\left(([j]_{q}-[i]_{q})\sum_{m=0}^{i+j-1}(-1)^{m}\left[\begin{array}[]{c}i+j-1\\\
m\end{array}\right]_{q}q^{-(i+j-1)(m+1)+\binom{m+1}{2}}\frac{(aq,cq;q)_{m}}{\left(abq^{2};q\right)_{m}}\right)_{i,j=0}^{2N-1}$
$\displaystyle=(-1)^{N}(ac)^{N(N-1)}q^{\frac{1}{6}N(4N^{2}-9N+11)}\prod^{N-1}_{k=0}\frac{(q;q)_{\infty}(q;q)_{2k+1}(a,b,c,abc^{-1}q;q)_{2k}}{(abq^{2};q)_{2(k+N)-2}}.$

######  Remark 3.2.

With $a=b=1$, big q-Jacobi polynomials reduce to the big q-Legendre
polynomials [23]. In particular, we have the following evaluation of
q-Catalan-Hankel Pfaffian

Pf
$\displaystyle\left(([j]_{q}-[i]_{q})\sum_{m=0}^{i+j-1}(-1)^{m}\left[\begin{array}[]{c}i+j-1\\\
m\end{array}\right]_{q}q^{-(i+j-1)(m+1)+\binom{m+1}{2}}\frac{(q,cq;q)_{m}}{\left(q^{2};q\right)_{m}}\right)_{i,j=0}^{2N-1}$
$\displaystyle=(-1)^{N}(c)^{N(N-1)}q^{\frac{1}{6}N(4N^{2}-9N+11)}\prod^{N-1}_{k=0}\frac{(q;q)_{\infty}(q;q)_{2k+1}(q,q,cq,c^{-1}q;q)_{2k}}{(q^{2};q)_{2(k+N)-2}}.$

##  4\. Further remarks

In this paper, we have developed a method based on the relation between
classical ($q$-)orthogonal and ($q$-)skew orthogonal polynomials to evaluate
certain $q$-Catalan-Hankel Pfaffians whose entries are composed of the moments
of classical orthogonal polynomials. Some examples are given to illustrate the
approach including the continuous ones (e.g. Hermite, Laguerre, Jacobi and
Cauchy) and discrete $q$-cases (e.g. Al-Salam $\&$ Carlitz I, little
$q$-Jacobi, Stieltjes-Wigert and big $q$-Jacobi polynomials). Amoung those
examples, the Al-Salam $\&$ Carlitz I and Little $q$-Jacobi case are compared
with the results obtained by M. Ishikawa and J. Zeng in [20] and we present
alternative proofs of [20, eq. (6.7) & Thm 5.2]. Besides, the examples in [20,
Conjecture 7.1] seem to be related to some discrete measure on the linear
lattice. However, as mentioned in [16], the skew moments defined in linear
lattices are of the form $\mu_{i,j}=(j-i)\mu_{i+j-1}+({j\choose 2}-{i\choose
2})\mu_{i+j-2}+\cdots$, it is unclear to us whether it could be written in the
Catalan-Hankel Pfaffian form.

As mentioned before, $q$-analogues of Selberg integral can be used to evaluate
$q$-Catalan-Hankel Pfaffians. For example, the formula (3.10) can be evaluated
by the following Askey-Habsieger-Kadell formula [20]

$\displaystyle\int_{[0,1]^{n}}\prod_{i<j}t_{i}^{2k}\left(q^{1-k}t_{j}/t_{i};q\right)_{2k}\prod_{i=1}^{n}t_{i}^{x-1}\frac{\left(t_{i}q;q\right)_{\infty}}{\left(t_{i}q^{y};q\right)_{\infty}}d_{q}\bm{t}$
$\displaystyle\quad=q^{kx\binom{n}{2}+2k^{2}\binom{n}{3}}\prod_{j=1}^{n}\frac{\Gamma_{q}(x+(j-1)k)\Gamma_{q}(y+(j-1)k)\Gamma_{q}(jk+1)}{\Gamma_{q}(x+y+(n+j-2)k)\Gamma_{q}(k+1)}.$

Conversely, in this paper, we give an explanation to this integral identity in
terms of the skew little $q$-Jacobi polynomials. It is still an open question
to construct classical skew orthogonal polynomials to make explanations of
some other discrete Selberg integrals such as [7, Prop. 5.1]

$\displaystyle\sum_{k_{1},\cdots,k_{r}=-n}^{n}\prod_{1\leq i<j\leq
r}|[k_{i}-k_{j}]_{q}|\prod_{i=1}^{r}q^{(k_{i}+n-r+i)^{2}/2}\left[\begin{array}[]{c}2n\\\
n+k_{i}\end{array}\right]$
$\displaystyle\quad=\frac{r!}{|r|_{q^{1/2}}!}\prod_{i=1}^{r}(-q^{1/2};q^{1/2})_{i}(-q^{i/2+1};q)_{2n-r}\prod_{i=1}^{r}\frac{\Gamma_{q}(1+i/2)}{\Gamma_{q}(3/2)}\frac{\Gamma_{q}(2n+1)\Gamma_{q}(2n-i+5/2)}{\Gamma_{q}(2n-i+2)\Gamma_{q}(2n-i/2+2)}.$

Moreover, combinatoric explanations of $q$-Catalan-Hankel Pfaffians are left
to be further investigated.

## References

  * [1] M. Adler, P. Forrester, T. Nagao and P. van Moerbeke.  _Classical skew orthogonal polynomials and random matrices_.  J. Stat. Phys. 99 (2000), 141–170.
  * [2] M. Adler, E. Horozov and P. van Moerbeke.  _The Pfaff lattice and skew-orthogonal polynomials_.  Internat. Math. Res. Notices 1999, 569–588.
  * [3] G. Andrews, R. Askey and R. Roy.  _Special functions_.  Encyclopedia Math. Appl., vol. 71, Cambridge University Press, Cambridge, 2000.
  * [4] G.E. Andrews and R. Askey.  _Classical orthogonal polynomials._ Polynômes Orthogonaux et Applications, Lecture Notes in Mathematics, 1985, V. 1171, 36–62. MR838970 (88c:33015b)
  * [5] R. Bacher.  _Determinants of matrices related to the Pascal triangle_.  J. Théor. Nombres Bordeaux, 14 (2002), 19-41.
  * [6] E. Bannai and T. Ito.  _Algebraic Combinatorics I: Association Schemes._ Benjamin & Cummings, Menlo Park, CA, 1984. MR882540 (87m:05001)
  * [7] R. Brent, C. Krattenthaler and O. Warnaar.  _Discrete analogues of Macdonald-Mehta integrals._ J. Comb. Theory Ser. A, 144 (2016), 80-126.
  * [8] X. Chang, Y. He, X. Hu and S. Li.  _Partial-skew-orthogonal polynomials and related integrable lattices with Pfaffian tau-functions_.  Commun. Math. Phys., 364 (2018), 1069-1119.
  * [9] T. Chihara.  _A characterization and a class of distribution functions for the Stieltjes-Wigert polynomials_.  Canad. Math. Bull., 13 (1970), 529-532.
  * [10] T.S. Chihara.  _An Introduction to Orthogonal Polynomials_.  Mathematics and its Applications, vol. 13, Gordon and Breach, 1978.
  * [11] J. Christiansen.  _The moment problem associated with the Stieltjes-Wigert polynomials._ J. Math. Anal. Appl., 277 (2003), 218-245.
  * [12] G. Han.  _Hankel continued fraction and its applications._ Adv. Math., 303 (2016), 295-321.
  * [13] P. Forrester.  _Log-gases and random matrices_.  London Mathematical Society Monographs Series, Vol. 34, Princeton University Press, Princeton, NJ, 2010.
  * [14] P. Forrester.  _Meet Andréief, Bordeaux 1886, and Andreev, Kharkov 1882–1883_.  Random Matrices Theory Appl. 8 (2019), 1930001, 9 pp.
  * [15] P. Forrester.  _Global and local scaling limit for the $\beta=2$ Stieltjes-Wigert random matrix ensemble_.  arXiv: 2011.11783.
  * [16] P. Forrester and S. Li.  _Classical discrete symplectic ensembles on the linear and exponential lattice: skew orthogonal polynomials and correlation functions_.  Trans. Amer. Math. Soc., 373 (2020), 665–698.
  * [17] P. Forrester and S. Warnaar.  _The importance of the Selberg integral_.  Bull. Amer. Math. Soc., 45 (2008), 489-534.
  * [18] G.H. Hardy and E.M. Wright.  _An Introduction to The Theory of Numbers._ 4th ed., Oxford University Press, Oxford, 1960.
  * [19] M. Ishikawa, H. Tagawa and J. Zeng.  _Pfaffian decomposition and a Pfaffian analogue of $q$-Catalan Hankel determinants_.  J. Comb. Theory. A, 120 (2013) 1263-1284.
  * [20] M. Ishikawa and J. Zeng.  _Hankel hyperpfaffian calculations and Selberg integrals_.  arXiv: 2008.09776v1.
  * [21] M. Ismail.  _Classical and quantum orthogonal polynomials in one variable_ , with two chapters by W. van Assche.  Encyclopedia of Mathematics and its Applications, Vol. 98, Cambridge University Press, Cambridge, England, 2005.
  * [22] R. Koekoek, P. Lesky and R. Swarttouw.  _Hypergeometric orthogonal polynomials and their $q$-analogues_.  Springer Monographs in Mathematics, Springer-Verlag Berlin Heidelberg 2010.
  * [23] S. Li and G. Yu.  _Christoffel transformations for (partial-)skew-orthogonal polynomials and applications_.  arXiv: 2008.00273, 2020.
  * [24] T. Masuda, K. Mimachi, Y. Kakagami, M. Noumi and K. Ueno.  _Representations of the quantum group $SU_{q}(2)$ and the little $q$-Jacobi polynomials_.  J. Func. Anal., 99 (1991), 357-386.
  * [25] M. Mehta and R. Wang.  _Calculation of a certain determinant_.  Commun. Math. Phys., 214 (2000) 227-232.
  * [26] A. Nikiforov and S. Suslov.  _Classical orthogonal polynomials of a discrete variable on nonuniform lattices_.  Lett. Math. Phys., 11 (1986), 27-34.
  * [27] P. N. Sadjang, W. Koepf and M. Foupouagnigni.  _On moments of classical orthogonal polynomials_.  J. Math. Anal. Appl. 424 (2015) 122–151.
  * [28] T.J. Stieltjes. _Recherches sur les fractions continues_. Ann. Fac. Sci. Toulouse 8 (1894) 1–122; Ann. Fac. Sci. Toulouse 9 (1895) 1–47; English version contained in T.J. Stieltjes, Collected Papers, G. van Dijk (Ed.), Vol. II, Springer, Berlin, 1993.
  * [29] R. Vein and P. Dale.  _Determinants and their applications in mathematical physics._ Applied Mathematical Sciences, Vol 134, Springer-Verlag, New York, 1999.
  * [30] S. Wigert.  _Sur les polynômes orthogonaux et l’approximation des fonctions continues_.  Arkiv för matematik, astronomi och fysik 17 (18) (1923).