# On First Integrals of Hamiltonian System with Holonomic Hamiltonian

Tomoyuki Iori This work was partly supported by JSPS KAKENHI Grant Number
JP21K21285.T. Iori is with the Department of Information and Physical
Sciences, Graduate School of Information Science and Technology, Osaka
University, 1-5 Yamadaoka, Suita, Osaka 565–0871, Japan
<EMAIL_ADDRESS>work has been submitted to the IEEE for possible
publication. Copyright may be transferred without notice, after which this
version may no longer be accessible.

###### Abstract

In this study, the solution of the Hamilton-Jacobi equation (HJE) with
holonomic Hamiltonian is investigated in terms of the first integrals of the
corresponding Hamiltonian system. Holonomic functions are related to a
specific type of partial differential equations called Pfaffian systems, whose
solution space can be regarded as a finite-dimensional real vector space. In
the finite-dimensional solution space, the existence of first integrals that
define a solution of the HJE is characterized by a finite number of algebraic
equations for finite-dimensional vectors, which can be easily solved and
verified. The derived characterization was illustrated through a numerical
example.

##  I Introduction

The Hamilton-Jacobi equation (HJE) is one of the most fundamental equations in
the analysis and control of nonlinear systems. It appears in various problems,
such as optimal control [1, 2], $H_{\infty}$ control [3], and balanced
realization [4] of nonlinear systems. For linear time-invariant systems, the
HJE is reduced to the algebraic Ricatti equation, a set of algebraic equations
for an unknown matrix, which can be solved using various solution methods.
However, for nonlinear systems, the HJE is formulated as a nonlinear partial
differential equation (PDE), which is difficult to solve analytically.
Numerical solution methods for the HJE have been studied based on various
mathematical techniques, such as series expansion [5], expansion with basis
functions [6], and data-driven approximation [7].

Hao et al. [8] proposed a numerical solution method for the HJE by generating
functions that consider the initial and terminal boundary conditions of
optimal control problems. The Taylor-series expansion of a generating function
up to any prescribed order is computed by solving first-order differential
equations. The approximated generating function is used to find a family of
optimal control for different boundary conditions. However, in these numerical
approaches, the computational complexity rapidly increases with the state
dimension, known as the curse of dimensionality.

The differential-geometric approach is an effective strategy for analyzing the
HJE [9]. From the perspective of symplectic geometry, the solutions of the HJE
are identified with Lagrangian manifolds in the cotangent bundle of the state
space. A Lagrangian manifold is characterized by a finite number of first
integrals of the corresponding Hamiltonian system as a set of their common
zeros. Hence, solving the HJE is reformulated as the problem of finding the
appropriate first integrals of the corresponding Hamiltonian system. This
characterization reduces the HJE to a set of ordinary differential equations
called the Lagrange-Charpit system, which is still difficult to solve
analytically.

Moreover, algebraic approaches based on commutative algebra are promising
techniques, especially for avoiding the curse of dimensionality. For time-
invariant Hamiltonian, Ohtsuka [10] characterized the solutions of the HJE
with algebraic gradients in terms of the existence of involutive zero-
dimensional ideals in a polynomial ring; Kawano and Ohtsuka [11] extended
these results to the case of time-varying Hamiltonian. It should be noted that
the first integrals that define a Lagrangian manifold are closely related to
involutive zero-dimensional ideals; Indeed, the existence of first integrals
in the polynomial ring implies the existence of an involutive zero-dimensional
ideal.

Considering these differential-geometric and algebraic approaches, in this
study, the HJE was analyzed from the perspective of the theory of $D$-modules,
which is a field of mathematics that studies the algebraic structures of PDEs
and their solution space [12]. In the theory of $D$-modules, the solution set
of PDEs is characterized by an ideal in the noncommutative ring of
differential operators, which can be computed using symbolic computation. The
symbolic computation of differential operators has been intensively studied
[13, 15, 16] and its applications can be found in statistics [17, 18] and
moment problems [19]. Moreover, for specific ideals, called zero-dimensional
ideals, the solutions are called holonomic functions and constitute a finite-
dimensional real vector space. This finite dimensional property is promising
for finding the first integrals because any first integral in the solution
space can be determined by a finite number of parameters.

In this paper, the first integrals of the Hamiltonian system are investigated
through the symbolic computation of differential operators. First, it is
assumed that all first integrals, including the Hamiltonian, that define a
Lagrangian manifold are included in the solution space of a zero-dimensional
ideal. This assumption allows us to associate the first integrals with vector-
valued functions that satisfy a specific type of first-order PDEs, called a
Pfaffian system. Moreover, any solution of the Pfaffian system is
characterized by only the boundary value at a single point where the Pfaffian
system is defined. This leads to the characterization of first integrals by a
finite number of algebraic equations for the boundary value. Through this
characterization, the first integrals can be determined by solving algebraic
equations, if such first integrals exist.

#### Notations

