# Poster: ARLIF-IDS — Attention augmented Real-Time Isolation Forest Intrusion
Detection System

Aman Priyanshu  Dept. of Information &
Communication Technology
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>Sarthak Shastri  Dept. of Information &
Communication Technology
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>Sai Sravan Medicherla  Dept. of
Information &
Communication Technology
Manipal Institute of Technology
Manipal, India
<EMAIL_ADDRESS>
###### Abstract

Distributed Denial of Service (DDoS) attack is a malicious attempt to disrupt
the normal traffic of a targeted server, service or network by overwhelming
the target or its surrounding infrastructure with a flood of Internet traffic.
Emerging technologies such as the Internet of Things and Software Defined
Networking leverage lightweight strategies for the early detection of DDoS
attacks. Previous literature demonstrates the utility of lower number of
significant features for intrusion detection [4, 2, 1]. Thus, it is essential
to have a fast and effective security identification model based on low number
of features.

In this work, a novel Attention-based Isolation Forest Intrusion Detection
System is proposed. The model considerably reduces training time and memory
consumption of the generated model. For performance assessment, the model is
assessed over two benchmark datasets, the NSL-KDD dataset & the KDDCUP’99
dataset. Experimental results demonstrate that the proposed attention
augmented model achieves a significant reduction in execution time, by 91.78%,
and an average detection F1-Score of 0.93 on the NSL-KDD and KDDCUP’99
dataset. The results of performance evaluation show that the proposed
methodology has low complexity and requires less processing time and
computational resources, outperforming other current IDS based on machine
learning algorithms.

###### Index Terms:

Attention, Isolation Forest, Intrusion Detection System, Memory Optimization,
Edge Deployment

##  I Introduction

The rapid advances in the communication technology have resulted in
transmission services handling huge network size and the corresponding data.
As a result, malicious entities have ample opportunity to attack and deny
service to target systems (DDoS). Despite existing traditional solutions,
Distributed Denial of Service (DDoS) attacks continue to be a prominent
cybercrime. Many attacks occur each day and pose security challenges for
client-networks. Therefore, the existence of accurate intrusion detection
systems becomes indispensable. An intrusion detection system (IDS) is a tool
that prevents the network from possible intrusions by inspecting the network
traffic, to ensure its confidentiality, integrity, and availability. Despite
enormous efforts by the researchers, IDS still faces challenges in improving
detection performance while reducing detection-time and increasing memory
optimization.

##  II Related Work

Existing surveys on anomaly detection techniques, only marginally consider the
utility of neural-network-based approaches, due to their intensive memory
consumption and high detection-time. Current research developments regarding
deep neural networks and LSTM architectures for anomaly detection are
oftentimes not scalable due to these issues. However, IDS which employ machine
learning methodology, such as classification-based, clustering-based and
information-theoretic approaches. Thus, techniques based on the isolation
forest (IF), random forest (RF), distance-based inference, entropy-based
inference constitute a major part of the investigations [4]. Theoretical
studies have also been conducted on LSTM and CNNs. An analysis is conducted
over these deep neural network methods.

##  III Methodology

Our proposed methodology, which we call ARLIF - Attention augmented Real-Time
Isolation Forest Intrusion Detection System, draws inspiration from the
transformer models [5]. However, unlike traditional attention-based model, we
employ learning on only the attention layer, whereas time-series datum is
computed using proba values predicted by Isolation Trees. The methodology
allows online learning, and real-time detection systems may employ it for
high-performance detection.

The architecture of this model relies solely on attention mechanisms to carry
out the seasonality learning, allowing the modulation and focus of predictions
to be placed on specific isolation trees, based on previous responses
(context).

###  III-A Attention Mechanism

Figure 1: Proposed methodology - Attention augmented Real-Time Isolation
Forest Intrusion Detection System

In order to utilize the attention mechanism in the network intrusion detection
domain and take advantage of its properties, we employ its contextualizing
attentive properties during forest response averaging. Taking each trees
predicted probability as its vector at $V_{k}^{i}$ for $tree^{i}$ at $time=k$.
Thereby, allowing the attention mechanism to create updated vector embeddings
for every tree, based on responses from others as well as, their time-series.

To formulate the algorithm presented in Fig 1, we limit previous time-series
predictions of each tree, to be previous $k$-responses. Therefore once a
forest of isolation trees are trained, the online learning will only be
executed through the attention mechanism, which would have:
$number\\_of\\_trainable\\_parameters=3*k*(k+1)$, which would be based on
$query,key,\&value$ layers. The model therefore, offers a memory optimal, time
efficient approach towards intrusion detection.

##  IV Results

TABLE I: Model Performance Inference on Datasets Model Name | F1-Score | Memory Acquired | Detection-Time
---|---|---|---
CNN | 0.95 | 105kB | 4.6s
LSTM | 0.96 | 238kB | 7.3s
ARLIF-IDS | 0.93 | 11kB | 0.6s
Isolation Forest | 0.87 | 6kB | 0.4s
Random Forest | 0.88 | 44kB | 0.9s
Entropy-based IDS | 0.83 | 3kB | 0.2s
Distance-based IDS | 0.84 | 4kB | 0.2s

We carry out a set of experimental evaluations in order to provide a deep
analysis of the performance of the proposed model. Then, we compare the
overall results against baseline models. The models are tested on KDDCUP’99
dataset and the NSL-KDD dataset as benchmark.

We present these results in Table I. As demonstrated the proposed methodology
is able to give high preliminary results for the memory-time-performance
tradeoff. The table presents neural network and tree based methods as the
highest performing algorithms, where as, entropy-based and distribution-based
methods provide low memory and execution-time utility. ARLIF-IDS gives a
satisfactory trade-off giving performance even higher than LSTM networks.

##  V Conclusion

This poster presents an intrusion detection model with attention mechanism for
isolation forest augmentation. We take into account the sparse correlation
majority traffic-features have with anomaly detection, and develop a strategy
to utilize models with low-feature space input further cutting down on
computational costs [4]. The proposed model was demonstrated to achieve
satisfactory performance of 0.93 F1-Score for the KDDCUP’99 dataset and NSL-
KDD dataset. We also show that it lowers average detection-time of 91.78% with
respect to LSTM anomaly detection model and 95.37% lower memory utilization
with respect to the same. For future work, we hope to utilize other
probabilistic methods instead of isolation trees.

## References

  * [1] Q. R. S. Fitni and K. Ramli, ”Implementation of Ensemble Learning and Feature Selection for Performance Improvements in Anomaly-Based Intrusion Detection Systems,” 2020 IEEE International Conference on Industry 4.0, Artificial Intelligence, and Communications Technology (IAICT), 2020, pp. 118-124, doi: 10.1109/IAICT50021.2020.9172014.
  * [2] A. A. Aburomman and M. Bin Ibne Reaz, ”Ensemble of binary SVM classifiers based on PCA and LDA feature extraction for intrusion detection,” 2016 IEEE Advanced Information Management, Communicates, Electronic and Automation Control Conference (IMCEC), 2016, pp. 636-640, doi: 10.1109/IMCEC.2016.7867287.
  * [3] C. Liu, Y. Liu, Y. Yan and J. Wang, ”An Intrusion Detection Model With Hierarchical Attention Mechanism,” in IEEE Access, vol. 8, pp. 67542-67554, 2020, doi: 10.1109/ACCESS.2020.2983568.
  * [4] Josy Elsa Varghese, Balachandra Muniyal, & Aman Priyanshu (2022). Finding an elite feature for (D)DoS fast detection—Mixed methods research. Computers & Electrical Engineering, 98, 107705.
  * [5] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, & Illia Polosukhin (2017). Attention Is All You Need. CoRR, abs/1706.03762.


# Tipping Cycles

Michael A.S. Thorne
British Antarctic Survey

Ecological systems are studied using many different approaches and
mathematical tools [3]. One approach, based on the Jacobian of Lotka-Volterra
type models [6], has been a staple of mathematical ecology for years, leading
to many ideas such as on questions of system stability [8]. Instability in
such systems is determined by the presence of an eigenvalue of the community
matrix lying in the right half plane. The coefficients of the characteristic
polynomial derived from community matrices contain information related to the
specific matrix elements that play a greater destabilising role. Yet the
destabilising circuits, or cycles, constructed by multiplying these elements
together, form only a subset of all the feedback loops comprising a given
system. This paper looks at the destabilising feedback loops in predator-prey,
mutualistic and competitive systems in terms of sets of the matrix elements to
explore how sign structure affects how the elements contribute to instability.
An analytic method is then described in which a variable acting on the self-
regulating terms of a community matrix enables one, depending on the field of
the maximal eigenvalue, to either determine the size of the key feedback that
leads to the system’s tipping point, or narrows down the search space.

Within mathematical ecology, one approach [3] to represent predator-prey,
mutualistic and competitive systems is through the community matrix, a real-
valued square matrix representing the linearisation of Lotka-Volterra type
dynamic equations [6]. Each of the forms is represented through differing sign
conventions as can be seen in the following $3\times 3$ examples:

Predator-Prey | Mutualistic | Competitive
---|---|---
$\begin{bmatrix}-a&b&c\\\ -d&-e&f\\\ -g&-h&-k\end{bmatrix}$ | $\begin{bmatrix}-a&b&c\\\ d&-e&f\\\ g&h&-k\end{bmatrix}$ |  $\begin{bmatrix}-a&-b&-c\\\ -d&-e&-f\\\ -g&-h&-k\end{bmatrix}$.

Stability of a community matrix is determined by whether all of the
eigenvalues lie in the left half of the complex plane. The characteristic
polynomial of an $n\times n$ community matrix is an $n^{th}$-order monic
polynomial (or made so by change of sign), whose roots are the eigenvalues of
the system. A necessary condition for the roots of a polynomial to all lie in
the left half plane, and therefore for the system to be stable, is for all of
its coefficients to be positive. However, this condition is not sufficient.
The Routh-Hurwitz stability criteria ([11],[5]) provides a further, and
sufficient, condition. The following equation is an example of a polynomial
satisfying the necessary condition of positive coefficients,

(1) $x^{4}+2x^{3}+3x^{2}+4x+5=0,$

yet whose Routh table (Table 1) reflects the existence of two roots in the
right half plane and is therefore not Routh-Hurwitz stable.

$x^{4}$ | 1 | 3 | 5
---|---|---|---
$x^{3}$ | 2 | 4 |
$x^{2}$ | 1 | 5 |
$x^{1}$ | -6 |  |
$x^{0}$ | 5 |  |
Table 1: Table 1 The Routh table for Eq. 1. The number of sign changes as one
moves down the first column indicate the number of roots in the right half
plane. In this case, two.

Polynomials of the example just described can only have roots in the right
half plane that are complex conjugate pairs. This was proven by Obrechkoff
[10], who showed that there are no solutions to polynomials with positive
coefficients that lie on the positive real axis (the right half plane roots of
Eq. 1 are $0.29\pm 1.42i$).

