Delay Differential Neural Networks

Srinivas Anumasa P. K. Srijith

Computer Science and Engineering Indian Institute of Technology Hyderabad
{cs16resch11004<EMAIL_ADDRESS>
###### Abstract

Neural ordinary differential equations (NODEs) treat computation of
intermediate feature vectors as trajectories of ordinary differential equation
parameterized by a neural network. In this paper, we propose a novel model,
delay differential neural networks (DDNN), inspired by delay differential
equations (DDEs). The proposed model considers the derivative of the hidden
feature vector as a function of the current feature vector and past feature
vectors (history). The function is modelled as a neural network and
consequently, it leads to continuous depth alternatives to many recent ResNet
variants. We propose two different DDNN architectures, depending on the way
current and past feature vectors are considered. For training DDNNs, we
provide a memory-efficient adjoint method for computing gradients and back-
propagate through the network. DDNN improves the data efficiency of NODE by
further reducing the number of parameters without affecting the generalization
performance. Experiments conducted on synthetic and real-world image
classification datasets such as Cifar10 and Cifar100 show the effectiveness of
the proposed models.

##  1 INTRODUCTION

The ability of deep learning models to capture rich representations of high
dimensional data has led to its successful application in computer vision
problems like image classification (Huang et al., 2017, He et al., 2016), and
image captioning (Vinyals et al., 2015). The backbone of many such tasks is a
deep learning model called Residual Networks (ResNets) (He et al., 2016). They
allowed deep learning to solve complex computer vision tasks by training deep
neural networks with more than 100 layers without suffering from vanishing
gradient problem. ResNets achieve this by using skip connections where input
at any layer is added to the output of that layer modelling the identity
mapping.

The feature mappings in ResNet can be considered as discretization of a
continuous solution modelled using ordinary differential equations (ODE) (Lu
et al., 2018). Based on this idea, a generalization of the ResNet was
introduced, called neural ordinary differential equations (NODE) (Chen et al.,
2018). NODE is an ordinary differential equation parameterized by a neural
network and can be seen to generalize ResNets to arbitrarily many layers.
Neural ODE has been shown to achieve a performance close to ResNet but with a
reduced number of parameters. Like in ResNets, the representations learned
through NODE blocks are finally mapped to the output through a fully connected
neural network (FCNN). Neural ODEs are also shown to be more robust than
convolutional neural networks (CNN) (Hanshu et al., 2019).

In this paper, we propose a novel continuous depth neural network model, delay
differential neural networks (DDNN), inspired by delay differential equations
(DDE). The proposed model assumes the derivative of the feature vector depends
not only on the current feature vector but also on the feature vectors
computed in the past 111A concurrent work (Anonymous, 2021) combines neural
networks and delay differential equations using a similar approach.. DDNN can
be considered as a generalization to NODE models that utilize past feature
vectors in addition to the present feature vector for modelling the feature
mappings. Similar to NODE, the depth is not fixed but the number of feature
vectors required from the past is fixed. The proposed DDNN model will improve
upon the modelling capabilities of NODE, and further reduces the number of
parameters in a NODE model. DDNN can be seen as the continuous depth
generalization to deep neural network models like DenseNet (Huang et al.,
2017) and MixNet (Wang et al., 2018). They showed that designing a neural
network architecture by utilizing the previously computed feature vectors can
improve generalization performance of ResNets. Inspired by these recent
advances in ResNet variants, we consider different architectures of the DDNN
model, which differs in the way historical feature vectors are combined and
demonstrate their usefulness in the synthetic and image classification data
sets.

Our contributions can be summarized as follows:

  * •

We propose delay differential neural networks (DDNN), a continuous depth deep
learning model based on delay differential equations parameterized by neural
networks.

  * •

We propose two different DDNN architectures, considering two different
approaches to combine current and past feature vectors in a DDNN model.

  * •

We provide a memory-efficient adjoint method for updating the parameters of
DDNN during backpropagation.

  * •

We conduct experiments on synthetic data sets and image classification
datasets such as Cifar10 and Cifar100 datasets demonstrating the effectiveness
of DDNN.

##  2 RELATED WORK

NODE (Chen et al., 2018) generalizes ResNet (He et al., 2016) as continuous
time representation of discrete feature vectors obtained from a Residual block
transformations. In (Chen et al., 2018), a memory efficient adjoint method was
proposed for learning parameters in a NODE model. Following this work, a
series of approaches were proposed (Gholami et al., 2019, Dupont et al., 2019,
Durkan et al., 2019) to address the learning issues in NODE. It was observed
that generalization performance of NODE got degraded when it was made more
deeper by concatenating multiple NODE blocks. The issue was caused due to
inconsistency among the feature vectors computed during forward propagation
and backpropagation which lead to incorrect computation of gradients.
Adapative checkpoint and modified adjoint approaches (Zhuang et al., 2020,
Gholami et al., 2019) were proposed to address the computation of accurate
gradients. On the other hand, Augmented NODE (Dupont et al., 2019) improves
the function modelling capability of the NODE by augmenting the space (adding
additional dimensions) in which NODE is solved. In contrast, we propose to
consider feature representations from previous time steps to improve the
function modelling capability in NODE.

##  3 BACKGROUND

Let $\mathcal{D}=\\{X,\mathbf{y}\\}=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{N}$ be
set of training datapoints with $\mathbf{x}_{i}\in\mathcal{R}^{D}$ and
$y_{i}\in\\{1,\ldots,C\\}$. We denote the test data point as
$(\mathbf{x}_{*},y_{*})$. The aim is to learn a function which maps from input
$\mathbf{x}$ to a class label $y$ from the data so that it will have good
generalization performance. We assume the function learnt using a neural
network model to be denoted by $f$. The hidden layers in a neural network are
denoted as $\mathbf{z}$. In this section, we will provide background
information required to understand the proposed model.

###  3.1 ResNets and Variants

In standard deep learning models, increasing the depth often hampers the
performance of the model, due to vanishing gradient problem (Glorot and
Bengio, 2010, Bengio et al., 1994, He et al., 2016).

Figure 1: A residual block in ResNet architecture (He et al., 2016) with skip
connection

ResNet (He et al., 2016) overcomes the vanishing gradient problem by allowing
skip connections into the architecture as shown in the Figure 1. Consequently,
the function mapping for a layer $t$ is given by

$\textbf{z}_{t+1}=\textbf{z}_{t}+f(\textbf{z}_{t},\boldsymbol{\theta}_{t})$
(1)

Where, $f$ is a neural network mapping parameterized by
$\boldsymbol{\theta}_{t}$. ResNet showed improvement in accuracy with increase
in depth, where the performance of other deep learning models got degraded
with increase in depth. ResNets although proven to be powerful, training time
and parameter complexity has become a major hurdle for training a deeper
ResNet. Stochastic depth ResNet Huang et al. (2016) randomly drops subset of
layers during training, resulting in reduced training time. Motivated by the
success of deep learning models built on the idea of maximum information
(gradients or features) flow (He et al., 2016, Larsson et al., 2017), DenseNet
architecture was proposed. DenseNet (Huang et al., 2017) utilizes all the
computed feature in the past by concatenating across channel dimension. As
there is no need to re-learn the redundant features, fewer number of
parameters are required compared to ResNet.

Figure 2: A 5-layer dense block (Huang et al., 2017), each layer takes all
preceding features as input.

ResNet (He et al., 2016) and DenseNet (Huang et al., 2017) differ in the way
past feature vector are considered. While Resnet adds the past feature vectors
to the current one, DenseNet concatenates them with the current one. Mixed
Link Network (MixNet) (Wang et al., 2018) exploits the advantages of both the
models and proposed an architecture to learn better representations. They
improve upon DenseNet and ResNet by using a more general architecture which
considers both ResNet(addition) and DenseNet(concatenation) operations on of
features vectors (Please see Figure 3). MixNets provided improved parameter
efficiency while achieving state-of-the-art accuracy on benchmark datasets.

Figure 3: Path topology of MixNet (Wang et al., 2018). The symbols “$+$” and
“$||$” denote element-wise addition and concatenation respectively.

###  3.2 Neural Ordinary Differential Equation(ODE)

Deep learning models learn a sequence of transformation through different
spatial domain to map input $\mathbf{x}_{i}$ to output $y_{i}$. In a ResNet
block computation of a hidden layer representation can be expressed using the
following transformation.

$\mathbf{z}(t+1)=\mathbf{z}(t)+f(t,\mathbf{z}(t),\boldsymbol{\theta}(t))$ (2)

where $\textbf{z}(t)$ is feature vector with $t\in\\{0...T\\}$ and $f$ is
neural network parameterized by parameters $\boldsymbol{\theta}(t)$. If we use
the same transformation at every step, the expression in (2) is equivalent to
computing the trajectory of an ordinary differential equation (ODE) using
Euler method with step size one. This forms the basis of neural ODEs (NODE)
(Chen et al., 2018) and it can be seen as solving an ordinary differential
equation

$\frac{dz(t)}{dt}=f(t,\textbf{z}(t),\boldsymbol{\theta})$ (3)

Given initial feature vector $\textbf{z}(0)$ the final feature vector
$\textbf{z}(T)$ can be computed by solving the ODE (3). On the final feature
vector $\textbf{z}(T)$, necessary transformations are applied using a fully
connected neural network (FCNN), involving multiple linear mapping and
activation to predict class probabilities. Parameters of the NODE model are
learnt using cross-entropy loss function and memory-efficient adjoint approach
to compute gradients.

###  3.3 Delay Differential Equations(DDE)

In ODEs, the derivative of the state vector with respect to time is modelled
as a function of the current state vector at that time.

$\frac{d\mathbf{z}(t)}{dt}=g(t,\mathbf{z}(t))$ (4)

where $\mathbf{z}(t)\in\mathbb{R}^{d}$, $g:\mathbb{R}^{d+1}\to\mathbb{R}$. The
solution to the ODE, i.e. the state vector value at time $t$ is obtained as a
solution to the initial value problems(IVP) (with initial state
$\mathbf{z}(0)$) using numerical techniques like Euler or Runga-kutta method.

In delay differential equations (DDE), the derivative of the state vector with
respect to time can be seen as a function of the current state vector and
previously computed state vectors. DDEs fall under the class of infinite-
dimensional dynamical systems, where the function evolves in time. A typical
delay differential equation is given by (5), where,
$\tau_{1},\tau_{2}...\tau_{n}$ are constant delay values and $\phi(t)$ is the
history function.

$\displaystyle\frac{d\mathbf{z}(t)}{dt}$ $\displaystyle=$ $\displaystyle
g(t,\mathbf{z}(t),\mathbf{z}(t-\tau_{1}),...,\mathbf{z}(t-\tau_{n})),\quad
t>t_{0}$ $\displaystyle\textbf{z}(t)$ $\displaystyle=$
$\displaystyle\phi(t),\quad t\leq t_{0}$ (5)

In this case, the change in state vector not only depends on the current state
vector but also on the previous $n$ state vectors
$\mathbf{z}(t-\tau_{1}),...,\mathbf{z}(t-\tau_{n})$. For a DDE, due to this
dependence initial value of the state $\mathbf{z}(t_{0})$ alone is not enough
to compute the trajectory. It also requires a history function which allows it
to compute the state values at time $t\leq t_{0}$. In general, the delays can
be a function of state vectors or can be any function of time. However, we
restrict to positive valued delays.

Figure 4: Solution to delay logistic equation as shown in (6), with a =1.4

Let us consider an example of a DDE with one constant delay.

$\displaystyle\frac{dz(t)}{dt}$ $\displaystyle=$ $\displaystyle
az(t)(1-z(t-1)),\quad t>0$ $\displaystyle z(t)$ $\displaystyle=$
$\displaystyle 0.1,\quad t\leq 0$ (6)

The above DDE is called delay logistic equation, famously known for its
ability to capture the dynamics of population (Kuang, 1993). Figure 4 shows
the plot of the solution to (6) computed using a DDE solver. DDEs can be
solved as a sequence of IVPs for ODEs (considering the history function) and
can be shown to have a unique solution in the interval of interest. DDE with
constant delays have a solution which behave differently from the ODEs. In
contrast to ODEs which produce smoothly varying trajectories, DDE trajectory
can exhibit abrupt changes and jumps. In fact, one can show that DDEs exhibit
a first order discontinuity at the initial point and higher order
discontinuities in the intervals which are multiple of time delays. Hence, the
abrupt changes in the DDE trajectory may get smoothen out eventually. This
behaviour could be quite useful in modelling the changes in feature
representations in neural networks, and motivates us to develop delay
differential neural networks which we explain in detail in the following
section.

###  3.4 Numerical Method for DDE

The numerical methods to solve DDE can be obtained by extending the numerical
methods to solve an ODE such as, Euler method and Runga-kutta (RK) methods
(Shampine, 1994) by considering the delay term. For instance, a $q^{th}$
order, fixed step size, explicit RK method to solve a DDE (Ismail et al.,
2002) is shown as in (3.4). Where, $h$ is the step size and values of the
parameters $\\{{b_{i},c_{i}}\\}_{i=1}^{q}$, $\\{a_{ij}\\}_{i=1,j=1}^{q,i-1}$
depend on the order $q$ of RK method.

$\displaystyle\mathbf{z}(t+h)=\mathbf{z}(t)+h\sum_{i=1}^{q}b_{i}k_{i}$
$\displaystyle
k_{i}=f(t+c_{i}h,\mathbf{z}(t)+h\sum_{j=1}^{i-1}a_{ij}k_{j},\mathbf{z}(t+c_{i}h-\tau))$

In this work, to solve a DDNN model we opted for RK12 adpative numerical
method, which is a $2$ stage numerical method as shown in (8). With values,
$b_{1}=0.5$, $b_{2}=0.5$, $a_{21}=1$, $c_{1}=0$, $c_{2}=1$.

$\displaystyle\mathbf{z}(t+h_{t})=\mathbf{z}(t)+h_{t}(0.5k_{1}+0.5k_{2})$ (8)
$\displaystyle k_{1}=f(t,\mathbf{z}(t),\mathbf{z}(t-\tau))$ $\displaystyle
k_{2}=f(t+h_{t},\mathbf{z}(t)+h_{t}k_{1},\mathbf{z}(t+h_{t}-\tau))$
$\displaystyle error=-h_{t}0.5k_{1}+h_{t}0.5k_{2}$

Where, $h_{t}$ is the step size accepted by the solver, which depends on the
computed error. For computing history term $\mathbf{z}(t+h_{t}-\tau)$, we need
to have access to the computed feature vectors in the past and their accepted
time values. In our work, we choose linear interpolation (9) to compute a
feature vector at any time in the past.

$\displaystyle\mathbf{z}(t_{j})=\mathbf{z}(t_{i})+\frac{\mathbf{z}(t_{k})-\mathbf{z}(t_{i})}{t_{k}-t_{i}}(t_{j}-t_{i})$
(9)

Let, the solver requested a feature vector from the past with time $t_{j}$. If
no feature vector is computed at time $t_{j}$, a pair of time points
$(t_{i},t_{k})$ satisfying $t_{i}<t_{j}<t_{k}$ are chosen for computing
$\mathbf{z}(t_{j})$ using linear interpolation (9).

##  4 DELAY DIFFERENTIAL NEURAL NETWORKS

Neural ODEs have been shown to be quite useful in many tasks which were
earlier solved using ResNets. They allowed the neural networks to grow
arbitrary deeper without increasing the number of parameters. However, their
generalization performance was not close to ResNets, and several approaches
were proposed to improve upon it. In this section, we propose an approach,
delay differential neural networks (DDNN), which would improve the function
modelling capability of NODE and further reduces the number of parameters in
the model. The proposed approach utilizes the past feature vectors to compute
the next feature vector, through the framework of delay differential
equations.

The proposed delay differential neural networks (DDNN) consider the following
dynamics to model the feature representations $\mathbf{z}(t)$

$\displaystyle\frac{d\mathbf{z}(t)}{dt}$ $\displaystyle=$ $\displaystyle
f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta}),\quad
t>t_{0},\tau>0$ $\displaystyle\textbf{z}(t)$ $\displaystyle=$
$\displaystyle\phi(t),\quad t\leq t_{0}$ (10)