For the field of real numbers ${\boldsymbol{\mathrm{R}}}$ and a vector of
indeterminates $y=[y_{1}\ \cdots\ y_{n}]^{\top}$,
${\boldsymbol{\mathrm{R}}}(y)$ denotes the field of rational functions in
$y_{1},\dots,y_{n}$ over ${\boldsymbol{\mathrm{R}}}$.
$\partial_{y}\coloneqq[\partial_{y_{1}}\ \cdots\ \partial_{y_{n}}]^{\top}$
denotes a vector of differential operators, where
$\partial_{y_{i}}={\partial{}/\partial{y_{i}}}$. Here, $\partial_{y}$ and
$\partial_{y_{i}}$ are abbreviated by $\partial$ and $\partial_{i}$,
respectively, if $y$ is clearly specified according to the context. For a
multi-index
$\alpha=(\alpha_{1},\dots,\alpha_{n})\in{\boldsymbol{\mathrm{Z}}}^{n}_{\geq
0}$, $y^{\alpha}$ and $\partial^{\alpha}$ denote the monomial
$y_{1}^{\alpha_{1}}\cdots y_{n}^{\alpha_{n}}$ and differential operator
$\partial_{1}^{\alpha_{1}}\cdots\partial_{n}^{\alpha_{n}}$, respectively. Let
${\boldsymbol{\mathrm{R}}}[y]$ and ${\boldsymbol{\mathrm{R}}}(y)$ be the set
of all polynomials and rational functions in $y_{1},\dots,y_{n}$,
respectively. The set of all $(m\times l)$ matrices and $(m\times m)$ skew
matrices with components in ${\boldsymbol{\mathrm{R}}}(y)$ is denoted by
${\boldsymbol{\mathrm{R}}}(y)^{m\times l}$ and
$\mathrm{Skew}_{m}({\boldsymbol{\mathrm{R}}}(y))$, respectively. The symbol
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$ denotes the
noncommutative ring of differential operators with coefficients in
${\boldsymbol{\mathrm{R}}}(y)$. The action of a differential operator
$\mathcal{P}$ on a sufficiently smooth function $F(y)=F(y_{1},\dots,y_{n})$ is
denoted by $\mathcal{P}\bullet F(y)$; for instance, $\partial_{i}\bullet
F(y)={\partial{F}/\partial{y_{i}}}(y)$. If $\mathcal{P}\bullet F=0$, a
differential operator $\mathcal{P}$ is said to _annihilate_ function $F$ and
$F$ is a _solution_ of $\mathcal{P}$. The left ideal in
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$ generated by a finite
number of differential operators $\mathcal{P}_{1},\dots,\mathcal{P}_{s}$ is
defined as
$\langle\mathcal{P}_{1},\dots,\mathcal{P}_{s}\rangle\coloneqq\\{\mathcal{Q}_{1}\mathcal{P}_{1}+\cdots+\mathcal{Q}_{s}\mathcal{P}_{s}\mid\mathcal{Q}_{1},\dots,\mathcal{Q}_{s}\in{\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle\\}$.
In this paper, the adjective “left” is omitted and it is referred to as an
“ideal” because there is no right ideal. If $\mathcal{P}\bullet F=0$ for all
$\mathcal{P}\in I$, the ideal
$I\subset{\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$ is said to
_annihilate_ function $F$ and $F$ is a _solution_ of $I$.

##  II First Integrals of Hamiltonian Systems and Hamilton-Jacobi Equation

For a scalar-valued function $h(x,p)$, we consider the following first-order
PDE for scalar-valued function $v(x)$:

$h(x,p)=0,\quad p=\nabla_{x}v(x),$ (1)

where $x\in{\boldsymbol{\mathrm{R}}}^{n}$. The PDE (1) is called the HJE for
the Hamiltonian $h$. The Hamiltonian system associated with $h$ is a dynamical
system with $2n$ state variables $z=(x,p)$ defined as follows.

$\begin{cases}\dot{x}&=\nabla_{p}h(x,p)\\\
\dot{p}&=-\nabla_{x}h(x,p).\end{cases}$ (2)

A first integral of the Hamiltonian system (2) is a function that is constant
along the trajectory of (2), i.e., a function $f$ that satisfies

$\dot{f}=(\nabla_{x}f)^{\top}\dot{x}+(\nabla_{p}f)^{\top}\dot{p}=\\{h,f\\}=0,$

where $\\{\cdot,\cdot\\}$ denotes the Poisson bracket, defined as

$\\{f,g\\}=(\nabla_{p}f)^{\top}\nabla_{x}g-(\nabla_{x}f)^{\top}\nabla_{p}g.$

The first integrals of the Hamiltonian system are closely related to the
solutions of the HJE. Let $\bar{z}=(\bar{x},\bar{p})$ be a fixed point with
$h(\bar{z})=0$ and $\pi$ be the projection to the first component, i.e.,
$\pi\colon z=(x,p)\mapsto x$. For a neighborhood
$U\subset{\boldsymbol{\mathrm{R}}}^{2n}$ of $\bar{z}$, a smooth function
$v(x)$ defined on $\pi(U)$ is a solution of the HJE if and only if $h(z)=0$
holds in the following subset

$\Lambda_{v}\coloneqq\\{z=(x,p)\in U\mid
p=\nabla_{x}v(x)\\}\subset{\boldsymbol{\mathrm{R}}}^{2n}.$

By symplectic geometry, $\Lambda_{V}$ can be written as

$\Lambda_{v}=\\{z\in U\mid f_{1}(z)=f_{2}(z)=\cdots=f_{n}(z)=0\\},$

where $f_{1}=h$ and $f_{2},\dots,f_{n}$ are the first integrals of (2) such
that the following conditions hold (see [9] for details):

$\displaystyle\\{f_{k}(z),f_{l}(z)\\}=0\quad(z\in U,\ 1\leq k<l\leq n),$ (3)
$\displaystyle\det[\nabla_{p}f_{1}(\bar{z})\ \cdots\
\nabla_{p}f_{n}(\bar{z})]\neq 0.$ (4)

Let us consider the problem to find analytic functions $f_{2},\dots,f_{n}$ at
$\bar{z}$ that satisfy the system of nonlinear PDEs (3) as well as the
inequality condition (4), which is difficult to solve analytically. The set of
all analytic functions at $\bar{z}$ comprises an infinite dimensional vector
space; this is reflected in the series expansion approach and Galerkin method
because they usually require an infinite number of parameters to exactly
determine a solution of the HJE, which cannot be accomplished in finite time.

To overcome this difficulty, the notion of Pfaffian system associated with
holonomic functions is used. The finite-dimensional solution space of a
Pfaffian system allows us to exactly reduce the problem of finding analytic
functions $f_{2},\dots,f_{n}$ to that of determining finite dimensional
vectors, which can be solved easily without approximations.

##  III Holonomic Functions and Pfaffian Systems

In this section, some notions related to the symbolic computation of
differential operators are introduced, referring to [13, 20, 16, 21] for most
of the definitions and lemmas.

First, let us consider the notion of holonomic functions, where all first
integrals are assumed to be included. Holonomic functions can be defined as
the solution of zero-dimensional ideals.

######  Definition 1 ([20])

For an ideal $I\subset{\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$, let
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle/I$ be the quotient space
as ${\boldsymbol{\mathrm{R}}}(y)$-linear spaces. The ideal $I$ is said to be
zero-dimensional if ${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle/I$ is
finite-dimensional.

######  Definition 2 ([20])

A function $f(y)$ is said to be holonomic if it is a solution of a zero-
dimensional ideal in ${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$.

In this paper, $h(x,p)$ is called the holonomic Hamiltonian if it is holonomic
as a function of $2n$ variables $z=(x,p)$. The following lemma can be used to
verify whether a given function is holonomic.

###### Lemma 1

A function $f(y)$ with $y=[y_{1}\ \cdots\ y_{m}]^{\top}$ is holonomic if it is
a solution of differential operators $\mathcal{P}_{i}$:

$\mathcal{P}_{i}\bullet f(y)=0\quad(i=1,\dots,m),$ (5)

where each $\mathcal{P}_{i}$ is defined as a finite sum of the following form:

$\mathcal{P}_{i}=\sum_{k}c_{ik}(y)\partial_{i}^{k}\quad\left(c_{ik}\in{\boldsymbol{\mathrm{R}}}(y)\right).$
(6)

###### Proof:

The proof is readily completed because the ideal
$\langle\mathcal{P}_{1},\dots,\mathcal{P}_{m}\rangle$ is zero-dimensional (see
Theorem 6.1.11 in [20]). ∎

The class of holonomic functions is closed under addition, multiplication, and
differentiation.

######  Lemma 2 ([21])

Suppose $f(y)$ and $g(y)$ are holonomic functions. Then, $f(y)+g(y)$,
$f(y)g(y)$, and $\partial_{i}f(y)\ (i=1,\dots,m)$ are also holonomic.

The zero-dimensional ideals that annihilate the sums and products of holonomic
functions can be computed through the symbolic computation of differential
operator; however, it requires the notion of holonomic ideals in the Weyl
algebra (see [15] for details).

A specific system of PDEs, called the Pfaffian system, can be derived from the
finite dimensional quotient space
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle/I$.

######  Definition 3 ([13])

The Pfaffian system defined by
$A_{1},\dots,A_{m}\allowbreak\in{\boldsymbol{\mathrm{R}}}(y)^{s\times s}$ is a
system of PDEs for an $s$-dimensional vector-valued function $q(y)$:

$\partial_{i}q(y)=A_{i}(y)q(y),\quad(i=1,\dots,m).$ (7)

The relationship between holonomic functions, zero-dimensional ideals, and
Pfaffian systems can be summarized as follows.

######  Lemma 3 ([13])

Let $I\subset{\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$ be a zero-
dimensional ideal such that the quotient space
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle/I$ is $s$-dimensional. Let
$\mathcal{B}$ be a vector of differential operators defined as

$\mathcal{B}\coloneqq[1\ \partial^{\alpha_{1}}\ \cdots\
\partial^{\alpha_{s-1}}]^{\top},$ (8)

where $1,\partial^{\alpha_{1}},\dots,\partial_{s-1}$ are the standard
monomials [13] of $I$. Then,
$A_{1}(y),\dots,A_{m}(y)\in{\boldsymbol{\mathrm{R}}}(y)^{s\times s}$ in (7)
can be computed such that a vector-valued function $Q(y)=\mathcal{B}\bullet
f(y)$ satisfies the Pfaffian system, where $f$ denotes any solution of $I$.
Moreover, $A_{1}(y),\dots,A_{m}(y)$ satisfy the integrability condition:

$\partial_{j}A_{i}+A_{i}A_{j}=\partial_{i}A_{j}+A_{j}A_{i}\quad(1\leq i<j\leq
m)$ (9)

The key property of the Pfaffian system is that its solution space is finite-
dimensional over ${\boldsymbol{\mathrm{R}}}$ as the space of analytic
functions.

######  Lemma 4 ([13])

With the same notations as Lemma 3, let $U$ denote any simply connected domain
in $\\{y\in{\boldsymbol{\mathrm{R}}}^{m}\mid D(y)\neq 0\\}$, where $D(y)$ is
the least common multiple of all denominators included in $A_{1},\dots,A_{m}$.
Let $\mathbb{S}(I)$ be the set of all analytic solutions of $I$ defined on
$U$. Let $\mathbb{VS}(A_{1},\dots,A_{m})$ be the set of all analytic solutions
of the Pfaffian system (7) defined on $U$. Then, $\mathbb{S}(I)$,
$\mathbb{VS}(A_{1},\dots,A_{m})$, and ${\boldsymbol{\mathrm{R}}}^{s}$ are
isomorphic to each other with isomorphisms
$\phi_{\mathcal{B}}\colon\mathbb{S}(I)\to\mathbb{VS}(A_{1},\dots,A_{m})$ and
$\psi_{\bar{y}}\colon\mathbb{VS}(A_{1},\dots,A_{m})\to{\boldsymbol{\mathrm{R}}}^{s}$
defined by

$\phi_{\mathcal{B}}(y)\coloneqq\mathcal{B}\bullet
f,\quad\psi_{\bar{y}}(q)\coloneqq q(\bar{y}),$

with any fixed $\bar{y}\in U$. The inverse $\phi_{\mathcal{B}}^{-1}(q)$ can be
defined as the first component of $q$; $\psi_{\bar{y}}^{-1}(\bar{q})$ can be
defined as the unique solution of (7) with boundary condition
$q(\bar{y})=\bar{q}$.

When $I=\langle\mathcal{P}_{1},\dots,\mathcal{P}_{m}\rangle$, $\mathcal{B}$
and $A_{1},\dots,A_{m}$ can be computed from
$\mathcal{P}_{1},\dots,\mathcal{P}_{m}$ using the Gröbner bases in
${\boldsymbol{\mathrm{R}}}(y)\langle\partial\rangle$ [13]. In this paper, the
solution space $\mathbb{S}(I)$ is denoted by
$\mathbb{S}(\mathcal{P}_{1},\dots,\mathcal{P}_{m})$ if
$I=\langle\mathcal{P}_{1},\dots,\mathcal{P}_{m}\rangle$.

##  IV First Integrals in Solution Space of Pfaffian System

For differential operators $\mathcal{P}_{1},\dots,\mathcal{P}_{2n}$ of the
form of (6) that annihilate $h(x,p)$, suppose that other first integrals
$f_{2},\dots,f_{n}$ that define a solution of the HJE are included in the
solution space $\mathbb{S}(\mathcal{P}_{1},\dots,\mathcal{P}_{2n})$. Now, a
condition for the existence of such first integrals is derived by considering
the Pfaffian system associated with $\mathcal{P}_{1},\dots,\mathcal{P}_{2n}$.
This condition is obtained as a set of algebraic equations for a finite number
of finite-dimensional vectors, which can be easily solved. Consequently, the
first integrals can be determined by solving the algebraic equations if they
have a sufficient number of solutions.

The condition can be obtained in two steps. First, the set of PDEs (3) is
converted into an infinite number of algebraic equations for some finite-
dimensional vectors. It should be noted that this first part is presented in
the preliminary form in the conference proceedings [14]. Subsequently, the
number of equations is reduced to a finite number. In this section,
$\mathcal{B}$ and $A_{1},\dots,A_{2n}\in{\boldsymbol{\mathrm{R}}}(z)^{d\times
d}$ denote the vector of differential operators and matrices of rational
functions derived from $\mathcal{P}_{1},\dots,\mathcal{P}_{2n}$ through Lemma
3. Moreover, $\mathbb{S}$ and $\mathbb{VS}$ denote
$\mathbb{S}(\mathcal{P}_{1},\dots,\mathcal{P}_{2n})$ and
$\mathbb{VS}(A_{1},\dots,A_{2n})$, respectively. In addition, a point
$\bar{z}\in U$ is fixed, where the domain $U$ is defined as in Lemma 3.

###  IV-A Reduction from PDEs to Infinite Number of Algebraic Equations

For the first step, suppose that all first integrals $f_{1},\dots,f_{n}$ lie
in $\mathbb{S}$. From Lemma 3, there exists $q_{k}\in\mathbb{VS}$ such that
$q_{k}=\phi_{\mathcal{B}}(f_{k})$ for each $f_{k}$.

From (8), we obtain

$\begin{cases}\nabla_{x}f_{k}(z)=B_{x}(z)q_{k}(z),\\\
\nabla_{p}f_{k}(z)=B_{p}(z)q_{k}(z),\end{cases}$ (10)

where $B_{x}\in{\boldsymbol{\mathrm{R}}}(z)^{n\times d}$ and
$B_{p}\in{\boldsymbol{\mathrm{R}}}(z)^{n\times d}$ consist of the first rows
of $A_{1},\dots,A_{n}$ and those of $A_{n+1},\dots,A_{2n}$, respectively. By
substituting (10) into (3) and (4), we obtain

$\displaystyle\begin{aligned} &\\{f_{k},f_{l}\\}(z)=(q_{k}^{\top}\Omega
q_{l})(z)=0,\\\ &\qquad\qquad\qquad(z\in U,\ 1\leq k<l\leq n)\end{aligned}$
(11) $\displaystyle\det\left\\{B_{p}(\bar{z})\left[q_{1}(\bar{z})\ \cdots\
q_{n}(\bar{z})\right]\right\\}\neq 0$ (12)

for a point $\bar{z}\in\mathbb{U}$ and $\Omega\coloneqq
B_{p}^{\top}B_{x}-B_{x}^{\top}B_{p}\in\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))$.

From Lemma 4, $q_{k}^{\top}\Omega q_{l}$ is analytic at any $\bar{z}\in U$;
hence, $q_{k}^{\top}\Omega q_{l}=0$ over $U$ if and only if

$\left.\partial^{\alpha}\bullet\left(q_{k}^{\top}\Omega
q_{l}\right)(z)\right|_{z=\bar{z}}=0$ (13)

holds for any $\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq 0}^{2n}$. For each
differential operator $\partial_{i}\ (i=1,\dots,2n)$, the Pfaffian system
yields

$\displaystyle\partial_{i}\bullet\left(q_{k}^{\top}\Omega q_{l}\right)$
$\displaystyle=(\partial_{i}q_{k})^{\top}\Omega
q_{l}+q_{k}^{\top}(\partial_{i}\Omega)q_{l}+q_{k}^{\top}\Omega(\partial_{i}q_{l})$
$\displaystyle=q_{k}^{\top}\left(A_{i}^{\top}\Omega+\partial_{i}\Omega+\Omega
A_{i}\right)q_{l}.$

By defining a mapping
$\mathcal{D}_{i}\colon\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))\to\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))$
for each $i=1,\dots,2n$ as

$\mathcal{D}_{i}\Omega\coloneqq A_{i}^{\top}\Omega+\Omega
A_{i}+\partial_{i}\Omega,$ (14)

we obtain

$\partial^{\alpha}\bullet(q_{k}^{\top}\Omega
q_{l})=q_{k}^{\top}(\mathcal{D}^{\alpha}\Omega)q_{l}.$ (15)

By substituting (15) into (13), an infinite set of algebraic equations for
vectors $\bar{q}_{k}$ and $\bar{q}_{l}$ is obtained:

$\bar{q}_{k}^{\top}\left\\{\left(\mathcal{D}^{\alpha}\Omega\right)(\bar{z})\right\\}\bar{q}_{l}=0\quad(\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n}),$ (16)

