# Sparse Inducing Points in Deep Gaussian Processes:
Enhancing Modeling with Denoising Diffusion Variational Inference

Jian Xu  Delu Zeng  John Paisley

###### Abstract

Deep Gaussian processes (DGPs) provide a robust paradigm for Bayesian deep
learning. In DGPs, a set of sparse integration locations called inducing
points are selected to approximate the posterior distribution of the model.
This is done to reduce computational complexity and improve model efficiency.
However, inferring the posterior distribution of inducing points is not
straightforward. Traditional variational inference approaches to posterior
approximation often lead to significant bias. To address this issue, we
propose an alternative method called Denoising Diffusion Variational Inference
(DDVI) that uses a denoising diffusion stochastic differential equation (SDE)
to generate posterior samples of inducing variables. We rely on score matching
methods for denoising diffusion model to approximate score functions with a
neural network. Furthermore, by combining classical mathematical theory of
SDEs with the minimization of KL divergence between the approximate and true
processes, we propose a novel explicit variational lower bound for the
marginal likelihood function of DGP. Through experiments on various datasets
and comparisons with baseline methods, we empirically demonstrate the
effectiveness of DDVI for posterior inference of inducing points for DGP
models.

Gaussian processes, denoising diffusions, variational inference

##  1 Introduction

Deep Gaussian Processes (DGPs) (Damianou & Lawrence, 2013) have emerged as a
robust framework for Bayesian deep learning (Fortuin, 2022) that allows for
flexible modeling of complex functions. DGPs extend the idea of Gaussian
Processes (GPs) (Rasmussen, 2003) to multiple layers, enabling the modeling of
hierarchical structures and capturing intricate dependencies within the data.
A crucial aspect of DGPs is the selection of inducing variables (Titsias,
2009; Snelson & Gharahmani, 2005; Quiñonero-Candela & Rasmussen, 2005), which
are sparse integration locations used to approximate the posterior
distribution of the model. By leveraging these inducing points, DGPs can
efficiently handle large datasets and reduce the computational burden.

Variational inference methods (Blei et al., 2017; Zhang et al., 2018) aim to
approximate the true posterior distribution with a parameterized variational
distribution by minimizing their KL divergence. In the context of DGPs,
traditional variational methods include mean-field Gaussian variational
inference (DSVI) (Salimbeni & Deisenroth, 2017) and Implicit Posterior
Variational Inference (IPVI) (Yu et al., 2019). However, both of these methods
have their limitations and can introduce significant bias when learning the
posterior distribution of inducing points.

DSVI approximates the posterior distribution of inducing points with a simple
Gaussian distribution. Although this approximation is analytically tractable,
it often leads to substantial bias when dealing with nonlinear likelihood
functions. The simplifying assumptions made in the mean-field approximation
can fail to capture the complex dependencies and correlations between the
inducing points, resulting in suboptimal results. On the other hand, IPVI uses
a neural network to parameterize the posterior distribution of inducing
points. Posterior inference is formulated as a Nash equilibrium (Awerbuch et
al., 2008) similar to that of generative adversarial networks (GANs)
(Goodfellow et al., 2014), requiring adversarial learning for the max-max
problem. However, optimizing this objective function can be challenging,
especially when dealing with non-convex neural networks, and lead to
instability during training and contribute to significant bias in the
posterior inference of inducing points (Jenni & Favaro, 2019).

These limitations of traditional variational methods for inference of DGPs
inspires the exploration of alternative approaches. Motivated by the success
of denoising diffusion models (Rombach et al., 2022) in deep learning, we
propose a Denoising Diffusion Variational Inference (DDVI) method that
utilizes the denoising diffusion SDE and incorporates principles similar to
the score matching method (Song et al., 2020) in order to construct the
objective function.

By employing the denoising diffusion SDE, we can accurately capture the
complex dependencies and correlations among the inducing points. Additionally,
similar to the score matching method, we can approximate the intricate score
functions required for accurate posterior inference using a neural network.
This combination finally allows us to explicitly derive a variational lower
bound for the marginal likelihood function by KL divergence minimization,
thereby addressing the bias introduced by traditional variational methods.
Furthermore, DDVI incorporates numerous unique insights, including the well-
developed mathematical theory of SDEs (Anderson, 1982; Haussmann & Pardoux,
1986), the bridge process trick, stochastic optimization techniques,
reparameterization techniques, and gradient backpropagation. These
collectively enable us to efficiently obtain posterior samples from the
denoising diffusion network. As a result, our approach improves not only the
computational efficiency but also ensures stable and reliable training in
DGPs.

In summary, our contributions can be outlined as follows:

  * •

We propose a novel parameterization approach for the posterior distribution of
inducing points in DGPs, utilizing a denoising diffusion process. This method
not only guarantees model efficiency by accurately capturing the complex
dependencies and correlations among the inducing points, but also facilitates
optimization and training.

  * •

We exploit the minimization of KL divergence between the approximate and true
processes to derive an explicit variational lower bound. To efficiently obtain
posterior samples, we employ stochastic optimization and reparameterization
techniques for gradient backpropagation within the denoising diffusion
network.

  * •

Through extensive experiments on various datasets and comparisons with
baseline methods, we empirically demonstrate the effectiveness of the DDVI
method in posterior inference of inducing points for DGP models.

##  2 Method

###  2.1 Model Review

####  2.1.1 Gaussian Process

Consider a random function $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ that maps
$N$ training inputs $\mathbf{X}\triangleq\\{\mathbf{x}_{n}\\}_{n=1}^{N}$ to a
set of noisy observed outputs $\mathbf{y}\triangleq\\{y_{n}\\}_{n=1}^{N}$.
Often, a zero mean Gaussian Process (GP) prior is assumed for the function,
$f\sim\mathcal{GP}(0,k)$, where $k$ denotes the covariance kernel function
$k:\mathbb{R}^{D}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$. Let
$\mathbf{f}\triangleq(f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{N}))^{\top}$
represent the latent function values at the inputs $\mathbf{X}$. The GP prior
assumption then induces a multivariate Gaussian prior over the function
values, expressed as
$p(\mathbf{f})=\mathcal{N}(\mathbf{f}|\mathbf{0},\mathbf{K}_{\mathbf{X}\mathbf{X}})$,
where the covariance matrix $\mathbf{K}_{\mathbf{XX}}$ is defined by
$\left[\mathbf{K}_{\mathbf{XX}}\right]_{ij}=k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$.
The observed outputs $\mathbf{y}$ are then assumed to be contaminated by
i.i.d. noise, modeled as
$p(\mathbf{y}|\mathbf{f})=\mathcal{N}(\mathbf{y}|\mathbf{f},\sigma^{2}\mathbf{I})$,
where $\sigma^{2}$ is the noise variance. The GP posterior distribution of the
latent output $p\left(\mathbf{f}|\mathbf{y}\right)$ has a closed-form
solution. However, the computational cost is $\mathcal{O}(N^{3})$ and the
storage requirement is $\mathcal{O}(N^{2})$, making it challenging to scale to
large datasets without introducing additional techniques.

Sparse methods have been developed that introduce _inducing points_
$\mathbf{Z}=\\{\mathbf{z}_{m}\\}_{m=1}^{M}$ from the input space, along with
corresponding _inducing variables_ :
$\mathbf{u}=\\{f(\mathbf{z}_{m})\\}_{m=1}^{M}$. These methods reduce the
computational complexity to $\mathcal{O}(NM^{2})$. In the _Sparse Gaussian
Processes_ (SGPs) framework, the inducing variables $\mathbf{u}$ and the
function values $\mathbf{f}$ share a joint multivariate Gaussian distribution,
expressed as $p(\mathbf{f},\mathbf{u})=p(\mathbf{f}|\mathbf{u})p(\mathbf{u})$,
with the conditional distribution given by

$p(\mathbf{f}|\mathbf{u})=\mathcal{N}(\mathbf{f}|\mathbf{K}_{\mathbf{X}\mathbf{Z}}\mathbf{K}_{\mathbf{Z}\mathbf{Z}}^{-1}\mathbf{u},\mathbf{K}_{\mathbf{X}\mathbf{X}}-\mathbf{K}_{\mathbf{X}\mathbf{Z}}\mathbf{K}_{\mathbf{Z}\mathbf{Z}}^{-1}\mathbf{K}_{\mathbf{Z}\mathbf{X}})$
(1)

and
$p\left(\mathbf{u}\right)=\mathcal{N}\left(\mathbf{u}|\mathbf{0},\mathbf{K}_{\mathbf{ZZ}}\right)$
is the prior over the outputs of the inducing points.

####  2.1.2 Deep Gaussian Processes

A multi-layer Deep Gaussian Process (DGP) model is a hierarchical composition
of Gaussian Process (GP) models constructed by stacking multiple-output Sparse
GPs (SGPs) together, as described in (Damianou & Lawrence, 2013). Consider a
DGP model with $L$ layers, where each layer $\ell=1,\dots,L$ consists of
$D_{\ell}$ independent random functions. The output of the
$(\ell-1)^{\mathrm{th}}$ layer, denoted as $\mathbf{F}_{\ell-1}$, serves as
the input to the $\ell^{\mathrm{th}}$ layer. Formally, the outputs of the
$\ell^{\mathrm{th}}$ layer are defined as
$\mathbf{F}_{\ell}\triangleq\\{f_{\ell,1}(\mathbf{F}_{\ell-1}),\cdots,f_{\ell,D_{\ell}}(\mathbf{F}_{\ell-1})\\},$
where $f_{\ell,d}\sim\mathcal{GP}(0,k_{\ell})$ for $d=1,\dots,D_{\ell}$, and
$\mathbf{F}_{0}\triangleq\mathbf{X}$. The inducing points and their
corresponding inducing variables for each layer are denoted by
$\boldsymbol{\mathcal{Z}}\triangleq\\{\mathbf{Z}_{\ell}\\}_{\ell=1}^{L}$ and
$\mathcal{U}\triangleq\\{\mathbf{U}_{\ell}\\}_{\ell=1}^{L}$, respectively.
Here,
$\mathbf{U}_{\ell}\triangleq\\{f_{\ell,1}(\mathbf{Z}_{\ell}),\cdots,f_{\ell,D_{\ell}}(\mathbf{Z}_{\ell})\\}$.
Let $\mathcal{F}\triangleq\\{\mathbf{F}_{\ell}\\}_{\ell=1}^{L}$. The design of
the DGP model leads to the following joint model density,

$p(\mathbf{y},\mathbf{F},\mathbf{U})=p\left(\mathbf{y}|\mathbf{F}_{L}\right)\prod_{\ell=1}^{L}{p}(\mathbf{F}_{\ell}|\mathbf{F}_{\ell-1},\mathbf{U}_{\ell})p\left(\mathbf{U}\right).$
(2)

Here we place independent GP priors within and across layers on $\mathbf{U}$,

$p(\mathbf{U})=\prod_{l=1}^{L}{p(\mathbf{U}_{l})}=\prod_{l=1}^{L}{\prod_{d=1}^{D}{\mathcal{N}\left(\mathbf{U}_{\ell,d}|0,\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}\right)}}$
(3)

and the condition similar to Eq. (1) is defined as follows,

$p\left(\mathbf{F}_{\ell}\mid\mathbf{F}_{\ell-1},\mathbf{U}_{\ell}\right)=\prod_{d=1}^{D_{\ell}}\mathcal{N}\left(\mathbf{F}_{\ell,d}\mid\mathbf{\mu}_{\ell,d},\mathbf{\Sigma}_{\ell,d}\right),$
(4)

where we define

$\displaystyle\mathbf{\mu}_{\ell,d}$ $\displaystyle=$
$\displaystyle\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{Z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{U}_{\ell,d}$
$\displaystyle\mathbf{\Sigma}_{\ell,d}$ $\displaystyle=$
$\displaystyle\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{F}_{\ell-1}}-\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{F}_{\ell-1}}.$

In the context of DGPs, traditional variational methods primarily include
mean-field Gaussian variational inference (DSVI) (Salimbeni & Deisenroth,
2017) and Implicit Posterior Variational Inference (IPVI) (Yu et al., 2019).
However, both of these methods have their limitations and can introduce
significant bias in inferring the posterior distribution of inducing points.

DSVI approximates the posterior by a mean-field Gaussian,
$q\left(\mathbf{U}_{\ell,1:D_{\ell}}\right)=\mathcal{N}\left(\mathbf{U}_{\ell,1:D_{\ell}}\mid\mathbf{m}_{\ell,1:D_{\ell}},\mathbf{S}_{\ell,1:D_{\ell}}\right)$,
where $\mathbf{m}_{\ell,1:D_{\ell}}$ and $\mathbf{S}_{\ell,1:D_{\ell}}$ are
variational parameters. However, this assumption is overly strict and may
limit the effectiveness and expressiveness of the model. The likelihood
$p(\mathbf{y}|\mathbf{U})$ is difficult to compute because the latent
functions $\mathbf{F}_{1},\cdots,\mathbf{F}_{L-1}$ are all non-linear kernel
functions.

On the other hand, IPVI utilizes a neural network $\xi$ to parameterize the
posterior distribution of inducing points. Posterior inference is formulated
as a Nash equilibrium (Awerbuch et al., 2008) similar to that of generative
adversarial networks (GANs) (Goodfellow et al., 2014), requiring adversarial
learning for the max-max problem,