Where $f$ is a neural network parameterized by $\theta$ which consider not
only the current feature vector but also the previously computed feature
vector at time $t-\tau$. Here we restrict ourselves to a single constant
delay, while the proposed framework can be generalized to multiple delays as
well. The initial representation $\mathbf{z}(0)$ is obtained through some
basic operation (e.g. downsampling) of the input $\mathbf{x}$. This is
transformed through DDNN defined in (10) to obtain the final representation
$\mathbf{z}(T)$ at some time $T$. The final representation $\mathbf{z}(T)$ is
transformed using a fully connected neural network to obtain the output $y$.
Figure 5 provides a pictorial representation of the transformations in a DDNN
model. The final representation $\mathbf{z}(T)$ is obtained by solving the
initial value problem defined by the DDNN (10). The solution to the DDNN model
can be obtained using numerical techniques such as the Runga-Kutta method
discussed in Section 3.4. Solving DDNN requires one to know the values of
$\mathbf{z}(t)$ for $t\leq t_{0}$, and is provided by the history function
$\phi(t)$. We consider the history function $\phi(t)=\mathbf{z}(0)$ for $t\leq
t_{0}$. The key component in the proposed DDNN model is the function $f$ which
consider the current and past feature vectors. We propose two variants of the
DDNN model, which differs in the way these feature vectors are used. In one
model, we consider a concatenation of the feature vectors while in the second
model we consider a combination of the feature vectors.

Figure 5: Feature transformations in a DDNN model considering either
concatenation or combination operations over past and present feature vectors

###  4.1 Concatenation of Feature Vectors

We first consider a DDNN model (DDNN_Cat) which concatenates past feature
vector $\mathbf{z}(t-\tau)$ with the current feature vector $\mathbf{z}(t)$,
and provides this concatenated feature vector as an input to the function $f$.
This is loosely inspired from the DenseNet architecture which does the
concatenation of feature vectors. For the DDNN_Cat model, the size of the
input and output feature vectors fed to the neural network transformation are
not the same. The neural network is designed in such a way, that the incoming
higher dimensional feature vector is transformed to a lower dimensional
feature vector. This process is continued until the time component $t$ reaches
$T$. In our experiments for image classification, a DDNN block is constructed
based on convolutional neural networks (CNN) which deal with tensors valued
features with multiple channels. These features are concatenated across
channel dimension and the CNN computes a tensor of reduced channels. An
advantage of the DDNN_Cat is that it can lead to a reduced number of
parameters with a competitive generalization performance compared to NODE and
Resnet. With DDNN_Cat, one can work with a half the output size compared to
NODE but without affecting the performance due to its consideration of
previous feature vectors. The dynamics of the DDNN_Cat can be written as

$\displaystyle\frac{d\mathbf{z}(t)}{dt}$ $\displaystyle=$ $\displaystyle
f(t,h(\mathbf{z}(t),\mathbf{z}(t-\tau)),\boldsymbol{\theta}),\quad t>t_{0}$
$\displaystyle\textbf{z}(t)$ $\displaystyle=$ $\displaystyle\phi(t),\quad
t\leq t_{0}$ (11)

where $h$ denotes an operation over the feature vectors. In the case of
DDNN_Cat, $h(\mathbf{z}(t),\mathbf{z}(t-\tau))=z(t)||z(t-\tau)$. If we use
Euler method to compute the trajectory of the DDE, it will take the following
form,

$\mathbf{z}(t+1)=\mathbf{z}(t)+f(t,h(\mathbf{z}(t),\mathbf{z}(t-\tau)),\boldsymbol{\theta})$
(12)

The parameter $\tau$ here is a hyperparameter which is determined through the
validation data.

Table 1: Comparing adjoint methods used for NODE and DDNN during Backpropagation. where $f_{\mathbf{v}}(t)$ is $\frac{\partial f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\mathbf{z}(t-\tau)}$, $f_{\mathbf{z}}(t)$ is $\frac{\partial f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\mathbf{z}(t)}$, $f_{\boldsymbol{\theta}}(t)$ is $\frac{\partial f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}$ and $\mathbf{z}_{\boldsymbol{\theta}}(t)$ is $\frac{\partial\mathbf{z}(t)}{\partial\boldsymbol{\theta}}$ NODE |  DDNN
---|---
$\alpha^{T}(t)=-\frac{\partial L(\mathbf{z}(t))}{\partial\mathbf{z}(t)}$ |  $\alpha^{T}(t)=-\frac{\partial L(\mathbf{z}(t))}{\partial\mathbf{z}(t)}$
$\frac{d\alpha^{T}(t)}{dt}=-\alpha^{T}(t)f_{\mathbf{z}}(t)$ |  $\frac{d\alpha^{T}(t)}{dt}=-\alpha^{T}(t)f_{\mathbf{z}}(t)-\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)$
$\frac{dL}{d\boldsymbol{\theta}}=\int_{0}^{T}\alpha^{T}(t)f_{\boldsymbol{\theta}}(t)dt$ |  $\frac{dL}{d\boldsymbol{\theta}}=\int_{0}^{T}\alpha^{T}(t)f_{\boldsymbol{\theta}}(t)dt-\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt$

###  4.2 Convex Combination of Feature Vectors

The second DDNN model we consider is DDNN_Conv, which considers a convex
combination of the current feature vector and the past feature vector in the
neural network transformation in the DDNN block. In particular, it considers
the operation over feature vectors as follows

$h(\mathbf{z}(t),\mathbf{z}(t-\tau))=\lambda\mathbf{z}(t)+(1-\lambda)\mathbf{z}(t-\tau)$
(13)

where, $\lambda$ and $\tau$ are hyperparameters which can be determined based
on validation data. Note that in this case, the input and output dimension of
the neural network function in the DDNN block is the same unlike the DDNN_Cat
model and the transformations take this into account.

##  5 BACKPROPAGATION USING ADJOINT METHOD

