# A general alternating-direction implicit Newton method for solving complex
continuous-time algebraic Riccati matrix equation ††thanks: The work was
supported in part by the National Natural Science Foundation of China
(12171412, 11771370), Natural Science Foundation for Distinguished Young
Scholars of Hunan Province (2021JJ10037), Hunan Youth Science and Technology
Innovation Talents Project (2021RC3110), the Key Project of Education
Department of Hunan Province (19A500, 21A0116).

Shifeng Li, Kai Jiang Juan Zhang Corresponding
author(zhangjuan@xtu.edu.cn).Key Laboratory of Intelligent Computing and
Information Processing of Ministry of Education, Hunan Key Laboratory for
Computation and Simulation in Science and Engineering, School of Mathematics
and Computational Science, Xiangtan University, Xiangtan, Hunan, 411105, P. R.
China.

Abstract In this paper, applying the Newton method, we transform the complex
continuous-time algebraic Riccati matrix equation into a Lyapunov equation.
Then, we introduce an efficient general alternating-direction implicit (GADI)
method to solve the Lyapunov equation. The inexact Newton-GADI method is
presented to save computational amount effectively. Moreover, we analyze the
convergence of the Newton-GADI method. The convergence rate of the Newton-GADI
and Newton-ADI methods is compared by analyzing their spectral radii.
Furthermore, we give a way to select the quasi-optimal parameter.
Corresponding numerical tests are shown to illustrate the effectiveness of the
proposed algorithms.

Keywords Complex algebraic Riccati equation; Newton method; alternating-
direction implicit method; convergence analysis.

##  1 Introduction

In this paper, we consider the complex continuous-time algebraic Riccati
matrix equation (CARE) defined by

$\displaystyle A^{*}X+XA-XKX+Q=0,$ (1)

where $A,~{}K,~{}Q\in\mathbb{C}^{n\times n}$, $K=K^{*}$, $Q=Q^{*}$, and $X$ is
an unknown matrix.

The CARE arises in many areas such as linear and non-linear optimal control
systems [1]-[3], Wiener-Hopf factorization of Markov chains [4], total least
squares problems [5], spectral factorization of rational matrix functions
[6]-[7], matrix sign functions [8]-[9], and transport theory [10]. For
example, consider the linear time-invariant system represented by

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)+B\bm{u}(t),~{}~{}\bm{p}(0)=\bm{p}_{0},$

where $B\in\mathbb{C}^{n\times m}$, $A$ is defined by (1),
$\bm{p}(t)\in\mathbb{C}^{n}$ and $\bm{u}(t)\in\mathbb{C}^{m}$ are the state
and control vectors of the system. We can represent
$u(t)=-R^{-1}B^{*}\bm{p}(t)$ for a vector $\bm{x}(t)$ such that the cost
functional

$\displaystyle
J(u)=\int^{t_{1}}_{t_{0}}[\bm{p}^{*}(t)Q\bm{p}(t)+\bm{u}^{*}(t)R\bm{u}(t)]\,dt+\bm{p}^{*}(t_{1})Z_{1}\bm{p}(t_{1})$
(2)

is minimized, where the cost penalty matrix $Q$ defined by (1) is positive
semi-definite, and the cost penalty matrices $R,~{}Z_{1}$ are positive
definite. Taken $K=BR^{-1}B^{*}$, defined by (1), the state equation becomes

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)-K\bm{x}(t).$

The cost function reduces to

$\displaystyle\hat{J}(u)=\int^{t_{1}}_{t_{0}}[\bm{p}^{*}(t)Q\bm{p}(t)+\bm{x}^{*}(t)K\bm{x}(t)]\,dt+\bm{p}^{*}(t_{1})Z_{1}\bm{p}(t_{1}).$

We wish to choose $\bm{x}(t)$ to minimize $\hat{J}$. This can be accomplished
by recognizing the fact input $\bm{x}(t)+\delta\hat{\bm{x}}(t)$, which is some
deviation from the desired stabilizing input. The corresponding state is
$\bm{p}(t)+\delta\hat{\bm{p}}(t)$ (by linearity). It can be shown that $J$ is
a minimum when the derivative of the adjusted cost function $\hat{J}$ with
respect to $\delta$ is zero [11], i.e.,

$\displaystyle\int^{t_{1}}_{t_{0}}[\bm{p}^{*}(t)Q\hat{\bm{p}}(t)+\bm{x}^{*}(t)A\hat{\bm{p}}(t)-\bm{x}^{*}(t)\dot{\hat{\bm{p}}}(t)]\,dt+\bm{p}^{*}(t_{1})Z_{1}\hat{\bm{p}}(t_{1})=0.$

Using integration by parts and $\hat{\bm{p}}(t_{0})=0$, we have

$\displaystyle\int^{t_{1}}_{t_{0}}[\bm{p}^{*}(t)Q+\bm{x}^{*}(t)A+\dot{\bm{x}}^{*}(t)]\hat{\bm{p}}(t)\,dt+[\bm{p}^{*}(t_{1})Z_{1}-\bm{x}^{*}(t_{1})]\hat{\bm{p}}(t_{1})=0.$

If the system is controllable, then a suitable choice of $\hat{\bm{x}}(t)$
gives any $\hat{\bm{p}}(t)$. This leads to the following requirements

$\displaystyle\begin{cases}\bm{p}^{*}(t)Q+\bm{x}^{*}(t)A+\dot{\bm{x}}^{*}(t)=0,\\\
\bm{p}^{*}(t_{1})Z_{1}-\bm{x}^{*}(t_{1})=0.\end{cases}$

Now the system becomes

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)-K\bm{x}(t).$

Since there exists a state transition matrix by linearity

$\displaystyle X(t,t_{1})=\begin{pmatrix}X_{11}(t,t_{1})&X_{12}(t,t_{1})\\\
X_{21}(t,t_{1})&X_{22}(t,t_{1})\end{pmatrix}$

such that

$\displaystyle\begin{pmatrix}\bm{p}(t)\\\
\bm{x}(t)\end{pmatrix}=\begin{pmatrix}X_{11}(t,t_{1})&X_{12}(t,t_{1})\\\
X_{21}(t,t_{1})&X_{22}(t,t_{1})\end{pmatrix}\begin{pmatrix}\bm{p}(t_{1})\\\
\bm{x}(t_{1})\end{pmatrix}.$

Then, we have $\bm{x}(t)=X(t)\bm{p}(t)$, where $X(t)$ is the solution of the
continuous time differential Riccati equation [11]

$\displaystyle-\frac{dX(t)}{dt}=A^{*}X(t)+X(t)A-X(t)KX(t)+Q.$ (3)

When $t_{1}\rightarrow\infty$, we get a steady state stabilizing solution
$\dot{X}(t)=0$. The continuous time differential Riccati equation (3) reduces
to the CARE (1). The optimal choice of input that minimizes the infinite
horizon cost function

$\displaystyle
J(u)=\int^{\infty}_{t_{0}}[\bm{p}^{*}(t)Q\bm{p}(t)+\bm{u}^{*}(t)R\bm{u}(t)]\,dt$

is $\bm{u}(t)=-R^{-1}B^{*}X\bm{p}(t)$.

Therefore, there are many scholars pay much attention to studying the CARE
(1). Many works are concerned with the existence of positive (semi-) definite
solution for this equation. In 1961, Kalman derived the existence condition of
the positive definite solution using observability and controllability [12].
Subsequently, the existence conditions of the positive (semi-) definite
solution have been considerably investigated under suitable assumptions in
[14]-[17].

######  Theorem 1.1.

[17] For the CARE (1), if the pair $(A,K)$ is stabilizable and $(A,Q)$ is
detectable, then this equation has a unique Hermitian positive semi-definite
solution $X$. $(A,K)$ is called stabilizable if there exists a matrix
$F\in\mathbb{C}^{n\times n}$ such that $A-KF$ is stable, i.e., all eigenvalues
are in the open left half-plane $\mathbb{C}_{-}$. The pair $(A,Q)$ is called
detectable if $(A^{*},Q^{*})$ is stabilizable.

Additionally, there are many works focusing on the numerical algorithms for
the CARE (1), such as the Schur method [18], the matrix sign function [19],
the structure-preserving doubling algorithm [20], and Krylov subspace
projection method [21]-[22]. Kleinman [23], Banks and Ito [24] applied the
Newton method, due to its quadratic convergence, to solve this equation.
Nevertheless, at each Newton iteration step, a Lyapunov equation needs to be
solved to get next iteration solution. Therefore, Navasca and Morris in
[25]-[26] combined the Newton method with a modified alternating-direction
implicit (ADI) method. Benner and his coauthors used the variant of Newton-ADI
algorithm to solve large-scale Riccati equations [27]-[28]. Moreover,
Feitzinger et al. in [30]-[31] proposed and analyzed the inexact Newton-ADI
method. Recently, [29] has proposed an iteration scheme and given its
theoretical analysis. It should be noticed that this algorithm still needs to
solve a Lyapunov equation twice at each iteration step. Based on this, we
concentrate on solving the CARE (1) by combing the Newton method with the
generalized alternating directions implicit (GADI).

The remaining part of this paper is organized as follows. In Section 2, we
propose two Newton-GADI methods, including the Newton-GADI and inexact Newton-
GADI algorithms. Moreover, we show the convergence analysis of the Newton-
GADI. A practical method is given to select the quasi-optimal parameter. In
Section 3, some numerical examples are devoted to showing the effectiveness of
the proposed algorithms. We draw some conclusion and remarks in the last
section.

Throughout this paper, let $\mathbb{C}^{n\times m}$ be the set of all $n\times
m$ complex matrices. $I_{n}$ is the identity matrix of order $n$. If
$A\in\mathbb{C}^{n\times n}$, the symbols $A^{\ast}$, $A^{-1}$ and $\|A\|_{2}$
express the conjugate transpose, the inverse, the spectral norm of $A$,
respectively. The eigenvalue and singular value sets of $A$ are denoted as
$\Lambda(A)=\\{\lambda_{i}(A),~{}i=1,2,\cdots,n\\},~{}~{}\Sigma(A)=\\{\sigma_{i}(A),~{}i=1,2,\cdots,n\\},$
where $\lambda_{i}(A)$ and $\sigma_{i}(A)$ represent the $i$-th component,
arranged in non-increasing order. $\rho(A)=\max\limits_{1\leq i\leq
n}\\{|\lambda_{i}(A)|\\}$ represents the spectral radius of $A$.

##  2 Newton-GADI algorithms

It is well known that one classical approach to solving the CARE (1) is to
tackle its nonlinearity with a Newton-type method. In this section, we apply
the Newton method to transform this equation into a Lyapunov equation. Then we
use the GADI method to solve the Lyapunov equation. The whole framework is
denoted by the Newton-GADI. To save computational amount, the inexact Newton-
GADI method is presented. Further, we give the convergence analysis of the
Newton-GADI algorithm.

###  2.1 Iteration schemes

Now we introduce the usual scheme of the Newton method for solving the CARE
(1).

Define the mapping $\mathcal{R}:\mathbb{C}^{n\times
n}\rightarrow\mathbb{C}^{n\times n}$:

$\displaystyle\mathcal{R}(X)=A^{*}X+XA-XKX+Q,~{}X\in\mathbb{C}^{n\times n},$

where $A,~{}K,~{}Q$ are defined in (1). The first Fr$\acute{\mbox{e}}$chet
derivative of $\mathcal{R}$ at a matrix $X$ is a linear map
$\mathcal{R}^{\prime}_{X}:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n\times
n}$ given by

$\displaystyle\mathcal{R}^{\prime}_{X}{(E)}=(A-KX)^{*}E+E(A-KX).$

The Newton method for the CARE (1) is

$\displaystyle
X_{k+1}=X_{k}-(\mathcal{R}^{\prime}_{X_{k}})^{-1}\mathcal{R}(X_{k}),~{}k=0,1,\cdots,$
(4)

given that the map $\mathcal{R}^{\prime}_{X_{k}}$ is invertible. The Newton
iteration scheme (4) is equivalent to

$\displaystyle(A-KX_{k})^{*}X_{k+1}+X_{k+1}(A-KX_{k})+X_{k}KX_{k}+Q=0.$ (5)

The convergence property for the Newton scheme (5) is illustrated by the
following result.

######  Theorem 2.1.

[23] For the CARE (1), assume that $(A,K)$ is stabilizable and $(A,Q)$ is
detectable, starting with any Hermitian matrix $X_{0}$ such that $A-KX_{0}$ is
stable, the sequence of Hermitian matrices $\\{X_{k}\\}_{k=1}^{\infty}$
determined by (5) quadratic converges to the Hermitian positive semi-definite
solution $X_{*}$ of (1). In other words, there exists a constant $\delta>0$
such that