where $\mathcal{D}^{\alpha}$ is the composition of mappings
$\mathcal{D}_{1}^{\alpha_{1}}\circ\mathcal{D}_{2}^{\alpha_{2}}\circ\cdots\circ\mathcal{D}_{2n}^{\alpha_{2n}}$.
Thus far, the discussion can be summarized as follows.

###### Theorem 1

Suppose the differential operators
$\mathcal{P}_{1},\dots,\allowbreak\mathcal{P}_{2n}$ are of the form of (6) and
annihilate the holonomic Hamiltonian $h$. For a fixed $\bar{z}\in U$ and
$\bar{q}_{1}=\psi_{\bar{z}}\circ\phi_{\mathcal{B}}(h)$, let
$\bar{q}_{2},\dots,\bar{q}_{n}\in{\boldsymbol{\mathrm{R}}}^{d}$ be the
constant vectors that satisfy (16) for $k,l\in\\{1,\dots,n\\}$. Then, the
functions
$f_{k}\coloneqq\phi_{\mathcal{B}}^{-1}\circ\psi_{\bar{z}}^{-1}(\bar{q}_{k})\in\mathbb{S}\
(k=2,\dots,n)$ satisfy conditions (3) and (4); in other words, they are the
first integrals that define a solution of the HJE (1).

Now, the condition (3) is converted into the set of algebraic equations (16)
for $n$ vectors
$\bar{q}_{1},\dots,\bar{q}_{n}\in{\boldsymbol{\mathrm{R}}}^{d}$. If
$\bar{q}_{2},\dots,\bar{q}_{n}$ satisfy (16) for $k,l\in\\{1,\dots,n\\}$ and
(12) with $\bar{q}_{1}$, the first integrals $f_{2},\dots,f_{n}$ that define a
solution of the HJE are obtained as
$f_{k}=\phi_{\mathcal{B}}^{-1}\circ\psi_{\bar{z}}^{-1}(\bar{q}_{k})\
(k=2,\dots,n)$. However, (16) consists of an infinite number of algebraic
equations, which cannot be solved or verified for some candidates of
$\bar{q}_{2},\dots,\bar{q}_{n}$ in finite time. The following section shows
that the number of equations can be reduced to be finite.

###  IV-B Reduction to Finite Set of Algebraic Equations

From the closure property of holonomic functions in Lemma 2, the Poisson
bracket $\\{f_{k},f_{l}\\}$ is a holonomic function of $z$ if $f_{k}$ and
$f_{l}$ are holonomic. Hence, a zero-dimensional ideal annihilating
$\\{f_{k},f_{l}\\}$ and the corresponding Pfaffian system exist. Furthermore,
it can be shown that there exists a Pfaffian system satisfied by all Poisson
brackets $\\{f,g\\}$ of any two functions $f,g\in\mathbb{S}$.

###### Lemma 5

A finite set of multi-indices
$\Gamma=\\{{\boldsymbol{\mathrm{0}}},\gamma_{1},\dots,\gamma_{t-1}\\}$ exists
such that for any pair of solutions $f,g\in\mathbb{S}$, a Pfaffian system

$\partial_{i}r(z)=T_{i}(z)r(z)\quad(i=1,\dots,2n)$ (17)

is satisfied by the vector-valued function $r=\mathcal{C}\bullet\\{f,g\\}$
with a vector of differential operators

$\mathcal{C}=[1\ \partial^{\gamma_{1}}\ \cdots\
\partial^{\gamma_{t-1}}]^{\top}.$ (18)

###### Proof:

First, the Pfaffian system (17) is derived from the infinite set of matrices
$\mathcal{D}^{\alpha}\Omega\ (\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n})$. For the ${\boldsymbol{\mathrm{R}}}(z)$-vector space
$\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))$, let
$\mathbb{M}\subseteq\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))$ be the
subspace spanned by all matrices $\mathcal{D}^{\alpha}\Omega$, i.e.,

$\mathbb{M}\coloneqq\mathrm{Span}_{{\boldsymbol{\mathrm{R}}}(z)}\\{\mathcal{D}^{\alpha}\Omega\mid\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n}\\}.$ (19)

The subspace $\mathbb{M}$ must be finite-dimensional because the dimension of
$\mathrm{Skew}_{d}({\boldsymbol{\mathrm{R}}}(z))$ is $(d(d-1)/2)$. Let
$t<(d(d-1)/2)$ be the dimension of $\mathbb{M}$ and the set of $t$ matrices
$\\{\mathcal{D}^{\gamma_{0}}\Omega,\dots,\mathcal{D}^{\gamma_{t-1}}\Omega\\}$
be a basis of $\mathbb{M}$. It should be noted that without loss of
generality, we can select
$\Gamma=\\{\gamma_{0},\gamma_{1},\dots,\gamma_{t-1}\\}$ such that
$\gamma_{0}={\boldsymbol{\mathrm{0}}}\in{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n}$ and
$\gamma_{0}\prec\gamma_{1}\prec\cdots\prec\gamma_{t-1}\prec\alpha$ for all
$\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq 0}^{2n}\setminus\Gamma$ with a
monomial order $\prec$ on ${\boldsymbol{\mathrm{Z}}}_{\geq 0}^{2n}$.

For any $k\in\\{0,\dots,t-1\\}$ and $i\in\\{1,\dots,2n\\}$, there exist
rational functions $c_{ikl}\in{\boldsymbol{\mathrm{R}}}(z)\ (l=0,\dots,t-1)$
such that the following holds:

$\mathcal{D}_{i}\mathcal{D}^{\gamma_{k}}\Omega=\sum_{l=0}^{t-1}c_{ikl}\mathcal{D}^{\gamma_{l}}\Omega.$
(20)

For every pair of functions $f,g\in\mathbb{S}$, (20) yields

$\displaystyle\partial_{i}\partial^{\gamma_{k}}\bullet\\{f,g\\}$
$\displaystyle=q_{f}^{\top}\\{\mathcal{D}_{i}\mathcal{D}^{\gamma_{k}}\Omega\\}q_{g}$
$\displaystyle=\sum_{l=0}^{t-1}c_{ikl}\left[q_{f}^{\top}\\{\mathcal{D}^{\gamma_{l}}\Omega\\}q_{g}\right]$
$\displaystyle=\sum_{l=0}^{t-1}c_{ikl}\left[\partial^{\gamma_{l}}\bullet\\{f,g\\}\right],$
(21)

where $q_{f}=\phi_{\mathcal{B}}(f)$ and $q_{g}=\phi_{\mathcal{B}}(g)$. By
combining (21) for $k=0,\dots,t-1$ into one, we obtain the Pfaffian system
(17) with matrices $T_{i}\ (i=1,\dots,2n)$ whose $(k,l)$-component is
$c_{ikl}$ and vector $\mathcal{C}$ defined by (18).

Next, the proof of matrices $T_{i}\ (i=1,\dots,2n)$ satisfying the
integrability condition (9) is presented. Let us consider (20). By applying
$\mathcal{D}_{j}$ with $j\neq i$ to both sides, we obtain

$\displaystyle\mathcal{D}_{j}\mathcal{D}_{i}\mathcal{D}^{\gamma_{k}}\Omega=$
$\displaystyle\mathcal{D}_{j}\left\\{\sum_{l=0}^{t-1}c_{ikl}\mathcal{D}^{\gamma_{l}}\Omega\right\\}$
$\displaystyle=$
$\displaystyle\sum_{l=0}^{t-1}\left\\{c_{ikl}\mathcal{D}_{j}\mathcal{D}^{\gamma_{l}}\Omega+(\partial_{j}c_{ikl})\mathcal{D}^{\gamma_{l}}\Omega\right\\}$
$\displaystyle=$
$\displaystyle\sum_{l=0}^{t-1}\left\\{c_{ikl}\sum_{s=0}^{t-1}c_{jls}\mathcal{D}^{\gamma_{s}}\Omega+(\partial_{j}c_{ikl})\mathcal{D}^{\gamma_{l}}\Omega\right\\}$
$\displaystyle=$
$\displaystyle\sum_{l=0}^{t-1}\left\\{\sum_{s=0}^{t-1}c_{iks}c_{jsl}+(\partial_{j}c_{ikl})\right\\}\mathcal{D}^{\gamma_{l}}\Omega.$
(22)

By replacing the subscript $i$ with $j$, we have

$\mathcal{D}_{i}\mathcal{D}_{j}\mathcal{D}^{\gamma_{k}}\Omega=\sum_{l=0}^{t-1}\left\\{\sum_{s=0}^{t-1}c_{jks}c_{isl}+(\partial_{i}c_{jkl})\right\\}\mathcal{D}^{\gamma_{l}}\Omega.$
(23)

The right-hand sides of (22) and (23) are equal to each other because
$\mathcal{D}_{i}\mathcal{D}_{j}\mathcal{D}^{\gamma_{k}}\Omega=\mathcal{D}_{j}\mathcal{D}_{i}\mathcal{D}^{\gamma_{k}}\Omega$
holds by considering (9) and (14). In particular, the coefficients of each
matrix are equal to each other because matrices
$\mathcal{D}^{\gamma_{0}}\Omega,\dots,\mathcal{D}^{\gamma_{t-1}}\Omega$ are
linearly independent, thereby yielding the integrability condition for
$T_{1},\dots,T_{2n}$:

$\partial_{j}T_{i}+T_{i}T_{j}=\partial_{i}T_{j}+T_{j}T_{i}\quad(1\leq i<j\leq
2n).$ (24)

This completes the proof. ∎

Thus, from Lemma 5, the main result is obtained as follows.

###### Theorem 2

With the same notations as Lemma 5, let $E(z)\in{\boldsymbol{\mathrm{R}}}[z]$
be the least common multiple of the denominators included in
$T_{1},\dots,T_{2n}$. For any fixed point $\bar{z}\in
U^{\prime}\coloneqq\\{z\in U\mid E(z)\neq 0\\}$, constant vectors
$\bar{q}_{k}\in{\boldsymbol{\mathrm{R}}}^{d}\ (k=1,\dots,n)$ satisfy (16) if
and only if they satisfy a finite number of algebraic equations:

$\bar{q}_{k}^{\top}\mathcal{D}^{\gamma}\Omega(\bar{z})\bar{q}_{l}=0\quad(\gamma\in\Gamma,\
1\leq k<l\leq n).$ (25)

###### Proof:

From Lemma 4 and the Pfaffian system (17), the Poisson bracket
$\\{f_{k},f_{l}\\}(z)$ is constantly equal to $0$ on $U^{\prime}$ if and only
if
$(\mathcal{C}\bullet\\{f_{k},f_{l}\\})(\bar{z})={\boldsymbol{\mathrm{0}}}\in{\boldsymbol{\mathrm{R}}}^{t}$
holds at any one point $\bar{z}\in U^{\prime}$. This implies that (16) holds
if and only if (25) is valid for
$\bar{q}_{k}=\psi_{\bar{z}}\circ\phi_{\mathcal{B}}(F_{k})\ (k=1,\dots,n)$. ∎

Theorem 2 ensures that we can determine $f_{2},\dots,f_{n}\in\mathbb{S}$ that
define a solution of the HJE with $h$ by finding
$\bar{q}_{2},\dots,\bar{q}_{n}\in{\boldsymbol{\mathrm{R}}}^{d}$ that satisfy
(25) and (12) with $\bar{q}_{1}=\psi_{\bar{z}}\circ\phi_{\mathcal{B}}(h)$.
More precisely, $f_{k}$ is specified as the analytic function that satisfies
$q_{k}=\mathcal{B}\bullet f_{k}$ for the solution $q_{k}(z)$ of the Pfaffian
system (7) defined by $A_{1},\dots,A_{2n}$ with boundary condition
$q_{k}(\bar{z})=\bar{q}_{k}$. Although it is difficult to analytically solve
the Pfaffian system, we can numerically evaluate $f_{k}$ and its derivatives
using the holonomic gradient method (HGM) [17]. Solving (25) is easier than
directly finding the first integrals that satisfy conditions (3) and (4) by
solving the Lagrange-Charpit system [9].