We discuss an approach to learn the parameters in the DDNN model. The
parameters associated with the neural networks of various DDNN blocks can be
learnt through back-propagation. To apply back-propagation standard loss
functions for regression and classification can be used. The output $\hat{y}$
obtained from the FCNN block is used to compute the loss function $L$. The
loss function can be taken as the least squares loss for regression problems
and cross-entropy loss for classification problems, for e.g.
$L(y,\hat{y})=-y\log p(\hat{y})-(1-y)\log 1-p(\hat{y})$. However, applying
standard back propagation to learn the DDNN parameters, requires computing and
storing the gradients with respect to every state and is not memory efficient.
We use the adjoint method, a memory-efficient way of computing the gradients,
for learning the parameters in the DDNN model. This is similar to the adjoint
method used in the NODE (Chen et al., 2018, Gholami et al., 2019, Zhuang et
al., 2020) models, but we derive it to the case with delay differential
equations.

Initialize $\alpha^{T}(T),\frac{dL}{d\boldsymbol{\theta}}=0$

For $N$ to 1:

(1) Compute the value of the function $f$ at $t_{i}$ and $t_{i}+\tau$ for
computing $f_{\mathbf{z}}(t_{i})$,$f_{\boldsymbol{\theta}}(t)$ and
$f_{\mathbf{v}}(t_{i}+\tau)$.

(2) local backward, update $\alpha(t_{i})$ and
$\frac{dL}{d\boldsymbol{\theta}}$ according to discretization of (15) and 16.

Algorithm 1 Adjoint method during backpropagation for DDNN model

(a) DDNN_Conv (b) ACA_NODE

Figure 6: Showing the learning and extrapolating capabilities of (a)DDNN_Conv
and (b)ACA_NODE, given samples from the solution of DDE (17). Plots in the
left half of both (a) and (b), shows learned and actual trajectories of
individual components with respect to time. Plots in the right half of both
(a) and (b) represent phase portraits of both learned and actual solution of
DDE (17). This is obtained by plotting the individual components against each
other over time $t$, where, $x=\mathbf{z}(t)[0]$ and $y=\mathbf{z}(t)[1]$.

For the DDNN model, we extend the adjoint method in (Calver and Enright, 2017)
for learning the parameters of a DDE. Algorithm 1 provides a snippet of
adjoint method used for the DDNN model and provide the detailed derivation in
the appendix. The adjoint $\alpha^{T}(t)$ models the derivative of the loss
function with respect to some feature vector $\mathbf{z}(t)$ at time $t$

$\alpha^{T}(t)=-\frac{\partial L(y,\hat{y})}{\partial\mathbf{z}(t)}$ (14)

The adjoint can be used to compute the gradient of parameters
$\boldsymbol{\theta}$ with respect to loss $L$. We can model the dynamics of
the adjoint $\alpha^{T}(t)$ over time and is given by the following delay
differential equation

$\frac{d\alpha^{T}(t)}{dt}=-\alpha^{T}(t)f_{\mathbf{z}}(t)-\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)$
(15)

where $f_{\mathbf{v}}(t)$ is $\frac{\partial
f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\mathbf{z}(t-\tau)}$,
$f_{\mathbf{z}}(t)$ is $\frac{\partial
f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\mathbf{z}(t)}$.
The value of the adjoint $\alpha^{T}(t)$ at time $t$ can be found by solving
the delay differential equation (15) backwards in time. The numerical methods
for DDE based on Runga-Kutta methods discussed in Section 3.4 can be used to
solve (15). Computing the gradients of the parameters $\boldsymbol{\theta}$
with respect to loss $L$ requires evaluating an integral, which depends on
$\mathbf{z}(t)$ and $\alpha^{T}(t)$.

$\frac{dL}{d\boldsymbol{\theta}}=\int_{0}^{T}\alpha^{T}(t)f_{\boldsymbol{\theta}}(t)dt-\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt$
(16)

Where $f_{\boldsymbol{\theta}}(t)$ is $\frac{\partial
f(t,\mathbf{z}(t),\mathbf{z}(t-\tau),\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}$
and $\mathbf{z}_{\boldsymbol{\theta}}(t)$ is
$\frac{\partial\mathbf{z}(t)}{\partial\boldsymbol{\theta}}$. Table 1 provides
a comparison of adjoint methods used for training NODE and the proposed DDNN
model.

Our approach uses adjoint method for computing gradients during
backpropagation and uses adaptive checkpoint (Zhuang et al., 2020) framework
for computing accurate gradients. The adaptive checkpoint approach requires
storing the feature vectors computed during forward propagation and is then
utilized during backpropagation. This is shown to be numerically more accurate
than naive adjoint method (Chen et al., 2018). Assuming $N$ to be the number
of states evaluated, the approach requires a memory complexity of $O(N)$ for
each DDNN block.

Figure 7: Validation loss of DDNN_Conv model with different delay $\tau$
values trained over synthetic data.

##  6 EXPERIMENTS

We conduct experiments to study the performance of the proposed models on
synthetic and real world image datasets. In the experiments, we compare
against the state-of-the-art NODE model (Zhuang et al., 2020) model. We
conduct experiments on a toy dataset to compare the learning capability of
DDNN over ACA_NODE (Zhuang et al., 2020). We also conduct experiments for
image classification on the real datasets such as Cifar10 and Cifar100, and
compare against ACA_NODE and ResNet18 models.

###  6.1 Synthetic Data experiments

In this section we demonstrate the performance of DDNN_Conv on a 2-D toy
dataset. The synthetic data is generated from a 2-dimensional DDE (17).

$\displaystyle\frac{d\mathbf{z}(t)}{dt}$ $\displaystyle=$ $\displaystyle
0.75\mathbf{z}(t)+0.25\mathbf{z}(t-2.5)\times A,t>t_{0}$ $\displaystyle A$
$\displaystyle=$ $\displaystyle\begin{bmatrix}-0.1&3.2\\\
-3.2&-0.1\end{bmatrix}$ $\displaystyle\textbf{z}(t)$ $\displaystyle=$
$\displaystyle[-0.2,0.1],\quad t\leq t_{0}$ (17)

We sampled a set of 1000 samples from the solution of the DDE using a scipy-
based DDE solver222https://pypi.org/project/ddeint/ available in python. The
1000 samples are sampled on a time scale $0$ to $10$ as shown in Figure 6.
Samples from time scale $0$ to $6$ were chosen as training samples, time scale
$6$ to $8$ as validation samples. Samples from time scale $8$ to $10$ as test
samples. A two layer neural network with hidden dimension 50 is used for both
the ACA_NODE and our proposed DDNN_Conv model with $\lambda=0.75$. Mean-square
error loss is computed between computed trajectory and the actual trajectory,
which will be used to update the parameters of the model.

While training DDNN_Conv we chose $\tau$ to be a hyperparameter, and varied
the value from $1.5$ to $3.5$ with an interval of $0.25$. Based on the
validation loss we choose the best $\tau$ value. Figure 7 shows the
performance of DDNN_Conv models using different delay($\tau$) values. One can
observe the delay with value $2.5$ is having the lowest validation loss
compared to other delay values which corroborates well with the delay used in
the DDE (17). Figure 6(a) shows the reconstruction capability of DDNN model
with a delay $2.5$ which is having the lowest validation loss among DDNN_Conv
models with other delay values. Figures showing the learning capabilities of
DDNN_Conv models with other delay values can be found in the appendix.
Although DDNN_Conv model has seen only the samples on a time scale $0$ to $6$,
it was able to generate well the samples on a time scale $6$ to $10$. Figure
6(b) shows the learning and reconstruction capability of ACA_NODE model, which
fails to learn. From this experiment, we conclude that DDNN_Conv model able to
learn the underlying dynamics provided samples from the solution of a DDE,
where ACA_NODE fails.

Table 2: Validation accuracy of DDNN models, under different operations, convex combination and concatenation with different delay values. All the models are trained for 90 epochs MODELS |  DATASET | DELAY
---|---|---
|  |  0.1 |  0.2 |  0.3 |  0.4 |  0.5 |  0.6 |  0.7 |  0.8 |  0.9 |  1.0
DDNN_Conv |  Cifar10 |  93.50 |  93.62 |  93.60 |  93.58 |  93.30 |  93.74 |  93.62 |  93.78 |  93.82 |  93.98
|  Cifar100 |  74.56 |  73.50 |  73.62 |  74.42 |  74.72 |  74.74 |  73.14 |  76.00 |  74.86 |  74.72
DDNN_Cat |  Cifar10 |  92.98 |  93.44 |  92.96 |  93.02 |  93.22 |  92.88 |  92.94 |  92.88 |  93.18 |  93.24
|  Cifar100 |  72.20 |  72.14 |  71.92 |  72.18 |  72.84 |  73.26 |  71.32 |  72.52 |  72.58 |  72.94
Table 3: Comparing test accuracy of ACA_NODE, ResNet18 models against DDNN_Cat and DDNN_Conv models. For the columns DDNN_Cat and DDNN_Conv delay values of the best models based on test accuracy are also given. DATASET |  ACA_NODE |  DDNN_Conv |  DDNN_Cat |  RESNET18
---|---|---|---|---
Cifar10 |  94.16 |  93.96(delay 1.0) |  94.30(delay 0.2) |  94.28
Cifar100 |  74.60 |  76.52(delay 0.8) |  74.34(delay 0.6) |  76.08

###  6.2 Image Classification

Under image classification, we conducted experiments on Cifar10 and Cifar100
datasets. Cross-entropy loss is used to compute gradients with respect to the
parameters. We split the dataset containing 60000 samples into 50000 as
training, 5000 as validation data, and 5000 as test data. For DDNN models,
validation data will be used to select the best delay hyperparameter. For
training the models we follow a similar setup as discussed in (Zhuang et al.,
2020), where all the models are trained for 90 epochs on both Cifar10 and
Cifar100. Scheduled decay learning rate is used with an initial learning rate
as $0.1$ with decay as $0.1$. The learning rate is rescheduled at epoch 30 and
60.

####  6.2.1 DDNN Convex Combination

We first study the performance of the DDNN model (DDNN_Conv) which consider a
convex combination of feature vectors. The parameter $\lambda$ was fixed to
$0.75$ after analysing the performance on the validation data. This
essentially means that we give a higher importance to the current feature
vector and a lower importance to the previous feature vector. As we stated
earlier, the delay($\tau$) value is treated as a hyperparameter, 10 different
DDNN_Conv models were trained with different $\tau$ values, ranging from $0.1$
to $1.0$ with an interval of $0.1$. The architecture design and the number of
parameters is same as of ACA_NODE.

In Table 2, the validation accuracies of DDNN_Conv models with different delay
($\tau$) values are provided. As can be observed from Table 2, the best
DDNN_Conv model for Cifar10 dataset is the one with delay $1.0$ and for
Cifar100 is one with the delay $0.8$. Thus, the DDNN_Conv model with the best
validation accuracy is found to be the one with higher delay value. The
possible explanation could be, re-utilizing the older features (higher delay
values) compared to re-utilizing the recently computed features (low delay
values) may be helping in better accuracy. From Table 3, one can observe that
the test accuracy of the best DDNN_Conv model for Cifar10 is $93.96\%$ and for
Cifar100 is $76.52\%$. On Cifar100 dataset, DDNN_Conv model achieved a
higher($76.52\%$) accuracy compared to ACA_NODE($74.60\%$) and
ResNet18($76.08\%$) models, showing the effectiveness of DDE based deep
leaning models.

####  6.2.2 Concatenation

The architecture design of DDNN_Cat is different when compared with DDNN_Conv
and ACA_NODE. In DDNN_Conv or ACA_NODE models, dimension of the expected input
feature vector (Tensor) to DDNN block or ACA_NODE block is same as the
computed feature vector (Tensor) dimension. But for DDNN_Cat, the input tensor
and output tensor differ in terms of number of channels. Let $C$ be expected
number of channels by both the models DDNN_Conv and DDNN_Cat. But, the
computed tensor is different in terms of number of channels. DDNN_Cat produces
a tensor with $C/2$ number of channels unlike $C$ for DDNN_Conv. This results
in a reduced number of parameters for DDNN_Cat ($7$ Million) compared to
DDNN_Conv and ACA_NODE models ($11$ Million).

We trained the DDNN_Cat model with different delay values and compared their
validation performance in Table 2. Test accuracy of the best performing
DDNN_Cat model trained on Cifar10 and Cifar100 datasets are given in Table 3.
We can find that the DDNN_Cat model We found that the DDNN_Cat model is very
competitive and even performed better than other models in the Cifar10
dataset, achieving an accuracy of $94.30$ for a delay value of $0.2$. Thus,
DDNN_Cat model gave a very good performance in spite of reduction in the
number of parameters.

##  7 CONCLUSION

In this work we proposed a novel continuous depth deep learning model, delay
differential neural networks (DDNN), based on the principles of delay
differential equations. We provided two different DDNN architectures, DDNN_cat
and DDNN_Conv based on the way DDNN block consider the present and past
feature vectors. We also discussed an adjoint method which provides a memory
efficient way to learn parameters in DDNN. Through experiments on toy-dataset
we showed that DDNN model can learn the trajectory of the solution to a DDE,
where ACA_NODE fails. We discussed the performance of DDNN models on image
classification datasets such as Cifar10 and Cifar100. We showed that under
concatenation operation, DDNN model with reduced number of parameters,
performs well without affecting the generalization performance. As a future
work, we develop more general DDNN models and apply to other problems in
computer vision.

## References

  * Anonymous (2021) Anonymous.  Neural delay differential equations.  In _Submitted to International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=Q1jmmQz72M2.  under review.
  * Bengio et al. (1994) Y. Bengio, P. Simard, and P. Frasconi.  Learning long-term dependencies with gradient descent is difficult.  _IEEE transactions on neural networks_ , 5(2):157–166, 1994.
  * Calver and Enright (2017) J. Calver and W. Enright.  Numerical methods for computing sensitivities for odes and ddes.  _Numerical Algorithms_ , 74(4):1101–1117, 2017\.
  * Chen et al. (2018) R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In _Advances in neural information processing systems_ , pages 6571–6583, 2018.
  * Dupont et al. (2019) E. Dupont, A. Doucet, and Y. W. Teh.  Augmented neural odes.  In _Advances in Neural Information Processing Systems_ , pages 3140–3150, 2019.
  * Durkan et al. (2019) C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios.  Neural spline flows.  In _Advances in Neural Information Processing Systems_ , pages 7511–7522, 2019.
  * Gholami et al. (2019) A. Gholami, K. Keutzer, and G. Biros.  Anode: unconditionally accurate memory-efficient gradients for neural odes.  In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , pages 730–736. AAAI Press, 2019.
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 249–256, 2010.
  * Hanshu et al. (2019) Y. Hanshu, D. Jiawei, T. Vincent, and F. Jiashi.  On robustness of neural ordinary differential equations.  In _International Conference on Learning Representations_ , 2019.
  * He et al. (2016) K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Huang et al. (2016) G. Huang, Y. Sun, Z. Liu, D. Sedra, and K. Q. Weinberger.  Deep networks with stochastic depth.  In _European conference on computer vision_ , pages 646–661. Springer, 2016.
  * Huang et al. (2017) G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4700–4708, 2017.
  * Ismail et al. (2002) F. Ismail, R. A. Al-Khasawneh, M. Suleiman, et al.  Numerical treatment of delay differential equations by runge-kutta method using hermite interpolation.  _MATEMATIKA: Malaysian Journal of Industrial and Applied Mathematics_ , 18:79–90, 2002.
  * Kuang (1993) Y. Kuang.  _Delay differential equations: with applications in population dynamics_.  Academic press, 1993.
  * Larsson et al. (2017) G. Larsson, M. Maire, and G. Shakhnarovich.  Fractalnet: Ultra-deep neural networks without residuals.  In _ICLR_ , 2017.
  * Lu et al. (2018) Y. Lu, A. Zhong, Q. Li, and B. Dong.  Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations.  In _International Conference on Machine Learning_ , pages 3276–3285. PMLR, 2018.
  * Shampine (1994) L. F. Shampine.  _Numerical solution of ordinary differential equations_.  CRC Press, 1994.
  * Vinyals et al. (2015) O. Vinyals, A. Toshev, S. Bengio, and D. Erhan.  Show and tell: A neural image caption generator.  _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , Jun 2015.
  * Wang et al. (2018) W. Wang, X. Li, T. Lu, and J. Yang.  Mixed link networks.  In _Proceedings of the 27th International Joint Conference on Artificial Intelligence_ , pages 2819–2825, 2018.
  * Zhuang et al. (2020) J. Zhuang, N. Dvornek, X. Li, S. Tatikonda, X. Papademetris, and J. Duncan.  Adaptive checkpoint adjoint method for gradient estimation in neural ODE.  In _Proceedings of the 37th International Conference on Machine Learning_ , Proceedings of Machine Learning Research, pages 11639–11649. PMLR, 2020.

##  Appendix A ADJOINT METHOD

In this section we give the derivation of the adjoint method based on Calver
and Enright (2017), which is used to update the parameters of the DDNN during
backpropagation. The derivations provided are for under one constant delay and
this can be easily extended to multiple constant delays. Given DDE with a
history function

$\displaystyle\frac{d\mathbf{z}(t)}{dt}=f(\mathbf{z}(t),\mathbf{z}(t-\tau),t,\boldsymbol{\theta}),t>t_{0}$
$\displaystyle\textbf{z}_{t}=\phi(t),t\leq t_{0}$ (18)

$f$ is represented by a neural network which is differentiable almost
everywhere. The forward pass is defined as:

$\hat{y}=\mathbf{z}(T)=\mathbf{z}(0)+\int_{0}^{T}f(\mathbf{z}(t),\mathbf{z}(t-\tau_{1}),t,\boldsymbol{\theta})dt$
(19)

The loss function is defined as:

$J(\hat{y},y)=J(\mathbf{z}(T),y))$ (20)