$\displaystyle\|X_{k}-X_{*}\|\leq\delta\|X_{k-1}-X_{*}\|^{2},~{}~{}k=0,1,\cdots,$

where $\|\cdot\|$ is any given matrix norm. Moreover, the iteration sequence
$\\{X_{k}\\}_{k=1}^{\infty}$ has a monotone convergence behavior, i.e.,

$\displaystyle 0\leq X_{*}\leq\cdots\leq X_{k+1}\leq X_{k}\leq\cdots\leq
X_{1}.$

Denote $A_{k}=KX_{k}-A$ and $F(X_{k})=X_{k}KX_{k}+Q$, where the mapping
$F:\mathbb{C}^{n\times n}\rightarrow\mathbb{C}^{n\times n}$. Actually, due to
$A-KX_{k}$ is stable, then $A_{k}$ is stable, i.e.,
$\mbox{Re}(\lambda_{i}(A_{k}))>0$. The iteration scheme can be written as

$\displaystyle A_{k}^{*}X_{k+1}+X_{k+1}A_{k}=F(X_{k}).$ (6)

It is required to solve the Lyapunov equation (6) in each Newton iteration
step. We recall an existence and uniqueness theorem for this equation.

######  Theorem 2.2.

The Lyapunov equation

$\displaystyle\hat{A}^{*}Y+Y\hat{A}+Z=0$

has a uniquely Hermitian positive semi-definite solution
$Y\in\mathbb{C}^{n\times n}$ if and only if $\hat{A}\in\mathbb{C}^{n\times n}$
is stable, for any Hermitian matrix $Z\in\mathbb{C}^{n\times n}$.

We apply the GADI framework in [32] to solve the Lyapunov equation (6). The
GADI scheme is

$\displaystyle\begin{cases}(\alpha_{k}I+A_{k}^{*})X_{k+1,\ell+\frac{1}{2}}=X_{k+1,\ell}(\alpha_{k}I-A_{k})+F(X_{k}),\\\
X_{k+1,\ell+1}(\alpha_{k}I+A_{k})=X_{k+1,\ell}[A_{k}-(1-\omega_{k})\alpha_{k}I]+(2-\omega_{k})\alpha_{k}X_{k+1,\ell+\frac{1}{2}},\end{cases}$
(7)

where $\ell=0,1,\cdots$, $X_{k+1,0}=X_{k}$, $\alpha_{k}>0$ and
$\omega_{k}\in[0,2)$.

Combining the Newton scheme with the GADI scheme, we have the Newton-GADI
algorithm, summarized in Algorithm 1 below.

Algorithm 1 (Newton-GADI Algorithm)

Given the matrices $A,~{}K,~{}Q$, the integer numbers $k_{\max}>1$ and
$\ell_{\max}>1$, the outer and inter iteration tolerance $\varepsilon_{out}$
and $\varepsilon_{inn}$, respectively.

Step 1. Compute $\beta=1+\|A\|_{\infty},~{}~{}B=A+\beta I$ and the initial
matrix $X_{0}=X^{-1}$ by solving Lyapunov equation $B^{*}X+XB-2Q=0$.

Step 2. Set $X_{k+1,0}=X_{k}$ and compute the matrices

$\displaystyle A_{k}=KX_{k}-A,~{}~{}F(X_{k})=X_{k}KX_{k}+Q.$

Step 3. Solve the GADI iteration scheme (7) to obtain
$X_{k+1}=X_{k+1,\ell_{k}}$ such that

$\displaystyle\|A^{*}X_{k+1,\ell_{k}}+X_{k+1,\ell_{k}}A-X_{k+1,\ell_{k}}KX_{k+1,\ell_{k}}+Q\|<\varepsilon_{inn}~{}~{}\mbox{or}~{}~{}\ell_{k}>\ell_{\max}.$

Step 4. Compute the normalized residual

$\displaystyle\mbox{NRes}(X_{k+1})=\frac{\|A^{*}X_{k+1}+X_{k+1}A-X_{k+1}KX_{k+1}+Q\|_{2}}{\|A^{*}X_{k+1}\|_{2}+\|X_{k+1}A\|_{2}+\|X_{k+1}KX_{k+1}\|_{2}+\|Q\|_{2}}.$

Step 5. If $\mbox{NRes}(X_{k+1})<\epsilon_{out}$ or $k>k_{max}$, then the
approximate solution of the CARE (1) is $\tilde{X}=X_{k+1}$; else set $k=k+1$
and return to Step 2.

######  Remark 2.3.

From Algorithm 1, it can be seen that when $\omega_{k}=0$, the GADI method
naturally reduces to the ADI method in [34]. We adopt the method given in [33]
to select the desired initial matrix $X_{0}$, see Step 1 in Algorithm 1.

For large-scale Lyapunov equations at each Newton step, it is important to
control the accuracy of the solution to gain efficiency, and wish to keep the
overall fast convergence property of the Newton method. The inexact Newton
method proposed in [30] has shown a rigorous guideline for the termination of
the inner iteration while the fast local rate of convergence is retained. The
Lyapunov equation can be approximately solved such that

$\displaystyle\|(A-KX_{k})^{*}X_{k+1}+X_{k+1}(A-KX_{k})+X_{k}KX_{k}+Q\|\leq\eta_{k}\|\mathcal{R}(X_{k})\|$
$\displaystyle=\eta_{k}\|A^{*}X_{k}+X_{k}A-X_{k}KX_{k}+Q\|.$

Formally, the iteration scheme is determined by solving

$\displaystyle(A-KX_{k})^{*}X_{k+1}+X_{k+1}(A-KX_{k})+X_{k}KX_{k}+Q=R_{k}.$
(8)

The following convergence property for the inexact Newton method is referenced
in [30].

######  Theorem 2.4.

[30] For the CARE (1), assume that $(A,K)$ is stabilizable and $(A,Q)$ is
detectable, and $X_{*}\in\mathbb{C}^{n\times n}$ is the Hermitian positive
semi-definite solution. If there exist $\tilde{\delta}>0$ and $\tilde{\eta}>0$
such that for an initial matrix $X_{0}\in\mathbb{C}^{n\times n}$ meeting
$\|X_{0}-X_{*}\|\leq\tilde{\delta}$, then the iteration sequence $X_{k}$
generated by (8) converges to $X_{*}$ if the residual $R_{k}$ satisfies

$\displaystyle\|R_{k}\|\leq\eta_{k}\|A^{*}X_{k}+X_{k}A-X_{k}KX_{k}+Q\|.$

The convergence rate is linear if $\eta_{k}\in(0,\tilde{\eta}]$. It is super-
linear if $\eta_{k}\rightarrow 0$ and is quadratic if

$\eta_{k}\leq K_{\eta}\|A^{*}X_{k}+X_{k}A-X_{k}KX_{k}+Q\|$

for some $K_{\eta}>0$.

As a consequence, the inexact Newton-GADI method is summarized in Algorithm 2.
Obviously, this algorithm consists of approximating the solution in the $k$-th
Newton step for a certain tolerance $\eta_{k}>0$. It decreases when $k$
increases, e.g., set $\eta_{k}=1/(k^{3}+1)$ (see Section 3 for details).

Algorithm 2 (Inexact Newton-GADI Algorithm)

Given the matrices $A,~{}K,~{}Q$, the integer numbers $k_{\max}>1$ and
$\ell_{\max}>1$, a sequence of positive numbers $\\{\eta_{k}\\}$ such that
$\eta_{k+1}<\eta_{k}$ for all $k$, the outer iteration tolerance
$\varepsilon_{out}$.

Step 1. Compute $\beta=1+\|A\|_{\infty},~{}~{}B=-(A+\beta I)$ and the initial
matrix $X_{0}=X^{-1}$ by solving the Lyapunov equation $B^{*}X+XB-2Q=0$.

Step 2. Set $X_{k+1,0}=X_{k}$ and compute the matrices

$\displaystyle A_{k}=KX_{k}-A,~{}~{}F(X_{k})=X_{k}KX_{k}+Q.$

Step 3. Solve the GADI iteration scheme (7) to obtain
$X_{k+1}=X_{k+1,\ell_{k}}$ such that

$\displaystyle\|R_{k}\|\leq\eta_{k}\|A^{*}X_{k}+X_{k}A-X_{k}KX_{k}+Q\|~{}~{}\mbox{or}~{}~{}\ell_{k}>\ell_{\max}.$

Step 4. Compute the normalized residual

$\displaystyle\mbox{NRes}(X_{k+1})=\frac{\|A^{*}X_{k+1}+X_{k+1}A-X_{k+1}KX_{k+1}+Q\|_{2}}{\|A^{*}X_{k+1}\|_{2}+\|X_{k+1}A\|_{2}+\|X_{k+1}KX_{k+1}\|_{2}+\|Q\|_{2}}.$

Step 5. If $\mbox{NRes}(X_{k+1})<\epsilon_{out}$ or $k>k_{max}$, then the
approximate solution of the CARE (1) is $\tilde{X}=X_{k+1}$; else set $k=k+1$
and return to Step 2.

###  2.2 Convergence analysis

In this subsection, we prove that the Newton-GADI scheme is convergent. The
convergence analysis of the Inexact Newton-GADI scheme is similar and hence is
omitted. We recall the definition and property of Kronecker product.

For $A_{1}=(a_{ij})\in\mathbb{C}^{n_{A}\times m_{A}}$ and
$B_{1}\in\mathbb{C}^{n_{B}\times m_{B}}$, the Kronecker product is defined as

$\displaystyle A_{1}\otimes
B_{1}=\begin{pmatrix}a_{11}B_{1}&a_{12}B_{1}&\cdots&a_{1m_{A}}B_{1}\\\
a_{21}B_{1}&a_{22}B_{1}&\cdots&a_{2m_{A}}B_{1}\\\ \vdots&\vdots&&\vdots\\\
a_{n_{A}1}B_{1}&a_{n_{A}2}B_{1}&\cdots&a_{n_{A}m_{A}}B_{1}\end{pmatrix}\in\mathbb{C}^{(n_{A}n_{B})\times(m_{A}m_{B})}.$

The vectorization operator (abbreviated as vec): $\mathbb{C}^{m\times
n}\rightarrow\mathbb{C}^{mn}$ is

$\displaystyle\mbox{vec}(X_{1})=(x^{T}_{1},x^{T}_{2},\cdots,x^{T}_{n})^{T}\in\mathbb{C}^{mn},~{}~{}\mbox{with}~{}~{}X_{1}=(x_{1},x_{2},\cdots,x_{n})\in\mathbb{C}^{m\times
n}.$

A brief review of some properties related to Kronecker products is required.

######  Proposition 2.5.

For $A_{1},~{}B_{1},~{}C,~{}D,~{}X_{1}$ with appropriate dimensions, we have

(i) $(A_{1}\otimes B_{1})(C\otimes D)=(A_{1}C)\otimes(B_{1}D)$;

(ii) $(A_{1}\otimes B_{1})^{-1}=A_{1}^{-1}\otimes B_{1}^{-1}$;

(iii) $\Lambda(A_{1}\otimes
B_{1})=\\{\lambda\mu:~{}\lambda\in\Lambda(A_{1})~{}\mbox{and}~{}\mu\in\Lambda(B_{1})\\}$;

(iv) $(A_{1}\otimes B_{1})x=\mbox{vec}(B_{1}^{T}X_{1}A_{1})$,
$\mbox{vec}(X_{1})=x$.

Thereby, we have the following convergence conclusion for the Newton-GADI
method.

######  Theorem 2.6.

For the CARE (1), assume that $(A,K)$ is stabilizable and $(A,Q)$ is
detectable, the parameters $\alpha_{k}>0$ and $0\leq\omega_{k}<2$. Then, for
any $k=0,1,\cdots$, the iteration sequence
$\\{X_{k+1,\ell}\\}^{\infty}_{\ell=0}$ defined by (7) converges to $X_{k+1}$.

######  Proof.

Taking the operation “vec” from the GADI iteration scheme (7) yields

$\displaystyle\begin{cases}(\alpha_{k}I+I\otimes
A_{k}^{*})\mbox{vec}(X_{k+1,\ell+\frac{1}{2}})=(\alpha_{k}I-A^{T}_{k}\otimes
I)\mbox{vec}(X_{k+1,\ell})+\mbox{vec}(F(X_{k})),\\\
(\alpha_{k}I+A^{T}_{k}\otimes I)\mbox{vec}(X_{k+1,\ell+1})=[A^{T}_{k}\otimes
I-(1-\omega_{k})\alpha_{k}I]\mbox{vec}(X_{k+1,\ell})+(2-\omega_{k})\alpha_{k}\cdot\mbox{vec}(X_{k+1,\ell+\frac{1}{2}}).\end{cases}$
(9)

Denote $x=\mbox{vec}(X)$ and $f(X_{k})=\mbox{vec}(F(X_{k}))$. The iteration
scheme (9) is equivalent to