Therefore, any systems that have a real maximal (largest real part) eigenvalue
(with imaginary part zero) have their tipping point between stability and
instability just where the last coefficient of their characteristic polynomial
becomes positive. This point is helpful in considering what contributes to the
destabilisation of a system, either in untangling the elements that make up
the feedbacks, as in the next section, or in isolating the size of the key
feedback or feedbacks that results in a system being pushed to its tipping
point, in the section following.

## Tipping Cycles

The idea of a tipping point is not very informative without its context,
without knowledge of the feedbacks that produce the instability. Ecological
systems that are modelled by community matrices have the benefit of the direct
relation through the coefficients of the characteristic polynomial which are
derived from the minors, or respectively the summations of the n-tuple sets of
the eigenvalues [1], of the matrix. These minors may be exhaustively described
by the terms consisting of sets of multiplied elements of the community
matrix. These terms are the feedback loops of the system, otherwise called
circuits, or cycles. Within each coefficient of a characteristic polynomial,
one can distinguish the cycles that contribute to it becoming more positive,
and therefore stabilising, or more negative and destabilising. For example,
consider the general characteristic polynomial of the $3\times 3$ predator-
prey community matrix described above,

$a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0},$

which expands, in terms of its matrix elements as,

$x^{3}+(a+e+k)x^{2}+(bd+ae+cg+fh+ak+ek)x+(ceg+bfg+afh+bdk+aek{\color[rgb]{1,0,0}-cdh}).$

Given that characteristic polynomials are monic, we know that the highest
order coefficient, $a_{3}$, is $1$ and therefore positive, $a_{2}$ is the
negative of the trace, and therefore positive, and as can be seen, $a_{1}$,
the summed combination of 2-tuples, is also positive. This leaves only $a_{0}$
able to switch between being positive or negative, and therefore able to
destabilise the system. For clarity of explanation, we have used the change of
sign of the coefficient for consideration of stability, even though this is
only a clear indication in systems with a maximal real-valued eigenvalue. But
the effect of increasing the positive or negative value of a coefficient on
its stability applies to all systems. Within $a_{0}$, highlighted above and
below in red, the only destabilising circuit, or tipping cycle, is $cdh$,

$\left[\begin{array}[]{rrr}-a&b&{\color[rgb]{1,0,0}c}\\\
{\color[rgb]{1,0,0}-d}&-e&f\\\
-g&{\color[rgb]{1,0,0}-h}&-k\end{array}\right].$

This means that however large the other elements are, they cannot force the
system to become unstable (if it is not already so), only the circuit of c, d
and h can do that. Given that there are five other cycles of size three that
make up the $a_{0}$ coefficient, it is helpful to consider a ratio, the
coefficient feedback sensitivity, $\tilde{a_{i}}$, which in this case is

$\tilde{a_{0}}=1/6.$

Table 2 shows the values of $\tilde{a_{i}}$ for predator-prey, mutualistic and
competitive systems up to $8\times 8$ sized community matrices.

| $n$ |  $a_{8}$$x^{8}$ |  $a_{7}$$x^{7}$ |  $a_{6}$$x^{6}$ |  $a_{5}$$x^{5}$ |  $a_{4}$$x^{4}$ |  $a_{3}$$x^{3}$ |  $a_{2}$$x^{2}$ |  $a_{1}$$x$ | $a_{0}$
---|---|---|---|---|---|---|---|---|---|---
P-P | $2$ |  |  |  |  |  |  | + | + | +
| $3$ |  |  |  |  |  | + | + | + | 1/6
| $4$ |  |  |  |  | + | + | + | 4/24 | 8/24
| $5$ |  |  |  | + | + | + | 10/60 | 40/120 | 52/120
| $6$ |  |  | + | + | + | 20/120 | 120/360 | 312/720 | 344/720
| $7$ |  | + | + | + | 35/210 | 280/840 | 1092/2520 | 2408/5040 | 2488/5040
| $8$ | + | + | + | 56/336 | 560/1680 | 2912/6720 | 9632/20160 | 19904/40320 | 20096/40320
M | $2$ |  |  |  |  |  |  | + | + | 1/2
| $3$ |  |  |  |  |  | + | + | 3/6 | 5/6
| $4$ |  |  |  |  | + | + | 6/12 | 20/24 | 20/24
| $5$ |  |  |  | + | + | 10/20 | 50/60 | 100/120 | 84/120
| $6$ |  |  | + | + | 15/30 | 100/120 | 300/360 | 504/720 | 424/720
| $7$ |  | + | + | 21/42 | 175/210 | 700/840 | 1764/2520 | 2968/5040 | 2680/5040
| $8$ | + | + | 28/56 | 280/336 | 1400/1680 | 4704/6720 | 11872/20160 | 21440/40320 | 20544/40320
C | $2$ |  |  |  |  |  |  | + | + | 1/2
| $3$ |  |  |  |  |  | + | + | 3/6 | 3/6
| $4$ |  |  |  |  | + | + | 6/12 | 12/24 | 12/24
| $5$ |  |  |  | + | + | 10/20 | 30/60 | 60/120 | 60/120
| $6$ |  |  | + | + | 15/30 | 60/120 | 180/360 | 360/720 | 360/720
| $7$ |  | + | + | 21/42 | 105/210 | 420/840 | 1260/2520 | 2520/5040 | 2520/5040
| $8$ | + | + | 28/56 | 168/336 | 840/1680 | 3360/6720 | 10080/20160 | 20160/40320 | 20160/40320
Table 2: Table 2 For the three forms under discussion $\tilde{a_{i}}$ is
presented for systems up $n=8$. By the inherent sign symmetries and
proportions it would be reasonable to suggest that, for any size $n$,
$\tilde{a_{i}}<1/2$ in predator-prey (P-P) systems (but as $n$ increases,
$a_{0}$ asymptotically tends to $1/2$), $\tilde{a_{i}}\geq 1/2$ in mutualistic
(M) systems (the case when $\tilde{a_{i}}=1/2$ only when the diagonal terms do
not play a role in the cycles contributing to destabilisation), and
$\tilde{a_{i}}=1/2$ in competitive (C) systems. Coefficients indicated by a
$+$ have no negative destabilising cycles and therefore $\tilde{a_{i}}=0$. The
total number of terms (denominator) for each $a_{i}$ is $n!/i!$. While the
numerators of the predator-prey and mutualistic $a_{i}$ do not show a coherent
sequence of growth for their tipping cycles for a given $n$ as $i$ increases,
the competitive systems follow the sequence of path polynomials of the
complete graph $K_{n}$ [12].

Consider the $4\times 4$ predator-prey community matrix,

$\left[\begin{array}[]{rrrr}-a&b&c&d\\\ -e&-f&g&h\\\ -k&-l&-m&p\\\
-q&-r&-s&-t\end{array}\right],$

then the coefficient $a_{1}$, consisting of 3-cycles ($n-i$) is,

cfk+bgk+agl+bem+afm+dfq+bhq+dmq+cpq+ahr+hmr+gpr+
aps+fps+bet+aft+ckt+glt+amt+fmt-cel-der-dks-hls.

We define the tipping cycle set of $a_{1}$ as
$\hat{a_{1}}=\\{cel,der,dks,hls\\}$. These four negative cycles
($\tilde{a_{1}}=4/24$), that if strong enough, could tip the $a_{1}$
coefficient from positive to negative, destabilising the system, consist of
only eight distinct matrix elements, $c,d,e,h,k,l,r$ and $s$, with each
element included in a differing number of cycles. That is, each element has a
specific weight. If we construct a weighted matrix of the elements of this
tipping cycle set ($[\hat{a_{1}}]$) one can clearly see the patterning of the
key elements that play a role in destabilisation,

$[\hat{a_{1}}]=\left[\begin{array}[]{rrrr}0&0&1&2\\\ 2&0&0&1\\\ 1&2&0&0\\\
0&1&2&0\end{array}\right].$

Applied to all the different forms, we can see that the sign structure of the
matrix has an effect on the different elements in their role in the
destabilising cycles,

$[\hat{a_{0}}]=$ Predator-Prey Mutualistic Competitive
$\begin{bmatrix}1&1&3&3\\\ 3&1&1&3\\\ 3&3&1&1\\\ 1&3&3&1\end{bmatrix}$
$\begin{bmatrix}5&5&5&5\\\ 5&5&5&5\\\ 5&5&5&5\\\ 5&5&5&5\end{bmatrix}$
$\begin{bmatrix}3&3&3&3\\\ 3&3&3&3\\\ 3&3&3&3\\\ 3&3&3&3\end{bmatrix}$

$[\hat{a_{1}}]=$ $\begin{bmatrix}0&0&1&2\\\ 2&0&0&1\\\ 1&2&0&0\\\
0&1&2&0\end{bmatrix}$ $\begin{bmatrix}3&4&4&4\\\ 4&3&4&4\\\ 4&4&3&4\\\
4&4&4&3\end{bmatrix}$ $\begin{bmatrix}3&2&2&2\\\ 2&3&2&2\\\ 2&2&3&2\\\
2&2&2&3\end{bmatrix}$

$[\hat{a_{2}}]=$ $\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}$ $\begin{bmatrix}0&1&1&1\\\ 1&0&1&1\\\ 1&1&0&1\\\
1&1&1&0\end{bmatrix}$ $\begin{bmatrix}0&1&1&1\\\ 1&0&1&1\\\ 1&1&0&1\\\
1&1&1&0\end{bmatrix}$

In addition, the total weighting of the elements in their role in
destabilisation across all the coefficients are a summation of the individual
coefficient tipping cycle set weightings,

$\sum_{i=0}^{n}[\hat{a_{i}}]=\begin{tabular}[]{l c c c}&Predator-
Prey&Mutualistic&Competitive\\\ &$\begin{bmatrix}1&1&4&5\\\ 5&1&1&4\\\
4&5&1&1\\\ 1&4&5&1\end{bmatrix}$&$\begin{bmatrix}8&10&10&10\\\ 10&8&10&10\\\
10&10&8&10\\\ 10&10&10&8\end{bmatrix}$&$\begin{bmatrix}6&6&6&6\\\ 6&6&6&6\\\
6&6&6&6\\\ 6&6&6&6\end{bmatrix}$ \vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\end{tabular}$

The symmetry seen in the competitive and mutualistic
$\sum_{i=0}^{n}[\hat{a_{i}}]$ are clearly a consequence of the symmetry in the
sign structure of these two forms, while for the predator-prey systems it
reflects its lack of sign symmetry. Even with sign symmetry, the mutualistic
weightings show the bias that comes in through the effect of moving away from
sign uniformity. Such bias in the weighting of certain individual elements
raises the question as to how much stability or instability is an artifact of
the structure. For example, in larger systems, positioning of one species over
another in the governing equations from which the community matrices are
constructed may well prove crucial for the outcome. That is, it posits the
question of whether the medium may well be the message.