To minimize the loss function we need to compute the
$\frac{dJ(\hat{y},y)}{d\boldsymbol{\theta}}$. Let $\alpha^{T}(t)$ be any
vector valued function defined for $t\in[0,T]$, then the lagrangian is,

$L=J(\mathbf{z}(T),\boldsymbol{\theta})+\int_{0}^{T}\alpha^{T}(t)[\frac{\mathbf{z}(t,\boldsymbol{\theta})}{dt}-f(t,\mathbf{z}(t,\boldsymbol{\theta}),\mathbf{z}(t-\tau,\boldsymbol{\theta}),\boldsymbol{\theta})]dt$
(21)

For ease of notation we write,

$\displaystyle\frac{\partial
f(t,\mathbf{z}(t,\boldsymbol{\theta}),\mathbf{z}(t-\tau,\boldsymbol{\theta}),\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\hskip
5.69046pt$ $\displaystyle\text{as}\hskip
5.69046ptf_{\boldsymbol{\theta}}(t),\frac{\partial
f(t,\mathbf{z}(t,\boldsymbol{\theta},\mathbf{z}(t-\tau,\boldsymbol{\theta}),\boldsymbol{\theta})}{\partial\mathbf{z}(T,\boldsymbol{\theta})}\hskip
5.69046pt\text{as}\hskip 5.69046ptf_{\mathbf{z}}(t),\frac{\partial
f(t,\mathbf{z}(t,\boldsymbol{\theta},\mathbf{z}(t-\tau,\boldsymbol{\theta}),\boldsymbol{\theta})}{\partial\mathbf{z}(t-\tau)}\hskip
5.69046pt\text{as}\hskip 5.69046ptf_{\mathbf{v}}(t),$
$\displaystyle\frac{\partial\mathbf{z}(t,\theta)}{\partial\boldsymbol{\theta}}\hskip
5.69046pt$ $\displaystyle\text{as}\hskip
5.69046pt\mathbf{z}_{\boldsymbol{\theta}}(t),\frac{\partial\mathbf{z}(t-\tau,\theta)}{\partial\boldsymbol{\theta}}\hskip
5.69046pt\text{as}\hskip 5.69046pt\mathbf{z}_{\boldsymbol{\theta}}(t-\tau)$
(22)

taking derivative w.r.t to $\boldsymbol{\theta}$

$\frac{dL}{d\boldsymbol{\theta}}=\frac{\partial
L(\mathbf{z}(T,\theta))}{\partial\mathbf{z}(T,\theta)}\mathbf{z}_{\boldsymbol{\theta}}(T)+\underbrace{\int_{0}^{T}\alpha^{T}(t)(\frac{d(\mathbf{z}_{\boldsymbol{\theta}}(t))}{dt})dt}_{\text{
applying integration by
parts}}-\int_{0}^{T}\alpha^{T}(t)[f_{\boldsymbol{\theta}}(t)+f_{\mathbf{z}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t)+f_{\mathbf{v}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t-\tau)]dt$
(23)

$\int_{0}^{T}\alpha^{T}(t)(\frac{d(\mathbf{z}_{\boldsymbol{\theta}}(t))}{dt})dt=[\alpha^{T}(t)\mathbf{z}_{\boldsymbol{\theta}}(t)]^{T}_{0}-\int_{0}^{T}\alpha^{T}(t)[\frac{d\alpha^{T}(t)}{dt}\mathbf{z}_{\boldsymbol{\theta}}(t)]$
(24)

$\mathbf{z}(0)$ is initial state doesnt depend on $\boldsymbol{\theta}$, so
$\mathbf{z}_{\boldsymbol{\theta}}(0)$ is $0$.

$=\underbrace{\frac{\partial
L(\mathbf{z}(T,\theta))}{\partial\mathbf{z}(T,\theta)}\mathbf{z}_{\boldsymbol{\theta}}(T)+\alpha^{T}(T)\mathbf{z}_{\boldsymbol{\theta}}(T)}_{\text{Equating
this part to
0}}-\int_{0}^{T}\alpha^{T}(t)[f_{\boldsymbol{\theta}}(t)+f_{\mathbf{z}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t)+f_{\mathbf{v}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t-\tau)]dt$
(25)

$\left[\frac{\partial
L(\mathbf{z}(T,\theta))}{\partial\mathbf{z}(T,\theta)}+\alpha^{T}(t)\right]\mathbf{z}_{\boldsymbol{\theta}}(T)=0$
(26)

$\boxed{\alpha^{T}(t)=-\frac{\partial
L(\mathbf{z}(T,\theta))}{\partial\mathbf{z}(T,\theta)}}$ (27)

This choice for $\alpha^{T}(T)$ eliminates the integral involving
$\mathbf{z}_{\boldsymbol{\theta}}(T)$

$=\int_{0}^{T}\frac{d\alpha^{T}(t)}{dt}\mathbf{z}_{\boldsymbol{\theta}}(t)dt-\int_{0}^{T}\alpha^{T}(t)[f_{\boldsymbol{\theta}}(t)+f_{\mathbf{z}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t)]-\underbrace{\int_{0}^{T}\alpha^{T}(t)f_{\mathbf{v}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t-\tau)dt}_{\text{Apply
change of variables}}$ (28)

$\int_{0}^{T}\alpha^{T}(t)f_{\mathbf{v}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t-\tau)dt=\int_{-\tau}^{T-\tau}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)]dt$
(29)

After change of variables

$=\int_{0}^{T}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt+\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt]$
(30)

After substituting Equation 30 in 28.

$\begin{split}\frac{\partial
L}{\partial\theta}=&\int_{0}^{T}\frac{d\alpha^{T}(t)}{dt}\mathbf{z}_{\boldsymbol{\theta}}(t)dt-\int_{0}^{T}\alpha^{T}(t)[f_{\boldsymbol{\theta}}(t)+f_{\mathbf{z}}(t)\mathbf{z}_{\boldsymbol{\theta}}(t)]\\\
&-\int_{0}^{T}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt-\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt\end{split}$
(31)

After rearranging the terms.

$\begin{split}&=\int_{0}^{T}\left[-\frac{d\alpha^{T}(t)}{dt}-\alpha^{T}(t)f_{\mathbf{z}}(t)-\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\right]\mathbf{z}_{\boldsymbol{\theta}}(t)dt-\int_{0}^{T}\alpha^{T}(t)f_{\boldsymbol{\theta}}(t)dt\\\
&\hskip
56.9055pt\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt\end{split}$
(32)

The dynamics of the adjoint $\alpha^{T}(t)$ given by the differential
equation.

$\boxed{\frac{d\alpha^{T}(t)}{dt}=-\alpha^{T}(t)f_{\mathbf{z}}(t)-\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)}$
(33)