$\displaystyle\begin{cases}(\alpha_{k}I+I\otimes
A_{k}^{*})x_{k+1,\ell+\frac{1}{2}}=(\alpha_{k}I-A^{T}_{k}\otimes
I)x_{k+1,\ell}+f(X_{k}),\\\ (\alpha_{k}I+A^{T}_{k}\otimes
I)x_{k+1,\ell+1}=[A^{T}_{k}\otimes
I-(1-\omega_{k})\alpha_{k}I]x_{k+1,\ell}+(2-\omega_{k})\alpha_{k}x_{k+1,\ell+\frac{1}{2}}.\end{cases}$
(10)

As $\alpha_{k}>0$ and $A_{k}$ is stable, then $\alpha_{k}I+I\otimes A_{k}^{*}$
and $\alpha_{k}I+A^{T}_{k}\otimes I$ are nonsingular. By (10), we get

$\displaystyle x_{k+1,\ell+1}$ $\displaystyle=(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}\Big{\\{}\Big{[}A^{T}_{k}\otimes
I-(1-\omega_{k})\alpha_{k}I\Big{]}+(2-\omega_{k})\alpha_{k}(\alpha_{k}I+I\otimes
A_{k}^{*})^{-1}(\alpha_{k}I-A^{T}_{k}\otimes I)\Big{\\}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}x_{k+1,\ell}+(2-\omega_{k})\alpha_{k}(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}f(X_{k})$
$\displaystyle=(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}(\alpha_{k}I+I\otimes
A_{k}^{*})^{-1}\Big{[}\alpha_{k}^{2}I+(I\otimes A_{k}^{*})(A^{T}_{k}\otimes
I)-(1-\omega_{k})\alpha_{k}C_{k}\Big{]}x_{k+1,\ell}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+(2-\omega_{k})\alpha_{k}(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}f(X_{k})$

where $C_{k}=I\otimes A_{k}^{*}+A^{T}_{k}\otimes I$. Denote

$\displaystyle T_{k}(\alpha_{k},\omega_{k})=(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}\Big{[}\alpha_{k}^{2}I+(I\otimes
A_{k}^{*})(A^{T}_{k}\otimes I)-(1-\omega_{k})\alpha_{k}C_{k}\Big{]},$
$\displaystyle
G_{k}(\alpha_{k},\omega_{k})=(2-\omega_{k})\alpha_{k}(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1},$

thus

$\displaystyle
x_{k+1,\ell+1}=T_{k}(\alpha_{k},\omega_{k})x_{k+1,\ell}+G_{k}(\alpha_{k},\omega_{k})f(X_{k}).$
(11)

Next, we prove that the iteration matrix
$\rho(T_{k}(\alpha_{k},\omega_{k}))<1$ for $0<\alpha_{k}$ and
$0\leq\omega_{k}<2$. Note that

$\displaystyle 2\alpha_{k}C_{k}=-(\alpha_{k}I-I\otimes
A_{k}^{*})(\alpha_{k}I-A^{T}_{k}\otimes I)+(\alpha_{k}I+I\otimes
A_{k}^{*})(\alpha_{k}I+A^{T}_{k}\otimes I),$

then

$\displaystyle T_{k}(\alpha_{k},\omega_{k})$
$\displaystyle=(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}(\alpha_{k}I+I\otimes
A_{k}^{*})^{-1}\Big{[}\alpha_{k}^{2}I+(I\otimes A_{k}^{*})(A^{T}_{k}\otimes
I)-(1-\omega_{k})\alpha_{k}C_{k}\Big{]}$
$\displaystyle=(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}(\alpha_{k}I+I\otimes
A_{k}^{*})^{-1}\Big{[}(\alpha_{k}I-I\otimes
A_{k}^{*})(\alpha_{k}I-A^{T}_{k}\otimes I)+\omega_{k}\alpha_{k}C_{k}\Big{]}$
$\displaystyle=\frac{1}{2}\Big{\\{}(2-\omega_{k})(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes
A_{k}^{*})(\alpha_{k}I-A^{T}_{k}\otimes I)+\omega_{k}I\Big{\\}}$
$\displaystyle=\frac{1}{2}[(2-\omega_{k})T_{k}(\alpha_{k})+\omega_{k}I]$ (12)

with

$\displaystyle T_{k}(\alpha_{k})=(\alpha_{k}I+A^{T}_{k}\otimes
I)^{-1}(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes
A_{k}^{*})(\alpha_{k}I-A^{T}_{k}\otimes I).$ (13)

Using (12), it is evident that

$\displaystyle\lambda_{i}(T_{k}(\alpha_{k},\omega_{k}))=\frac{1}{2}[(2-\omega_{k})\lambda_{i}(T_{k}(\alpha_{k}))+\omega_{k}](i=1,2,\cdots,n),$

and

$\displaystyle\rho(T_{k}(\alpha_{k},\omega_{k}))\leq\frac{1}{2}[(2-\omega_{k})\rho(T_{k}(\alpha_{k}))+\omega_{k}].$

Obviously, the matrix $T_{k}(\alpha_{k})$ defined in (13) is similar to

$\displaystyle\hat{T}_{k}(\alpha_{k})=(\alpha_{k}I+I\otimes
A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes A_{k}^{*})(\alpha_{k}I-A^{T}_{k}\otimes
I)(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}$

through the matrix $\alpha_{k}I+A^{T}_{k}\otimes I$. Thus we can obtain

$\displaystyle\rho(T_{k}(\alpha_{k}))$
$\displaystyle\leq\|(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes
A_{k}^{*})\|_{2}\cdot\|(\alpha_{k}I-A^{T}_{k}\otimes
I)(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}\|_{2}$
$\displaystyle\leq\|(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes
A_{k}^{*})\|_{2}\cdot\|(\alpha_{k}I-A^{T}_{k}\otimes
I)(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}\|_{2}$
$\displaystyle=\|A^{L}_{k}\|_{2}\|A^{R}_{k}\|_{2},$ (14)

where $A^{L}_{k}=(\alpha_{k}I+I\otimes A_{k}^{*})^{-1}(\alpha_{k}I-I\otimes
A_{k}^{*})$ and $A^{R}_{k}=(\alpha_{k}I-A^{T}_{k}\otimes
I)(\alpha_{k}I+A^{T}_{k}\otimes I)^{-1}$. Then, for
$\bm{x}\in\mathbb{C}^{n^{2}\times 1}$, we get

$\displaystyle\|A^{L}_{k}\|^{2}_{2}$
$\displaystyle=\max_{\|\bm{x}\|_{2}=1}\frac{\|(\alpha_{k}I-I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}}{\|(\alpha_{k}I+I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}}$
$\displaystyle=\max_{\|\bm{x}\|_{2}=1}\frac{\|(I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}-\alpha_{k}\bm{x}^{*}(I\otimes A_{k}+I\otimes
A_{k}^{*})\bm{x}+\alpha^{2}_{k}}{\|(I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}+\alpha_{k}\bm{x}^{*}(I\otimes A_{k}+I\otimes
A_{k}^{*})\bm{x}+\alpha^{2}_{k}}$
$\displaystyle\leq\max_{\|\bm{x}\|_{2}=1}\frac{\|(I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(I\otimes
A_{k}))+\alpha^{2}_{k}}{\|(I\otimes
A_{k}^{*})\bm{x}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(I\otimes
A_{k}))+\alpha^{2}_{k}}$ $\displaystyle\leq\frac{\|I\otimes
A_{k}^{*}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(I\otimes
A_{k}))+\alpha^{2}_{k}}{\|I\otimes
A_{k}^{*}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(I\otimes
A_{k}))+\alpha^{2}_{k}}$ $\displaystyle=\frac{\|I\otimes
A_{k}^{*}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}{\|I\otimes
A_{k}^{*}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}.$

Analogously, we obtain

$\displaystyle\|A^{R}_{k}\|^{2}_{2}\leq\frac{\|I\otimes
A_{k}^{T}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}{\|I\otimes
A_{k}^{T}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}.$

Since $A_{k}$ is stable and $\alpha_{k}>0$, by (14), it follows that
$\|A^{L}_{k}\|_{2}<1$ and $\|A^{R}_{k}\|_{2}<1$. Therefore,

$\displaystyle\rho(T_{k}(\alpha_{k}))<1,~{}~{}\mbox{and}~{}~{}\rho(T_{k}(\alpha_{k},\omega_{k}))\leq\frac{1}{2}[(2-\omega_{k})\rho(T_{k}(\alpha_{k}))+\omega_{k}]<1,$

which implies that the iteration scheme (11) is convergent. The proof is
completed.

###  2.3 Parameter selection

In this subsection, we first compare the convergence rate between the Newton-
GADI and Newton-ADI methods. That is, the relation between
$\rho(T_{k}(\alpha_{k},\omega_{k}))$ and $\rho(T_{k}(\alpha_{k}))$ defined in
Section 2.2 is discussed. Then we present a practical method to select the
quasi-optimal parameter.

######  Theorem 2.7.

Assume that the eigenvalues of $A_{k}$ have positive real parts. For
$j_{1},~{}j_{2}\in\\{1,\cdots,n^{2}\\}$, denote

$\rho(T_{k}(\alpha_{k},\omega_{k}))=|\zeta_{j_{1}}|=|a_{j_{1}}+i\,b_{j_{1}}|,~{}~{}\rho(T_{k}(\alpha_{k}))=|\eta_{j_{2}}|=|c_{j_{2}}+i\,d_{j_{2}}|.$

(i) When $|\eta_{j_{2}}|^{2}\leq c_{j_{2}}$, we have

$\displaystyle\rho(T_{k}(\alpha_{k}))<\rho(T_{k}(\alpha_{k},\omega_{k}))<1.$

(ii) When $|\eta_{j_{1}}|^{2}>c_{j_{1}}$ and
$0<\omega_{k}<\frac{4(|\eta_{j_{1}}|^{2}-c_{j_{1}})}{(1-c_{j_{1}})^{2}+d_{j_{1}}^{2}}<2,$
we have

$\displaystyle\rho(T_{k}(\alpha_{k},\omega_{k}))<\rho(T_{k}(\alpha_{k}))<1.$

######  Proof.

For $\zeta_{j}=a_{j}+ib_{j}\in\Lambda(T_{k}(\alpha_{k},\omega_{k}))$,
$\eta_{j}=c_{j}+id_{j}\in\Lambda(T_{k}(\alpha_{k}))(j=1,2,\cdots,n)$, then

$\displaystyle 2|\zeta_{j}|=|(2-\omega_{k})\eta_{j}+\omega_{k}|.$

Therefore,

$\displaystyle
4|\zeta_{j}|^{2}=|(2-\omega_{k})\eta_{j}+\omega_{k}|^{2}=[(2-\omega_{k})c_{j}+\omega_{k}]^{2}+[(2-\omega_{k})d_{j}]^{2}=(2-\omega_{k})^{2}|\eta_{j}|^{2}+\omega_{k}^{2}+2(2-\omega_{k})\omega_{k}c_{j},$

i.e.,

$\displaystyle 4|\zeta_{j}|^{2}-4|\eta_{j}|^{2}$
$\displaystyle=(\omega_{k}^{2}-4\omega_{k})|\eta_{j}|^{2}+\omega_{k}^{2}+2(2-\omega_{k})\omega_{k}c_{j}$
$\displaystyle=(1-2c_{j}+|\eta|^{2})\omega_{k}^{2}-4(|\eta_{j}|^{2}-c_{j})\omega_{k}$
$\displaystyle=[(1-c_{j})^{2}+d_{j}^{2}]\omega_{k}^{2}-4(|\eta_{j}|^{2}-c_{j})\omega_{k}.$

Note that

$\rho(T_{k}(\alpha_{k},\omega_{k}))=|\zeta_{j_{1}}|,~{}~{}\rho(T_{k}(\alpha_{k}))=|\eta_{j_{2}}|.$

(i) When $|\eta_{j_{2}}|^{2}\leq c_{j_{2}}$, we have

$\displaystyle\rho(T_{k}(\alpha_{k}))<\rho(T_{k}(\alpha_{k},\omega_{k}))<1.$

(ii) When $|\eta_{j_{1}}|^{2}>c_{j_{1}}$ and
$0<\omega_{k}<\frac{4(|\eta_{j_{1}}|^{2}-c_{j_{1}})}{(1-c_{j_{1}})^{2}+d_{j_{1}}^{2}}<2,$
we have

$\displaystyle|\zeta_{j_{1}}|<|\eta_{j_{1}}|,$

which means that
$\rho(T_{k}(\alpha_{k},\omega_{k}))<\rho(T_{k}(\alpha_{k}))<1$.

######  Remark 2.8.