It should be noted that for the solutions of algebraic equations (25) to
exist, the dimension $d$ of $\bar{q}_{k}$ must be sufficiently larger than the
number of equations $t$. There are more than one zero-dimensional ideals that
annihilate the holonomic Hamiltonian, and the dimensions $d$ and $t$ depend on
the choice of the ideal. It is a part of future work to clarify which zero-
dimensional ideal yields $d$ sufficiently larger than $t$.

##  V Derivation of Finite Equation Set

To find the algebraic equations (25), an appropriate finite set of multi-
indices
$\Gamma=\\{{\boldsymbol{\mathrm{0}}},\gamma_{1},\dots,\gamma_{t-1}\\}\subset{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n}$ must be determined. In the proof of Lemma 5, the multi-indices yield
a basis
$\Omega,\mathcal{D}^{\gamma_{1}}\Omega,\dots,\mathcal{D}^{\gamma_{t-1}}\Omega$
of the $t$-dimensional ${\boldsymbol{\mathrm{R}}}(z)$-vector space
$\mathbb{M}$. It should be noted that finitely many matrices in
$\\{\mathcal{D}^{\alpha}\Omega\mid\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq
0}^{2n}\\}$ can be recursively computed by applying $\mathcal{D}_{i}\
(i=1,\dots,2n)$ to $\Omega$. Moreover, for the finite set of matrices, the
computation of its linearly independent subset is a straightforward process.
However, we cannot determine whether the linearly independent set is
sufficient to span $\mathbb{M}$ because dimension $t$ is unknown. In this
section, a method is proposed to compute $\Gamma$, which provides a basis of
$\mathbb{M}$, from a finite number of matrices $\mathcal{D}^{\alpha}\Omega$.

For a finite number of matrices
$\mathcal{D}^{\alpha_{1}}\Omega,\dots,\allowbreak\mathcal{D}^{\alpha_{s}}\Omega$,
suppose that there exist several ${\boldsymbol{\mathrm{R}}}(z)$-linear
relations between them:

$c_{j,\alpha_{1}}\mathcal{D}^{\alpha_{1}}\Omega+\cdots+c_{j,\alpha_{s}}\mathcal{D}^{\alpha_{s}}\Omega=0\\\
(c_{j,\alpha_{i}}\in{\boldsymbol{\mathrm{R}}}(z);\ i=1,\dots,s;\
j=1,\dots,r).$ (26)

Equations (26) with (15) indicate that any Poisson bracket $\\{f,g\\}$ with
$f,g\in\mathbb{S}(\mathcal{P}_{1},\dots,\mathcal{P}_{2n})$ is annihilated by
differential operators $\mathcal{Q}_{j}\coloneqq
c_{j,\alpha_{1}}\partial^{\alpha_{1}}+\cdots+c_{j,\alpha_{s}}\partial^{\alpha_{s}}\
(j=1,\dots,r)$. If the ideal
$I=\langle\mathcal{Q}_{1},\dots,\mathcal{Q}_{r}\rangle$ is zero-dimensional,
the quotient space ${\boldsymbol{\mathrm{R}}}(z)\langle\partial\rangle/I$ is a
finite dimensional ${\boldsymbol{\mathrm{R}}}(z)$-vector space with a basis
$\\{1,\partial^{\beta_{1}},\dots,\partial^{\beta_{t^{\prime}-1}}\\}$
comprising the standard monomials [13] of $I$. Hence, any $\partial^{\alpha}\
(\alpha\in{\boldsymbol{\mathrm{Z}}}_{\geq 0}^{2n})$ can be expressed as a
unique ${\boldsymbol{\mathrm{R}}}(z)$-linear combination of the elements of
the basis. This implies that the set of matrices
$B\coloneqq\\{\Omega,\mathcal{D}^{\beta_{1}}\Omega,\dots,\mathcal{D}^{\beta_{t^{\prime}-1}}\Omega\\}$
spans $\mathbb{M}$. If $B$ is linearly independent, dimension $t$ is given as
$t=t^{\prime}$ and
$\Gamma=\\{{\boldsymbol{\mathrm{0}}},\beta_{1},\dots,\beta_{t^{\prime}-1}\\}$.
If not, the minimal linearly independent subset of $B$ is a basis of
$\mathbb{M}$. The procedure for computing $\Gamma$ is summarized in Algorithm
1.

Algorithm 1 Derivation of $\Gamma$ in (25)

1:Matrices of rational functions
$\Omega,A_{1},A_{2},\dots,A_{2n}\in{\boldsymbol{\mathrm{R}}}(z)^{d\times d}$

2:Finite set of multi-indices $\Gamma$

3:Compute finite number of matrices $\mathcal{D}^{\alpha}\Omega$

4:Find ${\boldsymbol{\mathrm{R}}}(z)$-linear relations between them and
corresponding differential operators $\mathcal{Q}_{j}\ (j=1,\dots,r)$

5:if $I=\langle\mathcal{Q}_{1},\dots,\mathcal{Q}_{r}\rangle$ is not zero-
dimensional then

6: Get back to line 3

7:Find a basis
$\\{1,\partial^{\beta_{1}},\dots,\partial^{\beta_{t^{\prime}-1}}\\}$ of
quotient space ${\boldsymbol{\mathrm{R}}}(z)/I$

8:if
$B=\\{\Omega,\mathcal{D}^{\beta_{1}}\Omega,\dots,\mathcal{D}^{\beta_{t^{\prime}-1}}\Omega\\}$
is linearly independent then

9: return $\Gamma\leftarrow\\{1,\beta_{1},\dots,\beta^{t^{\prime}-1}\\}$

10:else

11: Find minimal linearly independent subset $B^{\prime}\subset B$

12: return $\Gamma\leftarrow\\{\mbox{Multi-indices of matrices in
}B^{\prime}\\}$

##  VI Numerical Example

Consider the HJE for the following Hamiltonian.

$h(x,p)=-2p_{1}\sin(x_{1})+2x_{2}p_{2}-ap_{2}^{2}+bx_{1}^{4},$ (27)

where $x=[x_{1}\ x_{2}]^{\top}$, $p=[p_{1}\ p_{2}]^{\top}$, and
$a,b\in{\boldsymbol{\mathrm{R}}}$. To apply the proposed method to this HJE,
$h(x,p)$ should be verified to be holonomic. By virtue of Lemma 2, it is
sufficient to show that each term of $h(x,p)$ is a holonomic function. For
example, after several differentiations, it was found that the first term
$-2p_{1}\sin(x_{1})$ is annihilated by the following four differential
operators of the form (6): $\mathcal{P}_{1}=\partial_{x_{1}}^{2}+1$,
$\mathcal{P}_{2}=\partial_{x_{2}}$, $\mathcal{P}_{3}=\partial_{p_{1}}$, and
$\mathcal{P}_{4}=p_{2}\partial_{p_{2}}-1$. For the other terms, the
appropriate differential operators were found by differentiating them several
times, which guarantees that every term is holonomic and $h(x,p)$ is also
holonomic from Lemma 2. A zero-dimensional ideal that annihilates $h$ was
determined by computing the intersection of the ideals for all terms (see [15]
for details). It should be noted that this computation of intersection can be
performed independently of parameters $a$ and $b$.

Using the symbolic computation of differential operators, a Pfaffian system
satisfied by $h(x,p)$, i.e., the vector of differential operators
$\mathcal{B}$ and matrices of rational functions
$A_{i}(z)\in{\boldsymbol{\mathrm{R}}}(z)\ (i=1,\dots,4)$ are obtained as

$\displaystyle\mathcal{B}$ $\displaystyle=[1\ \partial_{p_{2}}\
\partial_{p_{1}}\ \partial_{x_{2}}\ \partial_{x_{1}}]^{\top},$ $\displaystyle
A_{x_{1}}$ $\displaystyle=\begin{bmatrix}0&0&0&0&1\\\ 0&0&0&0&0\\\
-\frac{4}{x_{1}p_{1}}&\frac{2p_{2}}{x_{1}p_{1}}&\frac{4}{x_{1}}&\frac{2x_{2}}{x_{1}p_{1}}&\frac{1}{p_{1}}\\\
0&0&0&0&0\\\
\frac{12}{x_{1}^{2}}&-\frac{6p_{2}}{x_{1}^{2}}&\frac{-12p_{1}}{x_{1}^{2}}-p_{1}&-\frac{6x_{2}}{x_{1}^{2}}&0\end{bmatrix},$

where $A_{2}$, $A_{3}$, and $A_{4}$ are omitted owing to space limitations.
From Lemma 4, the analytic solutions of the Pfaffian system defined on a
simply connected domain $U$ in $\\{z\in{\boldsymbol{\mathrm{R}}}^{4}\mid
x_{1}p_{1}p_{2}\neq 0\\}$ constitute the $5$-dimensional real vector space.

The matrices $B_{x},B_{p}\in{\boldsymbol{\mathrm{R}}}(z)^{2\times 5}$ and
$\Omega\in{\boldsymbol{\mathrm{R}}}(z)^{5\times 5}$ are derived from
$A_{1},\dots,A_{4}$ as constant matrices:

$B_{x}=\begin{bmatrix}0&0&0&0&1\\\
0&0&0&1&0\end{bmatrix},B_{p}=\begin{bmatrix}0&0&1&0&0\\\
0&1&0&0&0\end{bmatrix},$

and
$\Omega=B_{p}^{\top}b_{x}-B_{x}^{\top}B_{p}\in{\boldsymbol{\mathrm{R}}}(z)^{5\times
5}$. Using Algorithm 1, $\Gamma=\\{(0,0,0,0),(1,0,0,0),(0,0,1,0)\\}$ is
obtained, which indicates from Theorem 2 that (16) is valid if and only if

$\bar{q}_{1}^{\top}\Omega(\bar{z})\bar{q}_{2}=\bar{q}_{1}^{\top}\mathcal{D}_{1}\Omega(\bar{z})\bar{q}_{2}=\bar{q}_{1}^{\top}\mathcal{D}_{3}\Omega(\bar{z})\bar{q}_{2}=0$
(28)

holds. In this example, the condition (12) was obtained as

$\det\left\\{B_{p}[\bar{q}_{1}\
\bar{q}_{2}]\right\\}=\bar{q}_{1,3}\bar{q}_{2,2}-\bar{q}_{2,3}\bar{q}_{1,2}\neq
0.$ (29)

Consequently, any two vectors satisfying (28) and (29) yield the first
integrals of the Hamiltonian system associated with (27) that define a
solution of the HJE. It should be noted that all aforementioned computations
can be performed symbolically and include no approximations.

Based on the denominators of the components of $A_{1},\dots,A_{4}$, let
$U=\\{z\in{\boldsymbol{\mathrm{R}}}^{4}\mid x_{1}>0,p_{1}>0,p_{2}>0\\}$, and
let $\mathbb{VS}(A_{1},\dots,A_{4})$ and
$\mathbb{S}(\mathcal{P}_{1},\dots,\mathcal{P}_{4})$ be the solution spaces
defined on it. For a fixed point $\bar{z}=[\pi/6\ 1\ b(\pi/6)^{4}\
2/a]^{\top}$, which satisfies $h(\bar{z})=0$, we seek the other first integral
$f_{2}$ that defines a solution of the HJE. The constant vector $\bar{q}_{1}$
is computed as

$\bar{q}_{1}=(\mathcal{B}\bullet h)(\bar{z})=\left[0\ -2\ -1\ \frac{4}{a}\
-\frac{b\pi^{3}(\sqrt{3}\pi-24)}{1296}\right]^{\top}\in{\boldsymbol{\mathrm{R}}}^{5}$