$\displaystyle l_{\mathrm{IPVI}}(\xi)$
$\displaystyle=\mathbb{E}_{q_{\xi}(\mathbf{U})}\left[\log
p(\mathbf{y}|\mathbf{U})-D_{\psi^{\star}}(\mathbf{U})\right]$ (5)
$\displaystyle s.t.\quad\psi^{\star}$
$\displaystyle=\mathop{\mathrm{arg}\max}\limits_{\psi}\,\mathbb{E}_{p(\mathbf{U})}[\log\left(1-\sigma\left(D_{\psi}(\mathbf{U})\right)\right]$
$\displaystyle+\mathbb{E}_{q_{\xi}(\mathbf{U})}[\log\sigma\left(D_{\psi}(\mathbf{U})\right)],$

where $D_{\psi}$ is another discriminator network. However, optimizing such an
implict objective $l(\xi)$ can be challenging, especially when dealing with
the non-convex neural network $D_{\psi}$. This can lead to instability during
training and contribute to significant bias in the posterior inference of
inducing points (Jenni & Favaro, 2019).

To address this issue, we propose a novel parameterization approach for the
posterior distribution of inducing variables $\mathbf{U}$ that uses a
denoising diffusion process. This method not only ensures model efficiency by
accurately capturing the complex dependencies and correlations among the
inducing points, but also facilitates optimization and training.

###  2.2 Denoising Diffusion Variational Inference

####  2.2.1 Parameterizing Inducing Point Posteriors

Let $H=D\times M\times L$ denote the dimension of the inducing points. We aim
to sample from the true posterior distribution $q(\mathbf{U})$ in
$\mathbb{R}^{H}$, $q(\mathbf{U})=p(\mathbf{U}|\mathbf{y})$. Following a
similar setup to prior works (Tzen & Raginsky, 2019; Zhang & Chen, 2021;
Vargas et al., 2023), we start by sampling from a fixed distribution
$p_{\mathrm{fix}}$ and then follow a Markov process in which we consider a
sequential latent variable model with a joint distribution denoted as
$\mathcal{Q}\left(\mathbf{U}_{0},\ldots,\mathbf{U}_{T}\right)$ , for step
$t_{s}\in\\{0,..,T-1\\}$,

$\mathbf{U}_{t_{s}+1}\sim\mathcal{T}\left(\mathbf{U}_{t_{s}+1}\mid\mathbf{U}_{t_{s}}\right),\qquad\mathbf{U}_{0}\sim
p_{\mathrm{fix}}$ (6)

Here $\mathcal{T}\left(\mathbf{U}_{t_{s}+1}\mid\mathbf{U}_{t_{s}}\right)$
denotes a transition probability distribution. Through this sequence model, we
use the marginal distribution $\mathcal{Q}\left(\mathbf{U}_{T}\right)$ at the
terminal step $T$ to approximate the true posterior distribution
$q\left(\mathbf{U}_{T}\right)$.

####  2.2.2 Time-reversal Representation of Diffusion SDE

In this paper, we constrain the Markov process
$\mathcal{Q}\left(\mathbf{U}_{0},\ldots,\mathbf{U}_{T}\right)$ to be a time-
reversal process of the following forward noising diffusion stochastic
differential equation (SDE),

$\mathrm{d}\overrightarrow{\mathbf{U}}_{t}=\mathbf{h}(t,\overrightarrow{\mathbf{U}}_{t})\mathrm{d}t+g(t)\mathrm{d}B_{t},\qquad\overrightarrow{\mathbf{U}}_{0}\sim
q,$ (7)

where $\mathbf{h}(t,\cdot):\mathbb{R}^{H}\rightarrow\mathbb{R}^{H}$ is the
drift coefficient, $g(t)\in\mathbb{R}$ is the diffusion coefficient, and
$\left(B_{t}\right)_{t\in[0,T]}$ is an $H$-dimensional Brownian motion. This
diffusion induces the path measure $\mathcal{P}$ on the time interval $[0,T]$
and the marginal density of $\overrightarrow{\mathbf{U}}_{t}$ is denoted
$p_{t}$. Note that by definition we always have $p_{0}=q$ when using an SDE to
perturb this distribution. In DDPM (Ho et al., 2020; Song et al., 2020),
$p_{T}$ is an unstructured prior distribution that contains no information of
$p_{0}$, such as a Gaussian distribution with fixed mean and variance.

From (Anderson, 1982; Haussmann & Pardoux, 1986), the time-reversal
representation of Eq. (7),
$\overleftarrow{\mathbf{U}}_{t}=\overrightarrow{\mathbf{U}}_{T-t}$, where
equality is here in distribution, satisfies

$\displaystyle\mathrm{d}\overleftarrow{\mathbf{U}}_{t}$ $\displaystyle=$
$\displaystyle
g(T-t)^{2}\nabla\ln\big{(}p_{T-t}(\overleftarrow{\mathbf{U}}_{t})\big{)}\mathrm{~{}d}t-\mathbf{h}(T-t,\overleftarrow{\mathbf{U}}_{t})\mathrm{~{}d}t$
$\displaystyle+g(T-t)\mathrm{~{}d}W_{t},$
$\displaystyle\overleftarrow{\mathbf{U}}_{0}$ $\displaystyle\sim$
$\displaystyle p_{T},$ (8)

where $\left(W_{t}\right)_{t\in[0,T]}$ is another $H$-dimensional Brownian
motion. By definition, this time-reversal starts from

$\overleftarrow{\mathbf{U}}_{0}\sim p_{T}\approx p_{\mathrm{fix}}$

and is such that $\overleftarrow{\mathbf{U}}_{T}\sim q$. Since the
distribution of $\overleftarrow{\mathbf{U}}_{T}$ is consistent with the true
posterior $q$ , we can parameterize the transition probability
$\mathcal{T}\left(\mathbf{U}_{t_{s}+1}\mid\mathbf{U}_{t_{s}}\right)$ in the
Euler discretized form of Eq. (2.2.2).

####  2.2.3 Score Matching Technique

This suggests that if we could approximately simulate the diffusion of Eq.
(2.2.2), then we could obtain approximate samples from the target $q$.
However, putting this idea in practice requires being able to approximate the
intractable scores $\nabla\ln\big{(}p_{t}(\cdot)\big{)}$ for $t\in[0,T]$. To
achieve this, DDPM (Ho et al., 2020; Song et al., 2020) rely on score matching
techniques. Specially, to approximate $\mathcal{P}$ consider a path measure
$\mathcal{P}^{\phi}$ whose time-reversal is induced by

$\displaystyle\mathrm{d}\overleftarrow{\mathbf{U}}_{t}^{\phi}$
$\displaystyle=$ $\displaystyle
g(T-t)^{2}s_{\phi}(T-t,\overleftarrow{\mathbf{U}}_{t}^{\phi})\mathrm{d}t-\mathbf{h}(T-t,\overleftarrow{\mathbf{U}}_{t}^{\phi})\mathrm{d}t$
$\displaystyle+~{}g(T-t)\mathrm{~{}d}W_{t},$
$\displaystyle\overleftarrow{\mathbf{U}}_{0}^{\phi}$ $\displaystyle\sim$
$\displaystyle p_{\mathrm{fix}},$ (9)

so that the backward process
$\overleftarrow{\mathbf{U}}_{t}^{\phi}\sim\mathcal{Q}_{t}^{\phi}$, where
$p_{\mathrm{fix}}$ represents a fixed distribution. To obtain
$s_{\phi}(t,\cdot)\approx\nabla\ln\big{(}p_{t}(\cdot)\big{)}$, we parameterize
$s_{\phi}(t,\cdot)$ by a neural network whose parameters are obtained by
minimizing $\mathrm{KL}(\mathcal{P}||\mathcal{P}^{\phi})$. From the chain rule
for the KL divergence (Léonard, 2013) we have,

$\mathrm{KL}(\mathcal{P}||\mathcal{P}^{\phi})=\mathrm{KL}(p_{T}||p_{\mathrm{fix}})+\mathrm{KL}(\mathcal{P}\cdot|\overrightarrow{\mathbf{U}}_{T})||\mathcal{P}^{\phi}(\cdot|\overrightarrow{\mathbf{U}}_{T}^{\phi}))$
(10)

where by the well-known Girsanov Theorem (Oksendal, 2013) and the martingale
property of Itô integrals the second term on the RHS is

$\displaystyle\mathrm{KL}(\mathcal{P}(\cdot|\overrightarrow{\mathbf{U}}_{T})||\mathcal{P}^{\phi}(\cdot|\overrightarrow{\mathbf{U}}_{T}^{\phi}))=\hskip
86.72377pt$ (11)
$\displaystyle\frac{1}{2}\int_{0}^{T}\mathbb{E}_{\mathcal{P}_{t}}\left[g(t)^{2}\left\|\nabla\ln\big{(}p_{t}(\overrightarrow{\mathbf{U}}_{t})\big{)}-\boldsymbol{s}_{\boldsymbol{\phi}}(\overrightarrow{\mathbf{U}}_{t},t)\right\|_{2}^{2}\right]\mathrm{d}t$

From the denoising score matching derivation (Vincent, 2011), this can also be
written as

$\frac{1}{2}\int_{0}^{T}\mathbb{E}\left[g(t)^{2}\|\nabla\ln\big{(}p_{t}(\overrightarrow{\mathbf{U}}_{t}|\overrightarrow{\mathbf{U}}_{0})\big{)}-\boldsymbol{s}_{\boldsymbol{\phi}}(\overrightarrow{\mathbf{U}}_{t},t)\|_{2}^{2}\right]\mathrm{d}t$

plus a constant term, where the expectation is over the joint distribution
$p_{0}(\overrightarrow{\mathbf{U}}_{0})p_{t}(\overrightarrow{\mathbf{U}}_{t}|\overrightarrow{\mathbf{U}}_{0})$.

As the main loss function in DDPM, diffusion-based generative modeling
approaches typically rely on Eq. (10), which involves sampling from $p_{0}$,
the original data such as images, and then backpropagating to estimate the
parameters of the neural network $s_{\phi}$. However, unlike traditional score
matching techniques, this loss function is not applicable to our model since
our $p_{0}$ is the posterior probability $q$ and we only have access to the
joint likelihood, and cannot sample from it. To address this issue, we propose
an alternative approach by minimizing
$\mathrm{KL}(\mathcal{P}^{\phi}||\mathcal{P})$. Analogous to Eq. (10),
considering that we can only obtain samples from $\mathcal{Q}_{t}^{\phi}$, we
have,

$\displaystyle\mathrm{KL}(\mathcal{P}^{\phi}||\mathcal{P})$ $\displaystyle=$
$\displaystyle\mathrm{KL}(\mathcal{Q}^{\phi}||\mathcal{Q})$ $\displaystyle=$
$\displaystyle\mathrm{KL}(p_{\mathrm{fix}}||p_{T})+\mathrm{KL}(\mathcal{Q}^{\phi}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\phi})||\mathcal{Q}(\cdot|\overleftarrow{\mathbf{U}}_{0}))$

where

$\mathrm{KL}(\mathcal{Q}^{\phi}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\phi})||\mathcal{Q}(\cdot|\overleftarrow{\mathbf{U}}_{0}))=\frac{1}{2}\int_{0}^{T}\mathbb{E}\,g(T-t)^{2}\|\varphi(t)\|_{2}^{2}\mathrm{d}t$
(13)

where the expectation is over $\mathcal{Q}_{t}^{\phi}$ and we have defined

$\varphi(t)\triangleq\nabla\ln\big{(}p_{T-t}(\overleftarrow{\mathbf{U}}_{t}^{\phi})\big{)}-\boldsymbol{s}_{\boldsymbol{\phi}}(T-t,\overleftarrow{\mathbf{U}}_{t}^{\phi})$

However, the current challenge we face is that, although we can obtain samples
from $\mathcal{Q}_{t}^{\phi}$ by simulating the SDE (2.2.3), dealing with the
nonlinear drift function of SDE (2.2.3) makes it difficult to obtain
$\nabla\ln\big{(}p_{T-t}(\overleftarrow{\mathbf{U}}_{t}^{\phi})\big{)}$ in Eq.
(13).

####  2.2.4 Bridge Process Trick

Therefore, we propose an alternative approach by constructing a bridge process
$\mathcal{P}^{\mathrm{Bri}}$ to assist in measuring
$\mathrm{KL}(\mathcal{P}^{\phi}||\mathcal{P})$. First, we observe that

$\displaystyle\mathrm{KL(}\mathcal{P}^{\phi}||\mathcal{P})$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\phi}}{\mathrm{d}\mathcal{P}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\phi}}{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}+\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}{\mathrm{d}\mathcal{P}}$

where we represent the stochastic process $\mathrm{KL}$ with the Radon-Nikodym
derivative. Given the specific form in Eq. (2.2.4), we define the bridge
process $\mathcal{P}^{\mathrm{Bri}}$ to follow the diffusion formula as in Eq.
(7), but initialized at
$p_{0}^{\mathrm{Bri}}(\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}})=p_{\mathrm{fix}}$
instead of $q$, which aligns with the distribution of
$\overleftarrow{\mathbf{U}}_{0}$ in Eq. (2.2.3),

$\mathrm{d}\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}=\mathbf{h}(t,\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}})\mathrm{d}t+g(t)\mathrm{d}B_{t},\quad\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}}\sim
p_{\mathrm{fix}}.$ (15)

We typically assume $\mathbf{h}(\cdot,t)$ is affine,
$\mathbf{h}(x,t)=-\mathbf{\lambda}\left(t\right)x$ and
$p_{\mathrm{fix}}=\mathcal{N}(0,\sigma^{2}I)$. Then the transition kernel
$p_{t}(\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}|\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}})$
is always a Gaussian distribution $\mathcal{N}(l_{t},\Sigma_{t})$, where the
mean $l_{t}$ and variance $\Sigma_{t}$ are often known in closed-forms (Särkkä
& Solin, 2019) by

$\displaystyle\frac{\mathrm{d}l_{t}}{\mathrm{d}t}$
$\displaystyle=-\lambda(t)l_{t},\qquad\qquad\qquad\quad\,l_{0}=0$ (16)
$\displaystyle\frac{\mathrm{d}\Sigma_{t}}{\mathrm{d}t}$
$\displaystyle=-2\lambda(t)\Sigma_{t}+g(t)^{2}I,\qquad\Sigma_{0}=\sigma^{2}I$

By the calculations of ordinary differential equations (Hale & Lunel, 2013),
we obtain the following general solution to Eq. (16),

$\displaystyle l_{t}$
$\displaystyle=l_{0}e^{-\int_{0}^{t}\lambda(s)\mathrm{d}s},$ (17)
$\displaystyle\Sigma_{t}$
$\displaystyle=\left(\int_{0}^{t}g(r)^{2}e^{\int_{0}^{r}\lambda(s)\mathrm{d}s}\mathrm{d}rI+\Sigma_{0}\right)e^{-\int_{0}^{t}\lambda(s)\mathrm{d}s}$

According to Eq. (17) we can derive from the Gaussian linear transformation
principle that for any $t$, the distribution $p_{t}^{\mathrm{Bri}}$ of
$\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}$ is a zero-mean Gaussian
distribution,

$\displaystyle
p_{t}^{\mathrm{Bri}}(\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}})$
$\displaystyle=\int
p_{t}(\overrightarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}|\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}})p_{t}(\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}})d\overrightarrow{\mathbf{U}}_{0}^{\mathrm{Bri}}$
(18) $\displaystyle=\mathcal{N}(0,\kappa_{t}I)$

where we have defined the variance

$\kappa_{t}\triangleq\left(\int_{0}^{t}g(r)^{2}e^{\int_{0}^{r}\lambda(s)\mathrm{d}s}\mathrm{d}r+\sigma^{2}\right)e^{-\int_{0}^{t}\lambda(s)\mathrm{d}s}.$

We can write the SDE equation for the reverse process
$\mathcal{Q}^{\mathrm{Bri}}$ of $\mathcal{P}^{\mathrm{Bri}}$ as

$\displaystyle\mathrm{d}\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}$
$\displaystyle=$ $\displaystyle
g(T-t)^{2}\nabla\ln\big{(}p_{T-t}^{\mathrm{Bri}}\big{(}\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}\big{)}\big{)}\mathrm{d}t$
(19)
$\displaystyle-~{}\mathbf{h}(T-t,\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}})\mathrm{d}t$
$\displaystyle+~{}g(T-t)\mathrm{~{}d}W_{t},$
$\displaystyle\overleftarrow{\mathbf{U}}_{0}^{\mathrm{Bri}}$
$\displaystyle\sim$ $\displaystyle p_{T}^{\mathrm{Bri}}.$

According to Eq. (18), we can obtain an analytical expression for the
derivative of the log-likelihood function with respect to
$\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}$,