Computing the gradients with respect to the parameters $\boldsymbol{\theta}$
requires evaluating a integral, which depends on $\mathbf{z}(t)$ and
$\alpha^{T}(t)$.

$\boxed{\frac{dL}{d\boldsymbol{\theta}}=\int_{0}^{T}\alpha^{T}(t)f_{\boldsymbol{\theta}}(t)dt-\int_{-\tau}^{0}\alpha^{T}(t+\tau)f_{\mathbf{v}}(t+\tau)\mathbf{z}_{\boldsymbol{\theta}}(t)dt}$
(34)

##  Appendix B PERFORMANCE OF DDNN MODELS ON TOY-DATASET USING DIFFERENT DELAY
VALUES

(a) DDNN with delay 1.5 (b) DDNN with delay 1.75

(c) DDNN with delay 2.00 (d) DDNN with delay 2.25

(e) DDNN with delay 2.50 (f) DDNN with delay 2.75

(g) DDNN with delay 3.00 (h) DDNN with delay 3.25

(i) DDNN with delay 3.50

Figure 8: Showing the learning and extrapolating capabilities of DDNN models
with different delay values, training samples are from the solution of DDE
with delay 2.5.


# Harer–Zagier formulas for families of twisted hyperbolic knots

Andreani<EMAIL_ADDRESS>Okinawa Institude of Science and
Technology Graduate University, 1919-1 Tancha, Okinawa 904-0495, Japan.
Shinobu<EMAIL_ADDRESS>Okinawa Institude of Science and Technology
Graduate University, 1919-1 Tancha, Okinawa 904-0495, Japan.

###### Abstract

In an attempt to generalise knot matrix models for non-torus knots, which
currently remains an open problem, we derived formulas for the Harer–Zagier
transform of the HOMFLY–PT polynomial for some infinite families of twisted
hyperbolic knots. Among them, we found a family of Pretzel knots for which the
transform has a fully factorised form, while for the remaining families
considered it consists of sums of factorised terms. Their zeros have a
remarkable structure as the modulus of their product in all cases equals
unity.

#### Keywords

Knot matrix models, Superintegrability, Harer–Zagier transform, HOMFLY–PT
polynomial, Recursive formulas, Twisted hyperbolic knots, Pretzel knots

##  1 Introduction

Among the many uses of knots by humans since antiquity, their ability to store
information is remarkable. In particular, the ancient Chinese and Incan
civilisations used knotted strings as an alternative to writing. In the late
19th century, Lord Kelvin, through his vortex atom hypothesis [1], envisioned
to use knots to encode information about nature. Although his theory thwarted,
it gave birth to Knot theory, the mathematical study of knots. Among its major
achievements was the discovery of knot polynomial invariants, such as the
Alexander and Jones polynomials, or their 2 variable generalisation, called
the HOMFLY-PT polynomial. About a century later, a revolutionary work by E.
Witten [2] attributed a physical interpretation to such invariants, as
observables of Chern–Simons theory, hence reconnecting knots with physics and
resulting into a fruitful interchange. Indeed, more recently with the
development of matrix models, there has been an active effort to explore this
interrelation more deeply; and it is towards this goal that the present work
aims to contribute.

###  1.1 Chern–Simons theory and knot invariants

Chern–Simons (CS) theory is a Topological Quantum Field Theory on a
3-dimensional manifold that is invariant under the action of a gauge group
$G$. The Wilson loop operators $W_{\mathcal{K}}^{R}$ are the traces of
holonomies around a knot $\mathcal{K}$, evaluated in an irreducible
representation $R$ of $G$. The averages of these $\langle
W^{R}_{\mathcal{K}}\rangle$ are quantum, gauge invariant observables of the
theory. The special case when the manifold is $\mathbb{S}^{3}$, $G=SU(N)$ and
$R=\square$ (the fundamental representation) the observables yield the HOMFLY-
PT polynomial of the knot $\mathcal{K}$ (defined below in sec. 2) as
$\bar{H}_{\mathcal{K}}(q^{N},q)=\langle W_{\mathcal{K}}^{\square}\rangle$,
where $q$ depends on $N$ and $k$, the level (or coupling constant) of CS
theory [2, 3]. It can be ‘colored’ by different choices of irreducible
representations $R$, resulting in the colored HOMFLY (henceforth omitting –PT)
polynomial $\bar{H}_{\mathcal{K}}^{R}(q^{N},q)=\langle
W_{\mathcal{K}}^{R}\rangle$. It is a generalisation of both the Jones and
Alexander polynomials, which correspond to the particular cases $N=2$ and
$N=0$, respectively.

CS theory on $\mathbf{S}^{3}$ with gauge group $U(N)$ also admits a matrix
model formulation with measure for the average (up to constant factors) given
by

$\langle F\rangle_{CS}\sim\intop
F\prod_{i<j}^{N}\left(2\sinh\left(\frac{x_{i}-x_{j}}{2}\right)\right)^{2}\prod_{i=1}^{N}dx_{i}e^{-x_{i}^{2}/2g},$
(1)

where $\\{x_{i}\\}_{i=1}^{N}$ are the eigenvalues of an $N\times N$ Hermitian
matrix, $F$ is a function of the $\\{x_{i}\\}$, $g=\frac{2\pi i}{k+N}$ and the
factor in the bracket is known as the trigonometric Van-der-Monde function
[4].

###  1.2 Knot matrix models

More recently, Morozov et al. [5] conjectured a connection between knot
polynomial invariants and matrix models via the _superintegrability condition_

$\langle\chi^{R}\rangle_{\mathcal{K}}=\bar{H}^{R}_{\mathcal{K}}(q^{N},q).$ (2)

Superintegrability means that a complete set of averages are explicitly
calculable; and it is established that for (Hermitian) eigenvalue matrix
models the averages of characters are known to be again characters, i.e
$\langle\chi^{R}\rangle\sim\chi^{R}$ [6, 7]. Due to the dependence of Wilson
loop averages on representations, knot polynomial invariants can be thought of
as non-trivial generalisations of characters333In particular, the HOMFLY
polynomial for torus knots can be expressed in terms of Schur functions, see
[8] for details., hence allowing to use the condition $\left\langle
character\right\rangle=knot\;polynomial$ as the defining property (2).

Knot matrix models are, thus far, only consistently defined for the particular
case of torus knots444A torus knot (or link, when $(m,n)$ are not coprime) is
described algebraically as the intersection of the 3-sphere with a singular
complex curve
$V=\\{(\alpha,\beta)\in\mathbb{C}^{2}|\alpha^{m}-\beta^{n}=0\\}$, i.e.
$T(m,n)=T(n,m)=V\cap\mathbb{S}^{3}$. The integers $(m,n)$ give the number of
strands (toroidal windings) and number of leaves (poloidal windings),
respectively. $T(m,n)$, for which there exists an eigenvalue matrix model, the
TBEM model [4, 9], providing an explicit measure in the left hand side of (2),
given by

$\langle\chi^{R}\rangle_{T(m,n)}\sim\intop\chi^{R}\prod_{i<j}^{N}\sinh\left(\frac{x_{i}-x_{j}}{m}\right)\sinh\left(\frac{x_{i}-x_{j}}{n}\right)\prod_{i=1}^{N}dx_{i}e^{-x_{i}^{2}/2g}.$
(3)

Here $q=e^{\frac{g}{2mn}}$ and note that the trigonometric Van-der-Monde
function is $(m,n)$-deformed, but otherwise this expression is identical with
the one for the CS matrix model (1).

The  Harer–Zagier (HZ) transform, which is a discrete version of the Laplace
transform in $N$ explicitly given by

$Z_{\mathcal{K}}(q,\lambda)=\sum_{N=0}^{\infty}\bar{H}_{\mathcal{K}}(q^{N},q)\lambda^{N}$
(4)

provides an alternative manifestation of superintegrability:

$the\;HZ\;tranforms\;are\;completely\;factorised\;rational\;functions,$ (5)

i.e. they have zeroes and poles at positive and negative powers of $q$. This
is true for the case of torus knots, as shown in [5] using the quantum groups
technology and reconfirmed (via a different method) in the present work; and
it should be a minimum consistency requirement of any extension of the
definition (2) of knot matrix models to other families of knots. As a first
check, the HZ formula for the HOMFLY polynomial of the simplest hyperbolic
knot, the figure-8, was also computed in [5]. This turned out to be not
factorisable, and hence superintegrability fails in this case.

However, continuing this effort, in this article we derive the HZ formulas for
some infinite families of ‘twisted’ hyperbolic knots (which shall be described
in more detail below) and examine their factorisability properties (sec. 2.2),
their $q\rightarrow 1$ expansion (sec. 3.1), their poles (sec. 3.2) and zero
loci (sec. 3.3). The Appendix includes the HZ formulas and the $q\rightarrow
1$ expansion coefficients for some further families of twisted hyperbolic
knots.

##  2 The HOMFLY polynomial and its Harer–Zagier tranform

Figure 1: Resolving tree for $T(2,5)$

The HOMFLY polynomial $H_{\mathcal{K}}(v,z)$ of an oriented knot is a Laurent
polynomial in two variables, defined by the normalisation condition
$H_{\text{unknot}}=1$ and the _skein relation_

$v^{-1}H_{L_{+}}(v,z)-vH_{L_{-}}(v,z)=zH_{L_{0}}(v,z)$ (6)