By Theorem 2.7, if $|\eta_{j_{1}}|^{2}>c_{j_{1}}$ and
$0<\omega_{k}<\frac{4(|\eta_{j_{1}}|^{2}-c_{j_{1}})}{(1-c_{j_{1}})^{2}+d_{j_{1}}^{2}}<2,$
the Newton-GADI has a faster convergence rate than Newton-ADI proposed in
[34]. This is verified by numerical examples in Section 3.

Next, we provide a method for selecting theoretical quasi-parameters by
minimizing $\rho(T_{k}(\alpha_{k},\omega_{k}))$.

######  Theorem 2.9.

Assume that the eigenvalues of $A_{k}$ have positive real parts, set

$\displaystyle\nu_{n}=\max_{j}\sigma_{j}(A_{k}),~{}~{}\min_{j}Re(\lambda_{j}(A_{k}))=Re(\mu_{1})=Re(a_{\mu}+ib_{\mu})=a_{\mu}.$
(15)

Then the quasi-optimal parameter $\alpha_{k}$ can be obtained by

$\displaystyle\alpha_{k}^{*}=\mbox{arg}\min_{\alpha_{k}}\frac{\nu_{n}^{2}-2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}{\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}=\nu_{n}.$
(16)

######  Proof.

Using the bound defined by (14), we have

$\displaystyle\rho^{2}(T_{k}(\alpha_{k}))$
$\displaystyle\leq\|A_{L}\|^{2}_{2}\|A_{R}\|^{2}_{2}$
$\displaystyle\leq\frac{\|I\otimes
A_{k}^{*}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}{\|I\otimes
A_{k}^{*}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}\cdot\frac{\|I\otimes
A_{k}^{T}\|^{2}_{2}-2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}{\|I\otimes
A_{k}^{T}\|^{2}_{2}+2\alpha_{k}\min_{i}Re(\lambda_{i}(A_{k}))+\alpha^{2}_{k}}$
$\displaystyle=\frac{\nu_{n}^{2}-2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}{\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}\cdot\frac{\nu_{n}^{2}-2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}{\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha^{2}_{k}},$

i.e.,

$\displaystyle\rho(T_{k}(\alpha_{k}))\leq\frac{\nu_{n}^{2}-2\alpha_{k}a_{\mu}+\alpha^{2}_{k}}{\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha^{2}_{k}},$

where $\nu_{n}$ and $a_{\mu}$ are given in (15). Denote

$\displaystyle\Psi(\alpha_{k})=\frac{\nu_{n}^{2}-2\alpha_{k}a_{\mu}+\alpha_{k}^{2}}{\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha_{k}^{2}}.$

Then

$\displaystyle\frac{\partial\Psi(\alpha_{k})}{\partial\alpha_{k}}=\frac{4a_{\mu}(\alpha_{k}^{2}-\nu_{n}^{2})}{(\nu_{n}^{2}+2\alpha_{k}a_{\mu}+\alpha_{k}^{2})^{2}}.$

It is easy to show that $\Psi(\alpha_{k})$ gets to its minimum when
$\alpha_{k}=\nu_{n}$, i.e.,

$\alpha_{k}^{*}=\mbox{arg}\min\limits_{\alpha_{k}}\Psi(\alpha_{k})=\nu_{n},$

where $\nu_{n}$ depends on the matrix $A_{k}$.

######  Remark 2.10.

In terms of Theorem 2.9, we usually select $\omega_{k}=1$ and $\alpha_{k}^{*}$
given in (16) as the quasi-optimal parameters. It can be seen that this
parameter choice is reasonable using the residual analysis in Section 3.

##  3 Numerical experiments

In this section, we use several examples of the quadratic optimal control to
show the numerical feasibility and effectiveness of the Newton-GADI and
inexact Newton-GADI (InNewton-GADI) algorithms. The whole process is performed
on a computer with Intel Core 3.20GHz CPU, 4.00GB RAM and MATLAB R2017a.
IT(inn) and IT(out) represent the number of inter iteration steps and outer
iteration steps, respectively. Denote IT(cumul) by the sum of inter and outer
iteration steps, IT(ave) by the number of the GADI iterative steps required by
the average Newton step, and CPU by the computing time. The notation
$\mbox{Tri}(c,a,b)_{n\times n}$ demonstrates the $n\times n$ size tridiagonal
matrix

$\displaystyle\begin{pmatrix}a&b&&&\\\ c&a&b&&\\\ &c&\ddots&\ddots&\\\
&&\ddots&\ddots&b\\\ &&&c&a\end{pmatrix}.$

In our numerical experiment, set the outer iteration tolerance
$\varepsilon_{out}=1\times 10^{-8},~{}\varepsilon_{inn}=1\times 10^{-8}$,
$\ell_{\max}=1000,~{}k_{\max}=1000$. In the inexact Newton-GADI algorithm, set
$\eta_{k}=1/(k^{3}+1)$. Denote the numerical Hermitian semi-definite solution
as $\tilde{X}$, the finally normalized residual as

$\displaystyle\mbox{NRes}(\tilde{X})=\frac{\|A^{*}\tilde{X}+\tilde{X}A-\tilde{X}K\tilde{X}+Q\|_{2}}{\|A^{*}\tilde{X}\|_{2}+\|\tilde{X}A\|_{2}+\|\tilde{X}K\tilde{X}\|_{2}+\|Q\|_{2}},$

and the error of the inter iteration in $k$-th Newton iteration step as

$\displaystyle
R(X_{k})=A_{k}^{*}\tilde{X}^{in}_{k}+\tilde{X}^{in}_{k}A_{k}-X_{k}KX_{k}-Q,$

with $\tilde{X}^{in}_{k}$ is the numerical solution of Lyapunov equation in
$k$-th Newton iteration step.

When $\omega_{k}=0$, the Newton-GADI reduces to Newton-ADI proposed in [34].
It can be seen that the Newton-GADI converges faster than Newton-ADI in the
following examples.

######  Example 3.1.

Consider the linear time-invariant system of the form

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)+B\bm{u}(t)=\begin{pmatrix}-2+10\,i&0&-1\\\
0&-1+10\,i&0\\\ -1&-1&-2\,i\end{pmatrix}\bm{p}(t)+\begin{pmatrix}-2&0&-1\\\
0&-1&-1\\\ 1&0&-2\end{pmatrix}\bm{u}(t).$

We wish to minimize the cost function $J$ defined in (2). The cost penalty
matrices are

$\displaystyle R=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&4\end{pmatrix},~{}Q=\begin{pmatrix}0&0&0\\\ 0&1&0\\\ 0&0&5\end{pmatrix}.$

Then

$K=BR^{-1}B^{T}=\begin{pmatrix}17/4&1/4&-3/2\\\ 1/4&5/4&1/2\\\
-3/2&1/2&2\end{pmatrix}.$

We can verify that $K$ is positive definite and $Q$ is semi-positive definite.
The pair $(A,K)$ is stabilize as

$\displaystyle\Lambda(A-KI_{3})=\\{-1.9936-1.9179\,i,-6.2704+9.9858\,i,-2.2360+9.9321\,i\\}.$

The pair $(A,Q)$ is detectable as

$\displaystyle\Lambda(Q^{*}-A^{*}F)=\\{-2.0000+0.0000\,i,-3.5000-1.6583\,i,-3.5000+1.6583\,i\\}$

with

$\displaystyle F=\begin{pmatrix}-0.2785-0.5729\,i&0&0.2387-0.3660\,i\\\
-0.1651-0.0121\,i&-0.0297-0.2970\,i&0.2773-0.3460\,i\\\
0.2865-1.6393\,i&0&0.1830+3.1194\,i\end{pmatrix}.$

The numerical Hermitian semi-definite solution is

$\displaystyle\tilde{X}=\begin{pmatrix}-0.0020+0.0120i&-0.0034+0.0121i&-0.1269-0.0611i\\\
-0.0034+0.0121i&0.0001+0.0618i&-0.1251-0.0780i\\\
-0.1269-0.0611i&-0.1251-0.0780i&1.1794-0.9689i\end{pmatrix}.$

We know that the approximate solution of the Lyapunov equation can be obtained
by using the inexact algorithm at each Newton step. When the error generated
by the inexact algorithm is very large, it will affect the convergence rate in
next step. Figure 1 displays the error $\|R(X_{k})\|_{2}$ of the inexact inter
iterative algorithm in each Newton step. Moreover, due to the large error
$\|R(X_{2})\|_{2}$ caused by the GADI in the second iteration step, from
Figure 2, we know that the inner iteration step exceeds the maximum and is not
convergent in the third Newton step. Therefore, we reset
$\eta_{k}=\frac{1}{k^{4}+1}$ in the inexact Newton-ADI algorithm, see Figure
3.

Figure 1: Comparing the inner iteration accuracy of two inexact algorithm with
$\eta_{k}=\frac{1}{k^{3}+1}$. Figure 2: Comparing the number of inner
iteration steps of Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI at
each Newton step with $\eta_{k}=\frac{1}{k^{3}+1}$. Figure 3: Comparing two
inexact algorithms.

Figure 4 shows the convergence process of the four algorithms, including the
Newton-ADI ($\omega=0$), the Newton-GADI, the InNewton-ADI ($\omega=0$) and
the InNewton-GADI methods. It is not difficult to find that these algorithms
have almost the same convergence rate.

Figure 4: Comparing the convergence rate of the outer iteration for Newton-
ADI, Newton-GADI, InNewton-ADI and InNewton-GADI.

Figure 5 shows that the number of the GADI steps required for each Newton
step. It is easy to see that the inexact Newton-GADI method has the least
number of inter iteration steps among the four algorithms. The total number of
steps in the convergence process is the least compared other methods, see
Table 1 for details.

Figure 5: Comparing the number of inner iteration steps of Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI at each Newton step. Table 1: Numerical results of Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI. | $\omega_{k}^{*}$ | $NRes(\tilde{X})$ | IT(out) | IT(ave) | IT(cumul)
---|---|---|---|---|---
Newton-ADI | 0 | 8.1384e-09 | 19 | 48.21 | 916
Newton-GADI | 1 | 8.2030e-09 | 19 | 30.16 | 573
InNewton-ADI | 0 | 8.4516e-09 | 26 | 26.32 | 500
InNewton-GADI | 1 | 6.6070e-09 | 20 | 14.30 | 286

######  Example 3.2.

Consider the linear time-invariant system of the form

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)+B\bm{u}(t)=\begin{pmatrix}0&-1&0&0\\\
1&0&-1&0\\\ 0&1&0&-1\\\
0&0&1&0\end{pmatrix}\bm{p}(t)+10^{-3}\times\begin{pmatrix}3&-50&1&2\\\
1&-3&-2&1\\\ -3&1&3&4\\\ 3&-1&-4&3\end{pmatrix}\bm{u}(t).$

We wish to minimize the cost function $J$ defined in (2). The cost penalty
matrices are

$\displaystyle R=I_{4},~{}Q=\begin{pmatrix}0.0025&0&0&0\\\
0&0.0111&0.0025&0\\\ 0&0.0025&1.0006&0.0200\\\
0&0&0.0200&0.0004\end{pmatrix},~{}~{}\mbox{therefore}~{}~{}K=BB^{T}.$

It is easy see that $K$ is semi-positive definite and $Q$ is semi-positive
definite. The pair $(A,K)$ is stabilize as

$\displaystyle\Lambda(A-KI_{4})=\\{-0.0004+1.6180\,i,-0.0004-1.6180\,i,-0.0009+0.6180\,i,-0.0009-0.6180\,i\\}.$

The pair $(A,Q)$ is detectable as

$\displaystyle\Lambda(Q^{*}-A^{*}F)=\\{-0.0013+50.0608\,i,-0.0013-50.0608\,i,-0.0009+7.6862\,i,-0.0009-7.6862i\,\\}$

with

$\displaystyle F=\begin{pmatrix}17.5243&0.3931&-8.2505&-0.3925\\\
0.3931&25.8550&0.3901&-8.2707\\\ -8.2505&0.3901&25.8338&-0.0035\\\
-0.3925&-8.2707&-0.0035&17.5505\end{pmatrix}.$

The numerical Hermitian semi-definite solution is

$\displaystyle\tilde{X}=\begin{pmatrix}17.4818&0.3916&-8.2438&-0.3924\\\
0.3916&25.8038&0.3901&-8.2638\\\ -8.2438&0.3901&25.7818&-0.0035\\\
-0.3924&-8.2638&-0.0035&17.5055\end{pmatrix}.$

Figure 6 displays that the error of the inexact inter iterative algorithm with
the increase of Newton step. Figure 7 shows the convergence process of the
four algorithms and reports that these algorithms have almost the same
convergence rate.

Figure 6: Comparing the inner iteration accuracy of two inexact algorithm with
$\eta_{k}=\frac{1}{k^{4}+1}$. Figure 7: Comparing the convergence rate of the
outer iteration for Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI
with $\eta_{k}=\frac{1}{k^{4}+1}$ .