$\nabla\ln\big{(}p_{T-t}^{\mathrm{Bri}}\big{(}\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}\big{)}\big{)}=-\frac{\overleftarrow{\mathbf{U}}_{t}^{\mathrm{Bri}}}{\kappa_{T-t}}.$
(20)

Next we calculate the value of Eq. (2.2.4). For the first term, according to
the chain rule for KL and Girsanov Theorem (Oksendal, 2013), incorporating
Eqs. (2.2.3, 19, 20), we have

$\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\phi}}{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}=\mathrm{KL}\left(\mathcal{P}^{\phi}\|\mathcal{P}^{\mathrm{Bri}}\right)=\mathrm{KL}\left(\mathcal{Q}^{\phi}\|\mathcal{Q}^{\mathrm{Bri}}\right)$

which can be broken down into the sum of two terms,

$\operatorname{KL}\left(p_{\mathrm{fix}}\|p_{T}^{\mathrm{Bri}}\right)+\mathrm{KL}(\mathcal{Q}^{\phi}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\phi})||\mathcal{Q}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\mathrm{Bri}}))$
(21)

where

$\mathrm{KL}(\mathcal{Q}^{\phi}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\phi})||\mathcal{Q}(\cdot|\overleftarrow{\mathbf{U}}_{0}^{\mathrm{Bri}}))=\hskip
108.405pt$
$\frac{1}{2}\int_{0}^{T}\mathbb{E}_{\mathcal{Q}_{t}^{\phi}}~{}g(T-t)^{2}\Big{\|}\frac{\overleftarrow{\mathbf{U}}_{t}^{\mathrm{\phi}}}{\kappa_{T-t}}+\boldsymbol{s}_{\boldsymbol{\phi}}(T-t,\overleftarrow{\mathbf{U}}_{t}^{\phi})\Big{\|}_{2}^{2}~{}\mathrm{d}t$

At this point, we can simulate the SDE (2.2.3) to compute the first term in
Eq. (2.2.4). The integral term can be computed using either ODE solvers (Chen
et al., 2018) or by employing Riemann summation methods. For the second term,
$\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}{\mathrm{d}\mathcal{P}}$,
we can see from Eq. (7) and Eq. (15) that $\mathcal{P}$ and
$\mathcal{P}^{\mathrm{Bri}}$ have the same dynamic system $\tau$, except for
different initial values. Therefore, we have

$\displaystyle\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}{\mathrm{d}\mathcal{P}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathcal{P}^{\mathrm{Bri}}\left(\tau|\cdot\right)p_{0}^{\mathrm{Bri}}\left(\cdot\right)}{\mathcal{P}\left(\tau|\cdot\right)p_{0}\left(\cdot\right)}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log\frac{p_{0}^{\mathrm{Bri}}\left(\cdot\right)}{p_{0}\left(\cdot\right)}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log\frac{p_{\mathrm{fix}}}{q}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log\frac{p_{\mathrm{fix}}}{p(\mathbf{y}|\cdot)p(\cdot)}+\log
p(\mathbf{y})$

####  2.2.5 A New Evidence Lower Bound

Let
$l_{1}(\phi)=\mathbb{E}_{\mathcal{P}^{\phi}}\log\frac{\mathrm{d}\mathcal{P}^{\phi}}{\mathrm{d}\mathcal{P}^{\mathrm{Bri}}}$.
Combining Eqs. (2, 2.2.4, 21, 2.2.4), we obtain a new variational lower bound
$l(\phi)$ for the marginal likelihood $\log p(\mathbf{y})$ of our method,

$\displaystyle\log p(\mathbf{y})={}$
$\displaystyle\mathrm{KL(}\mathcal{P}^{\phi}||\mathcal{P})-l_{1}(\phi)-\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log\frac{p_{\mathrm{fix}}}{p(\mathbf{y}|\mathbf{\cdot})p(\mathbf{\cdot})}$
(23) $\displaystyle={}$
$\displaystyle\mathrm{KL(}\mathcal{P}^{\phi}||\mathcal{P})-l_{1}(\phi)-\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log{p_{\mathrm{fix}}}$
$\displaystyle+\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log
p(\cdot)+\mathbb{E}_{\mathcal{Q}_{T}^{\phi},\mathbf{F_{1}},...,\mathbf{F_{L}}}\log
p(\mathbf{y}|\mathbf{F_{L}})$ $\displaystyle\geqslant{}$
$\displaystyle\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log
p(\cdot)+\mathbb{E}_{\mathcal{Q}_{T}^{\phi},\mathbf{F_{1}},...,\mathbf{F_{L}}}\log
p(\mathbf{y}|\mathbf{F_{L}})$ $\displaystyle-
l_{1}(\phi)-\mathbb{E}_{\mathcal{Q}_{T}^{\phi}}\log{p_{\mathrm{fix}}}$
$\displaystyle={}$ $\displaystyle l(\phi)$

In our derivation, $p(\cdot)$ represents the prior function of $\mathbf{U}$.
By introducing a new variational lower bound for $\log p(\mathbf{y})$, our
proposed model, compared to the initial mean-field variational inference model
(DSVI) where $q$ is approximated by a Gaussian distribution, approximates the
posterior distribution through a denoising diffusion process. The flexibility
of the denoising neural network $\phi$ intuitively suggests that our model has
an advantage in approximating the posterior distribution. On the other hand,
compared to IPVI, DDVI provides an explicit evidence lower bound (ELBO), which
means it is easier to train and allows for efficient backpropagation.

###  2.3 Reparameterization Trick and SGD

Algorithm 1 Denoising Diffuison Variational Inference (DDVI) algorithm for
DGPs

Input: training data $\mathbf{X},\mathbf{y}$ mini-batch size $B$

Initialize diffusion coefficient $h,g$, all DGP hyperparameters $\gamma$,
denoising diffusion network parameters $\phi$ and set $l_{0}=0$

repeat

for $t_{s}=0$ to $T-1$ do

Draw $\boldsymbol{\epsilon}_{t_{s}}\sim\mathcal{N}(0,I)$ and set
$\mathbf{U}_{t_{s}+1}=\mathbf{U}_{t_{s}}-\mathbf{h}(\mathbf{U}_{t_{s}},T-t_{s})+g(T-t_{s})^{2}s_{\phi}\left(T-t_{s},\mathbf{U}_{t_{s}}\right)+g(T-t)\boldsymbol{\epsilon}_{t_{s}}$

Compute $\kappa_{T-(t_{s}+1)}$ by Eq. (18) and set
$l_{t_{s}+1}=l_{t_{s}}+g(T-(t_{s}+1))^{2}\|\frac{\mathbf{U}_{t_{s}+1}}{\kappa_{T-(t_{s}+1)}}+\boldsymbol{s}_{\boldsymbol{\phi}}(T-(t_{s}+1),\mathbf{U}_{t_{s}+1})\|_{2}^{2}$

end for

Sample mini-batch indices $I\subset\\{1,\ldots,N\\}\text{ with }|I|=B$ and set
$\\{\mathbf{U}_{\ell,1},...,\mathbf{U}_{\ell,D}\\}_{\ell=1}^{L}=\mathbf{U}_{T}$

for $\ell=1$ to $L$ do

Draw $\epsilon_{\ell,d}\sim\mathcal{N}(0,I)$ and calculate
$\mathbf{F}_{\ell,d}=\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{Z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{U}_{\ell,d}+\sqrt{\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{F}_{\ell-1}}-\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{F}_{\ell-1}}}\boldsymbol{\epsilon}_{\ell,d}$

end for

Set $l(\phi,\gamma)=-\log p_{\mathrm{fix}}\left(\mathbf{U}_{T}\right)+\log
p\left(\mathbf{U}_{T}\right)+\frac{N}{B}\log
p\left(\mathbf{y}_{I}\mid\mathbf{F}_{L}\right)-{\mathrm{KL}\left(p_{\mathrm{fix}}\parallel\mathcal{N}(0,\kappa_{T})\right)}-\frac{1}{2}l_{T}$

Make a gradient descent update of $l(\phi,\gamma)$

until $\phi,\gamma$ converge

For ease of sampling, we consider a reparameterization version of Eq. (23)
based on the approaximate transition probability
$\mathcal{T}_{\phi}\left(\mathbf{U}_{t_{s}+1}\mid\mathbf{U}_{t_{s}}\right)$
given by

$\displaystyle\begin{aligned}
\mathcal{T}_{\phi}(\mathbf{U}_{t_{s}+1})=~{}&\mathbf{U}_{t_{s}}-\mathbf{h}(\mathbf{U}_{t_{s}},T-t_{s})~{}+\\\
&g(T-t_{s})^{2}s_{\phi}\left(T-t_{s},\mathbf{U}_{t_{s}}\right)+g(T-t)\boldsymbol{\epsilon}_{t_{s}}\end{aligned}$
(24)

where $\boldsymbol{\epsilon}_{t_{s}}\sim\mathcal{N}(0,I)$. Given that
$\mathbf{U}_{t_{s}+1}=\mathcal{T}_{\phi}(\mathbf{U}_{t_{s}})$, we have a
representation of $\mathbf{U}_{t_{s}}$ by a stochastic flow,

$\displaystyle\mathbf{U}_{t_{s}+1}=\mathcal{T}_{\phi}\left(\mathbf{U}_{t_{s}}\right)=\mathcal{T}_{\phi}\circ\mathcal{T}_{\phi}\circ\cdots\mathcal{T}_{\phi}(\mathbf{U}_{0}).$
(25)

Moreover, for DGP models, we also have a reparameterization version (Salimbeni
& Deisenroth, 2017) of the conditional distribution in Eq. (4) of the form

$\displaystyle\begin{aligned}
\mathbf{F}_{\ell,d}=~{}&\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{Z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{U}_{\ell,d}\\\
&+\sqrt{\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{F}_{\ell-1}}-\mathbf{K}_{\mathbf{F}_{\ell-1}\mathbf{z}_{\ell}}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{Z}_{\ell}}^{-1}\mathbf{K}_{\mathbf{Z}_{\ell}\mathbf{F}_{\ell-1}}}\boldsymbol{\epsilon}_{\ell,d}\end{aligned}$
(26)

where $\boldsymbol{\epsilon}_{\ell,d}\in\mathbb{R}^{N}$ are standard Gaussian
random variables. In order to accelerate training and sampling in our
inference scheme, we propose a scalable variational bound that is tractable in
the large data regime based on stochastic variational inference (Kingma &
Welling, 2013; Hoffman & Blei, 2015; Salimbeni & Deisenroth, 2017; Naesseth et
al., 2020) and stochastic gradient descent (Welling & Teh, 2011; Chen et al.,
2014; Zou et al., 2019; Alexos et al., 2022). Specially, instead of computing
the full log likelihood, we use a stochastic variant to subsample datasets
into a mini-batches $\mathcal{D}_{I}$ with $\left|\mathbf{X}_{I}\ \right|=B$,
where $I\subset\\{1,2,..,N\\}$ is the index of a mini-batch. We present the
resulting stochastic inference for our Denoising Diffuison Variational
Inference algorithm for DGP models in Algorithm 1.

###  2.4 Predictive Distribution

To obtain the final layer density for making predictions, we first sample from
the optimized generator and transform the input locations $\mathbf{x}$ to the
test locations $\mathbf{x}^{\star}$ using Eq. (2). We subsequently compute the
function values at the test locations, which are represented as
$\mathbf{F}_{\ell}^{\star}$. Finally, we use the equation below to estimate
the density of the final layer, which enables us to make predictions for the
test data

$q(\mathbf{F}_{L}^{\star})=\int{{\prod_{d,\ell}{p(\mathbf{F}_{\ell,d}^{\star}|\mathbf{F}_{\ell-1}^{\star},\mathbf{U}_{\ell,d})\mathcal{Q}_{T}^{\phi}\left(\mathbf{U}_{\ell,d}\right)d\mathbf{F}_{\ell-1}^{\star}d\mathbf{U}_{\ell,d}}}}$

where $\mathcal{Q}_{T}^{\phi}$ represents the output of the denoising
diffusion process at time $T$ and the first term of the integral
$p(\mathbf{F}_{\ell,d}^{\star}|\mathbf{F}_{\ell-1}^{\star},\mathbf{U}_{\ell,d})$
is conditional Gaussian. We leverage this to draw samples from
$q\left(\mathbf{F}_{L}^{\star}\right)$ and further perform the sampling
according to the problem considered.

##  3 Experiments

Figure 1: Regression test RMSE results by our DDVI method (red), SGHMC (blue),
IPVI(green) and DSVI (black) for DGPs on ten UCI benchmark datasets. The
numbers 2, 3, 4, and 5 represent the layers of DGP methods. Lower is better.
The mean is shown with error bars of one standard error. The dimensions of the
data are displayed above each subgraph.

Figure 2: Regression test mean NLL results by our DDVI method (red), SGHMC
(blue), IPVI (green) and DSVI (black) for DGPs on ten UCI benchmark datasets.
The numbers 2, 3, 4, and 5 represent the layers of DGP methods. Lower is
better. The mean is shown with error bars of one standard error. The
dimensions of the data are displayed above each subgraph.

###  3.1 Baseline Models and Hyperparameter Settings

In order to evaluate the performance of our proposed method, we conducted
empirical evaluations on real-world datasets for both regression and
classification tasks, with both small and large datasets. We compare against
several other models, including Doubly Stochastic VI (DSVI) (Salimbeni &
Deisenroth, 2017), Implicit Posterior VI (IPVI) (Yu et al., 2019), and the
state-of-the-art SGHMC model (Havasi et al., 2018). All experiments were
conducted with the same hyper-parameters and initializations whenever possible
to obtain a fair comparison.

We constructed a random $0.9/0.1$ train/test split and normalized the features
of our datasets to the range $[-1,1]$. The depth $L$ of DGP models varied from
$2$ to $5$, with $128$ inducing points per layer, which were initialized by
sampling Gaussian random variables. The output dimension for each hidden layer
is set to $1$ for the final layer and the dimensionality of the data for all
others. We use the RBF kernel for all tasks. For all datasets, we have
optimized hyper-parameters and network parameters jointly and set learning
rate to $0.01$ using Adam optimizer (Kingma & Ba, 2014). We trained all models
on all datasets until convergence was achieved. In each experiment, we
repeated the process 10 times and reported the mean and standard deviation of
the metrics. The selection of the denoising diffusion networks are done
manually by the classical grid search approach for each experimental dataset.
It is worth mentioning that this work also benefits from the contributions of
the PyTorch platform, GPyTorch (Gardner et al., 2018), and related work on
neural SDE solvers (Li et al., 2020; Kidger et al., 2021). All our experiments
were conducted on an RTX 4090 GPU.

###  3.2 Regression Task

In our experiments, we evaluated the performance of the DDVI model on ten UCI
regression datasets, which varied in size from 308 to 2,055,733 data points.
We used the mean RMSE and mean NLL (Gneiting & Raftery, 2007) of the test data
as the performance metric, and the results are presented in Figure 1 and
Figure 2.