For any constant vector $\bar{q}_{2}$ that satisfies (28) and (12) with
$\bar{q}_{1}$, the other first integral that defines a solution of the HJE is
obtained as
$f_{2}=\phi_{\mathcal{B}}^{-1}\circ\psi_{\bar{z}}^{-1}(\bar{q}_{2})$. An
appropriate $\bar{q}_{2}\in{\boldsymbol{\mathrm{R}}}^{5}$ satisfying (28) and
(29) can be selected because the number of equations is $3$. For example,
$\bar{q}_{2}=[0\ a\ 0\ -2\ 0]^{\top}$ satisfies both conditions. Although it
is difficult to compute its explicit expression, $f_{2}$ is uniquely
characterized as the first component of the solution of the Pfaffian system
defined by $A_{1},\dots,A_{4}$ with boundary condition
$Q(\bar{z})=\bar{q}_{2}$ and can be evaluated using the HGM [17].

##  VII Conclusion

In this paper, the solutions of the HJE with holonomic Hamiltonian were
investigated. A solution of the HJE is characterized by a finite number of the
first integrals of the corresponding Hamiltonian system. Under the assumption
that the Hamiltonian is holonomic, the finite dimensional solution space of
the Pfaffian system associated with the holonomic Hamiltonian can be used to
obtain the condition for the other first integrals to characterize a solution
of the HJE. We obtained the condition as a finite number of algebraic
equations, which can be computed through the symbolic computation of
differential operators. Although the algebraic equations do not necessarily
have a solution, they can be easily solved if at least one solution exists.

Future work directions of this study include the investigations related to the
range of problems or conditions on the Hamiltonian to guarantee the existence
of solutions of the algebraic equations. When the HJE is considered in control
theory, a particular solution of the HJE called the stabilizing solution is
required. The characterization of the stabilizing solution in terms of
holonomic functions, zero-dimensional ideals, and Pfaffian systems will also
be considered in future work.

## References

  * [1] J. A. E. Bryson and Y.-C. Ho, Applied Optimal Control, 1st edition, John Wiley & Sons, 1975.
  * [2] B. D. O. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, 1979.
  * [3] A. J. van der Schaft, “L2-gain analysis of nonlinear systems and nonlinear state feedback Hinf control,” IEEE Trans. Automat. Contr., vol. 37, no. 6, pp. 770–784, 1992.
  * [4] J. M. Scherpen, “Balancing for nonlinear systems,” Syst. Control Lett., vol. 21, no. 2, pp. 143–153, 1993.
  * [5] D. L. Lukes, “Optimal regulation of nonlinear dynamical systems,” SIAM J. Control, vol. 7, no. 1, pp. 75–100, 1969\.
  * [6] R. W. Beard, G. N. Saridis, and J. T. Wen, “Galerkin approximations of the generalized Hamilton-Jacobi-Bellman equation,” Automatica, vol. 33, no. 12, pp. 2159–2177, 1997.
  * [7] Y. Ito, K. Fujimoto, and Y. Tadokoro, “Kernel-based Hamilton-Jacobi equations for data-driven optimal and H-infinity control,” IEEE Access, vol. 8, pp. 131047–131062, 2020.
  * [8] Z. Hao, K. Fujimoto, and Q. Zhang, “Approximate solutions to the Hamilton-Jacobi equations for generating functions,” J. Syst. Sci. Complex., vol. 33, no. 2, pp. 261–288, 2020.
  * [9] N. Sakamoto, “Analysis of the Hamilton-Jacobi equation in nonlinear control theory by symplectic geometry,” SIAM J. Control Optim., vol. 40, no. 6, pp. 1924–1937, 2002.
  * [10] T. Ohtsuka, “Solutions to the Hamilton-Jacobi equation with algebraic gradients,” IEEE Trans. Automat. Contr., vol. 56, no. 8, pp. 1874–1885, 2011.
  * [11] Y. Kawano and T. Ohtsuka, “Algebraic solutions to the Hamilton-Jacobi equation with the time-varying Hamiltonian,” SICE J. Control Meas. Syst. Integr., vol. 6, no. 1, pp. 28–37, 2013.
  * [12] S. C. Coutinho, A Primer of Algebraic D-Modules, Cambridge University Press, 1995.
  * [13] M. Saito, B. Sturmfels, and N. Takayama, Gröbner Deformations of Hypergeometric Differential Equations, Springer-Verlag, 2000.
  * [14] T. Iori, “Study on solution method to Hamilton-Jacobi equation with holonomic Hamiltonian,” in Proc. 64th Japan Joint Automat. Contr. Conf. (in Japanese), Kyoto, Japan, 2021, pp. 562–567.
  * [15] T. Oaku, Y. Shiraki, and N. Takayama, “Algebraic algorithms for D-Modules and numerical analysis,” Computer Mathematics (Proceedings of ASCM 2003), vol. 10, pp. 23–39, 2003.
  * [16] M. Kauers and P. Paule, The Concrete Tetrahedron, 1st edition, Springer-Verlag, 2011.
  * [17] H. Nakayama, K. Nishiyama, M. Noro, K. Ohara, T. Sei, N. Takayama, and A. Takemura, “Holonomic gradient descent and its application to the Fisher-Bingham integral,” Adv. Appl. Math., vol. 47, no. 3, pp. 639–658, 2011.
  * [18] A. Kume and T. Sei, “On the exact maximum likelihood inference of Fisher-Bingham distributions using an adjusted holonomic gradient method,” Stat. Comput., vol. 28, no. 4, pp. 835–847, 2018.
  * [19] F. Bréhard, M. Joldes, and J.-B. Lasserre, “On a moment problem with holonomic functions,” in Proc. 44th Int. Symp. Symb. Algebr. Comput., pp. 66–73, 2019.
  * [20] T. Hibi, ed., Gröbner Bases: Statistics and Software Systems, 1st edition, Springer Japan, 2013.
  * [21] M. Kauers, “The holonomic toolkit,” in C. Schneider and J. Bluemlein, eds., Computer Algebra in Quantum Field Theory, Integration, Summation, and Special Functions, Springer-Verlag, pp. 119–144, 2013.


# Effects of Auxiliary Knowledge on Continual Learning

Giovanni Bellitto, Matteo Pennisi,
Simone Palazzo, Concetto Spampinato  PeRCeiVe Lab,
University of Catania, Italy.
Email<EMAIL_ADDRESS>Lorenzo Bonicelli, Matteo Boschini,
Simone Calderara  AImageLab,
University of Modena and Reggio Emilia.
Email<EMAIL_ADDRESS>
###### Abstract

In Continual Learning (CL), a neural network is trained on a stream of data
whose distribution changes over time. In this context, the main problem is how
to learn new information without forgetting old knowledge (i.e., Catastrophic
Forgetting). Most existing CL approaches focus on finding solutions to
preserve acquired knowledge, so working on the “past” of the model. However,
we argue that as the model has to continually learn new tasks, it is also
important to put focus on the “present” knowledge that could improve following
tasks learning. In this paper we propose a new, simple, CL algorithm that
focuses on solving the current task in a way that might facilitate the
learning of the next ones. More specifically, our approach combines the main
data stream with a secondary, diverse and uncorrelated stream, from which the
network can draw auxiliary knowledge. This helps the model from different
perspectives, since auxiliary data may contain useful features for the current
and the next tasks and incoming task classes can be mapped onto auxiliary
classes. Furthermore, the addition of data to the current task is implicitly
making the classifier more robust as we are forcing the extraction of more
discriminative features. Our method can outperform existing state-of-the-art
models on the most common CL Image Classification benchmarks.

##  I Introduction

Human beings and animals are naturally able to memorize information presented
in a sequence [2]; on the contrary, Artificial Neural Networks (ANNs) learning
from a non-i.i.d. stream of data incur in Catastrophic Forgetting [3, 4].
Continual Learning (CL) [5, 6] aims at designing methods that compensate for
this issue and facilitate the retention of previous knowledge either by means
of regularization [7, 8], architectural designs [9, 10] or (pseudo-)replay of
past data [4, 11, 12].

The insurgence of catastrophic forgetting is ascribed to the tendency of
models to rewrite their hidden representations as they adjust their parameters
to best fit an input distribution that changes in time [13]. However, _McRae &
Hetherington_ highlight a meaningful difference in the way humans and ML
models learn from a sequence of data: whenever human subjects are evaluated on
their ability to memorize a sequence of concepts, they start out possessing an
already-large body of knowledge [14]. In other words, humans are generalists
that can anchor novel data in the context of previous knowledge, while ANNs
must specialize on a limited pool of data at each time without any additional
reference.

An obvious choice to bridge this gap is pre-training the models on a large
amount of available off-the-shelf i.i.d. data, leading to a better
initialization for the learning procedure [14, 15]. However, we observe that
pre-training is not always rewarding in a CL setting, especially in case of
small-size replay memories: the ever-changing stream of data entails large
changes in model parameters, leading to the forgetting of the pre-training.

We instead propose a learning strategy to limit catastrophic forgetting by
providing an additional data stream (uncorrelated from the target data), from
which the network can draw auxiliary knowledge. The role of this data stream
is to provide models with a more stable representation of the world that can
be re-used for incrementally learning new classes or categories leveraging the
already-learned low-level features. Indeed, it appears that the human brain
can adapt and rewire itself more easily when learning new things related to
familiar skills because pre-existing neuronal structure constrains what one
can learn [16]. We attempt to enforce this concept into CL through the
definition of _an associative rule_ that helps learning new classes by
measuring the simultaneous firing of neurons between past knowledge and the
current data stream. This is implemented through a simple yet effective
strategy named _MAH_ , that, during a new task, assigns new classes to model’s
corresponding most activated heads.

Experimental results carried out on standard CL settings, involving CIFAR-10
and (a subset of) Tiny-ImageNet benchmarks, demonstrate that using a separate
auxiliary data stream is mostly beneficial with limited size buffer leading to
a performance gain of several percent points w.r.t. state-of-the-art methods.
Analogously, the MAH strategy reveals to be more effective than the standard
class mapping procedure independently from the buffer size. We also
investigate the role of model pre-training as compared to sustained auxiliary
data employment highlighting that, for small-size buffers, auxiliary data is
to be preferred to pre-training.

Figure 1: During the training of the current Task $T_{i}$ the model is trained
combining the “current” data coming from the Task Stream and the “past” data
stored inside the buffer. In addition, the remaining heads are trained using
the auxiliary data stream. At the beginning of a new Task $T_{i+1}$ the MAH
procedure is conducted as follows: 1) Only the heads trained with the
auxiliary data are kept activated; 2) the $T_{i+1}$ task is forwarded to the
frozen model in order to store activation information about the heads; 3) for
each class in task $T_{i+1}$, each new class is assigned to the head that
activates the most, replacing the corresponding auxiliary data class.

In conclusion, our strategy is beneficial in Continual Learning from multiple
perspectives: the model avoids overfitting current examples, learns more
general features and – as auxiliary data-points stand in for future examples –
better prepares to learn future classes by suitably associating past knowledge
to the new acquired one. All these aspects are mainly observed with reduced
buffer size, thus contributing to the efforts that aim at generalizing CL
approaches to real-world scenarios.

##  II Related Work

The seminal study by McCloskey and Cohen first drew attention to the tendency
of ANNs to forget previously learned knowledge catastrophically [3]. In spite
of the outstanding results achieved by deep learning models in recent years
[17, 18], this problem still persists and prevents ANNs from learning flexibly
from non-i.i.d. data-streams. To tackle this issue, researchers and
practitioners design CL methods, i.e., strategies that make machine learning
models retain high accuracy on previously seen data when trained on an ever-
changing input distribution [5, 6]. While many distinct strategies have been
applied for this purpose, CL approaches can be broadly categorized into two
families: structuring or generalization.