Figure 8 reports the number of the GADI steps required for each Newton step.
It can be seen that the InNewton-GADI algorithm has a clear advantage over
other algorithms.

Figure 8: Comparing the number of inner iteration steps of Newton-ADI, Newton-
GADI, InNewton-ADI and InNewton-GADI at each Newton step.

Table 2 displays the numerical results of these four algorithms. It is not
difficult to find that this example is ill-conditioned. Some existing ADI
algorithms may be convergent slowly or not convergent, while the GADI scheme
proposed in this paper can accelerate the convergence process.

Table 2: Numerical results of Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI. | $\omega_{k}^{*}$ | $NRes(\tilde{X})$ | IT(out) | IT(ave) | IT(cumul)
---|---|---|---|---|---
Newton-ADI | 0 | 9.5234e-09 | 15 | 204.40 | 3066
Newton-GADI | 0.015 | 9.5846e-09 | 15 | 193.80 | 2907
InNewton-ADI | 0 | 1.2046e-10 | 15 | 175.87 | 2638
InNewton-GADI | 0.015 | 3.5208e-09 | 15 | 153.00 | 2295

######  Example 3.3.

Consider the linear time-invariant system of the form

$\displaystyle\dot{\bm{p}}(t)=A\bm{p}(t)+B\bm{u}(t)$

where

$\displaystyle A=\mbox{Tri}(-1-r,-4+8\,i,-1+r)_{n\times
n},~{}b=(1,0,\cdots,0)\in\mathbb{R}^{1\times n},~{}~{}B=(b^{T},I_{n}),$

with $r=\dfrac{1}{2n+2}$. To go further, we need to minimize the cost function
$J$ defined in (2). The cost penalty matrices are $R=I_{n}$, $K=BR^{-1}B^{T}$,
$Q=c^{T}c$ with $c=(\dfrac{1}{\sqrt{10}},0,\cdots,0)\in\mathbb{R}^{1\times
n}$.

When $n=128$, Figure 9 displays that the error of the inexact inter iterative
algorithm with the increase of Newton step.

Figure 9: Comparing the inner iteration accuracy of two inexact algorithm with
$\eta_{k}=\frac{1}{k^{3}+1}$ and $n=128$.

Figure 10 shows the convergence process of the four algorithms and reports
that these algorithms have almost the same convergence rate. Moreover, it is
easy to find that the convergence rate of the iteration is slow in the
beginning. From the fifth iteration step, these algorithms begin to converge
rapidly.

Figure 10: When $n=128$, comparing the convergence rate of the outer iteration
for Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI with
$\eta_{k}=\frac{1}{k^{3}+1}$ .

Figure 11 reports the number of the GADI steps required for each Newton step.
It can be seen that the InNewton-GADI algorithm has a clear advantage over
other algorithms.

Figure 11: When $n=128$, comparing the number of inner iteration steps of
Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI at each Newton step.

Table 3 displays the numerical results of these four algorithms when
$n=64,~{}n=128,~{}n=512$ and $n=1024$. We can see that as $n$ increases, the
number of iterative steps and the finally normalized residual
$NRes(\tilde{X})$ of each algorithm are almost the same, while the CPU time
increases gradually.

Table 3: Numerical results of Newton-ADI, Newton-GADI, InNewton-ADI and InNewton-GADI as $n$ increases. $n$ | Algorithm | $\omega_{k}^{*}$ | $NRes(\tilde{X})$ | IT(out) | IT(ave) | IT(cumul) | CPU(s)
---|---|---|---|---|---|---|---
64 | Newton-ADI | 0 | 6.3978e-09 | 11 | 12.73 | 140 | 0.2140
Newton-GADI | 1 | 9.4256e-09 | 7 | 14.71 | 103 | 0.1320
InNewton-ADI | 0 | 1.1340e-09 | 8 | 8.13 | 65 | 0.1070
InNewton-GADI | 1 | 4.2259e-09 | 8 | 4.50 | 36 | 0.0830
128 | Newton-ADI | 0 | 8.4557e-09 | 9 | 15.44 | 139 | 0.7540
Newton-GADI | 1 | 9.4644e-09 | 7 | 14.72 | 103 | 0.5400
InNewton-ADI | 0 | 1.1531e-09 | 8 | 8.13 | 65 | 0.5150
InNewton-GADI | 1 | 4.2424e-09 | 8 | 4.50 | 36 | 0.3680
512 | Newton-ADI | 0 | 8.4857e-09 | 9 | 15.44 | 139 | 22.1390
Newton-GADI | 1 | 9.4771e-09 | 7 | 14.72 | 103 | 16.9060
InNewton-ADI | 0 | 1.1596e-09 | 8 | 8.13 | 65 | 15.5520
InNewton-GADI | 1 | 4.2480e-09 | 8 | 4.50 | 36 | 10.8330
1024 | Newton-ADI | 0 | 8.4872e-09 | 9 | 15.44 | 139 | 256.4320
Newton-GADI | 1 | 9.4777e-09 | 7 | 14.72 | 103 | 201.9130
InNewton-ADI | 0 | 1.1600e-09 | 8 | 8.13 | 65 | 159.0060
InNewton-GADI | 1 | 4.2484e-09 | 8 | 4.50 | 36 | 113.2540

##  4 Conclusion

In this paper, we have proposed the Newton-GADI algorithm to solve the complex
CARE. Moreover, by inaccurately solving a Lyapunov equation in each Newton
step, the inexact Newton-GADI algorithm is presented. Finally, we have offered
numerical experiments to compare the effectiveness. It can be seen that the
Newton-GADI method has a faster inter iterative convergence rate than the
Newton-ADI method, while the inexact Newton-GADI method is more efficient.
Therefore, the Newton-GADI and inexact Newton-GADI methods could be valid and
attractive algorithms to solve the ill-conditioned Riccati equation. In the
future, we will further study parameter selection and solving the low-rank
large-scale CARE.

## References

  * [1] P. Lancaster, L. Rodman, Algebraic Riccati Equations, The Clarendon Press: Oxford, 1995.
  * [2] I. Petersen, Disturbance attenuation and $H^{\infty}$-optimization: a design method based on the algebraic Riccati equation, IEEE Transactions on Automatic Control, 32: 427-429, 1987.
  * [3] G. Tadmor, Worst-case design in the time domain: the maximum principle and the standard $H_{\infty}$ problem, Mathematics of Control, Signals, and Systems, 3: 301-324, 1990.
  * [4] D. Williams, A “potential-theoretic” note on the quadratic Wiener-Hopf equation for $Q$-matrices, In Seminar on Probability XVI, Lecture Notes in Mathematics, 920: 91-94, 1982.
  * [5] B. De Moor, J. David, Total linear least squares and the algebraic Riccati equation, Systems and Control Letters, 18: 329-337, 1992.
  * [6] K. Clancey, I. Gohberg, Factorization of matrix functions and singular integral operators, In Operator Theory: Advances and Applications, vol. 3. Birkh$\ddot{\mbox{a}}$user Verlag: Basel, 1981.
  * [7] I. Gohberg, M. Kaashoek, An inverse spectral problem for rational matrix functions and minimal divisibility, Integral Equations Operator Theory, 10: 437-465, 1987.
  * [8] J. Roberts, Linear model reduction and solution of the algebraic Riccati equation by use of the sign function, International Journal of Control, 32: 677-687, 1980.
  * [9] A. Laub, Invariant subspace methods for the numerical solution of Riccati equations, In The Riccati Equation, Bittanti S, Laub AJ, Willems JC (eds), Springer: Berlin, 1991.
  * [10] J. Juang, W. W. Lin, Nonsymmetric algebraic Riccati equations and Hamiltonian-like matrices, SIAM Journal on Matrix Analysis and Applications, 20: 228-243, 1999.
  * [11] J. Bellon, Riccati equations in optimal control theory, Thesis, Georgia State University, 2008.
  * [12] R. E. Kalman, New methods and results in linear prediction and estimation theory, RIAS Rep., 61: 99-119, 1961.
  * [13] M. L. Ni, A note on the maximum solutions of Riccati equations, Automatica, 27(6): 1059-1060, 1991.
  * [14] W. M. Wonham, On a matrix Riccati equation of stochastic control, SIAM Journal on Control and Optimization, 6: 681-697, 1968.
  * [15] V. Kucera, On nonnegative definite solutions to matrix quadratic equations, Automatica, 8: 413-423, 1972.
  * [16] H. Kano, Existence condition of positive-definite solutions for algebraic matrix Riccati equations, Automatica, 23: 393-397, 1987.
  * [17] D. A. Bini, B. Iannazzo, B. Meini, Numerical solution of algebraic Riccati equations, SIAM, 2012.
  * [18] R. Granat, B. Kagstrom, D. Kressner, A parallel Schur method for solving continuous-time algebraic Riccati equations, IEEE International Conference on Computer-aided Control Systems, IEEE, 2008.
  * [19] R. Byers, Solving the algebraic Riccati equation with the matrix sign function, Linear Algebra Appl., 85: 267-279, 1987.
  * [20] E. K. W. Chu, H. Y. Fan, W. W. Lin, C. S. Wang, A structure-preserving doubling algorithm for periodic discrete-time algebraic Riccati equations, International Journal of Control, 77: 767-788, 2004.
  * [21] P. Benner, J. Saak, Numerical solution of large and sparse continuous time algebraic matrix Riccati and Lyapunov equations: a state of the art survey, 36: 32-52, 2013.
  * [22] A. Bouhamidi, M. Hached, K. Jbilou, A preconditioned block Arnoldi method for large scale Lyapunov and algebraic Riccati equations, Kluwer Academic Publishers, 2016.
  * [23] D. Kleinman, On an iterative technique for Riccati equation computations, IEEE Transactions on Automatic Control, 13: 114-115, 1968.
  * [24] H. Banks, K. Ito, A numerical algorithm for optimal feedback gains in high dimensional linear quadratic regulator problems, SIAM Journal on Control and Optimization, 29: 499-515, 1991.
  * [25] K. Morris, C. Navasca, Solution of algebraic Riccati equations arising in control of partial differential equations, in Control and Boundary Analysis, J. P. Zolesio and J. Cagnol, eds., Lecture Notes in Pure Appl. Math. 240, CRC Press, 257-280, 2005.
  * [26] K. Morris, C. Navasca, Iterative solution of algebraic Riccati equations for damped systems, 2006 45th IEEE Conference on Decision and Control, 2436-2440, 2006.
  * [27] P. Benner, H. Mena, J. Saak, On the parameter selection problem in the Newton-ADI iteration for large-scale Riccati equations, Electronic Trans. Numer. Anal., 29: 136-149, 2008.
  * [28] P. Benner, J. Saak, A Galerkin-Newton-ADI method for solving large-scale algebraic Riccati equations, Technical Report SPP1253-090, DFG Priority Programme 1253 Optimization with Partial Differential Equations, Erlangen, Germany, January 2010.
  * [29] M. Dehghan, A. Shirilord, On the Hermitian and skew-Hermitian splitting-like iteration approach for solving complex continuous-time algebraic Riccati matrix equation, Applied Numerical Mathematics, 170: 109-127, 2021.
  * [30] F. Feitzinger, T. Hylla, E. W. Sachs, Inexact Kleinman-Newton method for Riccati equations, SIAM Journal on Matrix Analysis and Applications, 31: 272-288, 2009.
  * [31] P. Benner, M. Heinkenschloss, J. Saak, H. K. Weichelt, An inexact low-rank Newton-ADI method for large-scale algebraic Riccati equations, Applied Numerical Mathematics, 108: 125-142, 2016.
  * [32] K. Jiang, X. Su, J. Zhang, A general alternating-direction implicit framework with Gaussian process regression parameter prediction for large sparse linear systems, 2021, https://arxiv.org/abs/2109.12249.
  * [33] S. F. Xu. Matrix computation in control theory, Higher Education Press, 2011 (in Chinese).
  * [34] E. L. Wachspress, Iterative solution of the Lyapunov matrix equation, Applied Mathematics Letters, 1: 87-90, 1988.


11institutetext: Department of Radiology and Biomedical Research Imaging
Center (BRIC), University of North Carolina at Chapel Hill, NC 27599, USA
22institutetext: Department of Oral and Maxillofacial Surgery, Houston
Methodist Hospital, TX 77030, USA
33institutetext: Department of Surgery (Oral and Maxillofacial Surgery), Weill
Medical College, Cornell University, NY 10065, USA 44institutetext: Oral and
Craniomaxillofacial Surgery at Shanghai Ninth Hospital, Shanghai Jiaotong
University College of Medicine, Shanghai 200011, China 44email:
<EMAIL_ADDRESS>44email<EMAIL_ADDRESS>
# A Self-Supervised Deep Framework for Reference Bony Shape Estimation in
Orthognathic Surgical Planning