As shown in these two figures, our DDVI method consistently achieves
competitive results compared to three baselines on the majority of datasets.
This is attributed to the key advantages of our approach, which overcomes
limitations present in previous methods as discussed in the main text. Our
findings also suggest that deeper DGP models tend to perform better. It is
worth mentioning that the difference in performance may be attributed to the
nature of the datasets, such as their size or the presence of outliers or
singular values.

Using mini-batch algorithm and GPU acceleration, our method can also be
extended to larger datasets. Our evaluation of the performance of DDVI in
Figures 1 and 2 is also conducted on two real-world large-scale regression
datasets: the YearMSD dataset and the Airline dataset. The YearMSD dataset has
a large input dimension of $D=90$ and a data size of approximately 500,000.
The Airline dataset, on the other hand, has an input dimension of $D=8$ and a
large data size of approximately 2 million. For the YearMSD dataset, we split
the data into training and test sets, using the first 463,810 examples as
training data and the last 51,725 examples as test data. Similarly, for the
Airline dataset, we take the first 700K points for training and next 100K for
testing.

###  3.3 Image and Large-Scale Dataset Classification

We evaluate our method on multiclass classification tasks using the MNIST
(LeCun et al., 1998), Fashion-MNIST (Xiao et al., 2017), and CIFAR-10
(Krizhevsky et al., 2009) datasets. The first two datasets consist of
grayscale images of size $28\times 28$ pixels, while CIFAR-10 comprises
colored images of size $32\times 32$ pixels. The results are presented in
Table 1. We note that our method outperforms the other three methods on all
three datasets, with significantly less training time. Specially, for CIFAR-10
dataset, we utilize the convolutional layers of ResNet-20 (He et al., 2016) as
our feature extractor (Wilson et al., 2016) and achieve a remarkable accuracy
of 95.56 on the test set. Additionally, we evaluate our approach using two
large-scale classification datasets, the Higgs datase and the SUSY dataset,
which are presented in Table 2.

Table 1: Mean test accuracy ($\%$) and training details achieved by DSVI,
SGHMC, IPVI and DDVI (ours) DGP model for three image classification datasets.
Results are shown for 3 and 4 layers as indicated, and runtime is given per
iteration.

Data Set | Model | Time$3$ | Iter$3$ | Acc$3$ | Time4 | Iter4 | Acc4
---|---|---|---|---|---|---|---
| DSVI | 0.34s | 20K | 97.17 | 0.54s | 20K | 97.41
MNIST | IPVI | 0.49s | 20K | 97.58 | 0.62s | 20K | 97.80
| SGHMC | 1.14s | 20K | 97.25 | 1.22s | 20K | 97.55
| DDVI | 0.38s | 20K | 98.84 | 0.50s | 20K | 99.01
| DSVI | 0.34s | 20K | 87.45 | 0.50s | 20K | 87.99
Fashion | IPVI | 0.48s | 20K | 88.23 | 0.61s | 20K | 88.90
| SGHMC | 1.21s | 20K | 86.88 | 1.25s | 20K | 87.08
| DDVI | 0.40s | 20K | 90.36 | 0.55s | 20K | 90.85
| DSVI | 0.43s | 20K | 91.47 | 0.66s | 20K | 91.79
CIFAR-10 | IPVI | 0.62s | 20K | 92.79 | 0.78s | 20K | 93.52
| SGHMC | 8.04s | 20K | 92.62 | 8.61s | 20K | 92.94
| DDVI | 0.45s | 20K | 95.23 | 0.69s | 20K | 95.56

Table 2: Test AUC values for large-scale classification datasets. Uses random
90% / 10% training and test splits.

|  | SUSY | HIGGS
---|---|---|---
| $N$ | 5,500,000 | 11,000,000
| $D$ | 18 | 28
DSVI $M=128$ | $L=2$ | 0.876 | 0.830
$L=3$ | 0.877 | 0.837
$L=4$ | 0.878 | 0.841
$L=5$ | 0.878 | 0.846
IPVI $M=128$ | $L=2$ | 0.879 | 0.843
$L=3$ | 0.882 | 0.847
$L=4$ | 0.883 | 0.850
$L=5$ | 0.883 | 0.852
SGHMC $M=128$ | $L=2$ | 0.879 | 0.842
$L=3$ | 0.881 | 0.846
$L=4$ | 0.883 | 0.850
$L=5$ | 0.884 | 0.853
DDVI $M=128$ | $L=2$ | 0.883 | 0.849
$L=3$ | 0.885 | 0.852
$L=4$ | 0.887 | 0.856
$L=5$ | 0.886 | 0.857

Table 3: Mean RMSE and NLL achieved by DSVI, SGHMC, IPVI and DDVI (ours) GPLVM
model for data recovery task. Standard deviation is shown in parentheses.
Runtime is given per iteration.

Data Set | Model | Time | Iter | RMSE | NLL
---|---|---|---|---|---
| DSVI | 0.32s | 20K | 8.32 (0.2) | 1.49 (0.02)
Frey Faces | IPVI | 0.42s | 20K | 7.91 (0.4) | 1.33 (0.02)
| SGHMC | 1.13s | 20K | 7.95 (0.3) | 1.36 (0.03)
| DDVI | 0.36s | 20K |  7.64 (0.2) |  1.17 (0.01)

Figure 3: The Brendan faces reconstruction task with 75% missing pixels. The
top row represents the ground truth data and the bottom row showcases the
reconstructions from the 20-dimensional latent distribution.

###  3.4 Unsupervised Learning for Data Recovery Task

We conducted a reconstruction experiment on Frey Faces Data (Roweis & Saul,
2000), focusing on how models capture uncertainty when training with missing
data in structured inputs. We used the entire dataset with a latent variable
dimensionality of 20. The image data set contains 1965 images of a face taken
from sequential frames of a short video. Each image is of size 20$\times$28
yielding a 560 dimensional data space. In both cases, we chose 5% of the
training set as missing data samples and removed 75% of their pixels, seeking
to recover their original appearance. Figure 3 summarize the samples generated
from the learned latent distribution. This reconstruction experiment is
performed using the Gaussian Process Latent Variable Model (GPLVM) (Titsias &
Lawrence, 2010) and is similar to the related work by (Gal et al., 2014).

To demonstrate the effectiveness of our method in producing more accurate
likelihoods on image datasets, we present in Table 3 negative log-likelihood,
and RMSE for reconstructed images on the Frey Faces, comparing with baseline
methods. The results show that our method converges to higher likelihoods and
lower RMSE, indicating superior performance in high-dimensional and multi-
modal image data. This suggests that adding DDVI method can also improve the
convergence of the traditional GPLVM methods.

##  4 Conclusion

We have introduced Denoising Diffusion Variational Inference (DDVI) as an
alternative approach for inferring the posterior distribution of inducing
points in Deep Gaussian Processes (DGPs). By employing a denoising diffusion
stochastic differential equation (SDE) and utilizing the score matching
method, we are able to accurately approximate challenging score functions
using a neural network. Through extensive experiments and comparisons with
baseline methods on various datasets, we demonstrated the effectiveness of
DDVI in posterior inference of inducing points for DGP models. The DDVI method
addressed the limitations of traditional variational inference techniques,
reducing biases and improving accuracy in the posterior approximation. Our
proposed DDVI approach not only enhances computational efficiency, but also
provides a more robust framework for Bayesian deep learning with DGPs.

## Acknowledgements

The work is supported by the Fundamental Research Program of Guangdong, China,
under Grant 2023A1515011281; and in part by the National Natural Science
Foundation of China under Grant 61571005.

## Impact Statement

This paper aims to contribute to the advancement of the Machine Learning
field. Our work may have various potential societal implications, none of
which we believe need to be specifically emphasized here.