If the elements comprising tipping cycles are considered in terms of set
inclusion, in which for a given $\hat{a_{i}}$ set comparability is like an
antichain [2], and the set inclusion properties of smaller to larger sized
cycles is likewise determined by the collection of elements, then we can
describe the inclusion properties of $\hat{a_{i}}$ as $i$ decreases for a
given $n$. To see this, we again look at the $4\times 4$ predator-prey system,
where the 4-cycles are
$\hat{a_{0}}=\\{chlq,dgkr,demr,cepr,dfks,bhks,ahls,celt\\}$. In terms of the
constituent elements, the set of 3-cycles of $\hat{a_{1}}$ are subsets of some
of the 4-cycles of $\hat{a_{0}}$ (i.e. cel $\in$ celt):

$\hat{a_{1}}\subset\hat{a_{0}}.$

This property, in which the smaller tipping cycle sets are a subset of the
larger sized sets holds across all the forms of matrices for a given $n$:

$\hat{a_{i}}\subset\hat{a}_{j\ \forall j<i}.$

Yet the way the different forms (i.e. predator-prey, mutualistic and
competitive) result in different element weightings is also reflected in the
different ways their $\hat{a_{i}}$ flow through their $a_{i}$ as $i$
decreases.

For example, with the $4\times 4$ predator-prey system, there are four tipping
cycles in $\hat{a_{1}}$ and eight tipping cycles in $\hat{a_{0}}$. While the
four smaller sets of $\hat{a_{1}}$ are subsets of sets in $\hat{a_{0}}$, the
four remaining sets of $\hat{a_{0}}$ do not have such an embedding from sets
of $\hat{a_{1}}$, but are new configurations of elements from the $4\times 4$
system (i.e.
$\\{cel,der,dks,hls\\}\subset\\{celt,demr,dfks,ahls\\}+\\{chlq,dgkr,cepr,bhks\\}$).

In the $4\times 4$ competitive case, the six sets in $\hat{a_{2}}$ are subsets
in two sets each in $\hat{a_{1}}$, but in this case there are no new
configurations. Since all the $\hat{a_{1}}$ consist of embedded sets from
$\hat{a_{2}}$, we may call $\hat{a_{1}}$ a closure of $\hat{a_{2}}$. Six of
the twelve sets of $\hat{a_{0}}$ consist of two sets each from $\hat{a_{1}}$
(these are sets constructed through the possible combination of having two of
the four diagonal elements, $\binom{4}{2}$), with six new configurations:

$\hat{a_{2}}$ | $\hat{a_{1}}$ | $\hat{a_{0}}$
---|---|---
$\left\\{\begin{array}[]{c}be\\\ ck\\\ gl\\\ dq\\\ hr\\\ ps\end{array}\right\\}\subset$ | $\left\\{\begin{array}[]{c}bem\\\ bet\\\ cfk\\\ ckt\\\ agl\\\ glt\\\ dfq\\\ dmq\\\ ahr\\\ hmr\\\ aps\\\ fps\\\ \end{array}\right\\}\subset$ |  $\left\\{\begin{array}[]{c}bemt\\\ cfkt\\\ aglt\\\ dfmq\\\ ahmr\\\ afps\end{array}\right\\}+$ $\left\\{\begin{array}[]{c}chlq\\\ bgpq\\\ dgkr\\\ cepr\\\ bhks\\\ dels\end{array}\right\\}$

With the $4\times 4$ mutualistic system, there are six tipping cycles in
$\hat{a_{2}}$, and as with the competitive case, each of these are in two
tipping cycles in $\hat{a_{1}}$. However, unlike the competitive case,
$\hat{a_{1}}$ is not a closure of $\hat{a_{2}}$. There are eight new tipping
cycles that have none of the $\hat{a_{2}}$ as subsets. In the next tipping
cycle set, $\hat{a_{0}}$, it follows the compeitive case in having six sets,
each of which consists of two sets from $\hat{a_{1}}$, and six sets that are
new configurations (the same new configurations as in the competitive case).
In addition, there are eight cycles each of which has as a subset one of the
eight new configurations in $\hat{a_{1}}$. These eight cycles begin what might
be called a parallel inclusion pattern, complicating the flow yet further:

$\hat{a_{2}}$ | $\hat{a_{1}}$ | $\hat{a_{0}}$
---|---|---
$\left\\{\begin{array}[]{c}be\\\ ck\\\ gl\\\ dq\\\ hr\\\ ps\end{array}\right\\}\subset$ |  $\left\\{\begin{array}[]{c}bem\\\ bet\\\ cfk\\\ ckt\\\ agl\\\ glt\\\ dfq\\\ dmq\\\ ahr\\\ hmr\\\ aps\\\ fps\\\ \end{array}\right\\}+$ $\left\\{\begin{array}[]{c}bgk\\\ cel\\\ bhq\\\ cpq\\\ der\\\ gpr\\\ dks\\\ hls\end{array}\right\\}\subset$ |  $\left\\{\begin{array}[]{c}bemt\\\ cfkt\\\ aglt\\\ dfmq\\\ ahmr\\\ afps\end{array}\right\\}+$ $\left\\{\begin{array}[]{c}bgkt\\\ celt\\\ bhmq\\\ cfpq\\\ demr\\\ agpr\\\ dfks\\\ ahls\end{array}\right\\}+$ $\left\\{\begin{array}[]{c}chlq\\\ bgpq\\\ dgkr\\\ cepr\\\ bhks\\\ dels\end{array}\right\\}$

Thus, the differences in the specific sign structure of the community matrices
alters the flow of the inclusion properties of the tipping cycle sets as their
cycle size increases. A full characterisation of the flows in general should
yield the differences such as reflected in the $\tilde{a_{i}}$ of the
different forms, yet in a much more intricate way.

## Self-Regulating Cycles

The importance of self-regulation for stability has long been understood [7],
and mechanistically it plays a significant role in community matrices as the
self-regulating terms are those that lie on the diagonal (see [14] for a
discussion on the historical development of its role in analyses). By making
the diagonal input large enough one can force any system (the one restriction
is that the diagonal terms are negative, although see below) to become stable
through diagonal dominance [4]. In this section, we explore what information
about a system can be determined through applying a variable to the diagonal
of a community matrix.

Consider for example the $3\times 3$ predator-prey community matrix to which
we apply a complex variable $\sigma$ on the diagonal elements,

$\left[\begin{array}[]{rrr}-\sigma a&b&c\\\ -d&-\sigma e&f\\\ -g&-h&-\sigma
k\end{array}\right].$

The characteristic polynomial for this system is

$x^{3}+\sigma(a+e+k)x^{2}+(\sigma^{2}(ae+ak+ek)+bd+cg+fh)x+$
$\sigma^{3}aek+\sigma(ceg+afh+bdk)+bfg-cdh$,
where the coefficients of the characteristic polynomial in $x$ are themselves
polynomials in $\sigma$.

In general, for matrices of size $n\times n$ where the values $m_{i,j}$ are
real valued and taking any sign but where the diagonal $m_{i,i}$ are negative
(although this condition can be relaxed [13])

$\mathbf{M}=\begin{pmatrix}-m_{1,1}&m_{1,2}&\cdots&m_{1,n}\\\
m_{2,1}&-m_{2,2}&\cdots&m_{2,n}\\\ m_{3,1}&m_{3,2}&\cdots&m_{3,n}\\\
\vdots&\vdots&\ddots&\vdots\\\ m_{n,1}&m_{n,2}&\cdots&-m_{n,n}\end{pmatrix}.$

With a complex variable $\sigma$ applied to the diagonal,

$\mathbf{M_{\sigma}}=\begin{pmatrix}-\sigma m_{1,1}&m_{1,2}&\cdots&m_{1,n}\\\
m_{2,1}&-\sigma m_{2,2}&\cdots&m_{2,n}\\\ m_{3,1}&m_{3,2}&\cdots&m_{3,n}\\\
\vdots&\vdots&\ddots&\vdots\\\ m_{n,1}&m_{n,2}&\cdots&-\sigma
m_{n,n}\end{pmatrix},$

$\mathbf{M_{\sigma}}$ has the monic characteristic polynomial

(2)
$\sum_{i=0}^{n}p_{i}x^{i}=\sum_{i=0}^{n}(\sum_{j=0}^{n-i}q_{j}\sigma^{j})x^{i}$

with $p_{n}=1$ and $p_{i}=\sum_{j=0}^{n-i}q_{j}\sigma^{j}$ where the $q_{j}$
consist of terms of length $n-i$.

One could construct something similar to Table 2 listing the number of
positive and negative cycles affected by $\sigma$ that contribute to stability
or instability in the three forms previously discussed. The total number of
such positive and negative self-regulating cycles will consist of all cycles
in the system except those that comprise the coefficients of $\sigma^{0}$. But
the value of such counting would be limited as it would fail to capture the
non-linearity of $\sigma$ and therefore the possible importance of certain
cycles over others for the resulting maximal eigenvalue, and therefore the
stability.

When the maximal root, and therefore the maximal eigenvalue of the
corresponding system, of Eq. 2 is real-valued, the point of stability
coincides with the last coefficient becoming positive, as we know from
Obrechkoff [10]. We can solve each individual coefficient polynomial in
$\sigma$, and those roots that are real-valued indicate the values of $\sigma$
that when applied to the matrix $\mathbf{M_{\sigma}}$, are the points at which
the sign of the corresponding coefficient changes. The largest real-valued
root of each coefficient polynomial indicates the value of $\sigma$ at which
the coefficient becomes positive and remains positive for any larger value of
$\sigma$ applied to the matrix [13].

There is an interesting relation [9],[14],[13] between the matrix $\mathbf{M}$
and a new matrix formed by dividing the diagonal elements ($\mathbf{D}$) of
$\mathbf{M}$ and adding the identity matrix,