###  II-A Structuring approaches to Continual Learning

Methods in the first class aim at making interference between distinct
concepts less likely by endowing the stored knowledge with a disentangled
structure. [19] first pioneered the idea to reduce forgetting by
orthogonalizing feature representations. A similar approach was recently taken
by [20]. Alternatively, structuring can be pursued at an architectural level,
by explicitly allocating distinct subsets of model parameters to distinct
tasks [9, 21, 22], encouraging non-overlapping activation patterns for
different data [10, 23], or by simply applying dropout [24, 25]. Finally,
several approaches regularize back-propagation by projecting the gradient to
minimize the interference between tasks learned at different times [26, 27,
28]. While structuring approaches are usually characterized by a simpler
training procedure, they typically require the availability of a task-
identifier at test time.

###  II-B Generalization approaches to Continual Learning

At the opposite end of the spectrum, generalization methods prevent forgetting
by encouraging the model to compare and contrast input data encountered
throughout the sequence, thus recovering the i.i.d. property of training [14].
Most notably, rehearsal-based approaches do so by maintaining a working memory
of previously seen examples and interleaving them with the input data [4, 11,
29, 30, 31], while pseudo-rehearsal methods approximate this procedure with a
generative model [32, 12]. Other works prioritize the learning of high-level
representations either by adopting learning objectives designed not to disrupt
the performance on previous tasks [33, 34, 8, 35], or by making use of semi-
supervised learning techniques to learn general features [36, 37, 38]. The
generalization approach is taken to an extreme by [39], which shows
satisfactory results on CL benchmarks by training a model in an i.i.d. fashion
on samples gathered greedily from the input stream.

Generalization strategies naturally blend knowledge gathered at different
times to build a unified predictor, making them more reliable than their
structuring counterparts in the realistic settings where no task-identifier is
given at testing-time [40, 41]. The approach proposed in this paper aligns
with the former group of methods; indeed, we argue that generalization should
be extended beyond already-seen data and embrace yet-unseen knowledge as well.

TABLE I: Accuracy on Split-CIFAR-10 for several Continual Learning methods.
Best results in bold, second-best in italic.

Class-IL Buffer ER [42] GEM [26] A-GEM [27] iCaRL [33] FDR [43] GSS [30] HAL
[44] DER [29] DER++ [29] Ours 50 32.69 ± 0.39 22.10 ± 0.41 20.02 ± 0.08 55.51
± 1.64 28.32 ± 4.51 26.62 ± 1.36 25.26 ± 1.73 44.85 ± 2.71 49.28 ± 3.16 56.33
± 0.95 200 44.79 ± 1.86 25.54 ± 0.76 20.04 ± 0.34 49.02 ± 3.20 30.91 ± 2.74
39.07 ± 5.59 32.36 ± 2.70 61.93 ± 1.79 64.88 ± 1.17 70.86 ± 0.95 500 57.74 ±
0.27 26.20 ± 1.26 22.67 ± 0.57 47.55 ± 3.95 28.71 ± 3.23 49.73 ± 4.78 41.79 ±
4.46 70.51 ± 1.67 72.70 ± 1.36 75.07 ± 0.41 5120 82.47 ± 0.52 25.56 ± 3.46
21.99 ± 2.29 55.07 ± 1.55 19.70 ± 0.07 67.27 ± 4.27 59.12 ± 4.41 83.81 ± 0.33
85.24 ± 0.49 84.56 ± 0.55 Task-IL Buffer ER [42] GEM [26] A-GEM [27] iCaRL
[33] FDR [43] GSS [30] HAL [44] DER [29] DER++ [29] Ours 50 86.98 ± 1.19 81.36
± 1.43 81.09 ± 1.88 88.86 ± 2.51 85.23 ± 1.24 85.22 ± 1.03 78.73 ± 3.16 85.04
± 1.17 86.14 ± 2.56 89.57 ± 2.47 200 91.19 ± 0.94 90.44 ± 0.94 83.88 ± 1.49
88.99 ± 2.13 91.01 ± 0.68 88.80 ± 2.89 82.51 ± 3.20 91.40 ± 0.92 91.92 ± 0.60
93.30 ± 0.64 500 93.61 ± 0.27 92.16 ± 0.69 89.48 ± 1.45 88.22 ± 2.62 93.29 ±
0.59 91.02 ± 1.57 84.54 ± 2.36 93.40 ± 0.39 93.88 ± 0.50 93.62 ± 0.58 5120
96.98 ± 0.17 95.55 ± 0.02 90.10 ± 2.09 92.23 ± 0.84 94.32 ± 0.97 94.19 ± 1.15
88.51 ± 3.32 95.43 ± 0.33 96.12 ± 0.21 95.84 ± 0.42

##  III Method

Most CL methods use current and, if the method implies a rehearsal strategy,
past task classification heads during the training of the current task. Future
heads, that will be mapped to classes from following tasks, are not involved
in the process at all. This poses a potentially dangerous situation due to the
model minimizing its prediction scores for future heads, which results in a
high loss peak when these heads are used at the beginning of future tasks.

We propose to leverage an auxiliary data stream, not correlated with the main
task stream, in order to keep these future heads activated since the beginning
of training. The proposed strategy is also beneficial to learn more
distinguishing and reusable features, as the model cannot focus on simply
discriminating between the classes from the task at hand. Furthermore, since
auxiliary training leads future task heads to learn to recognize their own
specific patterns, we exploit this property to devise a “most activated heads”
(MAH) assignment strategy for future classes, that minimizes the loss peak
that the model typically incurs at the beginning of a new task. Hence, the use
of an auxiliary stream favors the current task and improves forward transfer
to future tasks. The proposed approach is illustrated in Fig. 1.

Formally, a typical CL classification problem requires solving several tasks
sequentially, where each task $T_{t}$, with $t\in\left\\{1,...,T\right\\}$ and
$T$ being the number of tasks, consists in learning to classify a set of
classes $\mathbf{C}_{t}$. In this work, we follow the common Class-IL and
Task-IL settings [41], which assume no overlap between classes from different
tasks.

Each task is associated with an i.i.d. distribution $D_{t}$ of
$(\mathbf{x},y)$ pairs of a data point with the corresponding class label from
$\mathbf{C}_{t}$. In practice, the distribution is approximated by a finite
set of samples, i.e.,
$\mathbf{D}_{t}=\left\\{\left(\mathbf{x}_{1},y_{1}\right),\left(\mathbf{x}_{2},y_{2}\right),\dots,\left(\mathbf{x}_{N_{t}},y_{N_{t}}\right)\right\\}$,
where $N_{t}$ is the number of examples for task $t$.

The objective of CL is to find a function $f_{\theta}$, depending on a set of
learnable parameters $\theta$, that minimizes a classification objective over
the entire task sequence, such as:

$\operatorname*{arg\,min}_{\theta}\sum_{t=1}^{T}\sum_{(\mathbf{x}_{i},y_{i})\in\mathbf{D}_{t}}\mathcal{L}_{\text{C}}\left(f_{\theta}(\mathbf{x}_{i}),y_{i}\right),$
(1)

where $\mathcal{L}_{\text{C}}$ is the classification loss (e.g., cross-
entropy).

While training for the current task, most recent CL approaches [7, 8, 11, 29]
attempt to reduce forgetting by adding an additional loss term that attempts
to retain accuracy on previously-seen tasks. The in-task objective at task $t$
then becomes:

$\operatorname*{arg\,min}_{\theta}\sum_{(\mathbf{x}_{i},y_{i})\in\mathbf{D}_{t}}\mathcal{L}_{\text{C}}\left(f_{\theta}(\mathbf{x}_{i}),y_{i}\right)+\mathcal{L}_{\text{CL}},$
(2)

where $\mathcal{L}_{\text{CL}}$ is a generic additional loss term that
implements countermeasures to catastrophic forgetting (e.g., experience replay
from a buffer of samples from previous tasks).

In the proposed scenario, an additional distribution of i.i.d. auxiliary data
$A$, where $A\neq D_{t}~{}\forall t$, is available to the model at training
time. Again, the distribution is represented by a set of sample/label pairs
$\mathbf{A}=\left\\{\left(\mathbf{x}_{1},y_{1}\right),\left(\mathbf{x}_{2},y_{2}\right),\dots,\left(\mathbf{x}_{N_{t}},y_{N_{t}}\right)\right\\}$,
where labels belong to the class set $\mathbf{C}_{A}$.

In the following, we explain the two key aspects of the proposed approach:
head pre-activation and “most activated heads” class mapping.

###  III-A Head pre-activation

To ensure that the model employs all of its classification heads from the
start, we use classes from the auxiliary dataset as “place-holders” for
classes from future tasks.

The basic requirement of an auxiliary dataset $\mathbf{A}$ is related to the
cardinality of its set of classes, $\left|\mathbf{C}_{A}\right|$, which should
satisfy the following condition:

$\left|\mathbf{C}_{A}\right|\geq\sum_{t=2}^{T}\left|\mathbf{C}_{t}\right|.$
(3)

In other words, the number of auxiliary classes should be at least equal to
the total number of classes in the sequence of continual learning tasks, minus
the number of classes from the first task. This guarantees that, when training
on the first task, the auxiliary dataset provides enough classes for the
classification heads reserved to future tasks.

_Before_ starting to train on the first task $t=1$, we randomly choose a
subset $\mathbf{C}_{A,t}\subseteq\mathbf{C}_{A}$, with cardinality
$\left|\mathbf{C}_{A,t}\right|=\sum_{t=2}^{T}\left|\mathbf{C}_{t}\right|$.
Samples from the selected classes are included in the auxiliary sub-dataset
$\mathbf{A}_{t}$ and class indexes from $\mathbf{C}_{A,t}$ are re-mapped to
the indexes of classes in $\cup_{t=2}^{T}\mathbf{C}_{t}$ corresponding to
future tasks.

At task $t$, we merge the corresponding dataset $\mathbf{D}_{t}$ and the
auxiliary sub-dataset $\mathbf{A}_{t}$ and train on the joint set of classes,
in order to minimize the following new in-task objective:

$\operatorname*{arg\,min}_{\theta}\sum_{(\mathbf{x}_{i},y_{i})\in\mathbf{D}_{t}\cup\mathbf{A}_{t}}\mathcal{L}_{\text{C}}\left(f_{\theta}(\mathbf{x}_{i}),y_{i}\right)+\mathcal{L}_{\text{CL}}.$
(4)

As a result, we ensure that all classification heads are employed, reducing
the risk of loss peaks on new tasks, and encourage the model to learn more
complex, discriminative and stable features.

###  III-B “Most activated heads” (MAH) class mapping

At the beginning of each task $t>1$, it is necessary to update the set of
auxiliary classes in $\mathbf{A}_{t}$, since $\left|\mathbf{C}_{t}\right|$
classes must be removed to make room for classes from the new task.

Moreover, in this scenario, it also makes sense to _assign_ the specific heads
that will correspond to classes in the new task, rather than simply
associating them to the next available heads. An appropriate class mapping can
make better (re)use of features learned by the model for classification of
auxiliary classes, and reduce high losses that may lead to forgetting
previously-learned features.

Our head assignment approach, named MAH from “most activated heads”, acts
before beginning to train on task $t>1$, by first computing the average logits
$\mathbf{l}_{c}$, i.e., pre-softmax head activations, for each task class
$c\in\mathbf{C}_{t}$: to this aim, we select the subset
$\mathbf{D}_{t,c}\subset\mathbf{D}_{t}$ which only contains elements of class
$c$, and average the corresponding logit vectors as returned by model
$f_{\theta}$:

$\mathbf{l}_{c}=\frac{1}{N_{t,c}}\sum_{\left(\mathbf{x}_{i},y_{i}\right)\in\mathbf{D}_{t,c}}f_{\theta}\left(\mathbf{x}_{i}\right),$
(5)

where $N_{t,c}$ is the number of elements of class $c$ in $\mathbf{D}_{t}$.

Then, each new class in $\mathbf{C}_{t}$ is simply associated to the
classification head that maximizes its predicted score, i.e.,
$\operatorname*{arg\,max}\mathbf{l}_{c}$. In case of index collisions, largest
values are given priority. Finally, the new auxiliary sub-dataset
$\mathbf{A}_{t}$ is updated by removing from $\mathbf{A}_{t-1}$ the set of
classes corresponding to selected indexes, i.e.,
$\left\\{\operatorname*{arg\,max}\mathbf{l}_{c}\right\\}_{c\in\mathbf{C}_{t}}$,
and training proceeds as previously described.

##  IV Experimental Settings

###  IV-A Datasets

We focus our experiments on two common evaluation protocols [41]: class-
incremental (Class-IL), where the model is asked to gradually solve the
complete problem but classes become available at different times; task-
incremental (Task-IL), where the model is guided by the task-identity and can
only focus to solve each task independently. Specifically, we leverage Split-
CIFAR-10 [45], a widely-used image classification dataset obtained by
splitting the 32$\times$32 images of CIFAR-10 into 5 binary tasks. For a more
comprehensive evaluation, we also test on the larger 64$\times$64 Split-Micro-
ImageNet: a novel benchmark composed of a 20-class subset of Tiny-ImageNet
[46], split into 5 tasks of 4 classes each.

As for the choice of the auxiliary data, we pair the original data with
similarly-sized datasets. In particular, the auxiliary dataset for Split-
CIFAR-10 consists of a subset of 10 super-classes from CIFAR-100, selected
among those which are not semantic-related to those contained in CIFAR-10. For
Split-Micro-ImageNet, we select a subset of 20 classes from ILSVRC-2012,
making sure that the chosen data is as unrelated as possible with the original
Tiny-ImageNet classes. In detail, we first remove Tiny-ImageNet classes from
the entire label set; then, we group the remaining 800 classes into 293 super-
classes, corresponding to synsets found at distance 8 from the _entity_ root
node. Finally, we apply Spectral Clustering to select the 20 classes which are
most representative of the super-classes.

###  IV-B Training procedure

We apply the approach described in Section III by adapting the DER++ [29]
method, a recent rehearsal-based approach inspired by knowledge distillation
principles. For a fair comparison among different models, in our experiments
we follow [29] and adopt the same training settings. As backbone, we use
ResNet-18 [47] (not pretrained). We optimize our model with SGD, for 50 epochs
on Split-CIFAR-10 and 100 on Split-Micro-ImageNet. During training, samples
from the current task and from auxiliary classes are combined so that each
mini-batch contains data from both domains. We apply random crops and
horizontal flips as data augmentation. All hyperparameters are as defined in
[29].

###  IV-C Results

To validate the effectiveness of our approach using auxiliary data during
training, we compare our method with other CL methods based on rehearsal
strategies: ER [42], GEM [26], A-GEM [27], iCaRL [33], FDR [43], GSS [30], HAL
[44], DER [29] and vanilla DER++ [29]. Performance for these methods is
reported from [29], except for the setup with buffer size equal to 50111In
this case, results were computed using the Mammoth framework for PyTorch:
https://github.com/aimagelab/mammoth.

As performance metrics, we report classification accuracy in the Class-IL
setting as the average of per-task classification accuracy computed after the
final task, and in the Task-IL setting as the average of per-task
classification accuracy computed at the end of each task.

TABLE II: Accuracy on Split-Micro-ImageNet for Experience Replay-based methods. Buffer | Method | Class-IL | Task-IL
---|---|---|---
50 | ER [42] | 22.58 ± 0.71 | 66.28 ± 1.83
DER [29] | 29.35 ± 2.16 | 70.88 ± 0.83
DER++ [29] | 31.92 ± 2.15 | 69.86 ± 1.96
Ours | 37.24 ± 2.76 | 71.84 ± 1.82
200 | ER [42] | 36.22 ± 1.06 | 77.82 ± 1.24
DER [29] | 46.18 ± 1.44 | 81.12 ± 1.59
DER++ [29] | 51.84 ± 1.32 | 82.66 ± 1.60
Ours | 52.83 ± 0.83 | 78.48 ± 1.42
500 | ER [42] | 49.70 ± 0.71 | 84.35 ± 0.79
DER [29] | 56.58 ± 2.44 | 84.56 ± 1.19
DER++ [29] | 60.28 ± 2.31 | 85.10 ± 0.93
Ours | 59.36 ± 0.81 | 79.95 ± 0.61
5120 | ER [42] | 70.40 ± 1.30 | 89.20 ± 0.01
DER [29] | 68.75 ± 0.25 | 89.35 ± 0.85
DER++ [29] | 74.98 ± 0.66 | 90.72 ± 0.65
Ours | 71.65 ± 0.82 | 85.93 ± 1.29

Table I and Table II report results on Split-CIFAR-10 and Split-Micro-
ImageNet, respectively. Our method yields the best Class-IL performance when
tested with small/medium buffer size. It is also noteworthy that as the
performance gain of our approach increases as the size of the buffer
decreases. When buffer size becomes significantly larger, vanilla DER++ still
achieves the best results, showing that retaining and replaying enough data
(5,120 samples represent more than 10% of the entire training set of CIFAR-10)
is still the best option to alleviate catastrophic forgetting, although this
goes in stark contrast to generalizing continual learning methods to real-
world problems. A similar behavior can be observed on the simpler Task-IL
setting, where our method obtains the highest performance or is on par with
existing methods under low–data availability regimes.

We also compare our approach to Co2L [36], that recently achieved state-the-
art performance in both settings222Co2L is not reported in Tables I and II, as
its training strategies are significantly different from the methods shown in
those tables.. Nevertheless, on Split-CIFAR-10, our method yields better
Class-IL performance than Co2L [36], that reaches 65.57$\pm$1.37 with buffer
size of 200 and 74.26$\pm$0.77 with buffer size of 500, respectively compared
70.86$\pm$0.95 and 75.07$\pm$0.41 by our method. The lower standard deviation
also shows that our approach tends to be more stable across tasks.

Finally, we monitor the loss over consecutive tasks in order to evaluate the
impact of auxiliary data during training. Fig. 2 shows the average training
loss for vanilla DER++ and our method on Split-Micro-ImageNet. It can be
observed that, as new tasks come in (every 100 epochs), the proposed approach
shows a smoother loss surface, conversely to the vanilla counterpart that,
instead, exhibits more noticeable peaks. Thus, the proposed strategy also
improves forward transfer and prevents disrupting gradient peaks when the
model switches to new tasks, resembling non-continual learning scenarios.

Figure 2: Training loss trend for our approach (red) and DER++ (blue) on
Split-Micro-ImageNet. When the model switches to a new task, MAH reduces loss
peaks by assigning new classes to the most suitable available heads trained on
auxiliary data.

###  IV-D Ablation study

In order to substantiate our design choices, we perform an ablation study to
quantify the contribution of a) using the auxiliary data stream and b) the MAH
strategy. The obtained results are reported in Table III and compared to
vanilla DER++ and DER++ with auxiliary data but without MAH. When the MAH
strategy is not used, the classification heads are selected in a sequential
order without making use of neural activation mapping between past (auxiliary)
and current classes.

Results show that training with auxiliary data yields significant performance
gains for all buffer sizes but 5,120, where replayed knowledge becomes
prevalent. In all cases, MAH outperforms sequential head mapping.

TABLE III: Ablation Study. Results obtained by the vanilla DER++ (first row of each block), DER++ with auxiliary data (second row), and the proposed method (third row), combining DER++ with auxiliary data and MAH strategy, for different buffer sizes. Buffer | Method | Class-IL | Task-IL
---|---|---|---
50 | DER++ | 49.28 ± 3.16 | 86.14 ± 2.56
$\hookrightarrow$ \+ AUX | 52.74 ± 1.02 | 88.51 ± 2.01
$\hookrightarrow$ \+ MAH | 56.33 ± 0.95 | 89.57 ± 2.47
200 | DER++ | 64.88 ± 1.17 | 91.92 ± 0.60
$\hookrightarrow$ \+ AUX | 69.91 ± 1.48 | 92.57 ± 1.02
$\hookrightarrow$ \+ MAH | 70.86 ± 0.95 | 93.30 ± 0.64
500 | DER++ | 72.70 ± 1.36 | 93.88 ± 0.50
$\hookrightarrow$ \+ AUX | 74.24 ± 0.61 | 93.93 ± 0.47
$\hookrightarrow$ \+ MAH | 75.07 ± 0.41 | 93.62 ± 0.58
5120 | DER++ | 85.24 ± 0.49 | 96.12 ± 0.21
$\hookrightarrow$ \+ AUX | 84.26 ± 0.22 | 95.58 ± 0.22
$\hookrightarrow$ \+ MAH | 84.56 ± 0.55 | 95.84 ± 0.42

###  IV-E Effect of Pre-training

We further investigate whether it is better to employ a backbone pre-trained
on auxiliary data or to train it from scratch using the proposed strategy. The
results of this analysis are reported in Table IV: pre-training on auxiliary
data appears to be always beneficial compared to training from scratch in
DER++ and as the buffer size increases. On the contrary, on small buffers,
using auxiliary data with our approach on a model trained from scratch yields
better performance than pre-training. Furthermore, with our method, pre-
training on auxiliary data leads instead to lower performance than training
from scratch, showing that pre-training is not always a reliable alternative
to continuously training with auxiliary data.

TABLE IV: Effect of Pre-training. Results obtained by the vanilla DER++ (first row of each block), DER++ pre-trained with auxiliary data (second row), the proposed method trained from scratch (third row), and the proposed model pre-trained with auxiliary data (fourth row), for different buffer sizes. Buffer | Method | Class-IL | Task-IL
---|---|---|---
50 | DER++ | 49.28 ± 3.16 | 86.14 ± 2.56
DER++ with pre-training | 50.82 ± 3.34 | 84.02 ± 2.98
Ours | 56.33 ± 0.95 | 89.57 ± 2.47
| Ours with pre-training | 53.52 ± 3.60 | 89.52 ± 0.88
200 | DER++ | 64.88 ± 1.17 | 91.92 ± 0.60
DER++ with pre-training | 68.71 ± 1.01 | 92.43 ± 0.53
Ours | 70.86 ± 0.95 | 93.30 ± 0.64
| Ours with pre-training | 65.17 ± 2.67 | 91.35 ± 1.71
500 | DER++ | 72.70 ± 1.36 | 93.88 ± 0.50
DER++ with pre-training | 75.91 ± 0.26 | 94.39 ± 0.29
Ours | 75.07 ± 0.41 | 93.62 ± 0.58
| Ours with pre-training | 71.39 ± 2.77 | 91.77 ± 0.76
5120 | DER++ | 85.24 ± 0.49 | 96.12 ± 0.21
DER++ with pre-training | 86.60 ± 0.42 | 96.29 ± 0.09
Ours | 84.56 ± 0.55 | 95.84 ± 0.42
| Ours with pre-training | 82.20 ± 1.37 | 94.65 ± 0.36

###  IV-F Generative Auxiliary Model

In the previous sections, we have consistently observed that using auxiliary
data helps retaining knowledge of previous tasks, especially with limited
buffer size. However, it is not efficient to maintain the auxiliary data in
memory, as in that case it is still preferable to simply use a larger buffer.