## References

  * Alexos et al. (2022) Alexos, A., Boyd, A. J., and Mandt, S.  Structured stochastic gradient mcmc.  In _International Conference on Machine Learning_ , pp. 414–434. PMLR, 2022.
  * Anderson (1982) Anderson, B. D.  Reverse-time diffusion equation models.  _Stochastic Processes and their Applications_ , 12(3):313–326, 1982.
  * Awerbuch et al. (2008) Awerbuch, B., Azar, Y., Epstein, A., Mirrokni, V. S., and Skopalik, A.  Fast convergence to nearly optimal solutions in potential games.  In _Proceedings of the 9th ACM conference on Electronic commerce_ , pp. 264–273, 2008.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., and McAuliffe, J. D.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Chen et al. (2018) Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018.
  * Chen et al. (2014) Chen, T., Fox, E., and Guestrin, C.  Stochastic gradient hamiltonian monte carlo.  In _International conference on machine learning_ , pp. 1683–1691. PMLR, 2014.
  * Cheng & Boots (2017) Cheng, C.-A. and Boots, B.  Variational inference for gaussian process models with linear complexity.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Damianou & Lawrence (2013) Damianou, A. and Lawrence, N.  Deep Gaussian processes.  In _Conference on Artificial Intelligence and Statistics_ , pp. 207–215, 2013.
  * Fortuin (2022) Fortuin, V.  Priors in bayesian deep learning: A review.  _International Statistical Review_ , 90(3):563–591, 2022.
  * Gal et al. (2014) Gal, Y., Van Der Wilk, M., and Rasmussen, C. E.  Distributed variational inference in sparse gaussian process regression and latent variable models.  _Advances in neural information processing systems_ , 27, 2014.
  * Gardner et al. (2018) Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G.  Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.  _Advances in neural information processing systems_ , 31, 2018.
  * Gneiting & Raftery (2007) Gneiting, T. and Raftery, A. E.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Goodfellow et al. (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27, 2014.
  * Hale & Lunel (2013) Hale, J. K. and Lunel, S. M. V.  _Introduction to functional differential equations_ , volume 99.  Springer Science & Business Media, 2013.
  * Haussmann & Pardoux (1986) Haussmann, U. G. and Pardoux, E.  Time reversal of diffusions.  _The Annals of Probability_ , pp. 1188–1205, 1986.
  * Havasi et al. (2018) Havasi, M., Hernández-Lobato, J. M., and Murillo-Fuentes, J. J.  Inference in deep Gaussian processes using stochastic gradient Hamiltonian Monte Carlo.  In _Conference on Neural Information Processing Systems_ , pp. 7517–7527, 2018.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proc. CVPR_ , pp. 770–778, 2016.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Hoffman & Blei (2015) Hoffman, M. D. and Blei, D. M.  Structured stochastic variational inference.  In _Artificial Intelligence and Statistics_ , pp. 361–369, 2015.
  * Jenni & Favaro (2019) Jenni, S. and Favaro, P.  On stabilizing generative adversarial training with noise.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 12145–12153, 2019.
  * Kidger et al. (2021) Kidger, P., Foster, J., Li, X. C., and Lyons, T.  Efficient and accurate gradients for neural sdes.  _Advances in Neural Information Processing Systems_ , 34:18747–18761, 2021.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma & Welling (2013) Kingma, D. P. and Welling, M.  Auto-encoding variational Bayes.  In _International Conference on Learning Representations_ , 2013.
  * Krizhevsky et al. (2009) Krizhevsky, A., Hinton, G., et al.  Learning multiple layers of features from tiny images.(2009), 2009.
  * Lalchand & Rasmussen (2020) Lalchand, V. and Rasmussen, C. E.  Approximate inference for fully bayesian gaussian process regression.  In _Symposium on Advances in Approximate Bayesian Inference_ , pp. 1–12. PMLR, 2020.
  * LeCun et al. (1998) LeCun, Y., Cortes, C., and Burges, C.  Mnist handwritten digit database, 1998.  _URL http://www. research. att. com/~ yann/ocr/mnist_ , 1998.
  * Léonard (2013) Léonard, C.  A survey of the schr$\backslash$” odinger problem and some of its connections with optimal transport.  _arXiv preprint arXiv:1308.0215_ , 2013.
  * Li et al. (2020) Li, X., Wong, T.-K. L., Chen, R. T., and Duvenaud, D.  Scalable gradients for stochastic differential equations.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 3870–3882. PMLR, 2020.
  * Ma et al. (2019) Ma, C., Li, Y., and Hernández-Lobato, J. M.  Variational implicit processes.  In _International Conference on Machine Learning_ , pp. 4222–4233. PMLR, 2019.
  * Mescheder et al. (2017) Mescheder, L., Nowozin, S., and Geiger, A.  Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks.  In _International conference on machine learning_ , pp. 2391–2400. PMLR, 2017.
  * Naesseth et al. (2020) Naesseth, C., Lindsten, F., and Blei, D.  Markovian score climbing: Variational inference with kl (p—— q).  _Advances in Neural Information Processing Systems_ , 33:15499–15510, 2020.
  * Oksendal (2013) Oksendal, B.  _Stochastic differential equations: an introduction with applications_.  Springer Science & Business Media, 2013.
  * Ortega et al. (2022) Ortega, L. A., Santana, S. R., and Hernández-Lobato, D.  Deep variational implicit processes.  _arXiv preprint arXiv:2206.06720_ , 2022.
  * Quiñonero-Candela & Rasmussen (2005) Quiñonero-Candela, J. and Rasmussen, C. E.  A unifying view of sparse approximate Gaussian process regression.  _Journal of Machine Learning Research_ , 6:1939–1959, 2005.
  * Rasmussen (2003) Rasmussen, C. E.  Gaussian processes in machine learning.  In _Summer school on machine learning_ , pp. 63–71. Springer, 2003.
  * Rombach et al. (2022) Rombach, R., Blattmann, A., Lorenz, D., Esser, P., and Ommer, B.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 10684–10695, 2022.
  * Rossi et al. (2021) Rossi, S., Heinonen, M., Bonilla, E., Shen, Z., and Filippone, M.  Sparse gaussian processes revisited: Bayesian approaches to inducing-variable approximations.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1837–1845. PMLR, 2021.
  * Roweis & Saul (2000) Roweis, S. T. and Saul, L. K.  Nonlinear dimensionality reduction by locally linear embedding.  _science_ , 290(5500):2323–2326, 2000.
  * Salimbeni & Deisenroth (2017) Salimbeni, H. and Deisenroth, M.  Doubly stochastic variational inference for deep Gaussian processes.  In _Conference on Neural Information Processing Systems_ , pp. 4588–4599, 2017.
  * Salimbeni et al. (2018) Salimbeni, H., Cheng, C.-A., Boots, B., and Deisenroth, M.  Orthogonally decoupled variational gaussian processes.  _Advances in neural information processing systems_ , 31, 2018.
  * Särkkä & Solin (2019) Särkkä, S. and Solin, A.  _Applied stochastic differential equations_ , volume 10.  Cambridge University Press, 2019.
  * Shi et al. (2020) Shi, J., Titsias, M., and Mnih, A.  Sparse orthogonal variational inference for gaussian processes.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1932–1942. PMLR, 2020.
  * Snelson & Gharahmani (2005) Snelson, E. L. and Gharahmani, Z.  Sparse Gaussian processes using pseudo-inputs.  In _Conference on Neural Information Processing Systems_ , pp. 1257–1264, 2005.
  * Song et al. (2020) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * Sun et al. (2019) Sun, S., Zhang, G., Shi, J., and Grosse, R.  Functional variational bayesian neural networks.  _arXiv preprint arXiv:1903.05779_ , 2019.
  * Sun et al. (2021) Sun, S., Shi, J., Wilson, A. G., and Grosse, R.  Scalable variational gaussian processes via harmonic kernel decomposition.  _arXiv preprint arXiv:2106.05992_ , 2021.
  * Titsias & Lawrence (2010) Titsias, M. and Lawrence, N. D.  Bayesian gaussian process latent variable model.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pp. 844–851. JMLR Workshop and Conference Proceedings, 2010.
  * Titsias (2009) Titsias, M. K.  Variational learning of inducing variables in sparse Gaussian processes.  In _Conference on Artificial Intelligence and Statistics_ , pp. 567–574, 2009.
  * Tzen & Raginsky (2019) Tzen, B. and Raginsky, M.  Theoretical guarantees for sampling and inference in generative models with latent diffusions.  In _Conference on Learning Theory_ , pp. 3084–3114. PMLR, 2019.
  * Vargas et al. (2023) Vargas, F., Grathwohl, W., and Doucet, A.  Denoising diffusion samplers.  _arXiv preprint arXiv:2302.13834_ , 2023.
  * Vincent (2011) Vincent, P.  A connection between score matching and denoising autoencoders.  _Neural computation_ , 23(7):1661–1674, 2011.
  * Welling & Teh (2011) Welling, M. and Teh, Y. W.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th international conference on machine learning (ICML-11)_ , pp. 681–688, 2011.
  * Wilson et al. (2016) Wilson, A. G., Hu, Z., Salakhutdinov, R. R., and Xing, E. P.  Stochastic variational deep kernel learning.  _Advances in neural information processing systems_ , 29, 2016.
  * Xiao et al. (2017) Xiao, H., Rasul, K., and Vollgraf, R.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  _arXiv preprint arXiv:1708.07747_ , 2017.
  * Yu et al. (2019) Yu, H., Chen, Y., Low, B. K. H., Jaillet, P., and Dai, Z.  Implicit posterior variational inference for deep gaussian processes.  _Advances in neural information processing systems_ , 32, 2019.
  * Zhang et al. (2018) Zhang, C., Bütepage, J., Kjellström, H., and Mandt, S.  Advances in variational inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(8):2008–2026, 2018.
  * Zhang & Chen (2021) Zhang, Q. and Chen, Y.  Path integral sampler: a stochastic control approach for sampling.  _arXiv preprint arXiv:2111.15141_ , 2021.
  * Zou et al. (2019) Zou, D., Xu, P., and Gu, Q.  Stochastic gradient hamiltonian monte carlo methods with recursive variance reduction.  _Advances in Neural Information Processing Systems_ , 32, 2019.

##  Appendix A Additional Related Works

A line of closely related work focusing on enhancing variational posteriors
with decoupled/orthogonal inducing points has garnered significant attention
in recent years (Cheng & Boots, 2017; Salimbeni et al., 2018; Shi et al.,
2020; Sun et al., 2021). These methods present variational Gaussian process
models that segregate the representation of mean and covariance functions
within the reproducing kernel Hilbert space. This novel parametrization
extends previous models and allows for solving the variational inference
problem using stochastic gradient ascent with linear time and space complexity
in the number of mean function parameters. In contrast to these approaches,
our work diverges in its emphasis on precise posterior inference for Gaussian
process models rather than complexity analysis concerning inducing points.

Another line of related work is on the idea of fully Bayesian Gaussian
processes (Lalchand & Rasmussen, 2020; Rossi et al., 2021), particularly
(Rossi et al., 2021) have applied MCMC-related methods and the ideas of fully
Bayesian approaches to deep GPs, achieving significant improvements. Our work,
on the other hand, focuses on improvements and advancements in the field of
variational inference.

Additionally, another line of related work is on implict stochastic processes
(Ma et al., 2019), that is not restricted to Gaussian predictive
distributions, which is closely related to function-space models (Sun et al.,
2019; Mescheder et al., 2017). Furthermore, (Ortega et al., 2022) further
extends the idea of DGP models to implict stochastic processes, enhancing the
flexibility of the models. We believe that this is a valuable complement to
our approach, and future work may involve integrating our method with these
advancements to explore more useful applications.


# Enumeration of linear codes with different hulls

Iliya Bouyukliev and Stefka Bouyuklieva I. Bouyukliev is with the Institute of
Mathematics and Informatics, Bulgarian Academy of Sciences, P.O. Box 323,
Veliko Tarnovo, BULGARIA. Email<EMAIL_ADDRESS>Bouyuklieva is with the
Faculty of Mathematics and Informatics, St. Cyril and St. Methodius University
of Veliko Tarnovo, BULGARIA. Email<EMAIL_ADDRESS>The research of this
author was supported, in part, by a Bulgarian NSF contract KP-06-N32/2-2019.

###### Abstract

The hull of a linear code $C$ is the intersection of $C$ with its dual code.
We present and analyze the number of linear $q$-ary codes of the same length
and dimension but with different dimensions for their hulls. We prove that for
given dimension $k$ and length $n\geq 2k$ the number of all $[n,k]_{q}$ linear
codes with hull dimension $l$ decreases as $l$ increases. We also present
classification results for binary and ternary linear codes with trivial hulls
(LCD and self-orthogonal) for some values of the length $n$ and dimension $k$,
comparing the obtained numbers with the number of all linear codes for the
given $n$ and $k$.

##  1 Introduction

Let $\mathbb{F}_{q}$ be a finite field with $q$ elements and
$\mathbb{F}_{q}^{n}$ be the $n$-dimensional vector space over
$\mathbb{F}_{q}$. Suppose that $\mathbb{F}_{q}^{n}$ is equipped with an inner
product. Any linear subspace $C$ of $\mathbb{F}_{q}^{n}$ is called a linear
$q$-ary code of length $n$, and its orthogonal complement $C^{\perp}$ (with
respect to the defined inner product) is its dual (or orthogonal) code.

The hull of a linear code $C$ was introduced by Assmus, Jr. and Key [4], and
it is the intersection of $C$ with its dual code $C^{\perp}$. The study of
hulls also has applications ranging from the construction of quantum error
correcting codes to post-quantum cryptography [16]. If the hull contains only
the zero vector, $C$ is called a linear complementary dual (LCD) code. LCD
codes were introduced by Massey [17] and gave an optimum linear coding
solution for the two user binary adder channel. The determination of the
largest minimum weight among all LCD $[n,k]$ codes and the classification of
all or only optimal LCD $[n,k]$ codes are fundamental problems in coding
theory. Recently, much work has been done concerning this fundamental problem
(see [7, 10, 11, 14]). It has been shown in [7] that if the considered inner
product is Euclidean and $q>3$, any linear code over $\mathbb{F}_{q}$ is
equivalent to an (Euclidean) LCD code, and if the inner product is Hermitian
and $q>4$, any $q$-ary linear code is equivalent to a Hermitian LCD code. Note
that a Hermitian inner product is considered only if $q$ is an even power of
the characteristic of the field.

The other extremal case is when the hull of a linear code $C$ coincides with
the code itself. Such a code is called self-orthogonal (SO) and it is a subset
of its dual code. If $C$ coincides with $C^{\perp}$, the code is self-dual.
Self-dual codes have been a topic of great interest since the 1960s, with
connections and applications to various branches of mathematics and coding
theory. Theoretically, self-dual codes are closely connected to other
mathematical structures such as designs, lattices, graphs and modular forms
[9]. Practically, many of the best-known codes are actually self-dual codes.

Self-dual and LCD codes are widely studied in the literature not only for
their mathematical attraction, but also for their applications. In this study,
we are also interested in linear codes that lie between these two extremal
cases. If $h(C)$ is the dimension of the hull of the linear code $C$, then

$0\leq h(C)\leq\min\\{\dim C,\dim C^{\perp}\\}.$

Sendrier in [20] gave a formula for the number of linear $q$-ary codes of a
given length $n$, dimension $k$ and hull dimension $l$. He proved that,
asymptotically, the proportion of $q$-ary linear codes whose hull has
dimension $l$ to all $q$-ary linear codes of the same length and dimension,
only depends on $l$ and $q$. It turns out that for a fixed length and
dimension and for LCD codes ($l=0$), this proportion is approximately $1-1/q$.
Hence most linear codes have hull dimension 0 (i.e. LCD) when $q$ is large.

By $A_{n,k,l}$ we denote the number of the linear $[n,k]_{q}$ codes whose hull
has dimension $l$ (for fixed $q$). The main result of the present study is
contained in the following theorem.

######  Theorem 1.

For fixed positive integers $k$ and $n\geq 2k$, the following inequalities
hold:

$A_{n,k,0}>A_{n,k,1}>\cdots>A_{n,k,k}.$ (1)

Note that $A_{n,k,0}$ and $A_{n,k,k}$ are the numbers of all LCD and SO
$q$-ary codes of length $n$ and dimension $k$, respectively. To prove the main
theorem, we use some known formulae which we present in Section 2. Section 3
is devoted to Theorem 1 and its proof. Note that $A_{n,k,l}$ is the number of
all different codes with the corresponding parameters, but the number of the
inequivalent among them is much smaller. Since we do not have a formula for
this number, we did some experiments and classified the binary and ternary
linear $[n,k]$ codes for some values of $n$ and $k$. We present and analyze
the obtained results in Section 4.

##  2 Preliminaries

In this section, we present the basic definitions, theorems and formulas that
we need for the proof of the main theorem, as well as for the computational
results in Section 4.

This study only applies to linear codes over a finite field. A linear
$[n,k,d]$ $q$-ary code $C$ is a $k$-dimensional subspace of the vector space
$\mathbb{F}_{q}^{n}$, and $d$ is the smallest weight among all non-zero
codewords of $C$, called the minimum weight (or minimum distance) of the code.
A code is called even if all its codewords have even weights. Let
$(u,v):\mathbb{F}_{q}^{n}\times\mathbb{F}_{q}^{n}\to\mathbb{F}_{q}$ be an
inner product in $\mathbb{F}_{q}^{n}$. The orthogonal complement of $C$
according to the defined inner product is called the dual code of $C$ and
denoted by $C^{\perp}$, i.e. $C^{\perp}=\\{u\in\mathbb{F}_{q}^{n}:(u,v)=0$ for
all $v\in C\\}$. Obviously, $C^{\perp}$ is a linear $[n,n-k]$ code. The
intersection $C\cap C^{\perp}$ is called the hull of the code and denoted by
$\mathcal{H}(C)$. As we mentioned in Introduction, the dimension $h(C)$ of the
hull can be at least 0 and at most $k$, as $h(C)=k$ if and only if the code is
self-orthogonal, and $h(C)=0$ if and only if $C$ is an LCD code.

Two linear $q$-ary codes $C_{1}$ and $C_{2}$ are equivalent if the codewords
of $C_{2}$ can be obtained from the codewords of $C_{1}$ via a sequence of
transformations of the following types:

  * (1)

permutation on the set of coordinate positions;

  * (2)

multiplication of the elements in a given position by a non-zero element of
$\mathbb{F}_{q}$;

  * (3)

application of a field automorphism to the elements in all coordinate
positions.

An automorphism of a linear code $C$ is a sequence of the transformations
$(1)-(3)$ which maps each codeword of $C$ onto a codeword of the same code.
The set of all automorphisms of $C$ forms a group called the automorphism
group $Aut({C})$ of the code. The presented formulas do not count the
equivalence but in Section 4 we enumerate only inequivalent codes.

If the considered inner product is Euclidean, i.e.
$(x,y)=\sum_{i=1}^{n}x_{i}y_{i}$ for any two vectors $x=(x_{1},\ldots,x_{n})$,
$y=(y_{1},\ldots,y_{n})$ from $\mathbb{F}_{q}^{n}$, the dimension of the hull
is an invariant of equivalent codes in the cases $q=2$ and $q=3$, but if
$q\geq 4$, any linear code over $\mathbb{F}_{q}$ is equivalent to an Euclidean
LCD code [7]. If $q=p^{2s}$, where $p$ is the characteristic of the field, we
can consider the Hermitian inner product over $\mathbb{F}_{q}$ defined by

$(x,y)=\sum_{i=1}^{n}x_{i}y_{i}^{\sqrt{q}},\ \forall
x=(x_{1},\ldots,x_{n}),y=(y_{1},\ldots,y_{n})\in\mathbb{F}_{q}^{n}.$

For this inner product, the dimension of the hull is an invariant in the case
of the quaternary codes ($q=4$), and if $q>4$, any linear code over
$\mathbb{F}_{q}$, $q=p^{2s}$, is equivalent to an Hermitian LCD code [7].

Formulae for counting the number $\sigma(n,k)$ of all different $q$-ary self-
orthogonal codes of length $n$ and dimension $k$ were proven in [19].

######  Theorem 2.

Let $m=\lfloor n/2\rfloor$ and
$\pi_{n,k}=\prod_{i=1}^{k}\frac{q^{2m-2i+2}-1}{q^{i}-1}$. For all $1\leq k\leq
m$, we have

$\sigma_{n,k}=\left\\{\begin{array}[]{rl}\pi_{n,k},&\mbox{if}\ n\ \text{is
odd},\\\ \frac{q^{n-k}-1}{q^{n}-1}\pi_{n,k},&\mbox{if}\ n\ \mbox{and}\ q\
\text{are even},\\\ \frac{q^{m-k}-1}{q^{m}-1}\pi_{n,k},&\mbox{if}\ n\equiv
2\pmod{4}\ \mbox{and}\ q\equiv 3\pmod{4},\\\
\frac{q^{m-k}+1}{q^{m}+1}\pi_{n,k},&\mbox{if}\ (n\equiv 0\pmod{4})\\\ &\
\mbox{or}\ (n\equiv 2\pmod{4}\ \mbox{and}\ q\equiv
1\pmod{4}).\end{array}\right.$

Note that $\sigma_{n,0}=1$. We need also the following corollary.

######  Corollary 3.

Let $m=\lfloor n/2\rfloor$ and $0\leq k\leq m-1$. Then

$\frac{\sigma_{n,k}}{\sigma_{n,k+1}}=\tau_{n,k}\frac{q^{k+1}-1}{q^{2m-2k}-1},$
(2)

where

$\tau_{n,k}=\left\\{\begin{array}[]{ll}1,&\mbox{if}\ n\ \text{is odd},\\\
\frac{q^{n-k}-1}{q^{n-k-1}-1},&\mbox{if}\ n\ \mbox{and}\ q\ \text{are
even},\\\ \frac{q^{m-k}+\epsilon}{q^{m-k-1}+\epsilon},&\mbox{if}\ n\ \text{is
even and}\ q\ \text{is odd},\end{array}\right.$ (3)

with $\epsilon=-1$ if $n\equiv 2\pmod{4}$ and $q\equiv 3\pmod{4}$, and
$\epsilon=1$ otherwise.

###### Proof.

The formula (2) follows directly from Theorem 2. $\Box$

The following theorem is the main tool in our investigations.

######  Theorem 4.

[20] Let $k$ and $n\geq 2k$ be positive integers. If $\sigma_{n,i}$ is the
number of all self-orthogonal $[n,i]_{q}$ codes, $i\leq k$, then the number of
all $[n,k]_{q}$ codes whose hull has dimension $l$, $l\leq k$, is equal to

$A_{n,k,l}=\sum_{i=l}^{k}{n-2i\brack k-i}{i\brack l}(-1)^{i-l}q^{i-l\choose
2}\sigma_{n,i}.$ (4)

In the above formula, ${n\brack k}$ is the Gaussian binomial coefficient,
defined by

${n\brack 0}=1,\ {n\brack
k}=\frac{(q^{n}-1)(q^{n-1}-1)\ldots(q^{n-k+1}-1)}{(q^{k}-1)(q^{k-1}-1)\ldots(q-1)},\mbox{if}\
1\leq k\leq n.$ (5)

We use the following properties of the Gaussian coefficients for $0\leq k\leq
n-1$:

$\displaystyle{n\brack n-k}$ $\displaystyle=$ $\displaystyle{n\brack k},$ (6)
$\displaystyle{n+1\brack k}$ $\displaystyle=$
$\displaystyle\frac{q^{n+1}-1}{q^{n-k+1}-1}{n\brack k},$ (7)
$\displaystyle{n\brack k+1}$ $\displaystyle=$
$\displaystyle\frac{q^{n-k}-1}{q^{k+1}-1}{n\brack k},$ (8)
$\displaystyle{n-2i\brack k-i}$ $\displaystyle=$
$\displaystyle\frac{(q^{n-2i}-1)(q^{n-2i-1}-1)}{(q^{n-k-i}-1)(q^{k-i}-1)}{n-2i-2\brack
k-i-1}.$ (9)

######  Remark 1.

In [8], the authors give more explicit formulas for the number of different
LCD codes with a given length and dimension in the cases when $q=2$ and
$q=p^{s}$ for an odd prime $p$.

##  3 Number of linear codes whose hull has dimension $l$ for different values
of $l$

For the calculations, we need to evaluate a few constants.

######  Lemma 5.

If $n=2m$ and $k\leq m$ are positive integers, then

$\tau_{n,k}\geq\left\\{\begin{array}[]{cl}q,&\mbox{if}\ q\ \text{is even, or}\
(n\equiv 2\pmod{4}\ \mbox{and}\ q\equiv 3\pmod{4})\\\
(q+1)/2,&\mbox{otherwise}.\end{array}\right.$

###### Proof.

We use the formula (3):

$\frac{q^{n-k}-1}{q^{n-k-1}-1}=q+\frac{q-1}{q^{n-k-1}-1}>q,\ \
\frac{q^{m-k}-1}{q^{m-k-1}-1}=q+\frac{q-1}{q^{m-k-1}-1}>q,$
$\frac{q^{m-i}+1}{q^{m-i-1}+1}=q-\frac{q-1}{q^{m-i-1}+1}\geq
q-\frac{q-1}{2}=\frac{q+1}{2}.$

This proves the lemma in all cases. $\Box$

######  Lemma 6.

Let $m=\lfloor n/2\rfloor$, $0\leq i<k\leq m$, and

$\lambda_{n,i}=\frac{(q^{n-2i}-1)(q^{n-2i-1}-1)}{(q^{2m-2i}-1)(q^{n-k-i}-1)(q^{k-i}-1)}=\left\\{\begin{array}[]{ll}\frac{q^{n-2i}-1}{(q^{n-k-i}-1)(q^{k-i}-1)}&\mbox{if}\
n=2m+1,\\\ \frac{q^{n-2i-1}-1}{(q^{n-k-i}-1)(q^{k-i}-1)}&\mbox{if}\ n=2m,\\\
\end{array}\right.$

Then $q^{2m+1-n}\lambda_{n,i}>1+\frac{1}{q^{k-i}-1}>1$.

###### Proof.

We have that

$\frac{q^{n-2i}-1}{(q^{n-k-i}-1)(q^{k-i}-1)}=1+\frac{1}{q^{n-k-i}-1}+\frac{1}{q^{k-i}-1}>1$

and

$\displaystyle\frac{q^{n-2i}-q}{(q^{n-k-i}-1)(q^{k-i}-1)}$
$\displaystyle=1+\frac{1}{q^{n-k-i}-1}+\frac{1}{q^{k-i}-1}-\frac{q-1}{(q^{n-k-i}-1)(q^{k-i}-1)}$
$\displaystyle\geq
1+\frac{1}{q^{n-k-i}-1}+\frac{1}{q^{k-i}-1}-\frac{q^{k-i}-1}{(q^{n-k-i}-1)(q^{k-i}-1)}$
$\displaystyle=1+\frac{1}{q^{k-i}-1}>1,$

therefore $\lambda_{n,i}>1$ if $n=2m+1$, and $q\lambda_{n,i}>1$ if $n=2m$.
$\Box$

We use the constant $\lambda_{n,i}$ in the following formula

${n-2i\brack k-i}=(q^{2m-2i}-1)\lambda_{n,i}{n-2i-2\brack k-i-1}.$ (10)

Consider the difference $A_{n,k,l}-A_{n,k,l+1}$ for $0\leq l\leq k-1$.

$\displaystyle A_{n,k,l}-A_{n,k,l+1}=$ $\displaystyle=$
$\displaystyle\sum_{i=l}^{k}{n-2i\brack k-i}{i\brack l}(-1)^{i-l}q^{i-l\choose
2}\sigma_{n,i}-\sum_{i=l+1}^{k}{n-2i\brack k-i}{i\brack
l+1}(-1)^{i-l-1}q^{i-l-1\choose 2}\sigma_{n,i}$ $\displaystyle=$
$\displaystyle{n-2l\brack
k-l}\sigma_{n,l}+\sum_{i=l+1}^{k}(-1)^{i-l}{n-2i\brack k-i}\left({i\brack
l}q^{i-l\choose 2}+{i\brack l+1}q^{i-l-1\choose 2}\right)\sigma_{n,i}$
$\displaystyle=$ $\displaystyle{n-2l\brack
k-l}\sigma_{n,l}+\sum_{i=l+1}^{k}(-1)^{i-l}a_{n,k,l,i},$

where

$\displaystyle a_{n,k,l,i}$ $\displaystyle={n-2i\brack k-i}\left({i\brack
l}q^{i-l\choose 2}+{i\brack l+1}q^{i-l-1\choose 2}\right)\sigma_{n,i}$
$\displaystyle={n-2i\brack k-i}{i\brack
l}(q^{i-l-1}+\frac{q^{i-l}-1}{q^{l+1}-1})q^{i-l-1\choose 2}\sigma_{n,i}$
$\displaystyle={n-2i\brack k-i}{i\brack
l}\frac{q^{i}-q^{i-l-1}+q^{i-l}-1}{q^{l+1}-1}q^{i-l-1\choose 2}\sigma_{n,i},$
$\displaystyle={n-2i\brack k-i}{i\brack
l}\frac{\alpha_{l,i}}{q^{l+1}-1}q^{i-l-1\choose 2}\sigma_{n,i},\ \ l+1\leq
i\leq k,$ $\alpha_{l,i}=q^{i-l-1}(q^{l+1}+q-1)-1\geq q^{l+1}+q-2\geq 2q-2.$

We take $a_{n,k,l,l}={n-2l\brack k-l}\sigma_{n,l}$. Obviously, the
coefficients $a_{n,k,l,i}$ are positive for all possible values of $n$, $k$,
$l$ and $i$. Firstly, we prove two lemmas about these coefficients.

######  Lemma 7.

If $0\leq l\leq k-1\leq m-1$, $m=\lfloor n/2\rfloor$, and $l+1\leq i\leq k-1$,
then $a_{n,k,l,i}>a_{n,k,l,i+1}$.

###### Proof.

We evaluate

$\displaystyle\frac{a_{n,k,l,i}}{a_{n,k,l,i+1}}$
$\displaystyle=\frac{{n-2i\brack k-i}{i\brack l}\alpha_{l,i}q^{i-l-1\choose
2}\sigma_{n,i}}{{n-2i-2\brack k-i-1}{i+1\brack l}\alpha_{l,i+1}q^{i-l\choose
2}\sigma_{n,i+1}}$
$\displaystyle=\frac{(q^{2m-2i}-1)\lambda_{n,i}(q^{i-l+1}-1)\alpha_{l,i}(q^{i+1}-1)\tau_{n,i}}{(q^{i+1}-1)\alpha_{l,i+1}q^{i-l-1}(q^{2m-2i}-1)}$
$\displaystyle=\frac{\lambda_{n,i}(q^{i-l+1}-1)\alpha_{l,i}\tau_{n,i}}{\alpha_{l,i+1}q^{i-l-1}}=\lambda_{n,i}\tau_{n,i}(q^{2}-\frac{1}{q^{i-l-1}})\frac{\alpha_{l,i}}{\alpha_{l,i+1}}.$

Using that

$\alpha_{l,i+1}=q^{i-l}(q^{l+1}+q-1)-1=q\alpha_{l,i}+q-1,$ (11)

we obtain

$\frac{\alpha_{l,i}}{\alpha_{l,i+1}}=\frac{\alpha_{l,i}}{q\alpha_{l,i}+q-1}=\frac{1}{q+\frac{q-1}{\alpha_{l,i}}}\geq\frac{1}{q+\frac{1}{2}}=\frac{2}{2q+1},$
(12)

and therefore

$\frac{a_{n,k,l,i}}{a_{n,k,l,i+1}}\geq\frac{2\lambda_{n,i}\tau_{n,i}}{2q+1}(q^{2}-\frac{1}{q^{i-l-1}})\geq\frac{2\lambda_{n,i}\tau_{n,i}(q^{2}-1)}{2q+1}.$

According to Lemma 6 and Lemma 5, we have that $\lambda_{n,i}\tau_{n,i}\geq 1$
in the cases when $n$ is odd, or $n$ and $q$ are even, or $n\equiv 2\pmod{4}$
and $q\equiv 3\pmod{4}$. Hence in these cases

$\frac{a_{n,k,l,i}}{a_{n,k,l,i+1}}\geq\frac{2(q^{2}-1)}{2q+1}>1.$

In the case when $n=2m$ is even and $q$ is odd ($q\geq 3$), such that (1)
$n\equiv 2\pmod{4}$ and $q\equiv 1\pmod{4}$, or (2) $n\equiv 0\pmod{4}$,
$\tau_{n,i}\geq(q+1)/2$, so

$\frac{a_{n,k,l,i}}{a_{n,k,l,i+1}}\geq\frac{2(q^{2}-1)(q+1)}{2q(2q+1)}>\frac{(q^{2}-1)}{2q}\geq\frac{8}{6}>1.$

Hence $a_{n,k,l,i}>a_{n,k,l,i+1}$ for $l+1\leq i\leq k-1$ in all cases. $\Box$

######  Lemma 8.

If $q\geq 3$ and $0\leq l\leq k-1$, then $a_{n,k,l,l}>a_{n,k,l,l+1}$.

###### Proof.

Now

$\displaystyle\frac{a_{n,k,l,l}}{a_{n,k,l,l+1}}$
$\displaystyle=\frac{{n-2l\brack k-l}(q-1)\sigma_{n,l}}{{n-2l-2\brack
k-l-1}\alpha_{l,l+1}\sigma_{n,l+1}}$
$\displaystyle=\frac{(q^{2m-2l}-1)\lambda_{n,l}(q-1)\tau_{n,l}(q^{l+1}-1)}{(q^{l+1}+q-2)(q^{2m-2l}-1)}$
$\displaystyle=\frac{\lambda_{n,l}\tau_{n,l}(q-1)(q^{l+1}-1)}{(q^{l+1}+q-2)}.$

As in the previous proof,

$\frac{a_{n,k,l,l}}{a_{n,k,l,l+1}}\geq\frac{(q-1)(q^{l+1}-1)}{q^{l+1}+q-2}=\frac{1}{\frac{1}{q-1}+\frac{1}{q^{l+1}-1}}\geq\frac{q-1}{2}\geq
1$

in the cases (1) $n=2m+1$, (2) $n=2m$ and $q=2^{s}$, and (3) $n=2m\equiv
2\pmod{4}$ and $q\equiv 3\pmod{4}$.

Let now $n\equiv 0\pmod{4}$ and $q\equiv 1\pmod{2}$, or $n\equiv 2\pmod{4}$
and $q\equiv 1\pmod{4}$. In this case

$\frac{a_{n,k,l,l}}{a_{n,k,l,l+1}}\geq\frac{(q+1)(q-1)(q^{l+1}-1)}{2(q^{l+1}+q-2)}\geq\frac{(q+1)(q-1)}{4}\geq\frac{q^{2}-1}{4}\geq
2.$

This gives us that $a_{n,k,l,l}>a_{n,k,l,l+1}$ if $q\geq 3$. $\Box$

Consider the binary case separately, i.e. $q=2$. Now

$a_{n,k,l,i}={n-2i\brack k-i}{i\brack
l}\frac{\alpha_{l,i}}{2^{l+1}-1}2^{i-l-1\choose 2}\sigma_{n,i},\ \ l+1\leq
i\leq k,\ \ a_{n,k,l,l}={n-2l\brack k-l}\sigma_{n,l},$

where $\alpha_{l,i}=2^{i-l-1}(2^{l+1}+1)-1$ for $l+1\leq i\leq k$. Using that

$\alpha_{l,l+1}=2^{l+1},\ \ \alpha_{l,l+2}=2^{l+2}+1,\ \
\alpha_{l,l+3}=2^{l+3}+3,$

we obtain

$a_{n,k,l,l+1}={n-2l-2\brack k-l-1}{l+1\brack
l}\frac{\alpha_{l,l+1}}{2^{l+1}-1}\sigma_{n,l+1}=2^{l+1}{n-2l-2\brack
k-l-1}\sigma_{n,l+1},$ $a_{n,k,l,l+2}={n-2l-4\brack k-l-2}{l+2\brack
l}\frac{\alpha_{l,l+2}}{2^{l+1}-1}\sigma_{n,l+2}=\frac{(2^{l+2}-1)(2^{l+2}+1)}{3}{n-2l-4\brack
k-l-2}\sigma_{n,l+2},$ $\displaystyle a_{n,k,l,l+3}$
$\displaystyle={n-2l-6\brack k-l-3}{l+3\brack
l}\frac{\alpha_{l,l+3}}{2^{l+1}-1}2^{2\choose 2}\sigma_{n,l+3}$
$\displaystyle={n-2l-6\brack
k-l-3}\frac{(2^{l+3}-1)(2^{l+2}-1)(2^{l+1}-1)(2^{l+3}+3)}{21(2^{l+1}-1)}2\sigma_{n,l+3}$
$\displaystyle=\frac{(2^{l+3}-1)(2^{l+2}-1)(2^{l+4}+6)}{21}{n-2l-6\brack
k-l-3}\sigma_{n,l+3}.$

Next, we consider separately the cases $k\leq 2l+1$ and $k>2l+1$.

######  Lemma 9.

If $q=2$ and $k\leq 2l+1$ then $a_{n,k,l,l}>a_{n,k,l,l+1}$.

###### Proof.

In this case

$a_{n,k,l,l}-a_{n,k,l,l+1}={n-2l\brack k-l}\sigma_{n,l}-2^{l+1}{n-2l-2\brack
k-l-1}\sigma_{n,l+1}$
$=[\lambda_{n,l}\tau_{n,l}(2^{l+1}-1)-2^{l+1}]{n-2l-2\brack
k-l-1}\sigma_{n,l+1}.$

According to Lemma 6,
$2^{2m+1-n}\lambda_{n,l}>1+\frac{1}{2^{k-l}-1}=\frac{2^{k-l}}{2^{k-l}-1}$. If
$n$ is odd then $\tau_{n,l}=1$ and therefore
$\lambda_{n,l}\tau_{n,l}>\frac{2^{k-l}}{2^{k-l}-1}$. We have the same
inequality in the opposite case, when $n$ is even, because then

$\lambda_{n,l}\tau_{n,l}>\frac{2^{k-l}(2^{n-l}-1)}{2(2^{k-l}-1)(2^{n-l-1}-1)}=\frac{2^{k-l}(2^{n-l}-1)}{(2^{k-l}-1)(2^{n-l}-2)}>\frac{2^{k-l}}{2^{k-l}-1}.$

Hence

$\displaystyle a_{n,k,l,l}-a_{n,k,l,l+1}$
$\displaystyle>[\frac{2^{k-l}(2^{l+1}-1)}{2^{k-l}-1}-2^{l+1}]{n-2l-2\brack
k-l-1}\sigma_{n,l+1}$
$\displaystyle=\frac{2^{l+1}-2^{k-l}}{2^{k-l}-1}{n-2l-2\brack
k-l-1}\sigma_{n,l+1}$
$\displaystyle=\frac{2^{2l+1}-2^{k}}{2^{k}-2^{l}}{n-2l-2\brack
k-l-1}\sigma_{n,l+1}\geq 0$

when $k\leq 2l+1$. $\Box$

We are now ready to prove the main theorem.

Proof of Theorem 1:

We use that

$A_{n,k,l}-A_{n,k,l+1}=\sum_{i=l}^{k}(-1)^{i-l}a_{n,k,l,i}=\sum_{j=0}^{k-l}(-1)^{j}a_{n,k,l,j+l}.$
(13)

We proved in Lemma 7, that $a_{n,k,l,i}>a_{n,k,l,i+1}$ if $l+1\leq i\leq k-1$
for all possible values of $l$, $k$, $n$ and $q$. From Lemma 8 we have that
$a_{n,k,l,l}>a_{n,k,l,l+1}$, if $q\geq 3$, and the same inequality holds in
the binary case if $k\leq 2l+1$. In these cases

  * •

if $k-l=2r+1$ then

$\sum_{j=0}^{k-l}(-1)^{j}a_{n,k,l,j+l}=\sum_{s=0}^{r}(a_{n,k,l,2s}-a_{n,k,l,2s+1})>0,$

  * •

and if $k-l=2r$ is even, we have

$A_{n,k,l}-A_{n,k,l+1}=\sum_{s=0}^{r-1}(a_{n,k,l,2s}-a_{n,k,l,2s+1})+a_{n,k,l,2r}>0.$

Consider now the case $q=2$ and $k>2l+1$. Then $k\geq l+3$ except when $l=0$
and $k=2$, but we will leave this special case for last. Let
$a=a_{n,k,l,l}-a_{n,k,l,l+1}+a_{n,k,l,l+2}-a_{n,k,l,l+3}$,
$k\geq\max\\{2l+2,l+3\\}$. We have

$a_{n,k,l,l+2}-a_{n,k,l,l+3}=\frac{(2^{l+2}-1)(2^{l+2}+1)}{3}{n-2l-4\brack
k-l-2}\sigma_{n,l+2}$
$-\frac{(2^{l+3}-1)(2^{l+2}-1)(2^{l+4}+6)}{21}{n-2l-6\brack
k-l-3}\sigma_{n,l+3}$
$=[\frac{(2^{2m-2l-4}-1)\lambda_{n,l+2}(2^{l+2}-1)(2^{l+2}+1)\tau_{n,l+2}(2^{l+3}-1)}{3(2^{2m-2l-4}-1)}$
$-\frac{(2^{l+3}-1)(2^{l+2}-1)(2^{l+4}+6)}{21}]{n-2l-6\brack
k-l-3}\sigma_{n,l+3}$
$=[\frac{\lambda_{n,l+2}(2^{l+2}+1)\tau_{n,l+2}}{3}-\frac{(2^{l+4}+6)}{21}](2^{l+2}-1)(2^{l+3}-1){n-2l-6\brack
k-l-3}\sigma_{n,l+3}$ $>(7(2^{l+2}+1)-2^{l+4}-6){n-2l-6\brack
k-l-3}\sigma_{n,l+3}$

For the last inequality we use that
$\lambda_{n,l+2}\tau_{n,l+2}>\frac{2^{k-l-2}}{2^{k-l-2}-1}>1$, and
$(2^{l+2}-1)(2^{l+3}-1)\geq 21$. It follows that

$a_{n,k,l,l+2}-a_{n,k,l,l+3}>(3.2^{l+2}+1){n-2l-6\brack k-l-3}\sigma_{n,l+3}.$

Let $\mu=\min\\{{n-2l-6\brack k-l-3}\sigma_{n,l+3},{n-2l-2\brack
k-l-1}\sigma_{n,l+1}\\}>0$. Then

$a>\left(3.2^{l+2}+1-\frac{2^{k}-2^{2l+1}}{2^{k}-2^{l}}\right)\mu=\left(3.2^{l+2}+\frac{2^{2l+1}-2^{l}}{2^{k}-2^{l}}\right)\mu>0.$

To prove that $A_{n,k,l}>A{n,k,l+1}$, we consider the same two cases depending
on the parity of $k-l$.

  * •

If $k-l=2r+1$ then

$\sum_{j=0}^{k-l}(-1)^{j}a_{n,k,l,j+l}=a+\sum_{s=2}^{r}(a_{n,k,l,2s}-a_{n,k,l,2s+1})>0.$

  * •

In the case when $k-l=2r$ is even, we have

$A_{n,k,l}-A_{n,k,l+1}=a+\sum_{s=2}^{r-1}(a_{n,k,l,2s}-a_{n,k,l,2s+1})+a_{n,k,l,2r}>0.$

In the last remaining case we have $k=2$ and $l=0$. Then

$\displaystyle A_{n,2,0}-A_{n,2,1}$ $\displaystyle={n\brack
2}\sigma_{n,0}-2{n-2\brack 1}\sigma_{n,1}+5{n-4\brack 0}\sigma_{n,2}$
$\displaystyle=\frac{(2^{n}-1)(2^{n-1}-1)}{3}-2(2^{n-2}-1)\sigma_{n,1}+5\sigma_{n,2}.$

After the calculations, we obtain

$A_{n,2,0}-A_{n,2,1}=\left\\{\begin{array}[]{cl}2^{n-3}(2^{n-1}-1)/3&\mbox{if}\
n=2m+1\\\ 2^{n-1}(2^{n-3}+1)/3&\mbox{if}\ n=2m.\end{array}\right.$

Hence $A_{n,k,l}>A_{n,k,l+1}$ over $\mathbb{F}_{q}$ in all cases. $\square$

##  4 Computational results

In this section, we present computational results with the number of
inequivalent linear $[n,k,\geq 2]$ codes of different types for given length
and dimension over $\mathbb{F}_{2}$ and $\mathbb{F}_{3}$, whose dual distance
is at least 2 (this means that their generator matrices have no zero columns).
Denote by $\mathcal{B}_{q}(n,k,l)$ the set of all inequivalent $q$-ary linear
codes whose hull has dimension $l$. In the binary case, the mass formula, that
helps to verify classification results for the linear $[n,k]_{2}$ codes $C$
with $h(C)=l$, is the following

$A_{n,k,l}=\sum_{C\in\mathcal{B}_{2}(n,k,l)}\frac{n!}{\mathrm{Aut}(C)}$ (14)

In the ternary case we have

$A_{n,k,l}=\sum_{C\in\mathcal{B}_{3}(n,k,l)}\frac{2^{n}n!}{\mathrm{Aut}(C)}.$
(15)

More information on these formulas for self-dual codes can be found in [15].
Araya and Harada in [1] used them in the classification of the binary LCD
codes of length $n\leq 13$ and ternary LCD codes of length $n\leq 10$. They
described clearly how to use these mass formulas for the verification of the
computational results for classification of the binary LCD $[6,3]$ codes with
all possible minimum distances $d\geq 1$ and dual distances $d^{\perp}\geq 1$.

We classify the binary and ternary linear, self-orthogonal and LCD codes of a
given dimension $3\leq k\leq 10$, length $k+3\leq n\leq 20$, minimum distance
$d\geq 2$ and dual distance $d^{\perp}\geq 2$. In the binary case, we classify
also the even codes with the corresponding length and dimension, i.e. all
linear codes whose codewords have only even weights. We do not count codes
with dual distance $d^{\perp}=1$ because if $C$ is an $[n,k,d]$ code with
$d^{\perp}=1$, then all its codewords share a common zero coordinate, so
$C=(0|C_{1})$ where $C_{1}$ is an $[n-1,k,d]$ code. In this case,
$C^{\perp}=(0|C_{1}^{\perp})\cup(1|C_{1}^{\perp})$,
$\mathcal{H}(C)=(0|\mathcal{H}(C_{1})$ and $h(C)=h(C_{1})$. If $C$ is an
$[n,k,1]$ code then $C\cong(0|C_{1})\cup(1|C_{1})$, where $C_{1}$ is an
$[n-1,k-1]$ code, and then $C^{\perp}\cong(0|C_{1}^{\perp})$. This gives us
that $\mathcal{H}(C)=(0|\mathcal{H}(C_{1})$ and $h(C)=h(C_{1})$. Therefore, if
we have the number $LCD^{*}(n,k)$ of all $[n,k,\geq 2]$ LCD codes with dual
distance $\geq 2$ for all length $\leq n$ and dimensions $\leq k$, we can
easily compute the number $LCD(n,k)$ of LCD $[n,k,\geq 1]$ codes with dual
distance $\geq 1$, using the following formula

$LCD(n,k)=\sum_{m=k+1}^{n}LCD^{*}(m,k)+LCD(n,k-1).$

We give a simple example that can be followed by hand.

######  Example 1.

Let $n=4$, $k=2$ and $q=2$. In this case

$A_{4,2,2}=\sigma_{4,2}=3,\ A_{4,2,1}=12,\ A_{4,2,0}=20,\
\mathcal{B}_{2}(4,2,2)=\mathcal{B}_{2}(4,2,1)=1,\ \mathcal{B}_{2}(4,2,0)=4.$

There are six inequivalent $[4,2]$ binary codes. The first one is obtained
from $\mathbb{F}_{2}^{2}$ by adding two zero columns and it is an LCD code.
There are two more binary $[4,2]$ inequivalent codes with zero columns, and
these are the codes $C_{2}=\\{0000,0110,0001,0111\\}$, $h(C_{2})=1$, and
$C_{3}=\\{0000,0110,0101,0011\\}$, $h(C_{3})=0$. The remaining three codes are
$C_{4}=\\{0000,1110,0001,1111\\}$, $h(C_{4})=0$,
$C_{5}=\\{0000,1110,0101,1011\\}$, $h(C_{5})=0$, and the self-dual $[4,2,2]$
code $C_{6}=\\{0000,1100,0011,1111\\}$, $h(C_{6})=2$.Only two of these codes
have minimum and dual distance $d=d^{\perp}=2$ and these are $C_{5}$ and
$C_{6}$.

We use the following properties of the considered types of binary and ternary
codes:

  * •

If $C$ is an LCD code, its dual code $C^{\perp}$ is also LCD. It follows that
$LCD^{*}(n,k)=LCD^{*}(n,n-k)$ and $LCD(n,k)=LCD(n,n-k)$ both in the binary and
the ternary case. The same holds for the number of all linear codes of length
$n$ and dimension $k$.

  * •

Self-orthogonal codes exist only when $k\leq n/2$.

  * •

All binary self-orthogonal codes are even.

  * •

Ternary self-dual codes exist only for lengths a multiple of 4 and only have
codewords of Hamming weight a multiple of 3.

We obtain the classification results by the program Generation of the software
package QExtNewEdition [5]. The computations were executed on a Windows 11 OS
in a single core of an Intel Xeon Gold 5118 CPU with a 2.30 GHz clock
frequency.

In the binary case, we classify all linear, even, self-orthogonal and LCD
codes of length $n\leq 20$ and dimension $k\leq 10$. The results are presented
in Table 1. For some values of $n$ and $k$, when all inequivalent codes are
too many (more than a milion), we classify only optimal $[n,k]$ codes. In this
cases, we put a $*$ after the number of codes. Consider for example $n=17$ and
$k=7$. The largest possible minimum distance for a binary $[17,7]$ code is
$d=6$ (see [12]). There are exactly 377 binary linear $[17,7,6]$ codes.
Exactly 329 of these codes are even, 7 are LCD, but none of them is self-
orthogonal. Furthermore, there are 497119 even, 58 self-orthogonal, and 14 734
654 LCD $[17,7,\geq 2]$ binary codes with dual distance $d^{\perp}\geq 2$.
Since we have the count of both all $[17,7,\geq 2]$ and optimal $[17,7,6]$
even codes without zero columns, we denote this in the table by 497119(329*).
The optimal binary self-orthogonal codes with the parameters presented in
Table 1 have been also classified in [13].

The results in the table confirm the fact that LCD codes are much more than
self-orthogonal and even more than even codes for a given length and size,
even when considering only inequivalent codes. However, this is not the case
if we consider only the optimal codes. In most cases, when the optimal code is
unique, it is not LCD. The optimal $[19,7,8]$ is self-orthogonal, as is the
optimal $[20,8,8]$ code. The optimal $[17,8,6]$ and $[17,9,5]$ codes are LCD,
but the optimal $[19,8,7]$, $[20,9,7]$ and $[18,9,6]$ are neither self-
orthogonal, nor LCD. We see interesting examples in the optimal codes of
dimension 7. Out of all 377 optimal $[17,7,6]$ codes, none is self-orthogonal,
but 7 are LCD codes. We have the opposite situation for length 20, namely out
of all 26 optimal $[20,7,8]$ codes, none is LCD, but four are self-orthogonal.
It is also worth noting the optimal $[20,10,6]$ codes, where out of all 1682
codes only one is odd-like (it contains codewords of odd weight), but none of
the 1681 even codes is self-orthogonal.

In the ternary case, we consider linear, self-orthogonal and LCD codes. We
classify all $[n,k,\geq 2]$ linear, LCD and SO codes without zero columns for
$k=3$ and $n\leq 20$, $k=4$ and $n\leq 15$, $k=5$ and $n\leq 14$, $6\leq k\leq
9$ and $n\leq 13$. For $k=7$ and $n=15,16,17$, $k=8$ and $n=16$, 17, 18, and
$k=9$, $n=19$ we classify only the self-orthogonal codes. The question mark
(?) in the table means that there are too many corresponding codes (more than
a million). Furthermore, we classify all optimal codes of these three types
with 21 parameters. The results are presented in Table 2. Let’s describe the
codes of dimension 5 in more detail. If we look at the number of codes with
length $n$, $10\leq n\leq 14$, we see that the LCD codes are more than half of
all inequivalent $[n,5,\geq 2]$ linear codes. The situation with self-
orthogonal codes is quite different - they occur much less often. For example,
for length 13, the linear codes are more than a million, but only 17 of them
are self-orthogonal. For the optimal codes, we have: (1) of four $[15,5,8]$
linear codes one is LCD (none is self orthogonal since 8 is not a multiple of
3), (2) the only $[16,5,9]$ code is self-orthogonal, (3) there are 1804 linear
$[17,5,9]$ codes, 35 of which are self-orthogonal and 400 are LCD, (4) none of
the seven $[18,5,10]$ and both $[19,5,11]$ codes is LCD, nor self-orthogonal,
(5) there are two linear $[20,5,12]$ codes and both are self-orthogonal.

We have interesting results with the optimal $[19,7,9]_{3}$ and $[20,8,9]_{3}$
codes. In these cases, all optimal linear codes are self-orthogonal.

Table 1: Classification of binary linear codes $(d^{\perp}\geq 2)$

$k=3$
---
$n=$ | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20
linear | 8 | 15 | 27 | 45 | 71 | 107 | 159 | 226 | 317 | 435 | 587 | 779 | 1024 | 1325 | 1699
even | 3 | 4 | 8 | 9 | 17 | 20 | 34 | 39 | 61 | 72 | 106 | 123 | 174 | 204 | 277
SO | 1 | 1 | 3 | 1 | 6 | 2 | 12 | 4 | 21 | 7 | 34 | 11 | 54 | 19 | 82
LCD | 2 | 5 | 7 | 17 | 20 | 42 | 47 | 91 | 98 | 180 | 189 | 328 | 340 | 565 | 580

$k=4$
---
$n=$ | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20
linear | 15 | 42 | 100 | 222 | 462 | 928 | 1782 | 3333 | 6058 | 10759 | 18694 | 31877 | 53357 | 87864
even | 4 | 10 | 18 | 37 | 63 | 122 | 202 | 366 | 602 | 1038 | 1671 | 2785 | 4411 | 7122
SO | - | 2 | 1 | 6 | 3 | 16 | 8 | 39 | 23 | 92 | 55 | 199 | 131 | 424
LCD | 5 | 16 | 30 | 82 | 139 | 345 | 568 | 1267 | 2040 | 4193 | 6631 | 12720 | 19734 | 35732

$k=5$
---
$n=$ | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20
linear | 27 | 100 | 331 | 1007 | 2936 | 8208 | 22326 | 59235 | 153711 | 390607 | 972726 | 2373644 | 5676542
even | 7 | 16 | 46 | 102 | 264 | 593 | 1448 | 3319 | 7886 | 18096 | 42193 | 96243 | 219712
SO | - | - | 2 | 2 | 11 | 8 | 38 | 33 | 134 | 123 | 442 | 462 | 1450
LCD | 7 | 30 | 84 | 297 | 816 | 2596 | 6908 | 20238 | 52248 | 142468 | 355083 | 908879 | 2177772

$k=6$
---
$n=$ | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | $[18,6,8]$ | $[19,6,8]$ | $[20,6,8]$
linear | 45 | 222 | 1007 | 4393 | 18621 | 78148 | 325815 | 1350439 | 5548052 | 2* | 28* | 1833*
even | 9 | 30 | 92 | 303 | 945 | 3166 | 10576 | 37017 | 131233 | 2* | 21* | 1418*
SO | - | - | - | 3 | 3 | 21 | 21 | 105 | 123 | 521(2*) | 746(2*) | 2758(23*)
LCD | 17 | 82 | 297 | 1418 | 5632 | 25954 | 108846 | 484648 | 2034711 | 8633817(0*) | 2* | 392*

$k=7$
---
$n=$ | 10 | 11 | 12 | 13 | 14 | 15 | 16 | $[17,7,6]$ | $[18,7,7]$ |  $[19,7,8]$) | $[20,7,8]$
linear | 71 | 462 | 2936 | 18621 | 121169 | 814087 | 5635181 | 377* | 2* | 1* | 26*
even | 13 | 46 | 194 | 774 | 3518 | 16714 | 87998 | 497119(329*) | 3010238(0*) | 1* | 21*
SO | - | - | - | - | 4 | 6 | 41 | 58(0*) | 300(0*) | 540(1*) | 2469(4*)
LCD | 20 | 139 | 816 | 5632 | 37166 | 272131 | 1968462 | 14734654(7*) | 0* | 0* | 0*

$k=8$
---
$n=$ | 11 | 12 | 13 | 14 | 15 | 16 | $[17,8,6]$ | $[18,8,6]$ | $[19,8,7]$ | $[20,8,8]$
linear | 107 | 928 | 8208 | 78148 | 814087 | 9273075 | 1* | 918* | 1* | 1*
even | 16 | 76 | 362 | 2020 | 12646 | 94136 | 818890(1*) | 907* | 0* | 1*
SO | - | - | - | - | - | 7 | 10(0*) | 86(0*) | 168(0*) | 1016(1*)
LCD | 42 | 345 | 2596 | 25954 | 272131 | 3315862 | 1* | 337* | 0* | 0*

$k=9$
---
$n=$ | 12 | 13 | 14 | 15 | 16 | $[17,9,5]$ | $[18,9,6]$ | $[19,9,6]$ | $[20,9,7]$
linear | 159 | 1782 | 22326 | 325815 | 5635181 | 1* | 1* | 1700* | 1*
even | 22 | 109 | 689 | 4973 | 46344 | 554238(0*) | 8547530(1*) | 1694* | 0*
SO | - | - | - | - | - | - | 9(0*) | 22(0*) | 194(0*)
LCD | 47 | 568 | 6908 | 108846 | 1968462 | 1* | 0* | 3* | 0*

$k=10$
---
$n=$ | 13 | 14 | 15 | 16 | $[17,10,4]$ | $[18,10,4]$ | $[19,10,5]$ | $[20,10,6]$
linear | 226 | 3333 | 59235 | 1350439 | 14390* | 11581361* | 31237* | 1682*
even | 26 | 165 | 1230 | 12257 | 169691(2614*) | 3433243(263147*) | 0* | 1681*
SO | - | - | - | - | - | - | - | 16(0*)
LCD | 91 | 1267 | 20238 | 484648 | 14734654(4550*) | 4535834* | 11554* | 601*

Table 2: Classification of ternary linear codes $(d^{\perp}\geq 2)$

$k=3$
---
$n=$ | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20
linear | 14 | 31 | 68 | 137 | 263 | 484 | 878 | 1538 | 2649 | 4474 | 7421 | 12093 | 19420 | 30680 | 47793
SO | 0 | 1 | 1 | 3 | 4 | 5 | 10 | 15 | 17 | 31 | 44 | 54 | 91 | 126 | 160
LCD | 7 | 15 | 33 | 67 | 132 | 253 | 471 | 839 | 1491 | 2560 | 4294 | 7142 | 11583 | 18423 | 29070

$k=4$
---
$n=$ | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | $[16,4,9]$ | $[17,4,10]$ | $[18,4,11]$
linear | 31 | 129 | 460 | 1638 | 5701 | 19996 | 69536 | 239681 | 809694 | 317* | 18* | 2*
SO | - | 1 | 1 | 3 | 5 | 16 | 26 | 52 | 121 | 255(13*) | 523 | 1267
LCD | 15 | 33 | 220 | 839 | 3077 | 11228 | 40668 | 144447 | 497679 | 124* | 8* | 0*

$k=5$
---
$n=$ | 8 | 9 | 10 | 11 | 12 | 13 | 14 | $[15,5,8]$ | $[16,5,9]$ | $[17,5,9]$ | $[18,5,10]$
linear | 68 | 460 | 3221 | 24342 | 202064 | 1767647 | 15604611 | 4* | 1* | 1804* | 7*
SO | - | - | 0 | 3 | 7 | 17 | 44 | 156 | 523(1*) | 1981(35*) | 9460
LCD | 33 | 220 | 1681 | 13537 | 118878 | 1080479 | 9737965 | 1* | 0* | 400* | 0*

$k=6$ | $k=7$
---|---
$n=$ | 9 | 10 | 11 | 12 | 13 | 10 | 11 | 12 | 13 | 15 | 16 | 17 | 20
linear | 137 | 1638 | 24342 | 474106 | 10956955 | 263 | 5701 | 202064 | 10956955 | ? | ? | ? | ?
SO | - | - | - | 3 | 4 | - | - | - | - | 12 | 50 | 249 | 2287775
LCD | 67 | 839 | 13537 | 283650 | 6807504 | 132 | 3077 | 118878 | 6807504 | ? | ? | ? | ?

$k=8$ | $k=9$
---|---
$n=$ | 11 | 12 | 13 | 16 | 17 | 18 | 12 | 13 | 14 | 19
linear | 484 | 19996 | 1767647 | ? | ? | ? | 878 | 69536 | 15604611 | ?
SO | - | - | - | 7 | 16 | 137 | - | - | - | 56
LCD | 253 | 11228 | 1080479 | ? | ? | ? | 471 | 40668 | 11835111 | ?

Optimal codes
---
| $[19,4,12]$ | $[20,4,12]$ | $[19,5,11]$ | $[20,5,12]$ | $[14,6,6]$ | $[15,6,7]$ | $[16,6,7]$ | $[17,6,8]$
linear | 1* | 84* | 2* | 2* | 47674* | 22* | $>10^{8}*$ | 2145181*
SO | 2867(1*) | 6893(32*) | 50618 | 294990(2*) | 15(4*) | 61 | 286 | 1504
LCD | 0* | 14* | 0* | 0* | 27776* | 10* | 53236943* | 807993*

Optimal codes
---
$n=$ | $[18,6,9]$ | $[19,6,9]$ | $[18,7,8]$ | $[19,7,9]$ | $[19,8,8]$ | $[20,8,9]$ | $[20,9,8]$
linear | 171* | ? | 827459* | 61* | 1508* | 23* | 32*
SO | 13831(105*) | 184980(18019*) | 2486 | 57551(61*) | 2281 | 112899(23*) | 1122
LCD | 4* | ? | 450403* | 0* | 363* | 0* | 2*

##  5 Conclusion

By a result of Sendrier [20], it is known that most linear codes are LCD when
$q$ is large. Moreover, the proportion of $q$-ary linear codes of length $n$,
dimension $k$ and specified hull dimension $l$ to all $q$-ary linear codes of
the same length and dimension is convergent when $n$ and $k$ goes to infinity.
Using the limit, Sendrier proved that the average dimension of the hull of a
$q$-ary linear code is asymptotically equal to $\sum_{i\geq
1}\frac{1}{q^{i}+1}$ [20].

In this paper, we obtain general results on hulls of linear codes, proving
that the number of all $q$-ary linear codes of a given length $n$ and
dimension $k$ decreases when the hull dimension increases, for all values of
$n$, $k$ and $q$.

In addition, we classify all binary linear, even, self-orthogonal and LCD
$[n\leq 20,k\leq 10,d\geq 2]$ codes and the ternary linear, self-orthogonal
and LCD $[n\leq 19,k\leq 10,d\geq 2]$ codes (with a few exceptions). For some
considered values of $n$ and $k$, when the number of all inequivalent linear
codes is huge, we classify only the optimal codes. The results are listed in
Tables 1 and 2.

## References

  * [1] Makoto Araya, Masaaki Harada, On the classification of linear complementary dual codes, Discrete Mathematics, Volume 342, Issue 1, 2019, 270–278.
  * [2] M. Araya, M. Harada, On the minimum weights of binary linear complementary dual codes, Cryptogr. Commun. 12, 285–300 (2020).
  * [3] M. Araya, M. Harada and K. Saito, Characterization and classification of optimal LCD codes, Des. Codes Cryptogr. 89, 617–640 (2021).
  * [4] Assmus, Jr. E. F., Key, J. D.: Affine and projective planes. Discrete Math. 83, 161–187 (1990).
  * [5] I. Bouyukliev, The Program Generation in the Software Package QextNewEdition. In: Bigatti A., Carette J., Davenport J., Joswig M., de Wolff T. (eds) Mathematical Software – ICMS 2020, Lecture Notes in Computer Science, vol 12097. Springer, Cham (2020).
  * [6] C. Carlet and S. Guilley, Complementary Dual Codes for Counter-measures to Side-Channel Attacks, Advances in Mathematics of Communications 10, 131–150 (2016).
  * [7] C. Carlet, S. Mesnager, C. Tang, Y. Qi and R. Pellikaan, Linear codes over $\mathbb{F}_{q}$ are equivalent to LCD codes for $q>3$, IEEE Trans. Inform. Theory 64, 3010–3017 (2018).
  * [8] C. Carlet, S. Mesnager, C. Tang, Y. Qi, New characterization and parametrization of LCD codes, IEEE Trans. Inf. Theory, 65(1), 39–49 (2019).
  * [9] Conway, J.H., Sloane, N.J.A., Sphere Packings, Lattices and Groups, third ed., Springer, New York, 1999.
  * [10] S.T. Dougherty, J.-L. Kim, B. Ozkaya, L. Sok and P. Solé, The combinatorics of LCD codes: linear programming bound and orthogonal matrices, Int. J. Inf. Coding Theory 4, 116–128 (2017).
  * [11] L. Galvez, J.-L. Kim, N. Lee, Y.G. Roe and B.-S. Won, Some bounds on binary LCD codes, Cryptogr. Commun. 10, 719–728 (2018).
  * [12] M. Grassl, Code Tables: Bounds on the parameters of various types of codes, http://www.codetables.de/.
  * [13] I. Bouyukliev, S. Bouyuklieva, T. A. Gulliver, and P. R. Östergård, “Classification of optimal binary self-orthogonal codes,” _Journal of Combinatorial Mathematics and Combinatorial Computing_ , vol. 59, pp. 33-87, 2006\.
  * [14] M. Harada, K. Saito, Binary linear complementary dual codes, Cryptogr. Commun. 11, 677–696 (2019).
  * [15] W. C. Huffman and V. Pless, Fundamentals of Error-Correcting Codes, Cambridge Univ. Press, 2003.
  * [16] Luo, G., Cao, X., Chen, X., MDS codes with hulls of arbitrary dimensions and their quantum error correction. IEEE Transactions on Information Theory, 65: 2944-2952, 2019.
  * [17] J.L. Massey, Linear codes with complementary duals, Discrete Math. 106/107, 337–342 (1992).
  * [18] McKay, B.: Isomorph-free exhaustive generation. J.Algorithms 26, 306–324 (1998).
  * [19] Vera Pless, The number of isotropic subspaces in a finite geometry, Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Nat. (8)39(1965), 418–421.
  * [20] Nicolas Sendrier: On the Dimension of the Hull. SIAM J. Discret. Math. 10(2): 282-293 (1997).