$\bar{\mathbf{M}}_{0}=-(\mathbf{D})^{-1}\mathbf{M}+\mathbf{I}=\begin{pmatrix}0&\frac{m_{1,2}}{m_{1,1}}&\cdots&\frac{m_{1,n}}{m_{1,1}}\\\
\frac{m_{2,1}}{m_{2,2}}&0&\cdots&\frac{m_{2,n}}{m_{2,2}}\\\
\frac{m_{3,1}}{m_{3,3}}&\frac{m_{3,2}}{m_{3,3}}&\cdots&\frac{m_{3,n}}{m_{3,3}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{m_{n,1}}{m_{n,n}}&\frac{m_{n,2}}{m_{n,n}}&\cdots&0\end{pmatrix}.$

If one uses an eigenvalue (ordinarily the maximal, $\lambda_{max}$) of
$\bar{\mathbf{M}}_{0}$ as $\sigma$ in $\mathbf{M_{\sigma}}$ then one of the
eigenvalues of $\mathbf{M_{\sigma=\lambda_{max}}}$ will be zero. This is due
to the characteristic polynomial of $\bar{\mathbf{M}}_{0}$ in $x$ being
equivalent to the determinant of $\mathbf{M_{\sigma}}$ in $\sigma$, and
therefore the same as the last coefficient in Eq. 2 [13]. For example, in the
$3\times 3$ predator-prey case, the scaled predator-prey matrix would be

$\left[\begin{array}[]{rrr}0&b/a&c/a\\\ -d/e&0&f/e\\\
-g/k&-h/k&0\end{array}\right],$

with the characteristic polynomial

(3) $x^{3}+x(ceg+afh+bdk)/aek+(bfg-cdh)/aek,$

which is the same as the $x^{0}$ term in $\sigma$ for the $3\times 3$
predator-prey case shown earlier,

(4) $\sigma^{3}aek+\sigma(ceg+afh+bdk)+bfg-cdh.$

Therefore, when a root of the polynomial in $x$ of Eq. 3 is chosen as $\sigma$
in Eq. 4, the latter also vanishes (becomes zero).

In general, when an eigenvalue, $\lambda$, of $\bar{\mathbf{M}}_{0}$ is chosen
as $\sigma$, whether real- or complex-valued, then the last coefficient of the
characteristic polynomial of $\mathbf{M_{\sigma=\lambda}}$ vanishes. With the
$p_{0}$ term in Eq. 2 vanishing, the remaining polynomial is

$x^{n}+p_{n-1}(\sigma)x^{n-1}+...+p_{2}(\sigma)x^{2}+p_{1}(\sigma)x=x(x^{n-1}+p_{n-1}(\sigma)x^{n-2}+...p_{2}(\sigma)x+p_{1}(\sigma)).$

Pulling an $x$ from the equation shows clearly that there is a root of zero,
with the remaining polynomial in $x$ of degree $n-1$.

If the maximal eigenvalue, $\lambda_{max}$, of $\bar{\mathbf{M}}_{0}$ is used
as $\sigma$ then if the maximal eigenvalue of
$\mathbf{M_{\sigma=\lambda_{max}}}$ is zero, it is clear that the $p_{0}$
coefficient, consisting of cycles of size $n$, contains the key cycle (for
simplicity, assume there is one key feedback) that tips the system from an
unstable to a stable state (and vice versa). If, however, the maximal
eigenvalue of $\mathbf{M_{\sigma=\lambda_{max}}}$ is non-zero, therefore with
the zero eigenvalue not maximal, it is necessary to consider the cases of
whether $\lambda_{max}$ of $\bar{\mathbf{M}}_{0}$ is real- or complex-valued
and whether the maximal eigenvalue of $\mathbf{M_{\sigma=\lambda_{max}}}$ is
real- or complex-valued (with positive real part).

If $\lambda_{max}$ of $\bar{\mathbf{M}}_{0}$ is real-valued and the maximal
eigenvalue of $\mathbf{M_{\sigma=\lambda_{max}}}$ is real-valued and positive,
then one (or more) of the coefficients of the remaining $(n-1)^{th}$-order
polynomial has a real-valued root that is larger than the largest real-valued
root of $p_{0}$. Any such coefficient in the $(n-1)^{th}$-order polynomial
with a larger real-valued root than in $p_{0}$ would be negative. This means
that the key tipping cycle is smaller in size than an $n$-cycle, and of the
size of the cycles comprising the negative coefficient with the largest real-
valued root. The value of $\sigma$ needed in $\mathbf{M_{\sigma}}$ to make the
system stable is this larger real-valued root.

If $\lambda_{max}$ of $\bar{\mathbf{M}}_{0}$ is real-valued and the maximal
eigenvalue of $\mathbf{M_{\sigma=\lambda_{max}}}$ is complex but with positive
real part, then the key tipping cycle is less than $n$. But the remaining
$(n-1)^{th}$-order polynomial may have coefficients that are all positive.
Increasing the value of $\sigma$, and therefore increasing the input of the
self-regulating cycles, will enable one to find the value of $\sigma$ at the
tipping point, but not necessarily pinpoint which sized cycle is key.

Finally, if $\lambda_{max}$ of $\bar{\mathbf{M}}_{0}$ is complex, then the
characteristic polynomial of $\mathbf{M_{\sigma=\lambda_{max}}}$ is a complex
polynomial. $p_{0}$ in Eq. 2 will have vanished but if the maximal eigenvalue
of $\mathbf{M_{\sigma=\lambda_{max}}}$ has positive real part (therefore the
size of the key tipping cycle is less than $n$), then as with the preceding
case, $\sigma$ may be adjusted to find the tipping point, but determining
which sized cycle less than $n$ is key will require another approach.

## Conclusion

The idea of a tipping point requires the context of which feedback is causing
a system’s instability. This paper explores community matrices that are widely
used to study ecological systems, detailing the feedbacks that push a system
to its tipping point at the level of the sign structure of the elements
comprising the matrices. Sign structure plays a significant role not only in
terms of determining how many cycles contribute to instability, but also in
biasing certain of the elements comprising the matrices. This bias extends to
the way the elements, if treated as sets, flow through the feedback loops in a
system as they increase in size.

Given a system, a variable can be applied to the self-regulating terms on the
diagonal. Not only can this be used under certain conditions to identify
analytically how much adjustment of self-regulation would be needed to bring a
system to its tipping point, but through the use of a scaling method, helps to
isolate the size of the key feedback producing instability.

The results and methods described in this paper are not limited to systems
ecology, but should find wider application in the stability of linear systems.

## References

  * [1] B. Brooks, The coefficients of the characteristic polynomial in terms of the eigenvalues and the elements of an n$\times$ n matrix, Applied mathematics letters, 19 (2006), pp. 511–515.
  * [2] K. Engel, Sperner Theory, Encyclopedia of Mathematics and its Applications, 65, Cambridge: Cambridge University Press, 1997.
  * [3] W. Geary, M. Bode, T. Doherty, E. Fulton, D. Nimmo, A. Tulloch, V. Tulloch, and E. Ritchie, A guide to ecosystem models and their environmental applications, Nature Ecology and Evolution, (2020).
  * [4] S. Gerschgorin, Über die abgrenzung der eigenwerte einer matrix, Izv. Akad. Nauk. USSR Otd. Fiz.-Mat., 6 (1931), pp. 749–754.
  * [5] A. Hurwitz, Über die bedingungen, unter welchen eine gleichung nur wurzeln mit negativen reellen theilen besitzt, Math. Ann., 46 (1895), pp. 273–284.
  * [6] A. Lotka, Elements of physical biology, Williams and Wilkins, Baltimore, Maryland (Reissued in 1956 by Dover), 1925.
  * [7] D. Martin, Self-regulation in living systems: Symposium at ottawa, Nature, 183 (1959), p. 370.
  * [8] R. May, Will a large complex system be stable?, Nature, 238 (1972), p. 413.
  * [9] A.-M. Neutel and M. Thorne, Interaction strengths in balanced carbon cycles and the absence of a relation between ecosystem complexity and stability, Ecology Letters, 17 (2014), pp. 651–661.
  * [10] N. Obrechkoff, Sur un problème de laguerre, C.R. Acd. Sci. Paris, 177 (1923), pp. 223–235.
  * [11] E. Routh, Treatise on the Stability of a Given State of Motion: Particularly Steady Motion, Macmillan and Co., London., 1877.
  * [12] N. Sloane, The on-line encyclopedia of integer sequences, https://oeis.org/.
  * [13] M. Thorne, $\sigma$-stable matrices, arxiv.org/abs/2012.04601, (2020).
  * [14] M. Thorne, E. Forgoston, L. Billings, and A.-M. Neutel, Matrix scaling and tipping points, In preparation, SIAM Journal of Applied Dynamical Systems.


11institutetext: The Chinese University of Hong Kong (Shenzhen), China
22institutetext: Shenzhen Research Institute of Big Data, Shenzhen, China
33institutetext: The Hong Kong University of Science and Technology
(Guangzhou), China 33email<EMAIL_ADDRESS>
44institutetext: Software School, Fudan University, Shanghai, China
55institutetext: Johns Hopkins University, Baltimore, USA

# MetaLR: Meta-tuning of Learning Rates for Transfer Learning in Medical
Imaging

Yixiong Chen 1122 Li Liu${}^{(\textrm{\Letter})}$ 33 Jingxian Li 44 Hua Jiang
1122
Chris Ding 11 Zongwei Zhou 55

###### Abstract

In medical image analysis, transfer learning is a powerful method for deep
neural networks (DNNs) to generalize well on limited medical data. Prior
efforts have focused on developing pre-training algorithms on domains such as
lung ultrasound, chest X-ray, and liver CT to bridge domain gaps. However, we
find that model fine-tuning also plays a crucial role in adapting medical
knowledge to target tasks. The common fine-tuning method is manually picking
transferable layers (e.g., the last few layers) to update, which is labor-
expensive. In this work, we propose a meta-learning-based LR tuner, named
MetaLR, to make different layers automatically co-adapt to downstream tasks
based on their transferabilities across domains. MetaLR learns appropriate LRs
for different layers in an online manner, preventing highly transferable
layers from forgetting their medical representation abilities and driving less
transferable layers to adapt actively to new domains. Extensive experiments on
various medical applications show that MetaLR outperforms previous state-of-
the-art (SOTA) fine-tuning strategies. Codes are released.

###### Keywords:

Medical image analysis Meta-learning Transfer learning.

##  1 Introduction

Transfer learning has become a standard practice in medical image analysis as
collecting and annotating data in clinical scenarios can be costly. The pre-
trained parameters endow better generalization to DNNs than the models trained
from scratch [21, 8]. A popular approach to enhancing model transferability is
by pre-training on domains similar to the targets [27, 26, 19, 9, 25].
However, utilizing specialized pre-training for all medical applications
becomes impractical due to the diversity between domains and tasks and privacy
concerns related to pre-training data. Consequently, recent work [13, 2, 20,
6] has focused on improving the generalization capabilities of existing pre-
trained DNN backbones through fine-tuning techniques.

Figure 1: The motivation of MetaLR. Previous works fix transferable layers in
pre-trained models to prevent them from catastrophic forgetting. It is
inflexible and labor-expensive for this method to find the optimal scheme.
MetaLR uses meta-learning to automatically optimize layer-wise LR for fine-
tuning.

Previous studies have shown that the transferability of shallower layers is
often higher than that of deeper layers [24]. Layer-wise fine-tuning [21],
which is a practical application of this finding, was introduced to preserve
the transferable low-level knowledge. Moreover, recent studys in medical fine-
tuning have revealed that the transferability of pre-trained knowledge
significantly varies across different downstream datasets [23] and layers
[22]. In some cases, transferability might even be irregular among layers for
domains far from pre-training data [7]. Given the diverse medical domains and
model architectures, there is currently no universal guideline to follow to
determine whether a particular layer should be retrained for a given target
domain.

To search for optimal layer combinations for fine-tuning, manually selecting
transferable layers [21, 2] can be a solution, but it requires a significant
amount of human labor and computational cost. In order to address this issue
and improve the flexibility of fine-tuning strategies, we propose controlling
the fine-tuning process with layer-wise learning rates (LRs), rather than
simply manually fixing or updating the layers (see Fig. 1). Our proposed
algorithm, Meta Learning Rate (MetaLR), is based on meta-learning [12] and
adaptively adjusts LRs for each layer according to transfer feedback. It
treats the layer-wise LRs as meta-knowledge and optimizes them to improve the
model generalization. Larger LRs indicate less transferability of
corresponding layers and require more updating, while smaller LRs preserve
transferable knowledge in the layers. Inspired by [18], we use an online
adaptation strategy of LRs with a time complexity of $O(n)$, instead of the
computationally-expensive bi-level $O(n^{2})$ meta-learning. We also enhance
the algorithm’s performance and stability with a proportional hyper-LR (LR for
LR) and a validation scheme on training data batches.

In summary, this work makes the following three contributions. 1) We introduce
MetaLR, a meta-learning-based LR tuner that can adaptively adjust layer-wise
LRs based on transfer learning feedback from various medical domains. 2) We
enhance MetaLR with a proportional hyper-LR and a validation scheme using
batched training data to improve the algorithm’s stability and efficacy. 3)
Extensive experiments on both lesion detection and tumor segmentation tasks
were conducted to demonstrate the superior efficiency and performance of
MetaLR compared to current SOTA medical fine-tuning techniques.