Deqiang Xiao 11 Hannah Deng 22 Tianshu Kuang 22 Lei Ma 11 Qin Liu 11 Xu Chen
11 Chunfeng Lian 11 Yankun Lang 11 Daeseung Kim 22 Jaime Gateno 2233 Steve
Guofang Shen 44 Dinggang Shen 11 Pew-Thian Yap, ✉ 11 James J. Xia, ✉ 2233

###### Abstract

Virtual orthognathic surgical planning involves simulating surgical
corrections of jaw deformities on 3D facial bony shape models. Due to the lack
of necessary guidance, the planning procedure is highly experience-dependent
and the planning results are often suboptimal. A reference facial bony shape
model representing normal anatomies can provide an objective guidance to
improve planning accuracy. Therefore, we propose a self-supervised deep
framework to automatically estimate reference facial bony shape models. Our
framework is an end-to-end trainable network, consisting of a simulator and a
corrector. In the training stage, the simulator maps jaw deformities of a
patient bone to a normal bone to generate a simulated deformed bone. The
corrector then restores the simulated deformed bone back to normal. In the
inference stage, the trained corrector is applied to generate a patient-
specific normal-looking reference bone from a real deformed bone. The proposed
framework was evaluated using a clinical dataset and compared with a state-of-
the-art method that is based on a supervised point-cloud network. Experimental
results show that the estimated shape models given by our approach are
clinically acceptable and significantly more accurate than that of the
competing method.

###### Keywords:

Orthognathic surgical planning Shape estimation Point-cloud network.

##  1 Introduction

Orthognathic surgery requires a detailed surgical plan to correct jaw
deformities. During computer-aided orthognathic surgical planning, a patient’s
head computed tomography (CT) or cone-beam computed tomography (CBCT) scan is
acquired to generate three-dimensional (3D) craniomaxillofacial (CMF) bony
shape models [1]. The deformed maxilla and mandible (i.e., “jaw” in Fig. 1(a))
are then osteotomized from the 3D models into multiple bony segments (Fig.
1(b)). Under the guidance of 3D cephalometric analysis [2], a surgeon
discretely moves each bony segment to a desired position, thus forming a new
normal-looking CMF skeletal anatomy. A group of surgical splints are then
designed and fabricated to correct the patient’s jaw according to the plan
during the surgery (Fig. 1(c)). Cephalometric analysis is based on a group of
specific distances and angles calculated from the anatomical landmarks on the
3D models, these measurements provide a limited guidance to the planning
process, forming an “ideal” surgical plan that still heavily relies on the
surgeon’s clinical experience and imagination of what the patient’s normal
bone should look like. Therefore, from the surgeon’s perspective, a reference
CMF bony shape model is highly desirable to accurately guide the movement of
the bony segments.

Figure 1: (a) The bony surface of a patient with jaw deformity, consisting of
a normal “midface” and a deformed “jaw” (in red). (b) The jaw is cut into
several segment pieces and moved to form a new normal-looking shape model. (c)
The patient’s facial bony shape model after treatments.

Wang et al. [3] proposed to estimate patient-specific reference bony shapes
based on sparse representation. Specifically, they represented a patient’s
bony shape with a set of landmarks, and split them into midface landmarks and
jaw landmarks. By representing the patient’s midface landmarks with a midface
dictionary from normal subjects, a set of sparse coefficients were obtained,
and applied to a normal jaw dictionary to estimate the patient’s normal jaw
landmarks. The reference bony shape model was then generated based on the
estimated landmarks. This method can work well in certain circumstances, but
is heavily dependent on a small group of pre-digitized landmarks, causing the
estimation performance to be sensitive to the locations of landmarks. To
alleviate this problem, a recent method [4] was introduced to first simulate
paired deformed-normal bony surfaces based on sparse representation, and then
train a supervised point-cloud network using the simulated data to generate
patient-specific reference bony surfaces. Although this network [4]
outperforms the method in [3], its generalizability is limited by the extent
deformity can be simulated.

We hypothesize that a fully end-to-end trainable network is capable of
estimating reference CMF bony shapes, and the network can be directly trained
using random pairs of deformed-normal bones. We propose a self-supervised deep
framework to first map jaw deformities of a patient bone to a normal bone to
generate a simulated deformed bone with a simulator network, and then restore
the simulated deformed bone back to normal with a corrector network. We train
the framework using a series of bony surfaces from patients and normal
subjects. The trained corrector is applied to patients to generate patient-
specific reference shape models. The proposed framework was evaluated using a
clinical dataset and showed superior performance over the method that is based
on supervised learning [4].

##  2 Methods

An overview of our framework is provided in Fig. 2. Given any pair of
deformed-normal bony surfaces, the simulator network takes as input the
coordinates of $N$ vertices from the two surfaces and outputs the vectors that
shift the vertices of the normal bony surface to generate the simulated
deformed bony surface with the jaw similar to the input deformed bone and the
midface same as the input normal bone. The corrector network learns a set of
displacement vectors from the vertex coordinates of the simulated deformed
bony surface and generates the corrected bony surface that matches the input
normal bony surface.

Figure 2: An overview of our self-supervised deep framework for estimating
reference CMF bony shape models.

###  2.0.1 Simulator Network:

The simulator network consists of a series of encoding, fusion, and decoding
layers to learn point features from the coordinates of vertices (Fig. 3).
There are two encoding branches with shared weights to extract point features
from the two input bony surfaces separately. Each encoding layer selects a
subset of $N_{\text{sub}}$ points from the input points via furthest point
sampling [5]. For each sampled point, its neighboring points in a 3D ball
region with radius $r$ are gathered from the input points. The point
convolution is then performed over point features of the gathered points via
PointConv [6], which is implemented based on shared multi-layer perceptron
(s-MLP) and max-pooling [7]. Following the encoding layers, the outputs of
corresponding layers from two encoding branches are fused via concatenation
and s-MLPs. The fused features are then decoded by upsampling via point
interpolation [5], grouping via 3D ball neighboring, and convolution via
PointConv. The decoding layer increases the number $N_{\text{sub}}$ of points
and reduces the dimension $C_{\text{feature}}$ of each feature vector. A set
of cascaded fusion and decoding layers are applied to generate a $N\times 3$
output vector $V_{\text{simulator}}$, which updates the positions of the
vertices of the input normal bone to obtain the simulated deformed bone.

Figure 3: The architecture of the simulator network, consisting of point-
feature encoding, fusion, and decoding layers.

We train the simulator network using loss function

$L_{\text{simulator}}=L_{\text{jaw}}+L_{\text{midface}}+\alpha
L_{\text{smooth}}+\beta L_{\text{simu\\_reg}},$ (1)

which encourages the jaw of the simulated deformed bone to be similar to the
input deformed bone based on the relative coordinates of vertices. Specially,
a set of landmarks are defined on each bony surface, the relative coordinate
vector $r_{\text{def\\_jaw}}(i,j)$ of the $i$-th jaw vertex relative to the
$j$-th landmark on the deformed bony surface is calculated as

$r_{\text{def\\_jaw}}(i,j)=c_{\text{def\\_jaw}}(i)-c_{\text{def\\_landmark}}(j),$
(2)

where $c_{\text{def\\_jaw}}(i)$ and $c_{\text{def\\_landmark}}(j)$ are
respectively the coordinate vectors of the $i$-th vertex and the $j$-th
landmark. Similarly, the relative coordinate vector for the $i$-th jaw vertex
of the simulated deformed bony surface is

$r_{\text{simu\\_jaw}}(i,j)=c_{\text{simu\\_jaw}}(i)-c_{\text{simu\\_landmark}}(j).$
(3)

We calculate $L_{\text{jaw}}$ by

$L_{\text{jaw}}=\frac{1}{N_{\text{jaw}}}\sum_{i=1}^{N_{\text{jaw}}}\sum_{j=1}^{K}\left\|r_{\text{def\\_jaw}}(i,j)-r_{\text{simu\\_jaw}}(i,j)\right\|_{2},$
(4)

where $N_{\text{jaw}}$ is the number of vertices that belong to the jaw, $K$
is the number of landmarks, and $\left\|\cdot\right\|_{2}$ is the
$\ell_{2}$-norm. $L_{\text{midface}}$ forces the network to fix the midface
during simulation. To keep the midface unchanged, i.e.,
$L_{\text{midface}}=0$, we set the displacement vectors corresponding to
midface vertices to zero in $V_{\text{simulator}}$. A smooth deformation field
is encouraged by calculating

$L_{\text{smooth}}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{\|\mathcal{N}(v_{i})\|}\sum_{v_{j}\in\mathcal{N}(v_{i})}\left\|u(i)-u(j)\right\|_{2},$
(5)

where $u(i)=c_{\text{simu}}(i)-c_{\text{norm}}(i)$ and
$u(j)=c_{\text{simu}}(j)-c_{\text{norm}}(j)$ are respectively the displacement
vectors of $v_{i}$ and $v_{j}$ on the simulated deformed bony surface.
$\mathcal{N}(v_{i})$ is the set of one-ring neighboring vertices of $v_{i}$.
$N$ is the total number of vertices on the surface. Finally, a $\ell_{2}$-norm
regularization $L_{\text{simu\\_reg}}$ of the network parameters is
incorporated to avoid overfitting.

###  2.0.2 Corrector Network:

The corrector is an encoding-decoding network (Fig. 4). Like the simulator
network, each encoding layer in the corrector performs sampling, grouping, and
convolution operations; each decoding layer performs interpolation, grouping,
and convolution. The skip connection between two corresponding layers from the
encoding and decoding streams is applied to learn point features integrating
local and global shape information. Following the decoding stream, the point
features are processed by a set of s-MLPs to obtain a $N\times 3$ vector,
which corrects the positions of vertices of the simulated deformed bone to get
a corrected bone.

Figure 4: The architecture of the corrector network with an encoding-decoding
structure.

The loss function for the corrector network is

$L_{\text{corrector}}=\frac{1}{N}\sum_{i=1}^{N}\left\|c_{\text{correct}}(i)-c_{\text{norm}}(i)\right\|_{2}+\lambda
L_{\text{correct\\_reg}},$ (6)

where $c_{\text{correct}}(i)$ and $c_{\text{norm}}(i)$ are the coordinate
vectors of the $i$-th vertex on the corrected and input normal bony surface,
respectively, $L_{\text{correct\\_reg}}$ enforces $\ell_{2}$-norm
regularization on the network parameters.

###  2.0.3 Network Training and Inference:

To train the proposed networks, vertex-wise correspondences of surfaces are
established by non-rigidly matching a template surface with all bony surfaces
in the training set. First, a group of corresponding landmarks localized on
the training surfaces are rigidly aligned. From the aligned surfaces, we then
choose the surface that is closest to the average landmarks as the template,
and warp it to each rest surface via non-rigid coherent point drift (CPD)
matching [8]. All warped surface templates are used for network training. We
reduce the number of vertices via surface simplification [9] to reduce the
computation cost. The vertex coordinates are min-max normalized. The loss
functions are minimized with the Adam [10] optimizer to determine the optimal
network parameters. The two networks are alternatively trained in turn on each
batch of samples. During the inference stage, the trained corrector is
directly applied on a patient bony surface to estimate a reference bony shape
model. The corrector network is invariant to the number and the order of
vertices on the testing surface [7].

##  3 Experimental Results

###  3.0.1 Materials and Settings:

We used CT scans of 67 normal subjects from a previous study [11] and CT scans
of 116 patients with jaw deformities to train our networks. The study was
approved by our Institutional Review Board (#Pro00009723). Following the
clinical routine [12], all CT data were segmented to reconstruct 3D bony shape
models. 51 landmarks were localized on each bony surface by experienced oral
surgeons for calculating the loss function during the training only. Following
Section 2.0.3, a bony surface template from the training surfaces was selected
based on the 51 landmarks to establish dense vertex correspondences.
Ultimately, a total of 7772 (67 $\times$ 116) random pairs of deformed-normal
bones were used to train the networks.

The simulator network was implemented with 4 encoding, 4 fusion, and 4
decoding layers. The corrector network was configured with 4 layers each for
encoding and decoding. The numbers of points for the encoding and decoding
layers were
$N_{\text{sub}}=\\{\frac{N}{4},\frac{N}{16},\frac{N}{64},\frac{N}{128},\frac{N}{64},\frac{N}{16},\frac{N}{4},N\\}$,
where $N=4724$ during training, and $N>10\,000$ during testing. The radius $r$
of the 3D neighboring ball was set to $\\{0.1,0.2,0.4,0.8\\}$ and
$\\{0.8,0.4,0.2,0.1\\}$ for the encoding and decoding layers, respectively.
The output point features from the 4 encoding and 4 decoding layers had
dimensions $C_{\text{feature}}=\\{64,128,256,512\\}$ and
$C_{\text{feature}}=\\{512,256,128,128\\}$, respectively. Empirically, we set
$\alpha=0.3$, $\beta=0.1$, and $\lambda=0.1$ in the loss functions (1) and
(6). The number $K=51$ of landmarks was set for calculating $L_{\text{jaw}}$
in (4). The networks were trained for 400 epochs with the learning rate of
$0.0001$.

For testing, we acquired the data from another 24 patients with paired pre-
and post-operative CT scans. The post-operative bones were used as ground
truth for performance evaluation. Following the method in [4], we remeshed
each post-operative bony surface according to its pre-operative bony surface
to construct the vertex-wise correspondences for calculating estimation
accuracy. Four evaluation metrics [4], i.e., vertex distance (VD), edge-length
distance (ED), surface coverage (SC), and landmark distance (LD), were
employed in the quantitative evaluation. We compared our method with DefNet
[4], which was trained using 6834 paired samples simulated by sparse
representation, 102 deformed bones were synthesized for each of 67 normal
bones from 116 patient bones. Qualitative evaluation was performed by an
experienced oral surgeon. Statistical significance of accuracy improvement was
evaluated using the paired t-test.

###  3.0.2 Results:

The trained framework was tested on the data of 24 patients. Figure 5 shows
the estimated bony surfaces and the vertex-wise distance heatmaps for three
randomly selected patients111A video illustration of the three estimation
examples is available as a supporting material.. All results from our method
inspected by the expert suggest that the estimated bones are clinically
acceptable. For quantitative evaluation, the estimation accuracy for the jaw
and the midface were calculated separately. Table 1 shows that our method is
statistically significantly more accurate in estimating the normal jaw
($p<0.05$) than DefNet based on the four metrics. The accuracy in maintaining
the midface is comparable between the two methods ($p>0.05$). Overall, the
proposed method outperforms DefNet in the term of accuracy.

Our method takes about 30 minutes for one training epoch, and 10 seconds for
testing on each patient data, evaluated based on an NVIDIA Titan Xp GPU.

Figure 5: Corrected bones estimated by the two methods for three randomly
selected patients. Heatmaps of surface vertex distance are calculated by
comparing the estimated bony surfaces with the corresponding post-operative
bony surfaces (ground truth). Table 1: Statistics for VD (mm), ED (mm), SC,
and LD (mm) based on evaluation using 24 patients.

| Method | VD | ED | SC | LD
---|---|---|---|---|---
Jaw | DefNet | 4.08 $\pm$ 1.02 | 0.31 $\pm$ 0.06 | 0.69 $\pm$ 0.05 | 4.26 $\pm$ 0.95
Ours | 3.49 $\pm$ 0.52 | 0.26 $\pm$ 0.05 | 0.72 $\pm$ 0.04 | 3.62 $\pm$ 0.53
Midface | DefNet | 1.36 $\pm$ 0.36 | 0.15 $\pm$ 0.04 | 0.96 $\pm$ 0.04 | 1.25 $\pm$ 0.31
Ours | 1.39 $\pm$ 0.37 | 0.15 $\pm$ 0.04 | 0.95 $\pm$ 0.04 | 1.32 $\pm$ 0.36

  * •

Note: A larger SC indicates better estimation performance.

##  4 Discussion and Conclusion

As sufficient paired deformed-normal bones are almost impossible to acquire
clinically, simulating paired data is necessary for training supervised deep
learning models for our task. Data simulation based on sparse representation
as done in DefNet [4] is limited in mimicking deformities, and requires
significant efforts from experts to confirm simulation quality. Our self-
supervised deep learning framework can be trained using unpaired data. The
embedded simulator network has a stronger non-linear representation ability
than sparse representation, and is able to accurately synthesize a broader
range of jaw deformities based on patient bones. Compared with DefNet, our
framework simplifies the training process and improves prediction accuracy.

## Acknowledgement

This work was supported in part by United States National Institutes of Health
(NIH) grants R01 DE022676, R01 DE027251, and R01 DE021863.

## References

  * [1] Xia, J., et al.: Algorithm for planning a double-jaw orthognathic surgery using a computer-aided surgical simulation (CASS) protocol. Part 1: planning sequence. Int. J. Oral Maxillofacial Surg. 44(12), 1431–1440 (2015)
  * [2] Xia, J., et al.: Algorithm for planning a double-jaw orthognathic surgery using a computer-aided surgical simulation (CASS) protocol. Part 2: three-dimensional cephalometry. Int. J. Oral Maxillofacial Surg. 44(12), 1441–1450 (2015)
  * [3] Wang, L., et al.: Estimating patient‐specific and anatomically correct reference model for craniomaxillofacial deformity via sparse representation. Med. Phys. 42(10), 5809-5816 (2015)
  * [4] Xiao, D., et al.: Estimating reference bony shape models for orthognathic surgical planning using 3D point-cloud deep learning. IEEE J. Biomed. Health Informat. (2021). 10.1109/JBHI.2021.3054494
  * [5] Qi, C.R., et al.: PointNet++: Deep hierarchical feature learning on point sets in a metric space. In: Proc. Adv. Neural Inf. Process. Syst., pp. 5099–5108. (2017)
  * [6] Wu, W., et al.: PointConv: Deep convolutional networks on 3D point clouds. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit., pp. 9621–9630. (2019)
  * [7] Qi, C.R., et al.: PointNet: Deep learning on point sets for 3D classification and segmentation. In: Proc. IEEE Conf. Comput. Vis. Pattern Recognit., pp. 652–660. (2017)
  * [8] Myronenko, A., Song, X.: Point set registration: Coherent point drift. IEEE Trans. Pattern Anal. Mach. Intell. 32(12), 2262–2275 (2010)
  * [9] Garland, M., Heckbert, P.: Surface simplification using quadric error metrics. In: Proc. SIGGRAPH 97, pp. 209–216. (1997)
  * [10] Kingma, D., Ba, J.: Adam: A method for stochastic optimization. In: Proc. Int. Conf. Learn. Representations, pp. 1–41. (2015)
  * [11] Yan, J., et al.: Three-dimensional CT measurement for the craniomaxillofacial structure of normal occlusion adults in Jiangsu, Zhejiang and Shanghai Area. China J. Oral Maxillofac. Surg. 8, 2–9 (2010)
  * [12] Yuan, P., et al.: Design, development and clinical validation of computer-aided surgical simulation system for streamlined orthognathic surgical planning. Int. J. Comput. Assist Radiol. Surg. 12(12), 2129–2143 (2017)


# BRAIN-INSPIRED FEATURE EXAGGERATION IN GENERATIVE REPLAY FOR CONTINUAL
LEARNING

###### Abstract

The catastrophic forgetting of previously learnt classes is one of the main
obstacles to the successful development of a reliable and accurate generative
continual learning model. When learning new classes, the internal
representation of previously learnt ones can often be overwritten, resulting
in the model’s ”memory” of earlier classes being lost over time. Recent
developments in neuroscience have uncovered a method through which the brain
avoids its own form of memory interference. Applying a targeted exaggeration
of the differences between features of similar, yet competing memories, the
brain can more easily distinguish and recall them. In this paper, the
application of such exaggeration, via the repulsion of replayed samples
belonging to competing classes, is explored. Through the development of a
’reconstruction repulsion’ loss, this paper presents a new state-of-the-art
performance on the classification of early classes in the class-incremental
learning dataset CIFAR100.

Index Terms—  Continual learning, catastrophic forgetting, generative replay,
VAE, computational neuroscience

##  1 Introduction

Artificial neural networks (ANNs) can provide state-of-the-art performance on
traditional image classification tasks. Yet, behind this success, there
remains the assumption that training inputs are independently and identically
distributed (IID). Typically, data is randomly shuffled prior to training,
reducing the dependence of consecutive examples. However, this is not always
possible; in the real world, samples tend to be temporally correlated, for
instance a robot learning new objects classes may first see many images of a
dog, then many of a cat. With ANNs, the result of such correlation is a strong
interference between learned parameters important to the prediction of current
and previously seen classes – causing the model to forget previously learned
classes almost entirely. An ideal remedy for such catastrophic forgetting (CF)
would involve the storage and interleaved replay of past examples, thus
removing any temporal correlation. In practice, however, the storage of an
every-growing set of samples is infeasible. Taking inspiration from the brain,
a more tractable solution involves the incorporation of a generative
model—such as VAE or GAN—into the continual learner, allowing for previously
seen classes to be ”replayed” alongside the current classes, in much the same
way memories are replayed in the mind. The current state-of-the-art generative
model, on which this paper builds [1], is able to obtain close to storage
based performance on a 10-class MNIST dataset, however, when applied to the
more challenging 100-class CIFAR100 dataset, classification precision compared
to the stored approach is greatly reduced—with earlier seen classes being
affected the most. This paper aims to address the reduction in early class
performance by imitating the newfound role of feature exaggeration for the
reduction of memory interference in the brain.

##  2 FEATURE EXAGGERATION IN THE BRAIN

Given the previous success of brain-inspired approaches aimed at mitigating CF
[1, 2, 3, 4] looking again to the brain—which has developed effective
processes to effortlessly tackle such recall—for inspiration seems inevitable.
Neuroscience studies have shown that participants trained on competing image-
name pairs (such as two similar ape images named Angus and Abner respectively)
showed greater repulsion in their neural representations [5]. What’s more,
greater repulsion is correlated with better recall of the image-name pairs,
indicating that the brain actively differentiates between highly similar,
competing memories to the improve their future recall. A more recent study by
Zhao et al [6] focused on identifying whether the distortion of specific
feature dimensions of memories (such as colour) plays a role in memory
repulsion. The study found that, given two images of the same object in
slightly different colours, participants tended to exaggerate the colour
difference between them, such that the pair were more distinguishable. As
above, a greater colour exaggeration was associated with a greater recall of
linked information. This suggests that the brain adaptively targets and
accentuates small differences in important discriminating features of similar
memories to reduce interference. Ultimately, this calls to attention to a new
approach to deal with interference in ANNs: instead of fighting against
interference, the model could be shaped by interference through an adaptive
feature exaggeration based approach.

##  3 METHODS

As described in Section 2, the phenomenon of feature-exaggeration in the brain
tends to occur in competing memories: where the combination of memories of
similar, yet distinct, objects can result in confusion during the recall of
either of item. This beneficial memory distortion focuses on exaggerating the
differences between the objects, resulting in a repulsion in the brain’s
internal representations of them which is linked to an improved rate of object
recall. In order to mimic this process in the case of continual learning, a
measure of competitive classes must first be defined, followed by a method
through which the feature-wise differences between these classes can be
exaggerated.

###  3.1 VAE architecture

Fig. 1: An overview of the BI-R model’s architecture [1]. Inputs,
$\mathbf{x}$, are fed into the encoder and FC layers, converting $\mathbf{x}$
to samples, $\mathbf{z}$, from the conditional distribution
$q_{\bm{\phi}}(\mathbf{z}|\mathbf{x})$. Depending on whether classification or
sample generation was occurring, the $\mathbf{z}$ samples would then be passed
through either the softmax (classification) or the split layer (generation).
At the split layer, the samples are re-parameterised into $\bm{\mu}$ and
$\log\bm{\sigma^{2}}$ before passing through two more FC layers and the
decoder. $\phi$, $\theta$ are the encoder, decoder weights respectively.

The model used in this work is built upon the brain-inspired replay (BI-R)
model proposed by van de Ven et al [1]. This combines a softmax classifier
with a VAE to generate samples based on previously seen classes and replay
them alongside the current inputs, as depicted in Fig. 1. The aim of a VAE’s
encoder is to distil the high-dimensional input space,
$\mathbf{x}\in\mathbb{R}^{n}$, into a lower-dimensional latent space,
$\mathbf{z}\in\mathbb{R}^{m}$ (where $m<n$), by approximating the intractable
posterior, $p(\mathbf{z}|\mathbf{x})$, as a multivariate Gaussian
distribution, $q_{\bm{\phi}}(\mathbf{z}|\mathbf{x})$, via variational
inference. This Gaussian is defined via its mean vector, $\bm{\mu}$, and the
logarithm of its variance, $\log\bm{\sigma^{2}}$, which are returned via a
split layer after the encoder and two FC layers, as depicted in Fig. 1. The
latent variables sampled from this distribution are then passed through two
more FC layers and the decoder, which acts as the conditional distribution,
$p_{\bm{\theta}}(\mathbf{x^{\prime}}|\mathbf{z})$, returning reconstructed
samples, $\mathbf{x^{\prime}}$. The VAE is then optimised over successive
batches by minimising the per-batch loss function:

$\displaystyle\begin{aligned}
&L(\mathbf{x},\mathbf{x}^{\prime},\mathbf{z};\bm{\phi},\bm{\theta},\bm{\lambda})=\\\
&\sum^{N_{s}}_{i=1}\left[\lambda_{r}\ell_{rec}(\bm{x}_{i},\bm{x}^{\prime}_{i};\bm{\phi},\bm{\theta})\right.\left.+\lambda_{d}\ell_{reg}(\bm{x}_{i};\bm{\phi},\bm{\lambda})\right],\end{aligned}$
(1)

where $N_{s}$ is the number of samples per batch; $\bm{\phi}$, $\bm{\theta}$
are the encoder/decoder weights and $\bm{\lambda}=(\bm{\mu},\bm{\sigma^{2}})$.
$\lambda_{r}$, $\lambda_{d}$ are the weights of the reconstruction and
regularisation losses (denoted as $\ell_{rec}$ and $\ell_{reg}$ respectively).
The reconstruction loss $\ell_{rec}(\cdot)$ of Eq. 1 is a feature-wise binary
cross-entropy between generated samples and their original counterparts (Eq.
2). This loss effectively ensures that the reconstructed samples appear as
similar to the original samples as possible.

$\displaystyle\begin{aligned}
&\ell_{rec}(\bm{x}_{i},\bm{x}^{\prime}_{i};\bm{\phi},\bm{\theta})=\\\
&\mathbb{E}_{\mathbf{\epsilon}\sim\mathbb{N}(\mathbf{0},\mathbf{I})}\left[\sum^{M_{f}}_{f=1}x_{i,f}\log(x^{\prime}_{i,f})+(1-x_{i,f})\log(1-x^{\prime}_{i,f})\right].\end{aligned}$
(2)

Where $M_{f}$ is the number of features. The regularisation loss
$\ell_{reg}(\cdot)$ aims to maintain as simplistic an approximating
distribution to the posterior as possible by minimising the difference between
the variational approximation, $q_{\bm{\phi}}(\mathbf{z}|\mathbf{x})$, and a
standard normal prior distribution, $p(\mathbf{x})$. To implement this, the
loss minimises the Kullback–Leibler (KL) divergence between the approximated
posterior and prior distributions.

$\displaystyle\begin{aligned}
\ell_{reg}(\bm{x}_{i};\bm{\phi},\bm{\lambda})=D_{KL}[q_{\bm{\phi}}(\mathbf{z}_{i}|\mathbf{x}_{i})||p(\mathbf{z}_{i})],\end{aligned}$
(3)

where $D_{KL}[\cdot]$ denotes the KL function that measures the distance
between two underlying distributions.

###  3.2 Competing classes

To determine the competing classes in practice, the batch of samples is passed
through the model’s own classifier: the softmax layer. This returns a discrete
probability distribution over all classes the model has so far seen, and can
be used to identify competing classes for a given sample according to the
model itself. Thus, competing classes are identified as those the model deems
most similar and vacillates between during classification. By selecting
similar classes this way, the newly proposed feature exaggeration inspired
loss can target the more frequently miss-classified classes.

To account for varying degrees of competitiveness between classes, an adaptive
process, involving the selection of a repulsion factor, $f$, where $f\geq 1$,
was used. This factor acts as a dynamic threshold that changes the selection
criteria of competing classes depending on the model’s confidence in its
original label. The highest softmax probability (corresponding to the sample
label) is divided by the repulsion factor, providing a threshold probability
for identifying competing classes for that sample. Any class with a
probability above that of the threshold is deemed competitive. To illustrate
this, Fig. 2 depicts the output of the softmax layer for a given sample and
process through which competitive classes are identified. If samples exist
that don’t have any classes other than the label above the competitive
threshold, no classes are selected for that sample and the repulsion loss is
not calculated over it.

Fig. 2: Competing class selection using a repulsion factor. The most likely
class (3) is taken as the sample label and its probability is divided by the
repulsion factor, $f$, to determine the repulsion threshold. Excluding the
label, any class with a probability greater than the threshold is selected as
a competing class to the label (2), any below this are ignored. In this
example $f$ is 3, making the repulsion threshold 0.2.

###  3.3 Reconstruction repulsion loss

To encourage the exaggeration of feature-wise differences between competing
classes, a new ”reconstruction repulsion” (RR) loss was developed. Taking
inspiration from the VAE’s reconstruction loss, this acts as a third
generative loss, alongside the reconstruction and regularisation terms. It
aims to maximise the binary cross-entropy (BCE) between a generated sample,
$\mathbf{x}_{i}^{\prime}$, and an original input belonging to a competing
class $\mathbf{x}_{i_{C}}$ to that of the sample, as shown in (Eq. 4). Through
the incorporation of this loss, it is hoped that VAE can be taught to not only
learn the characteristic feature vector of a generated class, but to also
exaggerate the differences between this class and other competing classes. An
important aspect to this, however, is that such repulsion between the features
of competing samples can only be done in moderation. Too much and the
generated samples may be unrecognisable from their intended classes.
Therefore, when combined with the original losses, the RR loss was weighted
with a factor, $\lambda_{RR}$, and the ratio between this factor and that of
the original reconstruction loss was explored. Furthermore, once the competing
classes to a given sample have been determined, a mean feature vector,
$\bm{\hat{x}}$, for each competing class is calculated via the mean of all
original samples in that class. Having the same shape and dimensionality as
the inputs, this mean feature vector can be viewed as a generalised sample
that can capture more of the characteristic feature patters of a given class.

$\displaystyle\begin{aligned}
&\ell_{RR}(\bm{x}^{\prime}_{i},\bm{\hat{x}}_{i_{C}};\bm{\phi},\bm{\theta})=\\\
&\mathbb{E}^{-1}_{\mathbf{\epsilon}\sim\mathbb{N}(\mathbf{0},\mathbf{I})}\left[\sum^{M_{f}}_{f=1}\hat{x}_{i_{C},f}\log(x^{\prime}_{i,f})+(1-\hat{x}_{i_{C},f})\log(1-x^{\prime}_{i,f})\right]\end{aligned}$
(4)

In addition, a similar reconstruction attraction (RA) loss was also used to
exaggerate the features of replayed samples belonging to the same class
towards each other, in an effort to attract their internal representations
(Eq. 5). The reasoning behind this relates again to Zhao et al’s [6]
observations of memory repulsion in the brain—namely that memories of related
objects tend to be exaggerated towards each other.

$\displaystyle\begin{aligned}
&\ell_{RA}(\bm{x}^{\prime}_{i},\bm{\hat{x}};\bm{\phi},\bm{\theta})=\\\
&\mathbb{E}_{\mathbf{\epsilon}\sim\mathbb{N}(\mathbf{0},\mathbf{I})}\left[\sum^{M_{f}}_{f=1}\hat{x}_{f}\log(x^{\prime}_{i,f})+(1-\hat{x}_{f})\log(1-x^{\prime}_{i,f})\right]\end{aligned}$
(5)

Due to the fact that the exaggeration observed in the brain only occurs
between memories themselves and not to the objects responsible for these
memories, the RR and RA losses were only applied to batches of replayed
samples, rather than new class batches.

##  4 RESULTS

$f$ | Improvement over baseline (%)
---|---
20% avg. |  50% avg. | Overall avg.
5 | 3.80/12.5 |  4.06/12.78 | 1.35/3.87
10 | 5.80/19.08 |  3.82/12.03 | 1.19/3.41
20 | 4.80/15.79 |  4.86/15.30 | 1.48/4.25
50 | 4.90/16.12 |  3.56/11.21 | 1.26/3.61
100 | 4.95/16.28 |  4.18/13.16 | 1.24/3.56
Table 1:  Table showing the absolute and relative change in precision above
that of the baseline for various repulsion factors, in the format: absolute
change/relative change. Comparisons are shown across the average (denoted as
avg.) of all test classes (Overall), the first 50% of classes seen and the
first 20% seen. Experiments were run using 5000 iterations with the pre-
determined optimal hyperparameters. Fig. 3: A comparison of the final
classification test set precision and convergence rates of the RR+RA
adaptation (RR) and baseline (BL). Presented are the overall precision
(average across all classes; a), precision across the first 50% of classes
encountered (b) and over the first 20% of classes encountered (c). This was
performed over current and replay batch sizes of 512, with $\lambda_{RR}$ and
$\lambda_{RA}$ weightings of $10^{-6}$ (the optimal values).

The final performance of the RR+RA model with optimised hyperparameters was
compared to that of the state-of-the-art BI-R ”baseline” model, using the
optimised hyperparameters as stated by van de Ven et al [1]. Given that CF
directly reduces the final classification precision of previously seen
classes, this was main metric employed as an indicator of continual learning
performance. All results stated relate to the average test precision across 5
runs of the experiment with optimal hyperparameters. Due to the nature of CF,
it will tend to have a greater detrimental effect on the learned
representations of earlier-seen classes, whose internal representations are
more likely to be ’over-written’. Thus, the average precision on the first 50%
of classes seen will tend to be lower than the average over the last 50% and
any successful method will likely lead to an improvement, not necessarily in
the overall average precision, but in the average precision of the first 50%
and most likely the first 20% of classes.

The final precision and convergence rate of the RR+RA adaptation was compared
to the baseline over a range of gradient descent iterations, as shown in Fig.
3. The average precision across all classes (Fig. 3a) showed a slight,
although not significant, improvement across the tested iterations. Comparing
the average precision of the first 50% of classes (Fig. 3b), however, the
inclusion of the RR+RA loss provided a significant increase over the baseline
for all tested iterations. For instance, the RR+RA model achieved an absolute
increase in percentage precision of 3.20% and a relative increase of 15.30%
when compared to baseline for 5000 iterations. The average performance over
the first 20% of classes seen was also compared to the baseline, showing
similar significant improvements to classification precision across all
iterations. A slightly larger absolute improvement of 4.30% and relative
increase of 15.79% over the baseline was obtained for 5000 iterations. The
best values for $f$ proved to be 10 and 20 (Table 1). Performance using the
repulsion factor, $f$, as described in 3.2 was also compared to that without
it without it, with the inclusion of $f$ providing a mean relative precision
improvement of 3.45% for the first 50 percent of classes.

##  5 DISCUSSIONS AND CONCLUSIONS

For all average metrics used, the RR+RA loss provided a significant
improvement over the baseline classification precision. This improvement was
amplified across the classification of the earlier-seen classes, with a mean
relative precision improvement of 15.30% for the first 50 percent of classes
and 15.79% for the first 20 percent of classes (absolute changes of 4.86% /
4.8%). These strong early class results act to confirm the hypothesis that a
targeted feature-based exaggeration of samples allows the model to more easily
learn the characteristic differences and similarities between samples of
competing and non-competing classes, resulting in a reduction of CF amongst
earlier-learnt classes. To further confirm the internal repulsion mechanism
through which this adaptation improves precision, future work could focus on
analysing the effect of the loss on the internal representations of competing
classes.

The RR+RA adaptation is certainly not perfect—relying on the assumption of
correct class labels being determined during the generation of replayed
batches. This assumption does not always hold true and could be a source of
interference in the optimisation of the RR loss. Further research could aim to
incorporate the model’s confidence in its predictions in the RR loss. This
alteration has the potential to lead to a further boost in classification
performance.

Taking a step back, this paper is centred around the task of continual image
classification, but modifications such as the RR+RA loss could also be applied
to other continual learning contexts with different data-types. For instance,
given that the loss only depends on the BCE between the feature-vectors of
competing samples, as long as the data can be translated into feature-vector
form this could be calculated regardless of the vector’s shape or dimensions.

## References

  * [1] Gido M van de Ven, Hava T Siegelmann, and Andreas S Tolias,  “Brain-inspired replay for continual learning with artificial neural networks,”  Nature communications, vol. 11, no. 1, pp. 1–14, 2020.
  * [2] James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, et al.,  “Overcoming catastrophic forgetting in neural networks,”  Proceedings of the national academy of sciences, vol. 114, no. 13, pp. 3521–3526, 2017.
  * [3] Friedemann Zenke, Ben Poole, and Surya Ganguli,  “Continual learning through synaptic intelligence,”  in International Conference on Machine Learning. PMLR, 2017, pp. 3987–3995.
  * [4] Hanul Shin, Jung Kwon Lee, Jaehong Kim, and Jiwon Kim,  “Continual learning with deep generative replay,”  arXiv preprint arXiv:1705.08690, 2017.
  * [5] JC Hulbert and KA Norman,  “Neural differentiation tracks improved recall of competing memories following interleaved study and retrieval practice,”  Cerebral Cortex, vol. 25, no. 10, pp. 3994–4008, 2015.
  * [6] Yufei Zhao, Avi JH Chanales, and Brice A Kuhl,  “Adaptive memory distortions are predicted by feature representations in parietal cortex,”  Journal of Neuroscience, vol. 41, no. 13, pp. 3014–3024, 2021.