where $L_{+}=$, $L_{-}=$, $L_{0}=$. For two disconnected knots
$\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ the product formula
$H_{\mathcal{K}_{1}\sqcup\mathcal{K}_{2}}=(v^{-1}-v)z^{-1}H_{\mathcal{K}_{1}}H_{\mathcal{K}_{2}}$
holds [10].

The unnormalised HOMFLY $\bar{H}_{\mathcal{K}}(v,z)$, is obtained by
multiplying with an overall factor $-(v^{-1}-v)z^{-1}=:\bar{H}_{unknot}(v,z)$;
and with the substitution $v=q^{N}$ and $z=q-q^{-1}$, where
$q=e^{\nicefrac{{\pi i}}{{(k+N)}}}$, we obtain
$\bar{H}_{\mathcal{K}}(q^{N},q)$ as arises from CS theory described
above555Due to a discrepancy in conventions between the mathematics and
physics literature, this holds up to some minus signs. For instance, an extra
overall minus sign is included in $\bar{H}_{unknot}$ in order to be in
agreement with the Wilson loop average of the circle in standard framing, as
derived in [2]. Such ambiguities, however, do not affect the essence of the
results in this article.. The HOMFLY polynomial can be combinatorially
computed using skein trees, as the one shown in the figure 1 for the example
of the torus knot $T(2,5)$, with the help of which we have derived recursive
or explicit formulas and computed their HZ transform for the following
families of knots.

###  2.1 Torus knots and links

The HOMFLY polynomial for 2–stranded torus knots and links666Whenever we refer
to torus links, we assume that all the components have parallel orientation,
as shown for example for $T(2,4)$ and $T(2,2)$ in fig. 1. The recursive
formula (7) is not valid for links with different relative orientation. can be
obtained by the following recursive relations with initial condition
$H_{T(2,2)}=\frac{v}{z}\left(1-v^{2}+z^{2}\right)$

$H_{T(2,n)}=v^{2}H_{T(2,n-2)}+zvH_{T(2,n-1)},\;\forall\;n\geq 3.$ (7)

For odd $n=2k+1$ with $k=1,2,3,...$ this can be restricted to knots only:

$H_{T(2,2k+1)}=v^{2}H_{T(2,2k-1)}+z^{2}\sum_{j=1}^{k}v^{2j}H_{T(2,2(k-j)+1)}+v^{2k}(1-v^{2})H_{T(2,1)}.$
(8)

For 3–stranded torus knots and links there are 3 different recursive formulas
corresponding to $n\mod 3=\\{0,1,2\\}$, with initial condition
$H_{T(3,3)}=v^{4}z^{2}\left(2-v^{2}+z^{2}\right)+v^{4}z^{-2}\left(1-v^{2}+z^{2}\right)\left(1-v^{2}+2z^{2}\right)$.
\- $\forall\;n\mod 3=2$, i.e. $n=2,5,8,...$

$H_{T(3,n)}=v^{2}H_{T(3,n-1)}+z^{2}\sum_{j=1}^{n-1}v^{2j}H_{T(3,n-j)}+v^{2(n-1)}(1-v^{2})H_{T(3,1)},$
(9)

\- $\;\forall\;n\mod 3=1\;(\geq 4)$, i.e. $n=4,7,10,...$

$H_{T(3,n)}=v^{4}H_{T(3,n-2)}+v^{2}z^{2}H_{T(3,n-1)}+2z^{2}\sum_{j=2}^{n-1}v^{2j}H_{T(3,n-j)}+2v^{2(n-1)}(1-v^{2})H_{T(3,1)},$
(10)

\- $\forall\;n\mod 3=0\;(\geq 6)$, i.e. $n=6,9,12,...$

$H_{T(3,n)}=v^{6}H_{T(3,n-3)}+v^{2}z^{2}H_{T(3,n-1)}+2v^{4}z^{2}H_{T(3,n-2)}+3z^{2}\sum_{j=3}^{n-1}v^{2j}H_{T(3,n-j)}+3v^{2(n-1)}(1-v^{2})H_{T(3,1)},$
(11)

where $T(m,1)$ is the unknot $\forall\;m$. To our knowledge, formulas (8)–(11)
are new results as they are nowhere to be found in the literature. Due to the
fact that skein trees are not unique and grow fast even at $m=4$, we were
unable to obtain a recursive relation for general777However, we did obtain the
general recursion formula
$V_{T(m,n)}(q)=q^{2(m-1)}V_{T(m,n-2)}(q)+(1-q^{2(m-1)})q^{(m+1)(n+1)}$ for the
single variable Jones polynomial corresponding to $N=2$, i.e.
$V_{\mathcal{K}}(q)=H_{\mathcal{K}}(q^{2},q-q^{-1})$. $(m,n)$. Hence, we used
instead the explicit formula for the HOMFLY polynomial of torus knots only
(i.e. for $m,n$ coprime) given in [11], which in our conventions reads

$\bar{H}_{T(m,n)}(q^{N},q)=\frac{q^{N}-q^{-N}}{q-q^{-1}}(q^{N}q)^{(m-1)(n-1)}\frac{1-q^{-2}}{1-q^{-2m}}{\displaystyle\sum_{\beta=0}^{m-1}q^{-2n\beta}\left(\prod_{i=1}^{\beta}\frac{q^{2N}q^{2i}-1}{q^{2i}-1}\right)\left(\prod_{j=1}^{m-1-\beta}\frac{q^{2N}-q^{2j}}{1-q^{2j}}\right)}.$

The corresponding HZ transform can be computed by applying the geometric
series to each $q^{N}$ power. Doing this calculation for sufficiently many
torus knots of fixed $m$ and arbitrary $n$, we inductively deduced the
following factorised formula

$Z_{T(m,n)}(q,\lambda)=\frac{\lambda\prod_{j=0}^{m-2}\left(1-\lambda
q^{(m+1)n+m-2-2j}\right)}{\prod_{j=0}^{m}\left(1-\lambda
q^{(m-1)n+m-2j}\right)}.$ (12)

Under $q\rightarrow q^{-1}$ and $\lambda\rightarrow q^{mn}\lambda$ (the latter
making the $m\leftrightarrow n$ symmetry of torus knots more explicit), this
indeed reproduces the result of [5], as claimed in the introduction.

###  2.2 Twisted hyperbolic knots

We repeated the above calculation for some families of ‘twisted’888Note the
difference with the standard knot theory jargon, in which _twist knots_ all
have unkotting number $1$ (as e.g. the families $\overline{2k}\;\overline{2}$,
$\overline{2k+1}\;\overline{2}$); and they don’t restrict to hyperbolic as
they also include 2-stranded torus knots $T(2,n)$. hyperbolic knots, obtained
as follows. Given a projection of a simple knot, which can be thought of as a
_generating knot_ , we choose a point where two strands are parallel
(adjacent) to each other and cut it open there, as in fig. 2(a)–2(d).
Introducing a number of whole twists at the place indicated with 3 dots in the
figure, and then reconnecting the strands, yields an infinite family of knots
with only even or odd number of crossings. The reason we avoid half twists is
because they sometimes result into more than one component, i.e. a link, which
we would like to omit in the current treatment. The twisted hyperbolic
families are labelled using Conway notation999If the reader is unfamiliar with
Conway notation, we refer to Chapter 2 of [12] for a concise and comprehensive
introduction., in which the juxtaposed numbers indicate the number of
crossings of the individual tangles used to compose the knot, while an over-
line (instead of the standard minus sign) is used to denote negative tangles.
The total number of crossings $n$ of the knot is equal to the sum of its
Conway numbers. Successive members of a family correspond to increasing
$k=1,2,3,...$, each having $k-1$ additional whole twists, which can be thought
of as $2k-1$ extra ‘bubbles’. The generating knot, corresponding to $k=1$, is
chosen in a way such that the bubbles consist of positive crossings (c.f.
$L_{+}$ in (6)) which amounts to sometimes using the mirror of the knots
listed in the Rolfsen table [12] (while we shall not be careful to explicitly
mention this whenever using Rolfsen notation, it should be clarified by the
respective Conway notation). Four such families, shown in fig. 2, along with
the obtained results for their HOMFLY polynomials and HZ transforms, are given
below. Some further examples are included in the Appendix.

(a) $\overline{2k}\;\overline{2}$

(b) $\overline{2k+1}\;\overline{2}$

(c) $\overline{2k+1}\;\overline{1}\;\overline{2}$

(d) $(2k+2)\;3$

Figure 2: Some families of twisted hyperbolic knots.

(a) The family $\overline{2k}\;\overline{2}$ is generated by the figure-8 knot
$\overline{2}\;\overline{2}$, or $4_{1}$, and includes the knots
$6_{1},\;8_{1},\;10_{1}$, in Rolfsen notation. Its (unnormalised) HOMFLY
polynomial is
$\bar{H}_{\overline{2k}\;\overline{2}}(v,z)=\frac{v-v^{-1}}{z}(v^{2k}(1-v^{-2})+v^{-2}-z^{2}\sum_{j=0}^{k-1}v^{2j}$),
while its Harer–Zagier transform in terms of the total number of crossings
$n=2+2k$ is

$Z_{\overline{2k}\;\overline{2}}(q,\lambda)=\frac{\lambda\left(1+\lambda
q^{-5}\right)\left(1-\lambda^{2}q^{3n-8}\right)}{\left(1-\lambda
q^{-1}\right)\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{n-5}\right)\left(1-\lambda q^{n-3}\right)\left(1-\lambda q^{n-1}\right)}\\\
-\frac{\lambda^{2}q^{n-3}\left(\left(q^{2}+1+q^{-2}\right)\left(1-\lambda
q^{n-7}\right)+q^{-3}\left(q^{n-3}+q^{-n+3}\right)\left(1-\lambda
q^{n-1}\right)-q^{n}\left(1-\lambda q^{-n-7}\right)\right)}{\left(1-\lambda
q^{-1}\right)\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{n-5}\right)\left(1-\lambda q^{n-3}\right)\left(1-\lambda q^{n-1}\right)}$

(b) The family $\overline{2k+1}\;\overline{2}$ is generated by
$\overline{3}\;\overline{2}$ or $5_{2}$; includes the knots $\;7_{2}$ and
$9_{2}$;
$\bar{H}_{\overline{2k+1}\;\overline{2}}(v,z)=\frac{v-v^{-1}}{z}(v^{2(k+1)}(1-v^{2})+v^{2}+z^{2}\sum_{j=1}^{k+1}v^{2j})$;
$n=3+2k$

$Z_{\overline{2k+1}\;\overline{2}}(q,\lambda)=\frac{\lambda\left(\left(1-\lambda
q^{3}\right)\left(1-\lambda q^{n-2}\right)\left(1-\lambda
q^{2n+3}\right)-\lambda\left(1-\lambda
q^{n+2}\right)\left(q^{n}-q^{5}\right)\left(1-q^{n-3}\right)\right)}{\left(1-\lambda
q\right)\left(1-\lambda q^{3}\right)\left(1-\lambda
q^{n-2}\right)\left(1-\lambda q^{n}\right)\left(1-\lambda q^{n+2}\right)}$