##  2 Method

This section provides a detailed description of the proposed MetaLR. It is a
meta-learning-based [12, 17] approach that determines the appropriate LR for
each layer based on its transfer feedback. It is important to note that fixing
transferable layers is a special case of this method, where fixed layers
always have zero LRs. First, we present the theoretical formulation of MetaLR.
Next, we discuss online adaptation for efficiently determining optimal LRs.
Finally, we demonstrate the use of a proportional hyper-LR and a validation
scheme with batched training data to enhance performance.

###  2.1 Formulation of Meta Learning Rate

Let $(x,y)$ denotes a sample-label pair, and $\\{(x_{i},y_{i}):i=1,...,N\\}$
be the training data. The validation dataset
$\\{(x^{v}_{i},y^{v}_{i}):i=1,...,M\\}$ is assumed to be independent and
identically distributed as the training dataset. Let $\hat{y}=\Phi(x,\theta)$
be the prediction for sample $x$ from deep model $\Phi$ with parameters
$\theta$. In standard training of DNNs, the aim is to minimize the expected
risk for the training set: $\frac{1}{N}\sum^{N}_{i=1}L(\hat{y}_{i},y_{i})$
with fixed training hyper-parameters, where $L(\hat{y},y)$ is the loss
function for the current task. The model generalization can be evaluated by
the validation loss $\frac{1}{M}\sum^{M}_{i=1}L(\hat{y}^{v}_{i},y^{v}_{i})$.
Based on the generalization, one can tune the hyper-parameters of the training
process to improve the model. The key idea of MetaLR is considering the layer-
wise LRs as self-adaptive hyper-parameters during the training and
automatically adjusting them to achieve better model generalization. We denote
the LR and model parameters for the layer $j$ at the iteration $t$ as
$\alpha_{j}^{t}$ and $\theta_{j}^{t}$. The LR scheduling scheme
$\alpha=\\{\alpha_{j}^{t}:j=1,...,d;~{}t=1,...,T\\}$ is what MetaLR wants to
learn, affecting which local optimal $\theta^{*}(\alpha)$ the model parameters
$\theta^{t}=\\{\theta_{j}^{t}:j=1,...,d\\}$ will converge to. The optimal
parameters $\theta^{*}(\alpha)$ are given by optimization on the training
data. At the same time, the best LR tuning scheme $\alpha^{*}$ can be
optimized based on the feedback for $\theta^{*}(\alpha)$ from the validation
loss. This problem can be formulated as the following bi-level optimization
problem:

$\displaystyle\min_{\alpha}\frac{1}{M}\sum^{M}_{i=1}L(\Phi(x^{v}_{i},\theta^{*}(\alpha)),y^{v}_{i}),$
(1) $\displaystyle
s.t.~{}\theta^{*}(\alpha)=\mathop{\arg\min}\limits_{\theta}\frac{1}{N}\sum^{N}_{i=1}L(\Phi(x_{i},\theta),y_{i}).$

MetaLR aims to use the validation set to optimize $\alpha$ through an
automatic process rather than a manual one. The optimal scheme $\alpha^{*}$
can be found by a nested optimization [12], but it is too computationally
expensive in practice. A faster and more lightweight method is needed to make
it practical.

###  2.2 Onilne Learning Rate Adaptation

Inspired by the online approximation [18], we propose efficiently adapting the
LRs and model parameters online. The motivation of the online LR adaptation is
updating the model parameters $\theta^{t}$ and LRs
$\\{\alpha^{t}_{j}:j=1,2,...,d\\}$ within the same loop. We first inspect the
descent direction of parameters $\theta^{t}_{j}$ on the training loss
landscape and adjust the $\alpha^{t}_{j}$ based on the transfer feedback.
Positive feedback (lower validation loss) means the LRs are encouraged to
increase.

Algorithm 1 Online Meta Learning Rate Algorithm

0: Training data $\mathcal{D}$, validation data $\mathcal{D}^{v}$, initial
model parameter $\\{\theta^{0}_{1},...,\theta^{0}_{d}\\}$, LRs
$\\{\alpha^{0}_{1},...,\alpha^{0}_{d}\\}$, batch size n, max iteration T;

0: Final model parameter $\theta^{T}=\\{\theta^{T}_{1},...,\theta^{T}_{d}\\}$;

1: for $t=0:T-1$ do

2: $\\{(x_{i},y_{i}):i=1,...,n\\}\leftarrow$ TrainDataLoader($\mathcal{D}$, n)
;

3: $\\{(x^{v}_{i},y^{v}_{i}):i=1,...,n\\}\leftarrow$
ValidDataLoader($\mathcal{D}^{v}$, n) ;

4: Step forward for one step to get
$\\{\hat{\theta_{1}^{t}}(\alpha_{1}^{t}),...,\hat{\theta_{d}^{t}}(\alpha_{d}^{t})\\}$
with Eq. (2);

5: Update $\\{\alpha^{t}_{1},...,\alpha^{t}_{d}\\}$ to become
$\\{\alpha^{t+1}_{1},...,\alpha^{t+1}_{d}\\}$ with Eq. (3);

6: Update $\\{\theta^{t}_{1},...,\theta^{t}_{d}\\}$ to become
$\\{\theta^{t+1}_{1},...,\theta^{t+1}_{d}\\}$ with Eq. (4);

7: end for

We adopt Stochastic Gradient Descent (SGD) as the optimizer to conduct the
meta-learning. The whole training process is summarized in Algorithm 1. At the
iteration $t$ of training, a training data batch
$\\{(x_{i},y_{i}),i=1,...,n\\}$ and a validation data batch
$\\{(x^{v}_{i},y^{v}_{i}):i=1,...,n\\}$ are sampled, where n is the size of
the batches. First, the parameters of each layer are updated once with the
current LR according to the descent direction on training batch.

$\hat{\theta_{j}^{t}}(\alpha_{j}^{t})=\theta_{j}^{t}-\alpha_{j}^{t}\nabla_{\theta_{j}}(\frac{1}{n}\sum_{i=1}^{n}L(\Phi(x_{i},\theta_{j}^{t}),y_{i})),~{}j=1,...,d.$
(2)

This step of updating aims to get feedback for LR of each layer. After taking
derivative of the validation loss w.r.t. $\alpha_{j}^{t}$, we can utilize the
gradient to know how the LR for each layer should be adjusted. So the second
step of MetaLR is to move the LRs along the meta objective gradient on the
validation data:

$\alpha_{j}^{t+1}=\alpha_{j}^{t}-\eta\nabla_{\alpha_{j}}(\frac{1}{n}\sum_{i=1}^{n}L(\Phi(x_{i}^{v},\hat{\theta_{j}^{t}}(\alpha_{j}^{t})),y_{i}^{v})),$
(3)

where $\eta$ is the hyper-LR. Finally, the updated LRs can be employed to
optimize the model parameters through gradient descent truly.

$\theta_{j}^{t+1}=\theta_{j}^{t}-\alpha_{j}^{t+1}\nabla_{\theta_{j}}(\frac{1}{n}\sum_{i=1}^{n}L(\Phi(x_{i},\theta_{j}^{t}),y_{i})).$
(4)

For practical use, we constrain the LR for each layer to be
$\alpha_{j}^{t}\in[10^{-6},10^{-2}]$.

###  2.3 Proportional Hyper Learning Rate

In practice, LRs are often tuned in an exponential style (e.g., 1e-3, 3e-3,
1e-2) and are always positive values. However, if a constant hyper-LR is used,
it will linearly update its corresponding LR regardless of numerical
constraints. This can lead to fluctuations in the LR or even the risk of the
LR becoming smaller than 0 and being truncated. To address this issue, we
propose using a proportional hyper-LR $\eta=\beta\times\alpha_{j}^{t}$, where
$\beta$ is a pre-defined hyper-parameter. This allows us to rewrite Eq. (3)
as:

$\alpha_{j}^{t+1}=\alpha_{j}^{t}(1-\beta\nabla_{\alpha_{j}}(\frac{1}{n}\sum_{i=1}^{n}L(\Phi(x_{i}^{v},\hat{\theta_{j}^{t}}(\alpha_{j}^{t})),y_{i}^{v}))).$
(5)

The exponential update of $\alpha_{j}^{t}$ guarantees its numerical stability.

###  2.4 Generalizability Validation on Training Data Batch

One limitation of MetaLR is that the LRs are updated using separate validation
data, which reduces the amount of data available for the training process.
This can be particularly problematic for medical transfer learning, where the
amount of downstream data has already been limited. In Eq. 2 and Eq. 3, the
update of model parameter $\theta_{j}^{t}$ and LR $\alpha_{j}^{t}$ is
performed using different datasets to ensure that the updated $\theta_{j}^{t}$
can be evaluated for generalization without being influenced by the seen data.
As an alternative, but weaker, approach, we explore using another batch of
training data for Eq. 3 to evaluate generalization. Since this batch was not
used in the update of Eq. 2, it may still perform well for validation in meta-
learning. The effect of this approach is verified in Sec. 3.2, and the
differences between the two methods are analyzed in Sec. 3.4.

##  3 Experiments and Analysis

###  3.1 Experimental Settings

We extensively evaluate MetaLR on four transfer learning tasks (as shown in
Tab. 1). To ensure the reproducibility of the results, all pre-trained models
(USCL [9], ImageNet [11], C2L [26], Models Genesis [27]) and target datasets
(POCUS [5], BUSI [1], Chest X-ray [16], LiTS [4]) are publicly available. In
our work, we consider models pre-trained on both natural and medical image
datasets, with three target modalities and three target organs, which makes
our experimental results more credible. For the lesion detection tasks, we
used ResNet-18 [14] with the Adam optimizer. The initial learning rate (LR)
and hyper-LR coefficient $\beta$ are set to $10^{-3}$ and $0.1$, respectively.
In addition, we use 25% of the training set as the validation set for meta-
learning. For the segmentation task, we use 3D U-Net [10] with the SGD
optimizer. The initial LR and hyper-LR coefficient $\beta$ are set to
$10^{-2}$ and $3\times 10^{-3}$, respectively. The validation set for the LiTS
segmentation dataset comprises 23 samples from the training set of size 111.
All experiments are implemented using PyTorch 1.10 on an Nvidia RTX A6000 GPU.
We report the mean values and standard deviations for each experiment with
five different random seeds. For more detailed information on the models and
hyper-parameters, please refer to our supplementary material.