A viable alternative would be to replace the auxiliary stream with a
generative replay model and use generated samples during task training. In
order to investigate the feasibility of this option, we use a generative
adversarial network (GAN) [48] to learn the distribution of auxiliary data,
and employ synthetic images in place of real ones in Eq. 4 in our method.

We then replicate the experiments carried out on Split-CIFAR-10 by pre-
training a BigGAN model [49] on super-classes of CIFAR-100 used as auxiliary
data, and compare the results with those obtained when using real images. As
it can be seen in Table V, performance achieved when using generated images
follow the same behavior observed with real data, i.e., performance increase
in settings with small buffer size, while the approach is less beneficial when
the replay memory increases.

TABLE V: Effect of replacing auxiliary data with a GAN. Results obtained on Split-CIFAR-10 when using no auxiliary data (vanilla DER++), real auxiliary data (the proposed method) and synthetic data. Buffer | Aux. data | Class-IL | Task-IL
---|---|---|---
50 | none | 49.28 ± 3.16 | 86.14 ± 2.56
real | 56.33 ± 0.95 | 89.57 ± 2.47
synthetic | 54.47 ± 3.05 | 89.23 ± 1.83
200 | none | 64.88 ± 1.17 | 91.92 ± 0.60
real | 70.86 ± 0.95 | 93.30 ± 0.64
synthetic | 68.84 ± 0.77 | 92.85 ± 0.24
500 | none | 72.70 ± 1.36 | 93.88 ± 0.50
real | 75.07 ± 0.41 | 93.62 ± 0.58
synthetic | 74.35 ± 1.04 | 93.42 ± 0.40
5120 | none | 85.24 ± 0.49 | 96.12 ± 0.21
real | 84.56 ± 0.55 | 95.84 ± 0.42
synthetic | 84.14 ± 0.45 | 95.65 ± 0.23

##  V Conclusion

In this work, we propose a novel approach for improving continual learning
accuracy by leveraging external data. Our experiments show that providing the
model with an additional auxiliary stream leads to an increase in performance,
especially when the employed memory buffer is small, and to more stable
training at the beginning of each task. We also observe that our approach
works consistently better, on small buffer settings, than alternative
knowledge transfer strategies such as direct pre-training on auxiliary data.
The approximation of the auxiliary data distribution through the use of
generative models also outperforms state-of-the-art models, thus indicating
the future direction of this work, i.e., a more effective modeling of previous
real-world knowledge as it seems to happen in the human hippocampus [50].

## Acknowledgment

This research has been partially supported by: _RehaStart_ project, Regione
Sicilia (PO FESR 2014/2020, Azione 1.1.5, N. 08ME6201000222, CUP
G79J18000610007); Piano della Ricerca di Ateneo 2020/2022, PIACERI Starting
Grant Università di Catania, Project: _BrAIn: Adaptive Brain-Derived
Artificial Intelligence Methods for Event Detection_ ; Piano della Ricerca di
Ateneo 2020/2022, Linea 2D, Università di Catania.

## References

  * [1]
  * [2] J. Cichon and W.-B. Gan, “Branch-specific dendritic ca 2+ spikes cause persistent synaptic plasticity,” _Nature_ , 2015.
  * [3] M. McCloskey and N. J. Cohen, “Catastrophic interference in connectionist networks: The sequential learning problem,” in _Psychol. Learn. Motiv._ , 1989.
  * [4] R. Ratcliff, “Connectionist models of recognition memory: constraints imposed by learning and forgetting functions.” _Psychol. Rev._ , 1990.
  * [5] G. I. Parisi, R. Kemker, J. L. Part, C. Kanan, and S. Wermter, “Continual lifelong learning with neural networks: A review,” _Neural Netw._ , 2019\.
  * [6] M. D. Lange, R. Aljundi, M. Masana, S. Parisot, X. Jia, A. Leonardis, G. G. Slabaugh, and T. Tuytelaars, “Continual learning: A comparative study on how to defy forgetting in classification tasks,” _IEEE TPAMI_ , 2021.
  * [7] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska _et al._ , “Overcoming catastrophic forgetting in neural networks,” _PNAS_ , 2017\.
  * [8] Z. Li and D. Hoiem, “Learning without forgetting,” _IEEE TPAMI_ , 2017.
  * [9] A. A. Rusu, N. C. Rabinowitz, G. Desjardins, H. Soyer, J. Kirkpatrick, K. Kavukcuoglu, R. Pascanu, and R. Hadsell, “Progressive neural networks,” _arXiv preprint arXiv:1606.04671_ , 2016.
  * [10] J. Serra, D. Suris, M. Miron, and A. Karatzoglou, “Overcoming catastrophic forgetting with hard attention to the task,” in _ICML_ , 2018.
  * [11] A. Chaudhry, M. Rohrbach, M. Elhoseiny, T. Ajanthan, P. K. Dokania, P. H. Torr, and M. Ranzato, “On tiny episodic memories in continual learning,” in _ICML Workshops_ , 2019.
  * [12] H. Shin, J. K. Lee, J. Kim, and J. Kim, “Continual learning with deep generative replay,” in _ANeurIPS_ , 2017.
  * [13] R. M. French, “Using semi-distributed representations to overcome catastrophic forgetting in connectionist networks,” in _Proc. Annu. Conf. Cogn. Sci. Soc._ , 1991.
  * [14] K. McRae and P. A. Hetherington, “Catastrophic interference is eliminated in pretrained networks,” in _Proc. Annu. Conf. Cogn. Sci. Soc._ , 1993.
  * [15] D. Hendrycks, K. Lee, and M. Mazeika, “Using pre-training can improve model robustness and uncertainty,” in _ICML_ , 2019.
  * [16] S. Sakellaridi, V. Christopoulos, T. Aflalo, K. Pejsa, E. Rosario, D. Ouellette, N. Pouratian, and R. Andersen, “Intrinsic variable learning for brain-machine interface control by human anterior intraparietal cortex,” _Neuron_ , 2019.
  * [17] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot _et al._ , “Mastering the game of go with deep neural networks and tree search,” _Nature_ , 2016.
  * [18] O. Vinyals, I. Babuschkin, W. M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D. H. Choi, R. Powell, T. Ewalds, P. Georgiev _et al._ , “Grandmaster level in starcraft ii using multi-agent reinforcement learning,” _Nature_ , 2019.
  * [19] S. Lewandowsky and S.-C. Li, “Catastrophic interference in neural networks: Causes, solutions, and data,” in _Interference and inhibition in cognition_. Elsevier, 1995.
  * [20] A. Chaudhry, N. Khan, P. K. Dokania, and P. H. Torr, “Continual learning in low-rank orthogonal subspaces,” in _ANeurIPS_ , 2020.
  * [21] A. Mallya and S. Lazebnik, “Packnet: Adding multiple tasks to a single network by iterative pruning,” in _CVPR_ , 2018.
  * [22] S. Jung, H. Ahn, S. Cha, and T. Moon, “Continual learning with node-importance based adaptive group sparse regularization,” in _ANeurIPS_ , 2020.
  * [23] J. Rajasegaran, M. Hayat, S. H. Khan, F. S. Khan, and L. Shao, “Random path selection for continual learning,” in _ANeurIPS_ , 2019.
  * [24] I. J. Goodfellow, M. Mirza, X. Da, A. C. Courville, and Y. Bengio, “An empirical investigation of catastrophic forgeting in gradient-based neural networks,” _ICLR_ , 2014.
  * [25] S. I. Mirzadeh, M. Farajtabar, R. Pascanu, and H. Ghasemzadeh, “Understanding the role of training regimes in continual learning,” in _ANeurIPS_ , 2020\.
  * [26] D. Lopez-Paz and M. Ranzato, “Gradient episodic memory for continual learning,” in _ANeurIPS_ , 2017.
  * [27] A. Chaudhry, M. Ranzato, M. Rohrbach, and M. Elhoseiny, “Efficient Lifelong Learning with A-GEM,” in _ICLR_ , 2019.
  * [28] M. Farajtabar, N. Azizan, A. Mott, and A. Li, “Orthogonal gradient descent for continual learning,” in _AISTATS_ , 2020.
  * [29] P. Buzzega, M. Boschini, A. Porrello, D. Abati, and S. Calderara, “Dark Experience for General Continual Learning: a Strong, Simple Baseline,” in _ANeurIPS_ , 2020.
  * [30] R. Aljundi, M. Lin, B. Goujaud, and Y. Bengio, “Gradient based sample selection for online continual learning,” in _ANeurIPS_ , 2019.
  * [31] R. Aljundi, E. Belilovsky, T. Tuytelaars, L. Charlin, M. Caccia, M. Lin, and L. Page-Caccia, “Online continual learning with maximal interfered retrieval,” in _ANeurIPS_ , 2019.
  * [32] M. Riemer, T. Klinger, D. Bouneffouf, and M. Franceschini, “Scalable recollections for continual lifelong learning,” in _AAAI_ , 2019.
  * [33] S.-A. Rebuffi, A. Kolesnikov, G. Sperl, and C. H. Lampert, “icarl: Incremental classifier and representation learning,” in _CVPR_ , 2017.
  * [34] S. Hou, X. Pan, C. C. Loy, Z. Wang, and D. Lin, “Learning a unified classifier incrementally via rebalancing,” in _CVPR_ , 2019.
  * [35] Y. Wu, Y. Chen, L. Wang, Y. Ye, Z. Liu, Y. Guo, and Y. Fu, “Large scale incremental learning,” in _CVPR_ , 2019.
  * [36] H. Cha, J. Lee, and J. Shin, “Co2l: Contrastive continual learning,” in _ICCV_ , 2021.
  * [37] Q. Pham, C. Liu, and S. Hoi, “Dualnet: Continual learning, fast and slow,” in _ANeurIPS_ , 2021.
  * [38] C. D. Kim, J. Jeong, S. Moon, and G. Kim, “Continual learning on noisy data streams via self-purified replay,” in _ICCV_ , 2021.
  * [39] A. Prabhu, P. H. Torr, and P. K. Dokania, “GDumb: A simple approach that questions our progress in continual learning,” in _ECCV_ , 2020.
  * [40] S. Farquhar and Y. Gal, “Towards Robust Evaluations of Continual Learning,” in _ICML Workshops_ , 2018.
  * [41] G. M. van de Ven and A. S. Tolias, “Three continual learning scenarios,” in _NeurIPS Workshops_ , 2018.
  * [42] M. Riemer, I. Cases, R. Ajemian, M. Liu, I. Rish, Y. Tu, and G. Tesauro, “Learning to Learn without Forgetting by Maximizing Transfer and Minimizing Interference,” in _ICLR_ , 2019.
  * [43] A. S. Benjamin, D. Rolnick, and K. Kording, “Measuring and regularizing networks in function space,” in _ICLR_ , 2019.
  * [44] A. Chaudhry, A. Gordo, P. K. Dokania, P. Torr, and D. Lopez-Paz, “Using hindsight to anchor past knowledge in continual learning,” in _AAAI_ , 2021\.
  * [45] F. Zenke, B. Poole, and S. Ganguli, “Continual learning through synaptic intelligence,” in _ICML_ , 2017.
  * [46] Stanford, “Tiny ImageNet Challenge (CS231n),” 2015, https://www.kaggle.com/c/tiny-imagenet.
  * [47] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _ICPR_ , 2016.
  * [48] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _ANeurIPS_ , 2014.
  * [49] A. Brock, J. Donahue, and K. Simonyan, “Large scale gan training for high fidelity natural image synthesis,” in _ICLR_ , 2019.
  * [50] G. Rothschild, E. Eban, and L. M. Frank, “A cortical-hippocampal-cortical loop of information processing during memory consolidation,” _Nat.Neurosci._ , 2017.