(c) The family $\overline{2k+1}\;\overline{1}\;\overline{2}$ is generated by
$\overline{3}\;\overline{1}\;\overline{2}$ or $6_{2}$; includes $\;8_{2}$ and
$10_{2}$;
$\bar{H}_{\overline{2k+1}\;\overline{1}\;\overline{2}}(v,z)=v^{-2}\bar{H}_{T(2,2k+1)}(v,z)-zv^{-1}\bar{H}_{T(2,2k+2)}(v,z)$;
$n=4+2k$

$Z_{\overline{2k+1}\;\overline{1}\;\overline{2}}(q,\lambda)=\frac{\lambda\left(\left(1+\lambda
q^{n-9}\right)\left(1+\lambda q^{3n-7}\right)-\lambda
q^{2n-8}\left(q^{-2}+q^{2}\right)\left(q^{-n+3}+q^{n-3}\right)\right)}{\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{n-5}\right)\left(1-\lambda
q^{n-3}\right)\left(1-\lambda q^{n-1}\right)}$

(d) The family $(2k+2)\;3$ is generated by $7_{3}$, but can also be thought of
as being generated by the $2\;3$ projection of $5_{2}$, corresponding to
$k=0$; includes $9_{3}$;
$\bar{H}_{(2k+2)\;3}(v,z)=v^{2}\bar{H}_{T(2,2k+3)}(v,z)+zv\bar{H}_{T(2,2k+2)}(v,z)$;
$n=5+2k$

$Z_{(2k+2)\;3}(q,\lambda)=\frac{\lambda\left(\left(1-\lambda
q^{n-2}\right)\left(1-\lambda q^{3n-2}\right)-\lambda
q^{2n-2}\left(q-q^{-1}\right)\left(q^{n-5}-q^{-n+5}\right)\right)}{\left(1-\lambda
q^{n-4}\right)\left(1-\lambda q^{n-2}\right)\left(1-\lambda
q^{n}\right)\left(1-\lambda q^{n+2}\right)}$

We deduce that the HZ transforms for these families of twisted hyperbolic
knots still have completely factorised denominators but the numerators now
consist of sums of two or more factorised terms. It is worth pointing out the
similarity of the recursive formulas in the latter two cases with the one for
2-strand torus knots (7), hence resulting in almost factorised HZ functions.
The only exception among these families is the case $5_{2}$ (i.e.
$\overline{3}\;\overline{2}$ or $2\;3$), which has a completely factorised HZ
function

$Z_{5_{2}}(q,\lambda)=\frac{\lambda\left(1-\lambda
q^{13}\right)}{\left(1-\lambda q\right)\left(1-\lambda
q^{5}\right)\left(1-\lambda q^{7}\right)}.$ (13)

Beyond these families, we computed the HZ transform for the HOMFLY polynomial
of all knots in the Rolfsen table with up to 8 crossings101010We have also
considered composite knots $\mathcal{K}_{1}\\#\mathcal{K}_{2}$, for which
$H_{\mathcal{K}_{1}\\#\mathcal{K}_{2}}=H_{\mathcal{K}_{1}}H_{\mathcal{K}_{2}}$,
but they seem to not have a factorised HZ transform even when
$\mathcal{K}_{1,2}$ are both torus knots.. Among them we found that, apart
from $5_{2}$, $8_{20}$ also has an HZ transform with a completely factorised
form. Subsequently, we realised that there is a whole family of twisted
hyperbolic knots generated by a 6 crossing projection of $5_{2}$.

Figure 3: $P(\overline{2},3,\overline{2k+1})$

These are the Pretzel knots $P(\overline{2},3,\overline{2k+1})$, shown in fig.
3, in which $5_{2}$ corresponds to $k=0$, while it includes the knots $8_{20}$
at $k=1$ and $10_{125}$ at $k=2$. Their HOMFLY polynomial and the
corresponding HZ transforms are

$\bar{H}_{P(\overline{2},3,\overline{2k+1})}=v^{-2}\bar{H}_{P(\overline{2},3,\overline{2k-1})}+z^{2}\sum_{j=1}^{k}v^{-2j}\bar{H}_{P(\overline{2},3,\overline{2(k-j)+1})}-v^{-2k}(1-v^{2}+z^{2})\bar{H}_{T(2,3)},$
$Z_{P(\overline{2},3,\overline{2k+1})}(q,\lambda)=\frac{\lambda\left(1-\lambda
q^{13-2k}\right)\left(1-\lambda q^{3(1-2k)}\right)}{\left(1-\lambda
q^{1-2k}\right)\left(1-\lambda q^{3-2k}\right)\left(1-\lambda
q^{5-2k}\right)\left(1-\lambda q^{7-2k}\right)}$ (14)

which agrees with eq. (13) at $k=0$. From this expression it is clear that the
family $P(\overline{2},3,\overline{2k+1})$ satisfies the property (5) and
hence it might be possible to derive an explicit measure for the average
$\langle...\rangle_{P(\overline{2},3,\overline{2k+1})}$, which would give the
first working definition of a knot matrix model for hyperbolic knots. This
will be the subject of future investigation.

##  3 Analysis of HZ functions

A few remarks about the results listed in the previous section are in order.
Remark 1 At $q=1$, all HZ formulas reduce to
$Z_{\mathcal{K}}(1,\lambda)=\frac{\lambda}{(1-\lambda)^{2}}$. Moreover, in the
limits $q\rightarrow\infty$ and $q\rightarrow 0$, only the formulas
$Z_{\overline{2k+1}\;\overline{2}}$, $Z_{(2k+2)\;3}$ (corresponding to knots
with odd number of crossings $n$) and
$Z_{\overline{2k+1}\;\overline{1}\;\overline{2}}$ for $k\geq 3$ (or $n\geq
10$) have a finite values, equal to $1/\lambda$ and $\lambda$, respectively.
These coincide with the hyperbolic families with a non-factorised HZ transform
that have no zeros on the negative real axis (c.f. sec. 3.3 below).
Remark 2 If $\lambda$ is set to $q$ or $q^{-1}$, the HZ transform of some
twisted hyperbolic knots becomes factorised. Examples are
$Z_{\overline{5}\;\overline{2}}(q,\lambda=q)=q(1-q^{16})/((1-q^{2})(1-q^{6})(1-q^{10}))$,
$Z_{4\;3}(q,\lambda=q^{-1})=-(1+q^{10})/(q(1-q^{2})(1-q^{6}))$ and
$Z_{6\;3}(q,\lambda=q)=-q(1+q^{14})/((1-q^{6})(1-q^{10}))$.
Remark 3 All of the above formulas are invariant under $q\mapsto q^{-1}$ and
$\lambda\mapsto\lambda^{-1}$, i.e.
$Z_{\mathcal{K}}(q,\lambda)=Z_{\mathcal{K}}(q^{-1},\lambda^{-1})$, while the
modular transformations $q\mapsto-q^{-1}$ and $\lambda\mapsto-\lambda^{-1}$
yield $Z_{\mathcal{K}}\mapsto-Z_{\mathcal{K}}$.

###  3.1 Expansion for $q$ close to $1$

The limit $q\rightarrow 1$ is equivalent to the limit of large $k$, which is
referred to as the weak coupling limit in the physics literature [2]. In this
regime, fixing $\lambda=1$, we can set $q=e^{x}$ for $|x|\ll 1$ and expand the
HZ formulas in powers of $x$. The expansions always take the form
$Z_{\mathcal{K}}(e^{x},1)=\sum_{i=-1}^{\infty}a_{2i}^{\mathcal{K}}x^{2i}$,
where $a_{2i}^{\mathcal{K}}\in\mathbb{Q}$ have denominators that are multiples
of a fixed odd number. We have explicitly computed $a_{-2}^{\mathcal{K}}$ for
the above twisted families of knots

$a_{-2}^{\overline{2k}\;\overline{2}}=-\frac{1}{3}+\frac{4}{(n-5)(n-3)(n-1)},\;\;a_{-2}^{\overline{2k+1}\;\overline{2}}=\frac{1}{3}+\frac{4}{(n-2)n(n+2)},$
$a_{-2}^{\overline{2k+1}\;\overline{1}\;\overline{2}}=\frac{3\left(n^{2}-6n+13\right)}{(n-7)(n-5)(n-3)(n-1)},\;\;a_{-2}^{(2k+2)\;3}=\frac{3\left(n^{2}-4n+8\right)}{(n-4)(n-2)n(n+2)},$
$a_{-2}^{P(\overline{2},3,\overline{2k+1})}=\frac{3(n-19)}{(n-13)(n-11)(n-9)}.$

###  3.2 Poles and holomorphicity

As can be easily seen from the fully factorised denominators of the HZ
formulas their $\lambda$ poles lie at positive and negative powers of $q$,
hence they lie on the unit circle (recall $q=e^{\pi i/(k+N)})$, while there is
an additional pole at $\lambda=\infty$. The sum of the residues over all the
finite $\lambda$ poles equals $1$. In fact, it is interesting to note that
this can be deduced by considering just the first (factorised) part of the HZ
formulas, as the sum of the residues of the $\lambda^{2}$ term always
vanishes. Finally, adding the residue at the pole at infinity, which always
equals $-1$, the total sum becomes $0$. Moreover, at fixed $\lambda=1$, the
$q$-poles of the HZ formulas lie at $0,\;1,\infty$ and at roots of unity.
Again the sum of all the residues, including infinity, equals $0$. Via Cauchy
theorem, this implies that the HZ formulas are holomorphic in the extended
complex $\lambda$ and $q$ planes.

###  3.3 Zero locus

It is also of interest to consider the zeros of the above derived HZ formulas.
In the figures 4–8 below we plot the vanishing sets
$\\{q\in\mathbb{C}|Z_{\mathcal{K}}(q,1)=0\\}$ with fixed $\lambda=1$, for a
few examples of both torus and twisted hyperbolic knots.

(a) $T(2,5)$

(b) $T(5,6)$

(c) $T(15,17)$

Figure 4: For torus knots $T(m,n)$ all zeros have norm equal to $1$, i.e. they
lie on the unit circle. As $(m,n)$ increase these become more dense, but none
seems to lie on the real axis.

(a) $P(\overline{2},3,\overline{1})$ or $5_{2}$

(b) $P(\overline{2},3,\overline{5})$ or $10_{125}$

(c) $P(\overline{2},3,\overline{85})$

Figure 5: For the Pretzel family $P(\overline{2},3,\overline{2k+1})$ all zeros
have norm equal to $1$, i.e. they lie on the unit circle. For $k=1,2$ their
density descreases but for $k\geq 3$ it increases. Again, none of the zeros
lies on the real axis.

(a) $\overline{2}\;\overline{2}$ or $4_{1}$

(b) $\overline{10}\;\overline{2}$

(c) $\overline{58}\;\overline{2}$

Figure 6: For hyperbolic knots in the $\overline{2k}\;\overline{2}$ family,
most zeros still lie on the unit circle, with exception a finite number of
pairs which, when multiplied have norm equal to $1$, i.e. they are related by
conformal inversion. One such pair lies on the negative real axis
$\forall\;n$.

(a) $\overline{5}\;\overline{2}$ or $7_{2}$

(b) $\overline{9}\;\overline{2}$

(c) $\overline{47}\;\overline{2}$

Figure 7: For hyperbolic knots in the $\overline{2k+1}\;\overline{2}$ family,
the zero structure for $n=3+2k\geq 7$ shows again a deviation from the unit
circle only in conformal pairs. As for torus knots, none of them lies on the
real axis.

(a) $\overline{3}\;\overline{1}\;\overline{2}$ or $6_{2}$

(b) $\overline{5}\;\overline{1}\;\overline{2}$ or $8_{2}$

(c) $\overline{7}\;\overline{1}\;\overline{2}$ or $10_{2}$

Figure 8: Zero locus for some hyperbolic knots in the
$\overline{2k+1}\;\overline{1}\;\overline{2}$ family. For $n=6,8$ there is a
pair of zeros on the real axis, as in the $\overline{2k}\;\overline{2}$ family
(c.f. fig. 6), but this stops being the case $\forall\;n\geq 10$.

From these plots we deduce that when the HZ formulas are factorised, as it is
the case for $P(\overline{2},3,\overline{2k+1})$ and torus knots, the zeros
have unit norm, i.e. they lie on the unit circle. When the HZ transform
consists of sums of factorised terms, as for the majority of the twisted
hyperbolic knots considered, deviations from the circle arise in conformal
pairs, i.e. there are zeros of the form $ae^{i\phi}$ and
$\frac{1}{a}e^{i\phi}$ with $|a|\neq 1$. The plots for the $(2k+2)\;3$ family
have similar traits as the ones in fig. 7 for $\overline{2k+1}\;\overline{2}$,
and hence are omitted. The resemblance of these results with the zeros of the
characteristic function for the exponents of a singular complex curve studied
in [13] is astounding. Moreover, there might be a relation of these zero
structures to the zeros of the Riemann-$\zeta$ function [14]. In fact, it is
remarkable that such zero structures appear in various areas of mathematics,
but they lack a systematic study. Hence, a more in depth exploration of these
connections deserves to be the subject of future research.

#### Acknowledgements

A.P. is grateful to Roland van der Veen for his feedback and suggestions at
the early stages of this project. We would also like to thank Reiko Toriumi
for participating in our discussions. This work is supported by JSPS Kakenhi
19H01813.

## References

  * [1] Lord Kelvin.  On vortex atoms.  In Proc. R. Soc. Edin, volume 6, pages 94–105, 1867.
  * [2] Edward Witten.  Quantum field theory and the jones polynomial.  Communications in Mathematical Physics, 121(3):351–399, 1989.
  * [3] J. M. F. Labastida and M. Mariño.  The homfly polynomial for torus links from chern–simons gauge theory.  International Journal of Modern Physics A, 10(07), 1995.
  * [4] Miguel Tierz.  Soft matrix models and Chern–Simons partition functions.  Modern Physics Letters A, 19(18):1365–1378, 2004.
  * [5] A. Morozov, A. Popolitov, and Sh Shakirov.  Harer–Zagier formulas for knot matrix models.  Physics Letters B, 818:136370, 2021.
  * [6] H. Itoyama, A. Mironov, A. Morozov, and A. Morozov.  Character expansion for HOMFLY polynomials III: All 3-strand braids in the first symmetric representation.  International Journal of Modern Physics A, 27(19):1250099, jul 2012\.
  * [7] Andrei Mironov and Alexei Morozov.  Superintegrability of kontsevich matrix model.  The European Physical Journal C, 81(3), mar 2021.
  * [8] P. Dunin-Barkowski, A. Mironov, A. Morozov, A. Sleptsov, and A. Smirnov.  Superpolynomials for torus knots from evolution induced by cut-and-join operators.  Journal of High Energy Physics, 2013(3):1–87, 2013.
  * [9] Andrea Brini, Marcos Mariño, and Bertrand Eynard.  Torus knots and mirror symmetry.  Annales Henri Poincaré, 13(8):1873–1910, mar 2012.
  * [10] Kunio Murasugi.  Knot theory and its applications.  Springer, 1996.
  * [11] Georgios Giasemidis and Miguel Tierz.  Torus knot polynomials and susy wilson loops.  Letters in Mathematical Physics, 104(12):1535–1556, 2014.
  * [12] Colin C Adams.  The knot book.  American Mathematical Soc., 1994.
  * [13] Kyoji Saito.  The zeroes of characteristic function $\chi_{f}$ for the exponents of a hypersurface isolated singular point.  Algebraic varieties and analytic varieties (S. Iitaka, ed.), Advanced Studies in Pure Math., 1, Kino- kuniya, Tokyo and North-Holland, Amsterdam, New York, Oxford, pages 195–217, 1983.
  * [14] Edouard Brézin and Shinobu Hikami.  Random matrix theory with an external source.  SpringerBriefs in Mathematical Physics Vol. 19, Springer, 2016.

## Appendix A

Here we include the computed HOMFLY polynomials and their HZ transforms for
some further families of twisted hyperbolic knots. The coefficient of the
poles in the $q\rightarrow 1$ expansion of sec. 3.1 are also given.

(a) $2k\;1\;1\;2$

(b) $2\;(2k-1)\;1\;2$

(c) $4\;(2k+2)$

(d) $(2k+2)\;1\;3$

(a)
$\bar{H}_{2k\;1\;1\;2}(v,z)=v^{-2}(1+z^{2})\bar{H}_{T(2,2k+1)}(v,z)-zv^{-3}\bar{H}_{T(2,2k)}(v,z)$;
includes $6_{3},8_{7},10_{5}$; $n=4+2k$

$Z_{2k\;1\;1\;2}(q,\lambda)=\frac{\lambda\left(\left(1-\lambda
q^{n-13}\right)\left(1-\lambda q^{3n-11}\right)-\lambda
q^{2n-12}\left(q-q^{-1}\right)\left(q^{2}+q^{-2}\right)\left(q^{-n+4}-q^{n-4}\right)\right)}{\left(1-\lambda
q^{n-9}\right)\left(1-\lambda q^{n-7}\right)\left(1-\lambda
q^{n-5}\right)\left(1-\lambda q^{n-3}\right)}$
$a_{-2}^{2k\;1\;1\;2}=\frac{3\left(n^{2}-14n+37\right)}{(n-9)(n-7)(n-5)(n-3)}$

(b)
$\bar{H}_{2\;(2k-1)\;1\;2}(v,z)=v^{-2}(1+z^{2})\bar{H}_{\overline{2k-1}\;\overline{2}}(v,z)-zv^{-3}\bar{H}_{T(2,2)}(v,z)$;
includes $6_{3},8_{8},10_{34}$; $n=4+2k$

$Z_{2\;(2k-1)\;1\;2}(q,\lambda)=\lambda\\{\left(1-\lambda
q^{-7}\right)\left(1-\lambda q\right)\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{2n-5}\right)\\\ +\lambda
q^{n-2}\left(1-\lambda
q^{-7}\right)\left(\left(q^{-n+5}+q^{n-5}\right)\left(1-\lambda
q^{n-7}\right)+q^{-1}\left(1+\lambda
q^{3}\right)\left(-1+q^{n-8}\right)\right)\\\ +\lambda q^{-5}\left(1-\lambda
q\right)\left(\left(1-\lambda
q^{3n-13}\right)\left(1-q^{2}\right)^{2}-q^{n}\left(1+\lambda
q\right)\left(1-q^{n-10}\right)\right)\\\ -\lambda
q^{2n-11}\left(q^{2}+q^{-2}\right)\left(1-\lambda q\right)\left(1-\lambda
q^{-n+3}\right)\\}\\\ \times\left(\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{-1}\right)\left(1-\lambda q\right)\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{n-5}\right)\left(1-\lambda
q^{n-3}\right)\right)^{-1}$
$a_{-2}^{2\;(2k-1)\;1\;2}=-\frac{7}{3}+\frac{4}{(n-7)(n-5)(n-3)}$