Table 1: Pre-training data, algorithms, and target tasks. Source | Pre-train Method | Target | Object | Task | Modality | Size
---|---|---|---|---|---|---
US-4 [9] |  USCL [9] |  POCUS [5] | Lung | COVID-19 detection | US | 2116 images
ImageNet [11] | supervised |  BUSI [1] | Breast | Tumor detection | US | 780 images
MIMIC-CXR [15] |  C2L [26] |  Chest X-ray [16] | Lung | Pneumonia detection | X-ray | 5856 images
LIDC-IDRI [3] |  Models Genesis [27] |  LiTS [4] | Liver | Liver segmentation | CT | 131 volumes

###  3.2 Ablation Study

In order to evaluate the effectiveness of our proposed method, we conduct an
ablation study w.r.t. the basic MetaLR algorithm, the proportional hyper-LR,
and batched-training-data validation (as shown in Tab. 2). When applying only
the basic MetaLR, we observe only marginal performance improvements for the
four downstream tasks. We conjecture that this is due to two reasons: Firstly,
the constant hyper-LR makes the training procedures less stable than direct
training, which is evident from the larger standard deviation of performance.
Secondly, part of the training data are split for validation, which can be
detrimental to the performance. After applying the proportional hyper-LR,
significant improvements are in both the performance and its stability.
Moreover, although the generalization validation on the training data batch
may introduce bias, providing sufficient training data ultimately benefits the
performance.

Table 2: Ablation study for MetaLR, hyper-LR, and validation data. The baseline is the direct tuning of all layers with constant LRs. The default setting for MetaLR is a constant hyper-LR of $10^{-3}$ and a separate validation set. MetaLR | Prop. hyper-LR | Val. on trainset | POCUS | BUSI | Chest X-ray | LiTS
---|---|---|---|---|---|---
|  |  |  91.6$\pm$0.8 |  84.4$\pm$0.7 |  94.8$\pm$0.3 |  93.1$\pm$0.4
✓ |  |  |  91.9$\pm$0.6 |  84.9$\pm$1.3 |  95.0$\pm$0.4 |  93.2$\pm$0.8
✓ | ✓ |  |  93.6$\pm$0.4 |  85.2$\pm$0.8 |  95.3$\pm$0.2 |  93.3$\pm$0.6
✓ |  | ✓ |  93.0$\pm$0.3 |  86.3$\pm$0.7 |  95.5$\pm$0.2 |  93.9$\pm$0.5
✓ | ✓ | ✓ |  93.9$\pm$0.4 |  86.7$\pm$0.7 |  95.8$\pm$0.3 |  94.2$\pm$0.5
$\bullet$ Final MetaLR outperforms baseline with p-values of
$0.0014,0.0016,0.0013,0.0054$.

###  3.3 Comparative Experiments

In our study, we compare MetaLR with several other fine-tuning schemes,
including tuning only the last layer / all layers with constant LRs, layer-
wise fine-tuning [21], bi-directional fine-tuning [7], and AutoLR [20]. The
U-Net fine-tuning scheme proposed by Mina et al.  [2] was also evaluated.

Table 3: Comparative experiments on lesion detection. We report sensitivities (%) of the abnormalities, overall accuracy (%), and training time on each task. Method | POCUS | BUSI | Chest X-ray
---|---|---|---
COVID |  Pneu. |  Acc |  Time |  Benign |  Malignant |  Acc |  Time |  Pneu. |  Acc |  Time
Last Layer |  77.9$\pm$2.1 |  84.0$\pm$1.3 |  84.1$\pm$0.2 | 15.8m |  83.5$\pm$0.4 |  47.6$\pm$4.4 |  66.8$\pm$0.5 | 4.4m |  99.7$\pm$1.3 |  87.8$\pm$0.6 | 12.7m
All Layers |  85.8$\pm$1.7 |  90.0$\pm$1.9 |  91.6$\pm$0.8 | 16.0m |  90.4$\pm$1.5 |  77.8$\pm$3.5 |  84.4$\pm$0.7 | 4.3m |  98.8$\pm$0.2 |  94.8$\pm$0.3 | 12.9m
Layer-wise |  87.5$\pm$1.0 |  92.3$\pm$1.3 |  92.1$\pm$0.3 | 2.4h |  90.8$\pm$1.2 |  75.7$\pm$2.6 |  85.6$\pm$0.4 | 39.0m |  97.9$\pm$0.3 |  95.2$\pm$0.2 | 1.9h
Bi-direc. |  90.1$\pm$1.2 |  92.5$\pm$1.5 |  93.6$\pm$0.2 | 12.0h |  92.2$\pm$1.0 |  77.1$\pm$3.5 |  86.5$\pm$0.5 | 3.2h |  98.4$\pm$0.3 |  95.4$\pm$0.1 | 9.7h
AutoLR |  89.8$\pm$1.6 |  89.7$\pm$1.5 |  90.4$\pm$0.8 | 17.5m |  90.4$\pm$1.8 |  76.2$\pm$3.2 |  84.9$\pm$0.8 | 4.9m |  95.4$\pm$0.5 |  93.0$\pm$0.8 | 13.3m
MetaLR |  94.8$\pm$1.2 |  93.1$\pm$1.5 |  93.9$\pm$0.4 | 24.8m |  92.2$\pm$0.7 |  75.6$\pm$3.6 |  86.7$\pm$0.7 | 6.0m |  97.4$\pm$0.4 |  95.8$\pm$0.3 | 26.3m

Results on Lesion Detection Tasks. MetaLR consistently shows the best
performance on all downstream tasks (Tab. 3). It shows 1% - 2.3% accuracy
improvements compared to direct training (i.e., tuning all layers) because it
takes into account the different transferabilities of different layers. While
manual picking methods, such as layer-wise and bi-directional fine-tuning,
also achieve higher performance, they require much more training time
(5$\times$ \- 50$\times$) for searching the best tuning scheme. On the other
hand, AutoLR is efficient, but its strong hypothesis harms its performance
sometimes. In contrast, MetaLR makes no hypothesis about transferability and
learns appropriate layer-wise LRs on different domains. Moreover, its
performance improvements are gained with only 1.5$\times$ \- 2.5$\times$
training time compared with direct training.

Table 4: Comparative experiments on LiTS liver segmentation task. Method |  PPV |  Sensitivity |  Dice |  Time
---|---|---|---|---
Last Layer |  26.1$\pm$5.5 |  71.5$\pm$4.2 |  33.5$\pm$3.4 | 2.5h
All Layers |  94.0$\pm$0.6 |  93.1$\pm$0.7 |  93.1$\pm$0.4 | 2.6h
Layer-wise |  92.1$\pm$1.3 |  96.4$\pm$0.4 |  93.7$\pm$0.3 | 41.6h
Bi-direc. |  92.4$\pm$1.1 |  96.1$\pm$0.2 |  93.8$\pm$0.1 | 171.2h
Mina et al. |  92.7$\pm$1.2 |  93.2$\pm$0.5 |  92.4$\pm$0.5 | 2.6h
MetaLR |  94.4$\pm$0.9 |  93.6$\pm$0.4 |  94.2$\pm$0.5 | 5.8h

Results on Segmentation Task. MetaLR achieves the best Dice performance on the
LiTS segmentation task (Tab. 4). Unlike ResNet for lesion detection, the U-Net
family has a more complex network topology. With skip connections, there are
two interpretations [2] of depths for layers: 1) the left-most layers are the
shallowest, and 2) the top layers of the “U” are the shallowest. This makes
the handpicking methods even more computationally expensive. However, MetaLR
updates the LR for each layer according to their validation gradients, and its
training efficiency is not affected by the complex model architecture.

###  3.4 Discussion and Findings

Figure 2: The LR curves for MetaLR on POCUS detection (a), on LiTS
segmentation (b), with constant hyper-LR (c), and with a separate validation
set (d).

The LRs learned by MetaLR change as the training progresses. For ResNet-18
(Fig. 2 (a)), the layer-wise LRs fluctuate drastically during the first 100
iterations. However, after iteration 100, all layers except the first layer
“Conv1” become stable at different levels. The first layer has a decreasing LR
(from $2.8\times 10^{-3}$ to $3\times 10^{-4}$) throughout the process,
reflecting its higher transferability. For 3D U-Net (Fig. 2 (b)), the middle
layers of the encoder “Down-128” and “Down-256” are the most transferable and
have the lowest LRs, which is difficult for previous fine-tuning schemes to
discover. As expected, the randomly initialized “FC” and “Out” layers have the
largest LRs since they are not transferable.

We also illustrate the LR curves with a constant hyper-LR instead of a
proportional one. The LR curves of “Block 3-1” and “Block 4-2” become much
more fluctuated (Fig. 2 (c)). This instability may be the key reason for the
instability of performance when using a constant hyper-LR. Furthermore, we
surprisingly find that the learned LRs are similar to the curves learned when
validated on the training set when using a separate validation set Fig. 2
(d)). With similar learned LR curves and more training data, it is reasonable
that batched training set validation can be an effective alternative to the
basic MetaLR.

##  4 Conclusion

In this work, we proposed a new fine-tuning scheme, MetaLR, for medical
transfer learning. It achieves significantly superior performance to the
previous SOTA fine-tuning algorithms. MetaLR alternatively optimizes model
parameters and layer-wise LRs in an online meta-learning fashion with a
proportional hyper-LR. It learns to assign lower LRs for the layers with
higher transferability and higher LRs for the less transferable layers. The
proposed algorithm is easy to implement and shows the potential to replace
manual layer-wise fine-tuning schemes. Future works include adapting MetaLR to
a wider variety of clinical tasks.

##  5 Acknowledgement

This work was supported by the National Natural Science Foundation of China
(No. 62101351) and the GuangDong Basic and Applied Basic Research Foundation
(No.2020A1515110376).

## References

  * [1] Al-Dhabyani, W., Gomaa, M., Khaled, H., Fahmy, A.: Dataset of breast ultrasound images. Data in brief 28, 104863 (2020)
  * [2] Amiri, M., Brooks, R., Rivaz, H.: Fine-tuning u-net for ultrasound image segmentation: different layers, different outcomes. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 67(12), 2510–2518 (2020)
  * [3] Armato III, S.G., McLennan, G., Bidaut, L., et al.: The lung image database consortium (lidc) and image database resource initiative (idri): a completed reference database of lung nodules on ct scans. Medical physics 38(2), 915–931 (2011)
  * [4] Bilic, P., Christ, P.F., Vorontsov, E., Chlebus, G., Chen, H., Dou, Q., Fu, C.W., Han, X., Heng, P.A., Hesser, J., et al.: The liver tumor segmentation benchmark (lits). arXiv preprint arXiv:1901.04056 (2019)
  * [5] Born, J., Wiedemann, N., Cossio, M., et al.: Accelerating detection of lung pathologies with explainable ultrasound image analysis. Applied Sciences 11(2), 672 (2021)
  * [6] Chambon, P., Cook, T.S., Langlotz, C.P.: Improved fine-tuning of in-domain transformer model for inferring covid-19 presence in multi-institutional radiology reports. Journal of Digital Imaging pp. 1–14 (2022)
  * [7] Chen, Y., Li, J., Ding, C., Liu, L.: Rethinking two consensuses of the transferability in deep learning. arXiv preprint arXiv:2212.00399 (2022)
  * [8] Chen, Y., Zhang, C., Ding, C.H., Liu, L.: Generating and weighting semantically consistent sample pairs for ultrasound contrastive learning. IEEE TMI (2022)
  * [9] Chen, Y., Zhang, C., Liu, L., Feng, C., Dong, C., Luo, Y., Wan, X.: Uscl: Pretraining deep ultrasound image diagnosis model through video contrastive representation learning. In: MICCAI. pp. 627–637. Springer (2021)
  * [10] Çiçek, Ö., Abdulkadir, A., Lienkamp, S.S., Brox, T., Ronneberger, O.: 3d u-net: learning dense volumetric segmentation from sparse annotation. In: MICCAI. pp. 424–432. Springer (2016)
  * [11] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: A large-scale hierarchical image database. In: CVPR. pp. 248–255 (2009)
  * [12] Franceschi, L., Frasconi, P., Salzo, S., Grazzi, R., Pontil, M.: Bilevel programming for hyperparameter optimization and meta-learning. In: ICML. pp. 1568–1577. PMLR (2018)
  * [13] Guo, Y., Shi, H., Kumar, A., Grauman, K., Rosing, T., Feris, R.: Spottune: transfer learning through adaptive fine-tuning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4805–4814 (2019)
  * [14] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: CVPR. pp. 770–778 (2016)
  * [15] Johnson, A.E., Pollard, T.J., Greenbaum, N.R., Lungren, M.P., Deng, C.y., Peng, Y., Lu, Z., Mark, R.G., Berkowitz, S.J., Horng, S.: Mimic-cxr-jpg, a large publicly available database of labeled chest radiographs. arXiv preprint arXiv:1901.07042 (2019)
  * [16] Kermany, D., Zhang, K., Goldbaum, M.: Large dataset of labeled optical coherence tomography (oct) and chest x-ray images. Mendeley Data 3, 10–17632 (2018)
  * [17] Li, Z., Zhou, F., Chen, F., Li, H.: Meta-sgd: Learning to learn quickly for few-shot learning. arXiv preprint arXiv:1707.09835 (2017)
  * [18] Ren, M., Zeng, W., Yang, B., Urtasun, R.: Learning to reweight examples for robust deep learning. In: ICML. pp. 4334–4343. PMLR (2018)
  * [19] Riasatian, A., Babaie, M., Maleki, D., Kalra, S., Valipour, M., Hemati, S., Zaveri, M., Safarpoor, A., Shafiei, S., Afshari, M., et al.: Fine-tuning and training of densenet for histopathology image representation using tcga diagnostic slides. Medical Image Analysis 70, 102032 (2021)
  * [20] Ro, Y., Choi, J.Y.: Autolr: Layer-wise pruning and auto-tuning of learning rates in fine-tuning of deep networks. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 35, pp. 2486–2494 (2021)
  * [21] Tajbakhsh, N., Shin, J.Y., Gurudu, S.R., Hurst, R.T., Kendall, C.B., Gotway, M.B., Liang, J.: Convolutional neural networks for medical image analysis: Full training or fine tuning? IEEE TMI 35(5), 1299–1312 (2016)
  * [22] Vrbančič, G., Podgorelec, V.: Transfer learning with adaptive fine-tuning. IEEE Access 8, 196197–196211 (2020)
  * [23] Wang, G., Li, W., Zuluaga, M.A., Pratt, R., Patel, P.A., Aertsen, M., Doel, T., David, A.L., Deprest, J., Ourselin, S., et al.: Interactive medical image segmentation using deep learning with image-specific fine tuning. IEEE TMI 37(7), 1562–1573 (2018)
  * [24] Yosinski, J., Clune, J., Bengio, Y., Lipson, H.: How transferable are features in deep neural networks? NeurIPS 27 (2014)
  * [25] Zhang, C., Chen, Y., Liu, L., Liu, Q., Zhou, X.: Hico: Hierarchical contrastive learning for ultrasound video model pretraining. In: ACCV. pp. 229–246 (2022)
  * [26] Zhou, H.Y., Yu, S., Bian, C., Hu, Y., Ma, K., Zheng, Y.: Comparing to learn: Surpassing imagenet pretraining on radiographs by comparing image representations. In: MICCAI. pp. 398–407 (2020)
  * [27] Zhou, Z., Sodha, V., Siddiquee, M.M.R., et al. : Models genesis: Generic autodidactic models for 3d medical image analysis. In: MICCAI. pp. 384–393 (2019)


longtable [a,b]Sofie Martins

# GPU-accelerated Higher Representations of Wilson Fermions with HiRep

Erik Kjellgren  Emiliano Molinaro  Claudio Pica  Antonio Rago

###### Abstract

We are improving one of the available lattice software packages HiRep by
adding GPU acceleration supporting highly-optimized simulations on both NVIDIA
and AMD GPUs. HiRep allows lattice simulations of theories with fermions in
higher representations and a variable number of colors in the gauge group. The
development is accompanied by an overall software quality improvement in the
build system, testing, and documentation, adding features for both CPUs and
GPUs. The software is available under https://github.com/claudiopica/HiRep.

##  1 Motivation

Predictions from models of physics beyond the standard model often depend on
input from strongly coupled gauge theories of variable numbers of colors and
fermion representations. Simulations on the lattice can quickly become
computationally expensive. Using Wilson Fermions constitutes, to this day, one
of the cheapest options to explore theories with different gauge groups. In
particular, the simulation of theories with higher dimensional fermion
representations at sufficient precision requires generally a large amount of
computational resources.

The lattice library HiRep allows simulations of higher representations of
Wilson fermions with general gauge groups. In the following we present results
for the porting of the library to state-of-the-art GPU accelerators. The
software is available at

https://github.com/claudiopica/HiRep

##  2 Software Quality

We improved the quality assurance of the software by introducing automated
unit tests, testing different numbers of colors and fermion representations,
MPI support, OpenMP support, GPU support, and clover improvement. for each
commit and code coverage reporting with codecov supplying the user with
detailed information about which lines of the code are tested. To ensure the
extensibility of the library, we generate documentation [6] including a user
manual, developer handbook and Doxygen function reference and improve code
readability using a code formatting check.

For a fast and simple compilation for different gauge groups and fermion
representations, we have updated the build system to ninja-build [1], a lower-
level build system that avoids large compilation times. The library is largely
independent of external modules or libraries, requiring gcc, perl and ninja-
build and optionally CUDA (GPU acceleration), MPI (parallelization) and hwloc
(hardware topology). Configuration of build setup and compilation variables is
supplied over a build help text.

New benchmarking code has been developed for gaining quick information on the
configuration of the software on a new cluster or supercomputer.

###  2.1 Memory Access Patterns

Lattice simulations are usually bound by the speed of the memory access. While
modern CPUs can access memory with high flexibility, access patterns are
crucial to optimal GPU performance. GPU manufacturers cite here a peak value
of global memory access, which can be reached by reading memory block after
block using sequential GPU threads. However, an ideal memory arrangement for
the GPU is not usually the most intuitive way to store lattice data.

Figure 1: Comparison of memory access patterns in CPU and GPU simulations

For simulations on the CPU, it is sufficient to store memory site-by-site as
in figure 1. However, this is a problem for reading memory from the GPU
because the kernels will only need to read the first component or first real
number of the site and then move on to the other components. Therefore, the
ideal memory access pattern is to store blocks of components as illustrated in
figure 1. Instead of components, we store the data in the smallest memory
units we read, single real numbers of single or double precision.

Unfortunately, this memory striding is incompatible with the legacy geometry
of HiRep. In the legacy geometry, we categorize sites of the local lattices
into sites that, during the hopping term application, do not depend on
communications (_bulk_), sites depending on communications (_boundary_), and
sites that the boundary depends on (_halo_ or _receive buffers_) and arrange
the sites of the four-dimensional lattice in such a way that the boundary and
receive buffers are contiguous in memory. This way, we can identify the
contiguous memory region and pass it to MPI.

On a four-dimensional lattice, this is rarely trivially possible. The simplest
solution to this problem is to collect the sites that need to be sent over MPI
into a separate send buffer, copying all sites on the boundary of the inner
lattice. The old geometry aims to minimize the number of copies by adding
duplicated sites in place and synchronizing them. Unfortunately, in
combination with the striding on the GPU, this is impossible because the
striding breaks the buffers into non-contiguous pieces.

Due to this, we implemented the new geometry using a complete copy to a send
buffer. This geometry has only two categories of sites: bulk and receive
buffers. Instead of waiting for communications before we compute the sites on
the boundary, we only mask the directions in the hopping term that can only be
executed after the communications are complete. This has the advantage that
more time is spent in the bulk computation, giving us more time to mask
communications. Find illustrations in figure 2, based on figures published in
[6].

Figure 2: New geometry (left) and legacy geometry (right), illustrated on 2D
lattices.

###  2.2 Send buffer synchronization

The new geometry has the disadvantage that the send and receive buffers are no
longer contiguous. Consequently, we need to synchronize the non-contiguous
pieces of information to a contiguous send buffer. This synchronization is
slightly more complicated than the synchronization on the legacy geometry;
however, the execution times of both operations are negligible compared to the
bulk calculations of the Dirac operator applications on typical local lattice
sizes.

Figure 3: Illustration of the send buffer synchronization in the new geometry

We are speeding up communications by partially applying the Dirac operator
before the synchronization to the send buffer. As a result, it is only
necessary to send a half-spinor, cutting the communication times
substantially, see for example [7].

##  3 Features

###  3.1 Wilson-Dirac Operator

The GPU-ported version of HiRep supports a GPU-accelerated Wilson-Dirac
operator in single and double precision for general higher representations of
SU($N_{c}$) for any number of colors. We are supporting even-odd
preconditioning, clover [17], exponentiated clover improvement [8], and the
Lüscher-Weisz gauge action [14]. Available are not only the HMC but also the
RHMC [13, 3], Hasenbusch acceleration [12] and a selection of integrators:
Leapfrog, 2nd and 4th order Omelyan integrators [16].

####  3.1.1 Strong and weak scaling

Figure 4: Weak scaling of time spent on a single application of the Wilson-
Dirac operator for an SU(2) gauge group with fermions in the fundamental
representation. Figure 5: Strong scaling of time spent on a single application
of the Wilson-Dirac operator for an SU(2) gauge group with fermions in the
fundamental representation.

We observe almost ideal weak scaling on larger clusters, even for multi-node
simulations. We show in figure 4 the weak scaling of the smallest possible
kernel: An SU(2) gauge group with fermions in the fundamental representation.
Both on LUMI-G and Tursa we see that the computations are masking the
communications for single and multi-node simulations.

The strong scaling in figure 5, however, is showing less ideal behavior. Both
on Tursa and LUMI-G the loss of efficiency is mainly caused by a strong
dependence of the kernel execution time on the total memory movement. In a
region where the local lattices become small, the execution time is dominated
by the kernel call overheads and ideal scaling is not expected.