(c)
$\bar{H}_{{4\;(2k+2)}}(v,z)=v^{-2}\bar{H}_{\overline{2(k+1)}\;\overline{2}}(v,z)-zv^{2k-1}\bar{H}_{T(2,2)}(v,z)-z(v-v^{-1})\sum_{j=0}^{k-1}v^{2j}$;
includes $8_{3}$, $10_{3}$; $n=6+2k$

$Z_{4\;(2k+2)}(q_{,}\lambda)=\lambda\\{\left(1+\lambda
q^{-9}\right)\left(1-\lambda q^{n-7}\right)\left(1+\lambda
q^{2n-7}\right)\left(1+\lambda^{2}q^{n-10}\right)\\\ -\lambda
q^{n-7}\left(1+\lambda q^{-9}\right)\left(q^{2}\left(1-\lambda
q^{-5}\right)\left(1+\lambda
q^{2n-9}\right)+q^{-2}\left(1-q^{n-2}\right)\left(1+\lambda^{2}q^{n-4}\right)\right)\\\
-\lambda q^{n-3}\left(1-\lambda
q^{-9}\right)\left(1+q^{n-10}\right)\left(1+\lambda^{2}q^{n-8}\right)\\\
-\lambda q^{-3}\left(2\left(1+\lambda
q^{-7}\right)\left(1-\lambda^{2}q^{4n-20}\right)+q^{2n-14}\left(1-\lambda
q^{-1}\right)^{2}\left(1-\lambda q^{3}\right)\right)\\\
+2\lambda^{2}q^{-8}\left(\left(1-\lambda
q^{n-3}\right)\left(1+q^{3n-14}\right)+q^{n+1}\left(q+q^{-1}\right)\left(1-\lambda
q^{2n-19}\right)\right)\\}\\\ \times\left(\left(1-\lambda
q^{-5}\right)\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{-1}\right)\left(1-\lambda q^{n-9}\right)\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{n-5}\right)\left(1-\lambda
q^{n-3}\right)\right)^{-1}$
$a_{-2}^{4\;(2k+2)}=\frac{1}{15}-\frac{8(n-6)}{(n-9)(n-7)(n-5)(n-3)}$

(d)
$\bar{H}_{(2k+2)\;1\;3}(v,z)=v^{-2}(1+z^{2})\bar{H}_{\overline{2(k+1)}\;\overline{2}}(v,z)-zv^{2k-3}\bar{H}_{T(2,2)}(v,z)-z(v-v^{-1})\sum_{j=-1}^{k-2}v^{2j}$;
includes $8_{4}$, $10_{4}$; $n=6+2k$

$Z_{(2k+2)\;1\;3}(q_{,}\lambda)=\frac{\lambda(\left(1+\lambda
q^{-11}\right)\left(1-\lambda q^{n-5}\right)\left(1+\lambda
q^{2n-13}\right))}{\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{-5}\right)\left(1-\lambda q^{n-9}\right)\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{n-5}\right)}\\\
+\frac{\lambda^{2}q^{-7}\left(\left(1-q^{n}\right)\left(1-q^{n-2}\right)(1-\lambda
q^{n-13})-q^{n-5}\left(q^{-2}+q^{2}\right)\left(q^{n-5}+q^{-n+5}\right)\left(1-\lambda
q^{n-5}\right)\right)}{\left(1-\lambda q^{-3}\right)\left(1-\lambda
q^{-5}\right)\left(1-\lambda q^{n-9}\right)\left(1-\lambda
q^{n-7}\right)\left(1-\lambda q^{n-5}\right)}$
$a_{-2}^{(2k+2)\;1\;3}=\frac{1}{15}-\frac{8}{(n-9)(n-7)(n-5)}$