####  3.1.2 Large-$N_{c}$ scaling

A crucial property of the library is the ability to use fermions in arbitrary
higher representations. For the CPU version the operation on the fields are
site-wise operations. As a result the computational throughput per thread or
process is increasing with the square of the dimension of the fermion
representation $N_{\mathrm{f}}^{2}$, the amount of computations in the matrix
multiplication of the represented gauge field with the spinor vector
components in the Dirac operator. Additionally, the amount of memory each
thread or process is reading increases, since the amount of data stored in a
single spinor and a single SU($N_{\mathrm{f}}$) matrix increases with $\sim
N_{\mathrm{f}}$ and $\sim N_{\mathrm{f}}^{2}$ respectively.

While on the CPU we pay a performance penalty proportional to the additional
computation and memory, the performance on the GPU can degrade very quickly
and become completely unfeasible due to register spilling: In contrast to CPUs
GPUs do not have an L3 cache and if too much memory is used in the local
thread, it ‘spills’ directly to global memory. Accessing global memory,
however, is substantially slower.

In order to allow scaling of the code to large fermion representations, we
have to parallelize in the fermion dimension. However, this parallelization is
non-trivial due to the matrix multiplication. Correspondingly, a
large-$N_{\mathrm{f}}$-improved kernel might not perform better than the naive
implementation for small representations such as the ones of QCD-like
theories. For this, we implemented two versions of this kernel, one is
parallelized. It can be activated by changing the compilation variables, which
is implemented by passing macro definitions to the preprocessor. The macro
LARGE_N switches to the parallelized kernel. However, on some cards the non-
parallelized kernel still performs better for QCD-like theories. As a result,
the default is the non-parallelized kernel.

Figure 6: Large-$N_{c}$ scaling of performance metrics of applications of the
Wilson-Dirac operator with an SU($N_{c}$) gauge group and fermions in the
fundamental representation.

Figure 6 shows the execution time, bandwidth, and computational throughput as
a function of the dimension of the fermion representation $N_{\mathrm{f}}$. To
avoid effects due to different memory sizes, we scaled the lattice size down
when increasing $N_{\mathrm{f}}$ to keep the total amount of memory moved
constant. The values given are interpolated since this was impossible due to
geometric constraints.

The figures show a steep increase in execution time for the unimproved kernel
for large gauge groups. While the scaling in time for the improved kernel is
not ideal, it is more stable. The computational throughput of the kernel is
largely constant. The drop in bandwidth can be explained by the fact that the
kernel becomes compute-bound for larger $N_{\mathrm{f}}$. Still, the peak
theoretical occupancy is limited because the kernel is not particularly
lightweight and needs many registers. Consequently, the card is exhausting its
capabilities in terms of computation for all $N_{\mathrm{f}}$. This
additionally causes a drop in the effective bandwidth.

###  3.2 Linear Algebra

An efficient inversion of the Wilson-Dirac operator further depends on
efficient linear algebra operations, particularly for the spinors. The single
and double-precision implementations of all fields available in the library
are easily extendable, simplify use, and reach 80-90% of the theoretical peak
performance on NVIDIA GPUs.

Figure 7: Bandwidths reached for linear algebra operations on different types
of NVIDIA GPUs. Tested for an SU(3) gauge group with fermions in the
fundamental representation.

In Figure 7, we test the achieved bandwidth of selected linear algebra
operations on a single device for different NVIDIA GPUs. We find that
approximately 90% of the peak performance is reached for simple operations
such as the identity operation (_Copy spinor field_) or a $\gamma_{5}$
application (_g5 application_) for sufficiently high memory movement. We see
that there is a dependence on the amount of data we moved since the kernels
have a call overhead and can only reach peak values if the GPUs have
sufficient work to do for this overhead to be small. More intensive operations
such as a multiply-assign with a complex number and application of
$\gamma_{5}$ (_g5 mulc add assign_) are also saturating capabilities of the
device.

Another critical operation type are _reductions_ , which map the field to a
single number, such as inner products and norms. These are implemented to take
the site-wise product or norms and reduce the result. Consequently, they
consist of multiple kernels in sequence, each needing to read a substantial
amount of data. It is expected that these operations do not saturate the
theoretical peak because it accounts only for a first single read. However,
the necessity of calling two kernels causes two stages where global memory
needs to be accessed.

While the site-wise operations are kernels that are part of the HiRep library,
the reductions are offloaded to cub, a library that has been part of the CUDA
SDK since CUDA 11.

###  3.3 Inverters

HiRep supports the Conjugate Gradient, stabilized Bi-Conjugate Gradient [19,
11], and QMR$\gamma_{5}$ inverters, which are ported to GPUs. In the future,
we are planning to implement domain decomposition methods, such as in [15, 2,
9] that reduce the amount of communications during the inversion additionally,
which can be useful for small local lattices and ill-conditioned operators.

###  3.4 HMC

Configuration generation with the HMC is supported for GPUs. The following
algorithmic checks are done following and checking consistency with [5] using
an SU(2) gauge group with two fermions in the adjoint representation, using
two RHMC monomials at $\beta=2.0$ and $m=-1.34$, both on an NVIDIA V100 GPU
(via UCloud DeiC Interactive HPC) and, in addition, for the HIP-ported version
tested on AMD MI250x GPUs (LUMI-G).

Figure 8: Average plaquette of ensembles Figure 9: Check that the average
plaquette value does not depend on the step size in the integrator (left).
Check that the simulations obey the Creutz equality (right).

The determined value of the average plaquette is independent of the step size
in the molecular dynamics integration. We further check the Creutz equality,
see [4], stating that

$\langle\mathrm{exp}(-\Delta H)\rangle=1$ (1)

While we do see, as expected, larger uncertainties for larger step sizes, the
ensembles show results consistent with this equality.

Figure 10: Distribution of Hamiltonian violations for selected step sizes

We further examine the distribution of Hamiltonian violations of the HMC
trajectories. For sufficiently small step sizes, the distribution is normal,
as expected, which is demonstrated both by the histograms on the right as well
as the quantile-quantile plots on the left, which show the values of
Hamiltonian violation on the y-axis and the expected quantile they should be
located in if they were Gaussian on the x-axis. A perfectly normal distributed
data vector forms a straight line in these plots. The slope of this line is
the standard deviation of the data.

Figure 11: Check analytic behavior of acceptance rate (bottom) and scaling of
the integrator with step size (top).

Lastly, we check that the Hamiltonian violations of the 2nd-order Omelyan
integrator [16] are increasing as

$\Delta H\sim\Delta\tau^{4},$ (2)

see [18], by fitting. Additionally, we check the asymptotic analytical
relation from [10]

$P_{\mathrm{acc}}\cong\mathrm{erfc}\left(\sqrt{\dfrac{\langle\Delta
H\rangle}{2}}\right)$ (3)

between the acceptance rate and Hamiltonian violations. We find that the
relations are fulfilled and that the algorithm is working.

##  4 Conclusion and outlook

The ported code is efficient and scales well. Substantial improvements in
overall software quality have been made. While configuration generation with
the HMC is efficient, performance improvements can still be achieved using
domain decomposition methods, as they minimize the need for communications.

##  5 Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement №813942. Testing, development, and benchmarking of this software was
possible using resources on LUMI-G provided by the Danish eInfrastructure
Consortium under grant application number DeiC-SDU-N5-2024055 and NVIDIA V100,
A100, and H100 nodes provided by the UCloud DeiC Interactive HPC system
managed by the eScience Center at the University of Southern Denmark. We are
grateful to Jacob Finkenrath and Martin Hansen for valuable discussions.

## References

  * [1] The ninja build system.  https://ninja-build.org/, 2013–2024.
  * [2] R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, and C. Rebbi.  Adaptive multigrid algorithm for the lattice Wilson-Dirac operator.  Phys. Rev. Lett., 105:201602, 2010.
  * [3] M. A. Clark and A. D. Kennedy.  The RHMC algorithm for two flavors of dynamical staggered fermions.  Nucl. Phys. B Proc. Suppl., 129:850–852, 2004.
  * [4] Michael Creutz.  Global monte carlo algorithms for many-fermion systems.  Phys. Rev. D, 38:1228–1238, Aug 1988.
  * [5] Luigi Del Debbio, Agostino Patella, and Claudio Pica.  Higher representations on the lattice: Numerical simulations. SU(2) with adjoint fermions.  Phys. Rev. D, 81:094503, 2010.
  * [6] Claudio Pica et. al.  Documentation of hirep.  https://claudiopica.github.io/HiRep/.
  * [7] Martin Lüscher et. al.  Openqcd.  https://luscher.web.cern.ch/luscher/openQCD/.
  * [8] Anthony Francis, Patrick Fritzsch, Martin Lüscher, and Antonio Rago.  Master-field simulations of O($a$)-improved lattice QCD: Algorithms, stability and exactness.  Comput. Phys. Commun., 255:107355, 2020.
  * [9] Andreas Frommer, Karsten Kahl, Stefan Krieg, Björn Leder, and Matthias Rottmann.  Adaptive Aggregation-Based Domain Decomposition Multigrid for the Lattice Wilson–Dirac Operator.  SIAM J. Sci. Comput., 36(4):A1581–A1608, 2014.
  * [10] Sourendu Gupta, A. Irback, F. Karsch, and B. Petersson.  The Acceptance Probability in the Hybrid Monte Carlo Method.  Phys. Lett. B, 242:437–443, 1990.
  * [11] Martin H. Gutknecht.  Variants of bicgstab for matrices with complex spectrum.  SIAM Journal on Scientific Computing, 14(5):1020–1033, 1993.
  * [12] Martin Hasenbusch.  Speeding up the hybrid Monte Carlo algorithm for dynamical fermions.  Phys. Lett. B, 519:177–182, 2001.
  * [13] A. D. Kennedy, Ivan Horvath, and Stefan Sint.  A New exact method for dynamical fermion computations with nonlocal actions.  Nucl. Phys. B Proc. Suppl., 73:834–836, 1999.
  * [14] M. Luscher and P. Weisz.  On-shell improved lattice gauge theories.  Commun. Math. Phys., 98(3):433, 1985.  [Erratum: Commun.Math.Phys. 98, 433 (1985)].
  * [15] Martin Luscher.  Solution of the Dirac equation in lattice QCD using a domain decomposition method.  Comput. Phys. Commun., 156:209–220, 2004.
  * [16] I.P. Omelyan, I.M. Mryglod, and R. Folk.  Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations.  Computer Physics Communications, 151(3):272–314, 2003.
  * [17] B. Sheikholeslami and R. Wohlert.  Improved Continuum Limit Lattice Action for QCD with Wilson Fermions.  Nucl. Phys. B, 259:572, 1985.
  * [18] Tetsuya Takaishi.  Choice of integrator in the hybrid monte carlo algorithm.  Computer Physics Communications, 133(1):6–17, 2000.
  * [19] H. A. van der Vorst.  Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 13(2):631–644, 1992.