#  $K\to\pi\pi$ matrix elements
from mixed action lattice QCD

Jack Laiho
Department of Physics and Astronomy, University of Glasgow, Glasgow, Scotland,
UK
E-mail<EMAIL_ADDRESS>Ruth S. Van de Water
Physics Department, Brookhaven National Laboratory, Upton, New York, USA
E-mail<EMAIL_ADDRESS>
###### Abstract:

We present a new method for determining $K\to\pi\pi$ matrix elements from
lattice simulations that is less costly than direct simulations of
$K\to\pi\pi$ at physical kinematics. It improves, however, upon the
traditional “indirect” approach of constructing the $K\to\pi\pi$ matrix
elements using NLO $SU(3)$ $\chi$PT, which can lead to large higher-order
chiral corrections. Using the explicit example of the $\Delta I=3/2$ $(27,1)$
operator to illustrate the method, we obtain a value for Re($A_{2}$) that
agrees with experiment and has a total uncertainty of $\mathrel{\hbox
to0.0pt{\lower 3.0pt\hbox{$\mathchar 536\relax$}\hss}\raise
2.0pt\hbox{$\mathchar 316\relax$}}$ 20%. Although our simulations use domain-
wall valence quarks on the MILC asqtad-improved gauge configurations, this
method is more general and can be applied to calculations with any fermion
formulation.

##  1 Motivation

Lattice calculations of $K\to\pi\pi$ matrix elements are important for
understanding the Standard Model and in constraining physics beyond the
Standard Model. For example, they are needed to explain the origin of the
$\Delta I=1/2$ rule and to compute the long-distance contributions to neutral
kaon mixing [1]. Because the lowest-order Standard Model contributions to
$\epsilon^{\prime}/\epsilon$ are from 1-loop electroweak penguin diagrams,
$K\to\pi\pi$ decay is sensitive to physics at very high scales. Many
extensions of the Standard Model lead to new particles that enter the loops,
and these contributions to $K\to\pi\pi$ may be sufficiently large that they
can be observed once the hadronic uncertainties in the weak matrix elements
are small enough.

A standard way of searching for new physics in the flavor sector is by
overconstraining the angles and sides of the CKM unitarity triangle [2]. This
requires precise experimental measurements and equally well-controlled
theoretical calculations of hadronic weak matrix elements using lattice QCD.
For many years improved measurements and calculations have simply confirmed
the Standard Model CKM framework at the few-percent level, but recent
$N_{f}=2+1$ flavor lattice calculations of $B_{K}$ with $\sim$ 4% precision
[3] have revealed a 2-3$\sigma$ tension in the CKM unitarity triangle [1, 4,
5]. This tension, which may be due to kaon or $B_{d}$-meson mixing, is
illustrated in Fig. 1. Almost all constraints on the CKM unitarity triangle,
however, come from the $B$-meson sector. Thus it is essential to place other
constraints on the unitarity triangle from the kaon sector in order to test
whether the amount of observed $CP$-violation in the $B$-meson sector is the
same as in the kaon sector. Once lattice QCD calculations of $K\to\pi\pi$
matrix elements are sufficiently precise, they can be combined with the
experimental measurement of $\epsilon_{K}^{\prime}/\epsilon_{K}$ to impose an
additional constraint on the apex of the CKM unitarity triangle.

Figure 1: Global fit of the CKM unitarity triangle [5].

##  2 New approach to $K\to\pi\pi$ matrix elements

The Maiani-Testa no-go theorem states that physical matrix elements cannot be
extracted from Euclidean correlation functions with multi-hadron states [6].
Therefore two general approaches have been developed to evade the Maiani-Testa
theorem and allow the determination of $K\to\pi\pi$ matrix elements using
lattice QCD. The “direct” Lellouch-Lüscher finite-volume method [7] is the
most straightforward to implement, but it is computationally demanding because
it requires a large ($\sim$ 6 fm) box and physical light-quark masses. The
“indirect” method constructs $K\to\pi\pi$ matrix elements using the low-energy
constants (LEC’s) of $SU(3)$ $\chi$PT obtained from calculating simpler
lattice quantities such as $K\to 0$ and $K\to\pi$. Although it was shown that
all LEC’s through next-to-leading order can be obtained from such “simple”
lattice quantities [8], this approach relies on the use of $SU(3)$ $\chi$PT at
the kaon mass, where the convergence of the chiral expansion is quite slow.

Li and Christ studied the extraction of $K\to\pi\pi$ matrix elements with
$N_{f}=2+1$ dynamical domain-wall lattice simulations using the “indirect”
method in Ref. [9]. They concluded that large uncertainties in the LO and NLO
$SU(3)$ LEC’s and the slow convergence of $SU(3)$ $\chi$PT at the scale of the
kaon mass lead to large errors that make the extraction of $K\to\pi\pi$ matrix
elements using the “indirect” method unreliable. Our procedure therefore
addresses these drawbacks of the traditional approach and improves upon it in
several ways. In the combined chiral-continuum extrapolation, we use the
physical pseudoscalar meson masses and decay constant. This leads to better
fits as measured by the correlated $\chi^{2}$/d.o.f.; nontrivial agreement
between the NLO mixed-action $\chi$PT prediction for the isovector scalar
correlator and lattice simulation data lends support to this approach [10].
When the fixed-order (NLO) fit is bad, we approximate higher order terms in
the chiral expansion by polynomials. This leads to larger leading-order terms
and hence suggests better convergence than was found in Ref. [9]. For example,
Fig. 2 shows the $SU(3)$ $\chi$PT fit of $B_{K}$ along with the sizes of the
various contributions [11]. The NLO corrections are only approximately 1/3 of
the LO terms even at $m_{s}/2$.

Figure 2: Chiral and continuum extrapolation of $B_{K}$ [11] (left plot) along
with convergence of the $SU(3)$ $\chi$PT fit (right plot). Circles (squares)
denote $a\sim 0.12$ fm ($a\sim 0.09$ fm) data. Only degenerate points are
shown, but the fit also includes non-degenerate data. The cyan band is the
degenerate quark mass full QCD curve ($m_{x}=m_{y}=m_{l}=m_{h}$) in the
continuum limit. The y-intercept of the band gives the LEC $B_{0}$, the value
of $B_{K}$ in the $SU(3)$ chiral limit. The right-most point on both plots
corresponds to $\sim m_{s}/2$.

Despite these findings, however, NLO $\chi$PT corrections can still be 50% or
more for some quantities. Therefore, to achieve the precision needed for
$K\to\pi\pi$ we do not rely on the “indirect” method alone. Rather, we combine
indirect and direct methods in a cost-effective way. We bypass the Maiani-
Testa theorem by simulating with both pions at rest. We fit the numerical data
to NLO mixed-action $\chi$PT plus higher-order analytic terms, extrapolate to
the continuum, and interpolate to the point at which
$m_{K}=m_{K}^{\textrm{phys.}}$ and $m_{\pi}=m_{K}^{\textrm{phys.}}/2$. Thus we
avoid relying upon $SU(3)$ $\chi$PT to extrapolate to the physical kaon mass,
where we expect higher-order corrections to be significant. We then correct
this unphysical kinematics point using fixed-order $SU(3)$ $\chi$PT. The low
energy constants needed for this correction can be obtained from simpler
quantities such as $f_{K}$, $K^{0}$-$\bar{K^{0}}$, and $K\to\pi$. Since the
kaon is tuned to its physical value, terms involving only kaons are correct to
all orders in the $SU(3)$ chiral expansion; we therefore expect higher-order
corrections to be small.

We can test this approach using the known quantities $f_{K}$ and $f_{\pi}$;
the results are shown in Fig. 3. The size of the NLO corrections are quite
small (below 10% for $f_{\pi}$ and below 5% for $f_{K}$), indicating that the
systematic uncertainty due to truncating the chiral expansion is under
control.

Figure 3: Demonstration of the method for $f_{\pi}$ (left plot) and $f_{K}$
(right plot). Errors on the circular points only show the statistical errors
in the numerator. Vertical error bands denote the total (statistical plus
systematic) uncertainties in $f_{\pi}$ and $f_{K}$. After the interpolation to
the unphysical kinematics point $m_{K}=m_{K}^{\textrm{phys.}}$ and
$m_{\pi}=m_{K}^{\textrm{phys.}}/2$, the size of the NLO $\chi$PT correction is
below 10% for $f_{\pi}$ and below 5% for $f_{K}$.

##  3 Preliminary determination of Re($A_{2}$)

We now use our approach to determine the $(27,1)$ $\Delta I=3/2$ $K\to\pi\pi$
matrix element, which can be combined with continuum Wilson coefficients [12]
to obtain Re($A_{2}$).

We compute the matrix element in unquenched lattice QCD using asqtad-improved
staggered sea quarks and domain-wall valence quarks. This mixed-action
approach shares the primary advantages of both staggered and domain-wall
lattice simulations. We use the publicly-available 2+1 flavor MILC gauge
configurations [13], and simulate with several valence and sea quark masses.
Although our preliminary analysis only uses two lattice spacings, we have
generated data at three lattice spacings ($a\sim 0.06$, 0.09 and 0.12 fm) and
will include all of it in a future publication. Our lightest taste-
pseudoscalar sea-sea pion has $m_{\pi,5}=240$ MeV while our lightest valence-
valence pion has $m_{\pi}=210$ MeV. On the $a\sim 0.06$ fm ensembles, the
heaviest (taste-singlet) sea-sea pion is also quite light $m_{\pi,I}=270$ MeV.
This gives us good control over our combined chiral-continuum extrapolation
using mixed-action $\chi$PT. The approximate chiral symmetry of the valence
sector ($m_{\textrm{res}}<3$ MeV on all lattice spacings) makes analysis of
the mixed-action simulation data simpler than the purely staggered case. Only
two additional parameters appear at 1-loop in the mixed action $\chi$PT
expressions for $m_{PS}$, $f_{PS}$, and $B_{K}$ as compared to the purely
domain-wall case [14], and they can both be obtained from spectrum
calculations. Furthermore, nonperturbative renormalization using the method of
Rome-Southampton [15] can be carried out in a straightforward manner. Finally,
the success of our earlier mixed-action lattice calculation of $B_{K}$ [11]
indicates that the mixed-action method is also a good way to determine
$K\to\pi\pi$ matrix elements.

Table 1: Data used for the preliminary determination of Re($A_{2}$). The
columns show the (i) approximate lattice spacings, (ii) lattice volumes,
(iii), nominal up/down ($m_{l}$) and strange quark ($m_{h}$) masses in the
sea, (iv) corresponding pseudoscalar taste pion mass, (v) partially quenched
valence quark masses ($m_{x}$), (vi) lightest available domain-wall pion mass,
and (vii) number of configurations analyzed on each ensemble.
| sea sector | valence sector |
---|---|---|---
$a$(fm) | $\left(\frac{L}{a}\right)^{3}\times\frac{T}{a}$ |  $am_{l}/am_{h}$ |  $am_{\pi,5}$ |  $am_{x}$ |  $am_{\pi}$ | $N_{\rm conf.}$
0.06 | $64^{3}\times 144$ | 0.0018/0.018 | 0.06678(3) | 0.0026, 0.0108, 0.033 | 0.06376(96) | 96
0.06 | $48^{3}\times 144$ | 0.0036/0.018 | 0.09353(7) | 0.0036, 0.0072, 0.0108, 0.033 | 0.07458(76) | 129
0.12 | $24^{3}\times 64$ | 0.005/0.05 | 0.15970(13) |  0.007, 0.02, 0.03, 0.05, 0.065 | 0.1718(11) | 218
0.12 | $20^{3}\times 64$ | 0.007/0.05 | 0.18887(8) |  0.01, 0.02, 0.03, 0.04, 0.05, 0.065 | 0.1968(08) | 279

Figure 4 shows the interpolation to the unphysical kinematics point
$m_{K}=m_{K}^{\textrm{phys.}}$ and $m_{\pi}=m_{K}^{\textrm{phys.}}/2$. Before
the interpolation, we adjust the data points by the known 1-loop finite volume
corrections, which only depend upon the valence quark masses [16]. The
interpolation currently uses LO $\chi$PT supplemented by NLO analytic terms,
including a term proportional to $a^{2}$ so that we can take the continuum
limit. We have finished calculating the mixed-action 1-loop chiral logs,
however, and are now working on incorporating them into the fit. We then
correct the unphysical kinematics “2$\pi$” point to physical kinematics using
$SU(3)$ $\chi$PT:

$\displaystyle\langle\pi^{+}\pi^{-}|{\mathcal{O}}_{i}|K^{0}\rangle_{\textrm{phys.}}$
$\displaystyle=$
$\displaystyle\langle\pi^{+}\pi^{-}|{\mathcal{O}}_{i}|K^{0}\rangle_{\textrm{$2\pi$}}\times(1+\delta_{\textrm{$\chi$PT}})\,,$
(1)

where the correction factor is given by

$\displaystyle\delta_{\textrm{$\chi$PT}}$ $\displaystyle=$
$\displaystyle\left(\langle\pi^{+}\pi^{-}|{\mathcal{O}}_{i}|K^{0}\rangle_{\textrm{phys.}}-\langle\pi^{+}\pi^{-}|{\mathcal{O}}_{i}|K^{0}\rangle_{\textrm{$2\pi$}}\right)/\langle\pi^{+}\pi^{-}|{\mathcal{O}}_{i}|K^{0}\rangle_{\textrm{$2\pi$}}\,.$
(2)

Because we have already interpolated to the physical kaon mass, terms in the
$\chi$PT expression that are only functions of the kaon mass cancel in the
numerator, and the correction is needed only for the short extrapolation from
$m_{K}/2$ to $m_{\pi}$. We therefore expect the convergence properties to be
better and the truncation errors to be smaller than if we were to extrapolate
up to the kaon mass.

At leading order, the expression for the $\Delta I=3/2$ $(27,1)$ $K\to\pi\pi$
matrix element is

$\displaystyle\langle\pi^{+}\pi^{-}|{\mathcal{O}}^{\Delta
I=3/2}_{(27,1)}|K^{0}\rangle_{\textrm{LO}}$ $\displaystyle=$ $\displaystyle
4iB_{0}f_{0}(m^{2}_{K}-m^{2}_{\pi})/3\,,$ (3)

where $f_{0}$ and $B_{0}$ are the pion decay constant and $B_{K}$ in the
chiral limits, respectively. The $\chi$PT correction factor is then

$\displaystyle\delta_{\textrm{$\chi$PT}}^{\textrm{LO}}=\left[(m_{K}/2)^{2}-m_{\pi}^{2}\right]/\left[m_{K}^{2}-(m_{K}/2)^{2}\right]$
(4)

and is only 23%, which is much smaller than the chiral correction factors
observed by Li and Christ using the standard direct approach [9]. Because the
low-energy constants $f_{0}$ and $B_{0}$ cancel in the ratio, there is no
ambiguity regarding the choice of $SU(3)$ LEC’s.

Figure 4: Left plot: interpolation of Re($A_{2}$) to the unphysical kinematics
point $m_{K}=m_{K}^{\textrm{phys.}}$ and $m_{\pi}=m_{K}^{\textrm{phys.}}/2$.
Circles (squares) denote $a\sim 0.12$ fm ($a\sim 0.06$ fm) data. Fit lines
correspond to the degenerate mass case and should pass through the filled
symbols. Right plot: correction of Re($A_{2}$) to physical kinematics using LO
$SU(3)$ $\chi$PT.

The renormalization factor for the $(27,1)$ $\Delta I=3/2$ operator is the
same as $B_{K}$, so we can use the result for $Z_{B_{K}}$ from Ref. [11]. We
obtain

$\displaystyle\textrm{Re}(A_{2})=1.568(86)\times 10^{-8},$ (5)

where the error is statistical only. This agrees with the experimental
measurement, $\textrm{Re}(A_{2})_{\textrm{exp}}=1.50\times
10^{-8}\textrm{GeV}$ [17]. Table 2 presents an estimated error budget for
$\textrm{Re}(A_{2})$. We assume that the truncation error due to leaving out
NLO corrections is half the size of the LO correction, and estimate that total
uncertainty in our preliminary result is below 20%. This should improve
further with the use of our full data set. The $\chi$PT truncation error and
error from the uncertainty in the LEC’s will likely also decrease with the use
of the NLO $\chi$PT correction factor. We restrict our lightest valence quark
mass to maintain $m_{\pi}L\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{$\mathchar
536\relax$}\hss}\raise 2.0pt\hbox{$\mathchar 318\relax$}}3.5$, and estimate
that our finite volume errors are a few percent using 1-loop FV$\chi$PT. We
will perform an explicit finite-volume study, however, before publication. Our
preliminary result is renormalized using lattice perturbation theory because
we have not yet completed the nonperturbative renormalization on the $a\sim
0.06$ fm ensembles. We expect these $Z$-factors to be reliable within
systematic uncertainties, however, because we find good agreement between
$Z_{B_{K}}$computed nonperturbatively and using lattice perturbation theory on
the $a\sim 0.12$ and $a\sim 0.09$ ensembles.

Table 2: Estimated total error budget for $\textrm{Re}(A_{2})$. Each source of uncertainty is given as a percentage. uncertainty | $\qquad\textrm{Re}(A_{2})$
---|---
statistics | 4.7%
$\chi$PT truncation error | 12%
uncertainty in leading-order LECs | 4%
discretization errors | 4%
finite volume errors | few percent
renormalization factor | 3.4%
scale and quark-mass uncertainties | 3%
Wilson coefficients | few percent
total | less than 20%

## References

  * [1] A. J. Buras and D. Guadagnoli, Phys. Rev. D 78, 033005 (2008).
  * [2] M. Antonelli, et al., Phys. Rept. 494, 197-414 (2010).
  * [3] See V. Lubicz, PoS LAT2009, 013 (2009) and Refs. therein.
  * [4] E. Lunghi and A. Soni, Phys. Lett. B 666, 162 (2008).
  * [5] J. Laiho, E. Lunghi and R. S. Van de Water, Phys. Rev. D 81, 034503 (2010); updates at www.latticeaverages.org.
  * [6] L. Maiani and M. Testa, Phys. Lett. B 245, 585 (1990).
  * [7] L. Lellouch and M. Lüscher, Commun. Math. Phys. 219, 31 (2001).
  * [8] J. Laiho and A. Soni, Phys. Rev. D 71, 014021 (2005).
  * [9] S. Li and N. H. Christ, PoS LATTICE2008, 272 (2008).
  * [10] C. Aubin, J. Laiho and R. S. Van de Water, Phys. Rev. D 77, 114501 (2008).
  * [11] C. Aubin, J. Laiho and R. S. Van de Water, Phys. Rev. D 81, 014507 (2010).
  * [12] M. Ciuchini, E. Franco, G. Martinelli, L. Reina and L. Silvestrini, Z. Phys. C 68, 239 (1995).
  * [13] C. Aubin et al. [MILC], Phys. Rev. D 70, 114501 (2004).
  * [14] O. Bar, C. Bernard, G. Rupak and N. Shoresh, Phys. Rev. D 72, 054502 (2005);
C. Aubin, J. Laiho and R. S. Van de Water, Phys. Rev. D 75, 034502 (2007).

  * [15] G. Martinell, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, Nucl. Phys. B 445, 81 (1995).
  * [16] C. J. D. Lin, G. Martinelli, C. T. Sachrajda et al., Nucl. Phys. B619, 467-498 (2001).
  * [17] A. J. Buras, arXiv:hep-ph/9806471.


#  Exploring Quantization and Mapping Synergy in Hardware-Aware Deep Neural
Network Accelerators

Jan Klhufek, Miroslav Safar, Vojtech Mrazek, Zdenek Vasicek, Lukas Sekanina
Brno University of Technology, Faculty of Information Technology, Brno, Czech
Republic
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>{mrazek, vasicek,
<EMAIL_ADDRESS>
###### Abstract

Energy efficiency and memory footprint of a convolutional neural network (CNN)
implemented on a CNN inference accelerator depend on many factors, including a
weight quantization strategy (i.e., data types and bit-widths) and mapping
(i.e., placement and scheduling of DNN elementary operations on hardware units
of the accelerator). We show that enabling rich mixed quantization schemes
during the implementation can open a previously hidden space of mappings that
utilize the hardware resources more effectively. CNNs utilizing quantized
weights and activations and suitable mappings can significantly improve trade-
offs among the accuracy, energy, and memory requirements compared to less
carefully optimized CNN implementations. To find, analyze, and exploit these
mappings, we: (i) extend a general-purpose state-of-the-art mapping tool
(Timeloop) to support mixed quantization, which is not currently available;
(ii) propose an efficient multi-objective optimization algorithm to find the
most suitable bit-widths and mapping for each DNN layer executed on the
accelerator; and (iii) conduct a detailed experimental evaluation to validate
the proposed method. On two CNNs (MobileNetV1 and MobileNetV2) and two
accelerators (Eyeriss and Simba) we show that for a given quality metric (such
as the accuracy on ImageNet), energy savings are up to 37% without any
accuracy drop.

††publicationid: pubid: 979-8-3503-5934-3/24/$31.00 ©2024 IEEE

##  I Introduction

Highly optimized hardware accelerators and memory subsystems have been
developed to meet the challenging energy, performance, and memory requirements
imposed on the fast and efficient processing of _convolutional neural
networks_ (CNNs) [1]. While CNNs are typically trained on GPUs that utilize
floating point (FP) data representation, the most energy-efficient CNN
hardware (inference) accelerators employ spatial architecture and advanced
mixed-precision quantization schemes operating with a few bits [2]. We will
primarily deal with _weight quantization_ , which can significantly reduce
memory size and transfers in addition to performance and energy efficiency
gains. Furthermore, we will also quantize the activations, which mainly affect
the smaller internal scratchpad memory. As quantization almost always
negatively impacts the CNN’s error, various automated quantization methods
enabling suitable trade-offs between the error and hardware parameters (such
as inference energy, latency, and weight memory size) have been proposed [3].

In this paper, we focus on one aspect of CNN quantization not addressed by the
literature: performance gains that can be obtained when the synergy of a
suitable _quantization_ and _mapping_ (i.e., placement and scheduling of CNN
elementary operations on hardware units) is exploited in CNN hardware
accelerators. Consider a mixed-precision quantization method guided by naïve
memory footprint estimation, e.g., a layer-wise optimization of weight bit-
widths minimizing the total number of bits in weight memory, i.e., the model
size. Fig. 1(a) shows that when the naïve estimation is applied in the
MobileNetV1 CNN, the correlation with the memory word count (after bit-
packing) is far from perfect (calculated for 1000 randomly generated quantized
versions of MobileNetV1). Moreover, only a weak correlation is visible for the
Energy-Delay-Product (EDP) evaluated for one inference on the Eyeriss
accelerator (Fig. 1(b)). These weak correlations exist because the
accelerator’s architecture, mapping, and memory subsystem, which significantly
influence hardware parameters, are not considered during quantization.

Figure 1: Correlation between Model Size (the number of bits for weights) and:
(a) Memory Word Count after bit-packing (b) and evaluated Energy-Delay Product
on Eyeriss for $1000$ unique configurations of MobileNetV1. The blue point
represents the original uniform setting.

The use of mixed-precision layer-wise quantization schemes is motivated by the
fact that different layers have different error resilience characteristics and
behave differently on the hardware [2]. This observation can be exploited by
various optimization and approximation strategies. However, the mixed-
precision quantization is not supported by state-of-the-art configurable
general-purpose models of hardware accelerators utilized to estimate the
hardware parameters of a given CNN implementation before its deployment. It
turns out that the automated quantization methods only consider either the CNN
model (not the hardware) [4] or exist for a particular hardware accelerator
[3], i.e., they are not usable across a reasonable set of CNN accelerators.
This is partly because determining hardware parameters using real hardware
accelerators is time-consuming.

To easily model, analyze, and investigate the impact of various mixed-
precision quantization schemes (particularly those based on computing the most
suitable bit-width for weights and activations for each layer) on the hardware
parameters of CNN hardware accelerators, we extended Timeloop to support
different mixed-precision quantization schemes. Timeloop [5] is a popular tool
searching for the most suitable mapping of the CNN model on the accelerator
resources, in which a rich set of accelerator architectures, configurations,
and mapping algorithms can be utilized. When coupled with Accelergy [6], an
early-stage energy and execution time estimation tool, detailed hardware
characteristics can be obtained quickly. Then, we show that enabling rich
mixed quantization schemes during the implementation of a CNN can open a
previously hidden space of mappings that utilize the hardware resources more
effectively than uniformly quantized layers accompanied by standard mappings.
CNNs utilizing quantized weights and activations (both input and output) and
suitable mappings can significantly improve trade-offs among the accuracy,
energy, and memory requirements compared to less carefully optimized CNN
implementations. We propose to determine the most suitable mixed-precision
quantization schemes with a multi-objective genetic algorithm NSGA-II [7]
(with extended Timeloop in the loop) that simultaneously minimizes the weight
memory size (reflecting the accelerator’s memory subsystems), inference
energy, and CNN error. The optimal mapping is obtained using Timeloop. To
stress the importance of mixed-precision quantization, we intentionally keep
the accelerator specification unchanged for all experiments.

In this paper, we make the following key contributions.

  1. 1.

To investigate the impact of mixed-precision quantization schemes on CNN
accelerator’s performance and energy consumption, we extend Timeloop to
support mixed-precision quantization. We show that adopting the mixed
precision quantization can increase the number of valid (CNN model-to-
hardware) mappings by one order of magnitude, enabling thus many optimization
strategies unreachable with a common approach. The extended version of
Timeloop is available as open-source code at: https://github.com/ehw-
fit/timeloop-with-quantization

  2. 2.

To analyze the impact of a particular CNN quantization on the CNN error, we
connect extended Timeloop with the PyTorch framework. The proposed extension
can be alternatively connected with the TensorFlow framework.

  3. 3.

To find a Pareto front containing CNNs that show high-quality tradeoffs
between the error and hardware parameters, we propose to adopt NSGA-II with
Timeloop in the loop. Timeloop coupled with Accelergy is much faster than
accelerator simulation or in-hardware measurement so that many candidate
solutions can be quickly generated and evaluated. The method is parallelized
and executed much faster than other approaches such as [3] and allows the
exploration of various accelerator architectures.

  4. 4.

On two CNNs (MobileNetV1 and MobileNetV2) and two CNN accelerators (Eyeriss
and Simba), we show that for a given quality metric (such as classification
accuracy on a subset of ImageNet), energy savings are up to 37% without any
accuracy drop.

##  II Related Work

_CNN hardware accelerators_ with spatial architecture employ an array of
locally communicating processing elements (PE). Each of them implements a
multiply-and-accumulate circuit, a small local memory (registers), and a
controller. PEs are usually organized as pipelined systolic arrays optimized
for fast execution of CNN operations. The term _mapping_ refers to the
dataflow strategy (i.e., computation order and parallelism strategy) coupled
with the tiling strategy (selecting the size of input data with respect to
available hardware resources such as buffers and PEs). A particular mapping of
a given CNN on hardware resources and the execution of this mapping are
implemented by the accelerator controller. The accelerator can use a common
dataflow strategy (e.g., input stationary, row stationary) [1] or highly
optimized dataflow determined using methods such as GAMMA (Genetic Algorithm-
based Mapper for ML Accelerators) [8]. In addition to optimizing the dataflow,
common approaches to reach energy-efficient processing and low memory
footprint are weight compression, CNN pruning, arithmetic operation
approximation, and quantization [2]. A survey of quantization methods for
efficient neural network inference was presented in [9]. The quantization can
be uniform, non-uniform, symmetric, asymmetric, static, and dynamic. It can be
applied layerwise, groupwise, channelwise, and sub-channelwise. Quantization-
aware training simulating quantization noise during F32 training by inserting
fake quantization to the model [10] and post-training quantization of the pre-
trained F32 model are two major quantization strategies. The mix-precision
quantization in which several quantization precisions are applied together,
can further improve the common uniform quantization. If the automated CNN
model design is performed together with the hardware optimization, we speak
about hardware-aware neural architecture search [11, 12]. A candidate CNN
implementation is evaluated in terms of hardware parameters (such as energy
and latency) using either simulation (which is usually slow), analytical
modeling (often not reliable), estimation tools (such as Timeloop [5],
Accelergy [6], BitFusion [13], MAESTRO [14], DNN chip predictor [15]),
surrogate ML models (fast but less reliable), or measurements on real hardware
[16].

Unfortunately, current tools developed for hardware parameter estimation of
CNN inference accelerators do not support mixed-precision quantization, which
is crucial for lowering energy consumption and memory utilization.

##  III Proposed framework

The proposed framework consists of three components (a) a mapping engine, (b)
a training engine, and (c) a search engine. The overall scheme is shown in
Fig. 2. The search engine is responsible for generating configurations based
on which a quantized version of the original CNN is obtained. The goal of the
search engine is to minimize CNN error and hardware resources. Given a
hardware accelerator, the mapping engine determines the best mapping for the
obtained quantized CNN. The mapper consumes a configuration file describing
the computational problem, parameters of the target HW accelerator, and
constraints. It outputs the best possible CNN execution plan and additional
information, such as the number of memory accesses, latency, and consumed
energy that are needed to run the search engine. The training engine is used
to perform retraining after quantization to recover from a quality loss caused
by introducing the quantized layers. To obtain the best configurations, the
search process is iterated for a predefined number of iterations. Both mapping
and training engines must support mixed-precision quantization.

Figure 2: Overall scheme of the proposed quantization and accelerator-mapping
optimization. Search engine (purple components), Training engine (blue),
Mapping engine (yellow).

###  III-A Mapping engine

The mapping engine is based on Timeloop [5] and Accelergy [6], which we have
extended to support mixed-precision mapping. Timeloop internally uses a mapper
that finds the best way to schedule operations and stage data on the specified
architecture. Various non-deterministic as well as deterministic heuristics
are implemented in the original software to explore the design space
effectively.

Our extension consists of the following modifications. First, the
specification is extended to permit the specification of bit-widths associated
with the data tensors of individual workloads. Second, the mapping algorithm
is modified to support so-called bit-packing [17]. This technique simply
allows as many data elements (i.e. weights or activations) to be stored in a
single memory word as the bit-width of the word allows. When properly
implemented, this technique not only allows for minimizing the amount of
memory used for storing weights and activations but also allows for maximizing
the memory bus utilization because multiple weights could be loaded in
parallel. This results in fewer memory accesses and lower energy consumption
for data transfer compared to a naïve approach, where one memory word is
occupied by one weight/activation. In addition to that, we modified a checker
which checks for mapping violations.

Timeloop is used to characterize the workload on a per-layer basis. That is,
scheduling is computed for each layer of the quantized network independently.
To fully exploit the impact of quantized activations on memory utilization,
the bit-width of both inputs and outputs are considered. The output bit-width
is determined by the bit-width of the following layer’s inputs, while constant
$8$ bits are set for the last layer’s outputs. The total energy is determined
as a sum of the energies required to compute every workload. The same is valid
also for total latency. Because we are using the mapping engine in a loop, we
have taken advantage of the per-layer approach, and we have implemented a
caching mechanism. Once a layer workload has been evaluated, the results are
stored in a cache. Subsequently, the cached results can be read and reused
when trying to find the best plan for the same workload, eliminating the need
for re-evaluation. This mechanism helps to accelerate substantially the design
space exploration because the candidate configurations typically contain many
similar parts.

###  III-B Training engine

The training engine is built on top of PyTorch framework and the principles of
Quantization Aware Training (QAT). We chose the PyTorch framework because of
its scalability, broad support on various computing platforms, and good
support of quantization.

The quantization is based on a per-tensor asymmetric approach, which is
natively supported in PyTorch. According to the tested configuration, the bit-
widths of the weights and the activations of the quantizable layers are set.
PyTorch can only handle 8-bit quantization; lower bit-widths are implemented
using specialized so-called observer modules that modify the allowed range of
values. Then, the fusing of layers to include batch normalization is applied,
and finally, fake-quantization observers are added. The modified model can
then be trained and validated.

When evaluating accuracy by running QAT, we typically start with a pre-trained
network at FP32 accuracy. We can quantize this network for a given bit-width
and perform training for $e$ epochs during the optimization loop. Since the
QAT is the most time-consuming part of the optimization and needs to be
performed as fast as possible, it is feasible to pre-quantize the input model
to some higher number of bits (e.g., 8) and only perform fine-tuning for a
smaller number of $e$ epochs in the loop.

###  III-C Search engine

In the mapping engine, only the memory path is considered in this study to
highlight the impact of mixed-precision weight quantization on energy and
memory usage. The computational MAC units remain untouched.

The accelerator configuration corresponding to the recipe for obtaining the
quantized network from the original network is modeled using a linear string
of tuples of integers. For MobileNetV1, for example, the string consists of 56
integers. Each tuple corresponds to a single layer and determines the bit-
width of the inputs and weights of the associated layer. The bit-width of the
outputs is determined by the bit-width of the inputs of the subsequent layer.

We use a multi-objective genetic algorithm (NSGA-II) [7] to obtain quantized
CNNs showing desired properties. The algorithm maintains a population of $|P|$
candidate configurations. The search starts from a population consisting of
configurations corresponding with uniformly quantized CNNs. The candidate
configurations are then iteratively optimized with respect to the accuracy of
quantized CNN and EDP required to perform one inference. For each candidate
solution, a corresponding PyTorch model and YAML configuration are created and
utilized in the training and mapping engine. Then, the energy and the accuracy
of the quantized CNN are evaluated.

The new population of candidate configurations is created from the current
parent population $P$ by the so-called uniform crossover followed by the
mutation operator. The number of offspring individuals is $|Q|$. The uniform
crossover randomly picks two configurations from $P$ and produces a single
offspring configuration where each integer is chosen with equal probability
from one of the two parents.

Then, with probability $p_{mutAcc}$ one randomly selected layer of the
offspring is set to the initial quantization 8/8, and with the probability
$p_{mut}$, one randomly selected integer is replaced by a randomly generated
but valid value. At the end of the optimization process, when a terminating
condition is met (typically, the maximum number of allowed iterations is
exceeded), the error of quantized CNNs is evaluated using the complete
training set. Solutions whose parameters are dominated by at least one other
solution are filtered out.

##  IV Experimental Setup

Two pre-trained CNN models are considered: MobileNetV1 and MobileNetV2. For
the evaluation, we classified images from the ImageNet dataset. To speed up
the evaluation, we have chosen a subset consisting of 100 classes. Both models
were trained using FP32 representation for 150 epochs (with top-1 accuracies
77.26% and 77.86% respectively) and then optionally quantized to 8 bits for 50
epochs. These quantized networks (FP32 and QAT-8 variants), as well as the
dataset, serve as input (see Original neural network and Dataset,
respectively, in Figure 2) to our tool. The proposed method is evaluated
independently for two HW accelerators supported by Timeloop/Accelergy: Eyeriss
and Simba. Eyeriss consists of 168 16-bit PEs, Simba employs 256 16-bit PEs.
The memory word size is 16. The characterization is done for 45nm technology.
The accelerators are provided to our tool in form of a text specification (see
Figure 2).

The number of epochs for QAT is $e=10$ for FP32 model and $e=5$ for QAT-8
model since it has already been adapted to the quantization. For the QAT-8
model we also considered $e=\\{10,20\\}$ number of epochs. The experiments
were executed on HPC nodes equipped with two AMD 7763 64-core CPUs and 8
NVIDIA A100 GPUs with 40 GB HBM2 RAM. The search is terminated after 48 hours.

The NSGA-II setup was determined after several test runs: the weights and
activations in each layer of CNN are allowed to be quantized to 2 – 8 bits;
$|P|=32$; $|Q|$ is chosen from $\\{8,16,32\\}$; $p_{mut}=10\%$;
$p_{mutAcc}=5\%$. It means that up to 32 Pareto-optimal solutions can be
generated at the end. Timeloop mapper is configured to use random search with
termination condition set to finding $2000$ valid mappings per workload.

##  V Results

Figure 3: Final Pareto fronts for (a) different Original neural network used
for QAT fine-tuning, (b) different settings for the number of offsprings
$|Q|$, and (c) different number of epochs $e$ during QAT. TABLE I: The total
number of exhaustively evaluated valid mappings for the second convolutional
layer of the MobileNet models on two HW accelerators for various quantization
settings. EDP of the least energy demanding mapping is given in
J$\cdot$Cycles.

Network | Bit-width Settings | Eyeriss | Simba
---|---|---|---
Layer | $q_{a},q_{w},q_{o}$ | Mappings | min. EDP | Mappings | min. EDP
MobileNet conv. layer #2 |  $16$, $16$, $16$ | $11,778$ | $638$ | $80,835$ | $204$
$8$, $8$, $8$ | $15,021$ | $388$ | $110,032$ | $159$
$8$, $4$, $8$ | $15,054$ | $372$ | $111,090$ | $159$
$8$, $2$, $8$ | $15,054$ | $364$ | $112,113$ | $158$
$4$, $4$, $4$ | $16,417$ | $281$ | $127,214$ | $137$
$2$, $2$, $2$ | $16,877$ | $228$ | $133,568$ | $126$

###  V-A Quantization-aware mapping

The purpose of the first experiment is to quantify the impact of the
modifications introduced in Timeloop on a) the size of the space of possible
valid mappings and b) energy savings. Table I shows the number of exhaustively
enumerated valid mappings for various quantization settings for one
convolutional layer. To conduct this analysis, we selected the second
convolutional layer (a depthwise convolutional layer) present in both analyzed
variants of MobileNet. Symbols $q_{a},q_{o}$ denote the activation and partial
sum bit-width; $q_{w}$ is the weight bit-width.

In the default setting, both accelerators employ 16-bit operands. For this
setting, the search space consists of 11,778 and 80,835 valid mappings for
Eyeriss and Simba. are considered, the number of valid mappings increases to
15,021 and 110,032, respectively. Keeping $q_{a}$ and $q_{o}$ at 8 bits and
further reducing $q_{w}$ will slightly increase the number of mappings. This
effect is visible, especially for the Simba accelerator, mainly due to the
fact that Eyeriss employs the row stationary dataflow. On the other hand,
decreasing both the activations and the weights bit-width results in a more
noticeable increase in the number of valid mappings. This is because the
execution of some mappings is bottle-necked by memories associated with
storing the activations. The results show that the proposed modification of
Timeloop allows for significantly extending the amount of valid mappings.

Figure 4 shows energy breakdowns for the whole MobileNetV1 and various uniform
quantization settings. The $x$b denotes the setting for which each layer was
independently calculated when $q_{a}=q_{w}=q_{o}=x$. The best mapping for each
setting was obtained by random search. As evident, the energy tied to the
memory transfers decreases with decreasing $x$, as it impacts the utilization
of memory subsystems that store the respective data. For example, in the case
of 4-bit quantization, when compared with the 8-bit setting, the total energy
is reduced by more than $32.5\%$, while the energy consumed in memory by
$54.5\%$. Moreover, note that for $x\geq 6$, the bit-packing yields no benefit
for the 16-bit word size. We can conclude that the proposed Timeloop extension
opens new mappings and allows us to exploit them to find more energy-efficient
memory utilization.

Figure 4: Energy breakdown of various versions of quantized MobileNetV1 on
Eyeriss accelerator.

###  V-B Automated quantization

The goal of NSGA-II is to find the most suitable trade-offs between CNN
accuracy and EDP for one inference. We selected Eyeriss and MobileNetV1
performing ImageNet classification to demonstrate the behavior of the proposed
automated quantization algorithm.

On one typical run ($e=10,|Q|=16$), plotted in Fig. 5, the optimization
progress is investigated. The Pareto front is improved significantly compared
to the initial uniform quantization. Most of the changes are done before the
11th generation.

Figure 5: Pareto fronts showing the best solutions found by NSGA-II across
various generations, applied to MobileNetV1.

As QAT is the most time-consuming part of the optimization process, we have to
determine the most suitable parameters of the proposed method to maximize its
benefits.

Firstly, we investigate the impact of a pre-trained model on the accuracy
achieved after QAT. For the FP32 model, we fine-tuned uniformly quantized
configurations for $e=10$, whereas for the QAT-8 model, we did the same for
$e=5$. The results in Fig. 3a show that better accuracies are obtained when
QAT-8 model is used, as it is already more accustomed to the effects of
quantization. Based on this outcome, we decided to use QAT-8 as the initial
model in all of the NSGA-II experiments.

Secondly, we ask what is the best offspring size if the number of evaluations
(i.e., the product $|Q|\times$ the number of generations) is constant. When
using the larger offspring population $|Q|$, it turns out that $|Q|=32$ is
unsuitable because only 8 generations are generated. There is no significant
difference in the resulting Pareto fronts for $|Q|=8$ or $32$, as shown in
Fig. 3b. Please note that the size of the parent population ($|P|$) has only a
minimal effect on the performance because these individuals are stored from
the previous generations, and their evaluation is already known.

Thirdly, we investigate the best trade-offs between the number of generations
and the number of epochs $e$ allowed for training candidate quantized CNNs.
Fig. 3c shows that better trade-offs between accuracy and EDP are reached when
$e$ is higher. However, the search is faster for lower $e$ because the number
of generations is 28 and 14 for $e$ = 10 and 20, respectively. We prefer a
larger number of epochs to achieve better accuracy-EDP trade-off.

###  V-C Final results

Figure 6: Tradeoff between accuracy and EDP evaluated on Eyeriss accelerator
running MobileNetV1 inference. The parameters are relative to the uniform
8-bit implementation.

For the final comparisons and based on the previous experiments, we used the
following settings of the proposed method: $e=20,|P|=32,|Q|=16$.

The proposed method utilizing mixed-precision quantization for the Eyeriss
accelerator is compared with the setting employing (a) the uniform
quantization (see Uniform in Fig. 6), (b) naïve quantization (see Naïve in
Fig. 6), and a different accelerator (see Proposed for Simba in Fig. 6). Fig.
6 reveals that if the optimization algorithm exploits the properties of the
target Eyeriss accelerator, it can reach better results than for an unseen
accelerator (represented by Simba in our case). A naïve approach considering
only memory aspects reached Pareto front incomparable with the proposed one.
Neither the uniform quantization is able to deliver better results than our
approach.

Table II shows the energy savings in the memory subsystem for two different
CNNs and two accelerators in the context of the accuracy drops on the ImageNet
subset. The uniform approach represents the SoA solutions that do not explore
the quantization of individual layers. It can be seen that this approach has
managed to find a small number of relatively good solutions, but worse than
when mixed-quantization is considered. The naïve approach, which is compared
against approaches that do not exploit accelerator properties, delivered good
results, but the proposed mixed-quantization method is better. The best trade-
offs are reported for the Eyeriss accelerator because its memory subsystem has
a higher impact on the overall energy.

Our method can also help in designing new hardware accelerators for CNN
because it can cheaply estimate the impact of complex quantization schemes on
the resulting performance and other parameters without the need to implement
the accelerator.

TABLE II: Reduction in memory energy $\Delta_{e}{{}_{m}}$ and relative
accuracy loss $\Delta_{acc}$ achieved by different types of automated
quantization algorithms.

Architecture | Eyeriss | Simba
---|---|---
Network | MobileNetV1 | MobileNetV2 | MobileNetV1 | MobileNetV2
|  | $\Delta_{e}{{}_{m}}$ | $\Delta_{acc}$
---|---
| $\Delta_{e}{{}_{m}}$ | $\Delta_{acc}$
---|---
| $\Delta_{e}{{}_{m}}$ | $\Delta_{acc}$
---|---
| $\Delta_{e}{{}_{m}}$ | $\Delta_{acc}$
---|---
Uniform 1) |  | -34.9% | -3.0%
---|---
-54.3% | -8.8%
| -33.0% | -0.7%
---|---
-47.9% | -7.9%
| -34.3% | -3.0%
---|---
-51.5% | -8.8%
| -27.3% | -0.7%
---|---
-41.0% | -7.9%
Naïve 2) |  | -57.2% | +0.2%
---|---
-63.4% | -0.5%
| -30.5% | +0.9%
---|---
-41.3% | +0.8%
-47.6% | +0.5%
-49.2% | +0.4%
| -48.7% | +0.3%
---|---
-58.0% | -0.5%
| -25.9% | +0.9%
---|---
-34.5% | +0.8%
-44.2% | +0.7%
-49.7% | +0.3%
Proposed |  | -45.2% | +0.8%
---|---
-56.8% | +0.4%
-61.9% | +0.1%
-63.1% | -0.4%
| -44.2% | +1.3%
---|---
-46.4% | +0.8%
-48.0% | +0.6%
-53.6% | +0.5%
| -51.7% | +0.4%
---|---
-53.9% | -0.0%
-58.0% | -0.5%
-58.6% | -0.5%
| -43.5% | +0.9%
---|---
-44.2% | +0.8%
-47.5% | -0.2%
-48.4% | -1.3%
1) Uniform quantization has been used e.g. in [4, 18].
2) A general automated quantization approach without exploring target
accelerator parameters e.g., [19].

##  VI Conclusions

We show that the introduction of mixed-precision quantization of weights and
activations in hardware accelerators of CNNs increases the number of available
CNN-to-hardware mappings and allows for better trade-offs between energy and
CNN errors than common quantization schemes. In contrast to other methods,
e.g. [3], these trade-offs can be investigated and quantified for general
accelerator architectures in software (i.e., in a design phase of a HW
accelerator) because the proposed extension of the Timeloop estimation tool
has this capability.

Acknowledgement This work was supported by the Czech Science Foundation
project 24-10990S and by the Ministry of Education, Youth and Sports of the
Czech Republic through the e-INFRA CZ (ID:90254).

## References

  * [1] V. Sze, Y. Chen _et al._ , “Efficient processing of deep neural networks: A tutorial and survey,” _Proc. IEEE_ , vol. 105, no. 12, pp. 2295–2329, 2017.
  * [2] G. Armeniakos, G. Zervakis _et al._ , “Hardware approximate techniques for deep neural network accelerators: A survey,” _ACM Comput. Surv._ , vol. 55, no. 4, pp. 83:1–83:36, 2023.
  * [3] K. Wang, Z. Liu _et al._ , “HAQ: Hardware-aware automated quantization with mixed precision,” in _CVPR ’19’_ , 2019, pp. 8604–8612.
  * [4] P. Gysel, J. Pimentel _et al._ , “Ristretto: A framework for empirical study of resource-efficient inference in convolutional neural networks,” _IEEE Trans. Neural Networks Learning Systems_ , vol. 29, no. 11, 2018.
  * [5] A. Parashar _et al._ , “Timeloop: A systematic approach to dnn accelerator evaluation,” in _ISPASS ’19_ , 2019, pp. 304–315.
  * [6] Y. N. Wu, J. S. Emer, and V. Sze, “Accelergy: An architecture-level energy estimation methodology for accelerator designs,” in _ICCAD ’19_ , 2019, pp. 1–8.
  * [7] K. Deb, A. Pratap _et al._ , “A fast and elitist multiobjective genetic algorithm: NSGA-II,” _IEEE Trans. Evolutionary Computation_ , vol. 6, no. 2, pp. 182–197, April 2002.
  * [8] S.-C. Kao and T. Krishna, “Gamma: Automating the hw mapping of dnn models on accelerators via genetic algorithm,” in _ICCAD ’20_ , 2020.
  * [9] A. Gholami, S. Kim _et al._ , “A survey of quantization methods for efficient neural network inference,” _CoRR_ , vol. abs/2103.13630, 2021.
  * [10] B. Jacob, S. Kligys _et al._ , “Quantization and training of neural networks for efficient integer-arithmetic-only inference,” _CoRR_ , vol. abs/1712.05877, 2017.
  * [11] T. Wang, K. Wang _et al._ , “APQ: Joint search for network architecture, pruning and quantization policy,” in _CVPR ’20_ , 2020, pp. 2075–2084.
  * [12] Y. Lin, M. Yang, and S. Han, “NAAS: Neural accelerator architecture search,” in _DAC ’21_ , 2021.
  * [13] H. Sharma, J. Park _et al._ , “Bit Fusion: Bit-level dynamically composable architecture for accelerating deep neural network,” in _ISCA ’18_ , 2018.
  * [14] H. Kwon, P. Chatarasi _et al._ , “Maestro: A data-centric approach to understand reuse, performance, and hardware cost of dnn mappings,” _IEEE Micro_ , vol. 40, no. 3, pp. 20–29, 2020.
  * [15] Y. Zhao, C. Li _et al._ , “Dnn-chip predictor: An analytical performance predictor for DNN accelerators with various dataflows and hardware architectures,” _CoRR_ , vol. abs/2002.11270, 2020.
  * [16] V. Sze, Y.-H. Chen _et al._ , “How to evaluate deep neural network processors: TOPS/W (alone) considered harmful,” _IEEE Solid-State Circuits Magazine_ , vol. 12, no. 3, pp. 28–41, 2020.
  * [17] Y. Hu, J. Zhai _et al._ , “Bitflow: Exploiting vector parallelism for binary neural networks on cpu,” in _2018 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2018, pp. 244–253.
  * [18] Y.-H. Chen, T. Krishna _et al._ , “Eyeriss: An energy-efficient reconfigurable accelerator for deep convolutional neural networks,” _IEEE Journal of Solid-State Circuits_ , vol. 52, no. 1, pp. 127–138, 2017.
  * [19] J. Choi, Z. Wang _et al._ , “PACT: parameterized clipping activation for quantized neural networks,” _CoRR_ , vol. abs/1805.06085, 2018.


\{ X, Y \}  := \mathrm{ad}(X) Y = T_e (\Theta) (X) Y = \frac{d^2}{dt dt} |_{t, s =0} \mathrm{Ad}_{  \mathcal{A} (\mathrm{exp}^{tX}) }  \mathrm{exp}^{sY} ~~ (\text{since } L_x^{\diamond_\mathcal{A}} = L^\diamond_{\mathcal{A} (x)}).
\end{align}

It follows from the expressions (<ref>) and (<ref>) that one obtains the commutative diagram

\begin{align}\label{first-diag}
\xymatrix{
\text{Lie group} \ar[d]_{\text{Diff}} \ar[rr]^{\text{pointed avg operator}} & & \text{pointed Lie rack} \ar[d]^{\text{Diff}} \\
\text{Lie algebra} \ar[rr]_{\text{avg operator}} & & \text{Leibniz algebra.}
\end{align}

Next, we consider pointed Lie group-racks and show that the differentiation at the identity element gives rise to the structure of a Lie-Leibniz algebra. We begin with the algebraic notion first.

A Lie-Leibniz algebra is a triple $(\mathfrak{g} , [~,~], \{ ~, ~ \})$ in which $(\mathfrak{g} , [~,~])$ is a Lie algebra and $(\mathfrak{g}, \{ ~, ~ \})$ is a Leibniz algebra satisfying the additional compatibility

\begin{align}\label{ll-iden}
\{ x, [y, z ] \} = [ \{ x, y \}, z] + [y, \{ x, z \}], \text{ for } x, y, z \in \mathfrak{g}.
\end{align}

Note that the compatibility condition (<ref>) says that all the left translations $\{ x, - \}$ by the Leibniz bracket is a derivation for the Lie bracket on $\mathfrak{g}$.

Any Lie algebra $(\mathfrak{g} , [~,~])$ can be regarded as a Lie-Leibniz algebra by taking $\{ ~, ~ \} = [~,~]$. In the following, we give two important classes of examples, one arising from representations and another arising from averaging operators on Lie algebras.

Let $(\mathfrak{g}, [~, ~])$ be a Lie algebra and $(V, \rho)$ be a representation (i.e. $\rho : \mathfrak{g} \times V \rightarrow V$ is a bilinear map satisfying $\rho ([x, y], v) = \rho (x , \rho (y, v)) - \rho (y, \rho (x, v))$, for $x, y \in \mathfrak{g}$ and $v \in V$). Then the direct sum $\mathfrak{g} \oplus V$ inherits the semidirect product Lie algebra and the hemi-semidirect product Leibniz algebra structure with the brackets

\begin{align*}
&\text{semidirect product Lie} \qquad \qquad \qquad [(x, u), (y, v)]_\ltimes := ([x, y], \rho (x, v) - \rho (y, u)),\\
&\text{hemi-semidirect product Leibniz} \qquad \{ (x, u), (y, v)\} := ([x, y], \rho (x, v) ),
\end{align*}

for $(x, u), (y, v) \in \mathfrak{g} \oplus V$. Then it is easy to see that $(\mathfrak{g} \oplus V, [~,~]_\ltimes , \{ ~, ~ \})$ is a Lie-Leibniz algebra.

Let $(\mathfrak{g}, [~,~])$ be a Lie algebra and $\mathcal{A} : \mathfrak{g} \rightarrow \mathfrak{g}$ be an averaging operator on $\mathfrak{g}$. Then $(\mathfrak{g}, \{ ~,~ \}_\mathcal{A})$ is a Leibniz algebra (called the descendent Leibniz algebra), where $\{ x, y \}_\mathcal{A}:= [ \mathcal{A} (x), y]$, for $x, y \in \mathfrak{g}$. Moreover, the triple $(\mathfrak{g} , [~,~], \{ ~,~ \}_\mathcal{A})$ is a Lie-Leibniz algebra.

For any $x, y, z \in \mathfrak{g}$, we have

\begin{align*}
\{ x, \{ y, z \}_\mathcal{A} \}_\mathcal{A} = [\mathcal{A}(x), [\mathcal{A} (y), z]] =~& [[ \mathcal{A}(x), \mathcal{A}(y)], z ] + [ \mathcal{A}(y), [\mathcal{A}(x), z]] \\
=~& \{ \{ x, y \}_\mathcal{A} , z \}_\mathcal{A} + \{ y, \{ x, z \}_\mathcal{A} \}_\mathcal{A}.
\end{align*}

Hence $(\mathfrak{g}, \{ ~,~ \}_\mathcal{A})$ is a Leibniz algebra. Moreover, we have

\begin{align*}
\{ x, [y, z ] \}_\mathcal{A} = [\mathcal{A}(x) , [y, z] ] =~& [[\mathcal{A} (x), y] , z ] + [y, [\mathcal{A}(x), z ]] \\
=~& [\{ x, y \}_\mathcal{A}, z] + [y, \{ x, z \}_\mathcal{A}].
\end{align*}

This shows that $(\mathfrak{g} , [~,~], \{ ~,~ \}_\mathcal{A})$ is a Lie-Leibniz algebra.

The next result gives the partial converse of the above proposition.

Every Lie-Leibniz algebra can be embedded into an averaging Lie algebra.

Let $(\mathfrak{g}, [~,~], \{ ~, ~ \})$ be a Lie-Leibniz algebra. Consider the Lie algebra $\mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie}$ associated to the Leibniz algebra $(\mathfrak{g}, \{ ~, ~ \})$. That is, $\mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} = \mathfrak{g} / \{ \mathfrak{g}, \mathfrak{g} \} $ equipped with the Lie bracket

\begin{align*}
\llbracket \langle x \rangle, \langle y \rangle \rrbracket := \langle \{ x, y \} \rangle, \text{ for } \langle x \rangle, \langle y \rangle \in \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie}.
\end{align*}

The Lie algebra $( \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie}, \llbracket ~, ~ \rrbracket )$ has a representation on the vector space $\mathfrak{g}$ with the action map $\rho : \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \times \mathfrak{g} \rightarrow \mathfrak{g}$ given by $\rho (\langle x \rangle, y) := \{ x, y \}$, for $\langle x \rangle \in \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie}$ and $y \in \mathfrak{g}$. It follows from (<ref>) that the map $\rho$ additionally satisfies

\begin{align*}
\rho (\langle x \rangle,  [y, z]) = [\rho (\langle x \rangle, y) , z] + [y , \rho (\langle x \rangle, z)], \text{ for } y, z \in \mathfrak{g}.
\end{align*}

In other words, the Lie algebra $\mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie}$ acts on the Lie algebra $(\mathfrak{g}, [~,~ ])$ by derivations. As a consequence, the space $\mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}$ inherits a Lie algebra structure with the Lie bracket

\begin{align*}
\llbracket (\langle x \rangle, y), (\langle x' \rangle, y') \rrbracket_\ltimes := \big( \llbracket \langle x \rangle, \langle x' \rangle \rrbracket, \rho (\langle x \rangle, y') - \rho (\langle x' \rangle,  y) + [y, y'] \big),
\end{align*}

for $(\langle x \rangle, y), (\langle x' \rangle, y') \in \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}.$ Moreover, the map $P : \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g} \rightarrow \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}$ defined by $P (\langle x \rangle, y) := (\langle y \rangle, 0)$ is an averaging operator on the Lie algebra $\big( \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}, \llbracket ~, ~ \rrbracket_\ltimes \big)$. That is, $\big( \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}, P  \big)$ is an averaging Lie algebra. Finally, it is easy to verify that the inclusion map $i : \mathfrak{g} \hookrightarrow  \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}, ~x \mapsto (0, x)$ is a morphism of Lie-Leibniz algebras, where the Lie-Leibniz algebra structure on $ \mathfrak{g}^{ \{ ~, ~ \} }_\mathrm{Lie} \oplus \mathfrak{g}$ is induced by the averaging operator $P$. Hence the result follows.

A pointed Lie group-rack is a triple $(G, \cdot, \smallblackdiamond)$ in which $(G, \cdot)$ is a Lie group with the identity element $e$ and $(G, \smallblackdiamond, e)$ is a pointed Lie rack satisfying $x \smallblackdiamond (yz) = (x \smallblackdiamond y) (x \smallblackdiamond z)$, for $x, y, z \in G$.

It follows that a pointed Lie group-rack is a pointed group-rack in which the underlying set is a smooth manifold such that the group operations (multiplication and inverse) and the rack operation are all smooth.

Let $G, H$ be two Lie groups and $\Phi : G \rightarrow \mathrm{Aut}(H)$ be a Lie group homomorphism. Then the cartesian product $G \times H$ carries the semidirect product Lie group structure with the multiplication

\begin{align*}
(x,h) \cdot (y, k) = (xy, h \Phi_x k), \text{ for } (x, h), (y, k) \in G \times H.
\end{align*}

On the cartesian product $G \times H$, there is also a Lie rack structure with the operation

\begin{align*}
(x,h) \diamond (y,k) = (xyx^{-1}, \Phi_x k), \text{ for } (x, h), (y, k) \in G \times H.
\end{align*}

Then it is easy to verify that $(G \times H, \cdot , \diamond)$ is a pointed Lie group-rack.

Let $(G, \cdot)$ be a Lie group and $\mathcal{A} : G \rightarrow G$ be a smooth pointed averaging operator. Then the triple $(G, \cdot , \diamond_\mathcal{A})$ is a pointed Lie group-rack, where $x \diamond_\mathcal{A} y := \mathcal{A} (x) y \mathcal{A} (x)^{-1}$ is the descendent rack structure.

Let $(G, \cdot, \smallblackdiamond)$ be a pointed Lie group-rack. Since $(G, \cdot)$ is a Lie group, it yields a Lie algebra structure on the distinguished tangent space $ T_e G$. On the other hand, $(G, \smallblackdiamond, e)$ is a pointed Lie rack implies that $T_e G$ inherits a Leibniz algebra structure with the bracket $\{ X, Y \} := \mathrm{ad}(X) Y$, for $X, Y \in T_e G$. Note that the compatibility condition $x \smallblackdiamond (yz) = (x \smallblackdiamond y) (x \smallblackdiamond z)$  simply means that all left multiplications $L_x^\smallblackdiamond$'s are group automorphisms. Hence for each $x \in G$, the map $T_e (L_x^\smallblackdiamond) : T_eG \rightarrow T_eG$ is a Lie algebra automorphism. Thus, in this case, we obtain a smooth map $\Theta : G \rightarrow \mathrm{Aut} (T_eG), x \mapsto T_e (L_x^\smallblackdiamond)$. By differentiating at the identity element $e$, we precisely get the map $\mathrm{ad} : T_e G \rightarrow \mathrm{Der} (T_e G)$. Hence for any $X, Y, Z \in T_e G$, we have

\begin{align*}
\{ X, [Y, Z] \} = \mathrm{ad} (X)[Y, Z]
=~& [\mathrm{ad} (X) Y, Z] + [Y, \mathrm{ad} (X) Z] \\
=~& [\{ X, Y \}, Z] + [Y, \{ X, Z \}].
\end{align*}

Therefore, we obtain a Lie-Leibniz algebra $(T_e G, [~,~], \{ ~, ~ \})$, called the tangent Lie-Leibniz algebra of the pointed Lie group-rack $(G, \cdot, \smallblackdiamond)$.

Let $(G, \cdot)$ be a Lie group and $\mathcal{A} : G \rightarrow G$ be a smooth pointed averaging operator. Then the Lie group $(G, \cdot)$ and the descendent pointed Lie rack $(G, \diamond_\mathcal{A}, e)$ together form a pointed Lie group-rack. On the infinitesimal level, the map $\mathcal{A}$ induces an averaging operator $\mathcal{A}_{* e} : \mathfrak{g} \rightarrow \mathfrak{g}$ on the corresponding Lie algebra $(\mathfrak{g} = T_e G, [~,~])$. The Lie algebra $(\mathfrak{g}, [~,~])$ and the descendent Leibniz algebra $(\mathfrak{g}, \{ ~, ~ \}_{\mathcal{A}_{* e}})$ forms a Lie-Leibniz algebra. Since the differentiation of the pointed Lie rack $(G, \diamond_\mathcal{A}, e)$ is the Leibniz algebra $(\mathfrak{g}, \{ ~, ~ \}_{\mathcal{A}_{* e}})$ and we have the commutative diagram (<ref>), we get another commutative diagram

\[
\xymatrix{
\text{Lie group} \ar[d]_{\text{Diff}} \ar[rr]^{\text{pointed avg operator}} & & \text{Lie group-rack} \ar[d]^{\text{Diff}} \\
\text{Lie algebra} \ar[rr]_{\text{avg operator}} & & \text{Lie-Leibniz algebra.}
\]

Finally, we consider smooth pointed averaging operators on a pointed Lie rack and show that their differentiation gives rise to averaging operators on the corresponding Leibniz algebra.

Let $(Q, \diamond, e)$ be a pointed Lie rack with the corresponding Leibniz algebra $(T_e Q, \{ ~, ~\})$. Let $\mathcal{A} : Q \rightarrow Q$ be a smooth pointed averaging operator on $(Q, \diamond, e)$.

(i) Then the differentiation $\mathcal{A}_{*e} : T_e Q \rightarrow T_e Q$ is a left averaging operator on the Leibniz algebra $(T_e Q, \{ ~, ~\})$.

(ii) The Leibniz algebra corresponding to the descendent pointed Lie rack $(Q, \diamond_\mathcal{A}, e)$ is precisely given by $(T_e Q, \{ ~, ~ \}_{\mathcal{A}_{*e}}).$

(i) Since $\mathcal{A} : Q \rightarrow Q$ is a smooth map with $\mathcal{A} (e) =e$, the differentiation $\mathcal{A}_{* e}$ is a linear map from $T_e Q$ to itself. On the other hand, $\mathcal{A}$ is an averaging operaor implies that $L_{\mathcal{A} (x)}^\diamond \circ \mathcal{A} = \mathcal{A} \circ L_{\mathcal{A} (x)}^\diamond$, for all $x \in Q$. By taking the differentiation at the distinguished element $e$, one obtains

\begin{align}\label{leib-r-a}
T_e (L^\diamond_{\mathcal{A} (x)}) \circ \mathcal{A}_{* e} = \mathcal{A}_{* e} \circ T_e (L^\diamond_{\mathcal{A} (x)}), \text{ for all } x \in Q.
\end{align}

Next, consider the smooth map $\Theta : Q \rightarrow GL (T_e Q)$ given by $\Theta (x) = T_e (L_x^\diamond)$, for $x \in Q$. Then it follows from (<ref>) that $\Theta (\mathcal{A} (x) ) \circ \mathcal{A}_{* e} = \mathcal{A}_{* e} \circ \Theta (\mathcal{A} (x) ) $, for all $x \in Q$. Hence for any $Y \in T_e Q$, we have

\begin{align*}
\Theta (\mathcal{A}(x)) (\mathcal{A}_{*e} (Y)) = \mathcal{A}_{* e} \big(  \Theta ( \mathcal{A} (x)) (Y)  \big).
\end{align*}

Finally, by differentiating the above identity with respect to $x$ at the distinguished point $e$, one get $\{ \mathcal{A}_{* e} (X) ,  \mathcal{A}_{* e} (Y) \} =  \mathcal{A}_{* e} \{  \mathcal{A}_{* e} (X), Y \}$, for any $X, Y \in T_e Q$. This proves that $\mathcal{A}_{* e} : T_e Q \rightarrow T_e Q$ is a left averaging operator on the Leibniz algebra $(T_e Q, \{ ~, ~ \}).$

(ii) Let $\Theta_\mathcal{A} : Q \rightarrow GL (T_e Q)$ be the smooth map given by $\Theta_\mathcal{A} (x) := T_e (L_x^{\diamond_\mathcal{A}})$, for $x \in Q$. Since $L_x^{\diamond_\mathcal{A}} = L_{\mathcal{A} (x)}^\diamond$, we have $\Theta_\mathcal{A} (x) = T_e ( L^\diamond_{\mathcal{A} (x)}) = \Theta (\mathcal{A}(x))$, for all $x \in Q$. Hence we get $\Theta_\mathcal{A} = \Theta \circ \mathcal{A}$. If $(T_e Q, \{ ~, ~ \}_{\diamond_\mathcal{A}})$ is the Leibniz algebra corresponding to the pointed Lie rack $(Q, \diamond_\mathcal{A}, e)$ then for any $X, Y \in T_e Q$, we have

\begin{align*}
\{ X, Y \}_{\diamond_\mathcal{A}} := T_e (\Theta_\mathcal{A}) (X) Y
=~& T_e (\Theta) (\mathcal{A}_{* e} (X)) Y  \quad (\text{as } \Theta_\mathcal{A} = \Theta \circ \mathcal{A}) \\
=~& \mathrm{ad} ( \mathcal{A}_{* e} (X) ) Y \quad (\because \mathrm{ad} = T_e (\Theta) ) \\
=~& \{ \mathcal{A}_{* e} (X), Y \} = \{ X, Y \}_{\mathcal{A}_{* e }}.
\end{align*}

This proves the desired result.

The last part of the above theorem yields the first commutative diagram of (<ref>). By taking the higher powers of the averaging operator $\mathcal{A}$, one obtain the second commutative diagram of (<ref>),

\begin{align}\label{two-comm}
\xymatrix{
\substack{\text{pointed Lie rack} \\ (Q, \diamond, e)} \ar[d]_{\text{Diff}} \ar[rr]^{ \substack{\text{pointed avg operator} \\ \mathcal{A}}} & & \substack{\text{pointed Lie rack} \\ (Q, \diamond_\mathcal{A}, e)}   \ar[d]^{\text{Diff}} &  \substack{\text{pointed Lie rack} \\ (Q, \diamond, e)} \ar[d]_{\text{Diff}} \ar[rr]^{ \substack{\text{pointed avg operator} \\ \mathcal{A}^k}} & & \substack{\text{pointed Lie rack} \\ (Q, \diamond_{\mathcal{A}^k}, e)}   \ar[d]^{\text{Diff}} \\
\substack{ \text{Leibniz algebra} \\ (T_e Q, \{ ~, ~ \})} \ar[rr]_{ \substack{\text{left avg operator} \\ \mathcal{A}_{*e} } } & & \substack{ \text{Leibniz algebra} \\ (T_e Q, \{ ~, ~ \}_{\mathcal{A}_{*e}})} &   \substack{ \text{Leibniz algebra} \\ (T_e Q, \{ ~, ~ \})} \ar[rr]_{ \substack{\text{left avg operator} \\ {\mathcal{A}^k_{*e} }}  } & & \substack{ \text{Leibniz algebra} \\ (T_e Q, \{ ~, ~ \}_{\mathcal{A}^k_{*e}}).}
\end{align}

One may also define pointed Lie rack-pairings and pointed Lie di-racks, and find their possible connections to di-Leibniz algebras. We will consider these topics in some other work.

§ SOME GENERALIZATIONS AND FUTURE WORKS

In this section, we will give a quick sketch of some generalizations of the topics covered in this paper and mention some future works.

Our main objective of the present paper was to study averaging operators on various algebraic structures and figure out the corresponding induced structures. Many of the structures related to racks are intimately connected to the corresponding structures related to groups. For instance, rack-pairings are related to two-sided skew braces by conjugation. Recently, [5] Bardakov and Gubarev introduced the notion of a skew $k$-brace in the study of Rota-Baxter operators on groups.

A (left) skew $k$-brace is a tuple $(G, \cdot_0, \cdot_1, \ldots, \cdot_k)$ consisting of a nonempty set $G$ with binary operations $\cdot_0, \cdot_1, \ldots, \cdot_k : G \times G \rightarrow G$ such that

- $(G, \cdot_i)$ is a group, for each $0 \leq i \leq k$,

- for each $1 \leq i \leq k$, the triple $(G, \cdot_{i-1}, \cdot_i)$ is a skew brace, that is,

\begin{align*}
x \cdot_i (y \cdot_{i-1} z) = (x \cdot_i y) \cdot_{i-1} x^{-1 (\cdot_{i-1})} \cdot_{i-1} (x \cdot_i z), \text{ for } x, y, z \in G.
\end{align*}

Here $x^{-1 (\cdot_{i-1})}$ is the inverse of the element $x$ in the group $(G, \cdot_{i-1})$.

One may also define the notion of a two-sided skew $k$-brace as a skew $k$-brace $(G, \cdot_0, \cdot_1, \ldots, \cdot_k)$ such that $(G, \cdot_{i-1}, \cdot_i)$ is a two-sided skew brace, for any $1 \leq i \leq k$.

Generalizing Definition <ref> in the context of racks, one gets the following notion.

Let $k$ be a natural number. A rack $k$-pairing is a tuple $(Q, \diamond_0, \diamond_1, \ldots, \diamond_k)$ consisting of a nonempty set $Q$ equipped with $(k+1)$ binary operations $\diamond_0, \diamond_1, \ldots, \diamond_k : Q \times Q \rightarrow Q$ such that the following conditions are hold:

- $(Q, \diamond_i)$ is a rack, for each $0 \leq i \leq k$,

- for all $x, y, z \in Q$ and $1 \leq i \leq k$,

\begin{align*}
x \diamond_i (y \diamond_{i-1} z) = (x \diamond_i y) \diamond_{i-1} (x \diamond_i z).
\end{align*}

It follows from the above definition that a rack $1$-pairing is simply a rack-pairing. More generally, a rack $k$-pairing is a tuple $(Q, \diamond_0, \diamond_1, \ldots, \diamond_k)$ such that $(Q, \diamond_{i-1} , \diamond_i)$ is a rack-pairing, for each $1 \leq i \leq k$. The following result is a generalization of Proposition <ref>.

Let $(G, \cdot_0, \cdot_1, \ldots, \cdot_k)$ be a two-sided skew $k$-brace. Then the tuple $(G, \diamond_0, \diamond_1, \ldots, \diamond_k)$ is a rack $k$-pairing, where

\begin{align*}
x \diamond_i y := x \cdot_i y \cdot_i x^{-1 (\cdot_i)}, \text{ for any } x, y \in G \text{ and } 0 \leq i \leq k.
\end{align*}

In Proposition <ref>, we have shown that any averaging operator on a rack yields the structure of a rack-pairing. In the following result, we show that an averaging operator also induces a rack $k$-pairing.

Let $(Q, \diamond)$ be a rack and $\mathcal{A}: Q \rightarrow Q$ be an averaging operator on $Q$. Define binary operations $\diamond_1, \ldots, \diamond_k : Q \times Q \rightarrow Q$ by

\begin{align*}
x \diamond_i y = \mathcal{A}^i (x) \diamond y, \text{ for all } x, y \in Q.
\end{align*}

Then $(Q, \diamond_0 = \diamond, \diamond_1, \ldots, \diamond_k)$ is a rack $k$-pairing.

It is not hard to define pointed rack $k$-pairings and related structures. Is a pointed rack $k$-pairing always embedded into an averaging rack? On the other hand, one can define pointed Lie rack $k$-pairings by considering smooth structures. It is interesting to find their infinitesimal version. Another interesting question could be the use of rack $k$-pairings in knot theory.

The notion of $n$-valued groups (also called multivalued groups) are generalization of usual abstract groups (see, for instance, [13]). The theory of $n$-valued groups generalizes many results from classical group theory.

In a subsequent paper, we aim to define Rota-Baxter operators and averaging operators on $n$-valued groups and the concept of $n$-valued skew braces. We are also looking forward to defining $n$-valued Lie algebras and finding possible connections to Lie $n$-valued groups.

Another higher generalization of racks was introduced by Biyogmam [10] who considered $n$-racks. Note that $n$-valued racks and $n$-racks are two different notions. In the former one, the defining operation takes values in the $n$ cartesian product of the underlying set, while in the latter one, the defining operation is defined on the $n$ cartesian product of the underlying set. The author in [10] also considered Lie $n$-racks and showed that the differentiation of a pointed Lie $n$-rack yields an $n$-Leibniz algebra. In a forthcoming article, we will define the notion of an averaging operator on $n$-racks and show that such an operator induces an averaging operator on the corresponding $n$-Leibniz algebra. We also aim to formulate and study the $n$-ary generalization of the set-theoretical solutions of the Yang-Baxter equation and show that an $n$-rack naturally gives rise to a solution of such an equation.

Acknowledgements. The author would like to thank the Department of Mathematics, IIT Kharagpur for providing the beautiful academic atmosphere where the research has been conducted.

Data Availability Statement. Data sharing does not apply to this article as no new data were created or analyzed in this study.

[1] M. Aguiar, Pre-Poisson algebras, Lett. Math. Phys. 54 (2000), 263-277.

[2] C. Bai, O. Bellier, L. Guo and X. Ni, Splitting of operations, Manin products and Rota-Baxter operators, Int.
Math. Res. Not. 2013, no. 3 (2013), 485-524.

[3] C. Bai, L. Guo, Y. Sheng and R. Tang, Post-groups, (Lie-)Butcher groups and the Yang–Baxter equation, Math. Ann. to appear, https://doi.org/10.1007/s00208-023-02592-z

[4] V. G. Bardakov and V. A. Bovdi, Rota-Baxter and averaging operators on racks and rack algebras, arXiv preprint, arXiv:2402.11660.

[5] V. G. Bardakov and V. Gubarev, Rota-Baxter groups, skew left braces, and the Yang-Baxter equation, J. Algebra 596 (2022), 328-351.

[9] J. Baez and A. S. Crans, Higher-dimensional algebra VI. Lie $2$-algebras, Theor. Appl. Categ. 12 (2004), 492-538.

[10] G. R. Biyogmam, Lie $n$-racks, Comptes Rendus. Mathématique, 349, no. 17-18 (2011), 957-960.

[11] R. Bonezzi and O. Hohm, Leibniz gauge theories and infinity structures, Commun. Math. Phys. 377 (2020), 2027–2077.

[12] B. Brainerd, On the structure of averaging operators, J. Math. Anal. 5 (1962), 347-377.

[13] V. M. Buchstaber, $n$-valued groups: theory and applications, Moscow Math. J. 6,
no. 1 (2006), 57-84.

[14] W. Cao, An algebraic study of averaging operators, Ph.D. thesis, Rutgers University at Newark (2000).

[16] A. Das, Controlling structures, deformations and homotopy theory for averaging algebras, https://arxiv.org/abs/2303.17798.

[18] A. Das and N. Rathee, Extensions and automorphisms of Rota-Baxter groups, J. Algebra 636 (2023), 626-665.

[19] A. Das and S. Sen, $2$-term averaging $L_\infty$-algebras and non-abelian extensions of averaging Lie algebras, J. Algebra 644 (2024), 126-151.

[21] M. Eisermann, Yang-Baxter deformations of quandles and racks, Alg. Geom. Topology 5 (2005), 537-562.

[22] P. Etingof, T. Schedler and A. Soloviev, Set-theoretical solutions to the quantum Yang-Baxter equation, Duke Math. J. 100 (2) (1999), 169-209.

[23] M. Goncharov, Rota-Baxter operators on cocommutative Hopf algebras, J. Algebra 582 (2021), 39-56.

[24] L. Guarnieri and L. Vendramin, Skew braces and the Yang-Baxter equation, Mathematics of Computation 86, no. 307 (2017), 2519-2534.

[25] L. Guo, H. Lang and Y. Sheng, Integration and geometrization of Rota-Baxter Lie algebras, Adv. Math. 387 (2021), 107834.

[27] J. Kampde Fériet, Létat actuel du probléme de la turbulaence (I and II), La Sci. Aérienne 3 (1934), 9-34 &
4 (1935), 12-52.

[28]  J. L. Kelley, Averging operators on $C_\infty (X)$, Illinois J. Math. 2 (1958), 214-223.

[29] M. K. Kinyon, Leibniz algebras, Lie racks, and digroups, J. Lie Theory 17 (2007), 99-114.

[30]  A. Kotov and T. Strobl, The embedding tensor, Leibniz–Loday algebras, and their higher Gauge theories,
Commun. Math. Phys. 376 (2020), 235–258.

[32] Y. Li, Y. Sheng and R. Tang, Relative Rota-Baxter operator of weight $0$ on groups, pre-groups, braces, the Yang-Baxter equation and $T$-structures, arXiv preprint, arXiv:2312.01337.

[33]  J.-L. Loday, Dialgebras, Dialgebras and related operads, 7-66, Lecture Notes in Math., 1763, Springer, Berlin, 2001.

[34]  J. B. Miller, Averaging and Reynolds operators on Banach algebra I, Representation by derivation and
antiderivations, J. Math. Anal. Appl. 14 (1966), 527-548.

[35]  S.-T. C. Moy, Characterizations of conditional expectation as a transformation on function spaces, Pacific
J. Math. 4 (1954), 47-63.

[36] D. Quillen, Rational homotopy theory, Ann. of Math. 90 (1969), 205-295.

[37] J. Pei, C. Bai, L. Guo and X. Ni, Replicators, Manin white product of binary operads and average operators, New
Trends in Algebras and Combinatorics, (2020), pp. 317-353.

[38] J. Pei and L. Guo, Averaging algebras, Schröder numbers, rooted trees and operads, J. Algebraic Combinatorics 42 (2015), 73-109.

[39]  O. Reynolds, On the dynamic theory of incompressible viscous fluids and the determination of the criterion, Phil. Trans. Roy. Soc. A 136 (1895), 123-164.

[40] W. Rump, Braces, radical rings, and the quantum Yang–Baxter equation, J. Algebra 307 (2007), 153-170.

[41] Y. Sheng, R. Tang and C. Zhu, The controlling $L_\infty$-algebra, cohomology and homotopy of embedding tensors and Lie–Leibniz triples, Commun. Math. Phys. 386 (2021), 269-304.

[42] A. Soloviev, Non-unitary set-theoretical solutions to the quantum Yang–Baxter equation, Math. Res. Lett. 7 (2000), no. 5-6, 577–596.

[43] S. Trappeniers, On two-sided skew braces, J. Algebra 631 (2023) 267-286.

[44] V. Turaev, Multi-quandles and topological pairs, https://arxiv.org/abs/2205.00951v2.


# Multitask Learning for Improved Late Mechanical Activation Detection of
Heart From Cine DENSE MRI

###### Abstract

The selection of an optimal pacing site, which is ideally scar-free and late
activated, is critical to the response of cardiac resynchronization therapy
(CRT). Despite the success of current approaches formulating the detection of
such late mechanical activation (LMA) regions as a problem of activation time
regression, their accuracy remains unsatisfactory, particularly in cases where
myocardial scar exists. To address this issue, this paper introduces a multi-
task deep learning framework that simultaneously estimates LMA amount and
classify the scar-free LMA regions based on cine displacement encoding with
stimulated echoes (DENSE) magnetic resonance imaging (MRI). With a newly
introduced auxiliary LMA region classification sub-network, our proposed model
shows more robustness to the complex pattern cause by myocardial scar,
significantly eliminates their negative effects in LMA detection, and in turn
improves the performance of scar classification. To evaluate the effectiveness
of our method, we tests our model on real cardiac MR images and compare the
predicted LMA with the state-of-the-art approaches. It shows that our approach
achieves substantially increased accuracy. In addition, we employ the
gradient-weighted class activation mapping (Grad-CAM) to visualize the feature
maps learned by all methods. Experimental results suggest that our proposed
model better recognizes the LMA region pattern.

##  1 Introduction

Cardiac resynchronization therapy (CRT) is widely used to treat cardiac
conduction system disorders, such as left bundle branch block and intrinsic
myocardial diseases [1, 2]. Standard CRT, however, suffers from a low response
rate that is approximately $40\%$ [3, 4]. It has been shown that the choice of
an optimal pacing site has considerable influence on non-response rate, which
could be effectively decreased by implanting the CRT left ventricular (LV)
lead at areas with late activation [5, 6]. Developing accurate LV activation
time measure methods is critical to improve CRT response rate.

A recent work [7] developed a semi-automatic LV late activation detection
method based on the circumferential myocardium strain that is computed from
cine DENSE displacement field. The authors first build strain matrices that
contain spatial-temporal strain information, and apply active contour
algorithm [8] to estimate the activation time of each location. However, the
applicability of this method is limited due to the requirement of a case-by-
case model parameter tuning and low optimization efficiency because of an
iterative search scheme. To alleviate these problems, a deep learning based
late activation prediction model was first introduced in [9]. Despite the
success of achieving significantly higher accuracy and speed, the predicted
LMA remains unsatisfactory in challenging cases, i.e., when myocardial scars
occur. The existence of scar tissues severely affects and produces activation
features that can mislead the LMA detection.

In this paper, we introduce a multi-task learning (MTL) network that
simultaneously detects LMA and classify scar-free LMA regions based on cine
DENSE MR images. In contrast to the previously proposed deep networks that
detect LMA without considering the negative scar effects [9], we jointly train
an LMA estimation regression network, as well as an auxiliary LMA region
classification network that guides the LMA to avoid non-LMA regions. To
evaluate the effectiveness of our method, we run tests using real cardiac MR
images and compare the predicted LMA regions with the state-of-the-art deep
learning-based approaches [9]. Experimental results show that our approach
achieves substantially increased LMA estimation accuracy. Additionally, the
attention heatmaps generated by the gradient-weighted class activation mapping
(Grad-CAM) [10] indicates that our proposed model focus better on the real LMA
features.

##  2 Background

Strain analysis. The time to the onset of circumferential shortening (TOS) has
been proven to have a close linear relationship with electrical activation
time and the ability to detect LMA regions [11]. Following the same
principles, we use the TOS of segmental myocardial regions to measure the
cardiac mechanical activation. In order to estimate the TOS, we use
circumferential myocardium strain (computed on the displacement field of cine
DENSE images) [12] that shows several advantages of: (i) representing
myocardial contraction in a relatively low-dimensional space that reduces the
cost of network training and risk of overfitting; (ii) robustness to motion
artifacts; and (iii) the effectiveness of dyssynchrony quantification and high
reproducibility for validation [12].

Building Strain Matrix. Consider a cine DENSE image with $T$ time frames. For
each time frame, we compute a $N$-dimensional strain vector from a number of
$N$ myocardial sectors, beginning from the the middle of the intersection
points of the left and right ventricle and following counter-clockwise order.
A $N\times T$ strain matrix that includes information from all time frames is
further built by concatenating all strain vectors across time. A curve of TOS
values of all myocardial sectors can be overlaid on the strain matrix (see
examples in Fig. 1(d) and (e)).

Fig. 1: A flowchart of strain analysis of cine DENSE MRIs. Here, $t_{1}$,
$t_{2}$ and $t_{3}$ represent different time frames. (a) LV segmental
myocardial strains; and (b) generating a strain matrix of all time frames.

Challenging strain pattern introduced by myocardial scarring. Strain in the
scarring region would show stretching in the early frames due to the
compression from neighbor sectors, but show no or little contraction in the
late frames, which may fool the existing algorithm to wrongly recognize it as
LMA. Two strain matrices with TOS curves are shown in Fig. 2, where in (a)
there is no scar region and the predicted TOS curves from all methods overlap
well with the ground truth. In (b), however, both the previously proposed
algorithms are misled by the strong early stretching in the scar region and
make wrong prediction, while our method gives much more reliable result.

##  3 Our Method

Enhancing algorithm robustness to myocardial scarring by multi-task learning.
To solve the problem of over-estimating the LMA regions, an auxiliary task of
determining the LMA region is introduced to solve the problem. When the model
is trained only to reduce the difference between ground-truth and predicted
TOS curves, the model may focus much more on predicting accurate TOS values
than the range of LMA region. For example, in Fig. 2(b), while visually the
previously proposed models make obvious mistake in the scar region, this error
actually contribute little to the total error. On the contrary, the auxiliary
LMA region classification task will result in much larger penalty for
incorrect LMA region prediction, and thus guide the model to recognize the
misleading regions and give more accurate TOS prediction. Besides, adding the
auxiliary has been theoretically proven to also adding extra regularization
and increase the generalization ability [13, 14].

Fig. 2: Examples of strain matrices with ground truth and predicted TOS
curves.

Proposed MTL framework. As is shown in Fig. 3, three main components are
included in the framework: the joint convolutional sub-network, the regression
sub-network, and the classification sub-network. Given the input strain
matrix, the joint convolutional sub-network first extracts deep features that
are used by both following sub-networks. Taking the deep features as input,
the regression sub-network predicts the TOS values of all sectors, while the
classification sub-network determines whether each sector has LMA.

Fig. 3: An overview of our proposed multi-task deep learning framework.

The joint convolutional sub-network contains a sequence of convolutional, max-
pooling, batch normalization and ReLU layers. Both the classification sub-
network and regression sub-network contain the same types of layers as the
joint network and multiple fully connected layers. A shifted leaky ReLU layer,
which is proposed in [9], is adopted as the final activation function for the
regression network. Compared with other common activation functions, the
shifted leaky ReLU layer prevents negative TOS predictions as the TOS value
should be always positive, and at the same time avoids the gradient vanishing
problem of ReLU [15].

The loss is the weighted sum of the regression loss, the classification loss
and the regularization term:

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{reg}+\lambda_{cls}\cdot\mathcal{L}_{cls}+r\cdot\mathcal{R}(w)$
$\displaystyle=\|\hat{\mathbf{t}}-\mathbf{t}\|_{2}+\lambda_{cls}\cdot\text{CE}(\hat{\mathbf{l}},\mathbf{l})+r\cdot|\mathbf{w}|,$

where $\mathbf{t}$ and $\hat{\mathbf{t}}$ are the ground-truth and predicted
TOS curves, respectively. $\mathbf{l}$ and $\hat{\mathbf{l}}$ are the ground-
truth and predicted one-hot label vectors, which indicate whether each sector
has late activation or not. $\mathbf{w}$ is the trainable parameter of the
whole network. $|\cdot|$ refers to the $L_{1}$ norm, and
$\mathcal{R}(\mathbf{w})$ works as model regularizer with weighting parameter
$r$. $\text{CE}(\cdot,\cdot)$ denotes the cross entropy, which is defined as
$\text{CE}(\hat{\mathbf{l}},\mathbf{l})=-\frac{1}{N_{s}}\sum_{i}^{N_{s}}\sum_{j}^{N_{c}}\hat{\mathbf{l}}_{ij}\log{(\mathbf{l}_{ij})}$.

Interpreting the models with Grad-CAM. To better check whether our model
recognizes the scar region, we adapts Gradient-weighted Class Activation
Mapping (Grad-CAM) algorithm [10]. Grad-CAM highlights the regions which
affects the model prediction most. As a good LMA prediction model should make
prediction based on both pre-stretching and late activation patterns, the
corresponding Grad-CAM result should highlight both regions. Since the
original Grad-CAM works for only a single dimension of the output, we
visualize the Grad-CAM of the dimension corresponding to the central sector of
the LMA region.

##  4 Experiential Results

###  4.1 Experimental Dettings

Data acquisition. Data were acquired using a 1.5T MR scanner (Avanto, Siemens,
Erlangen, Germany) with a four-channel phased-array radiofrequency coil. Cine
DENSE was performed in 4 short-axis planes at basal, two mid-ventricular, and
apical levels. Cine DENSE parameters included a temporal resolution of 17 ms,
pixel size of $2.65\times 2.65\text{ mm}^{2}$ and slice thickness = $8$ mm.
Displacement was encoded in two orthogonal directions and a spiral k-space
trajectory was used with 6 interleaves per image. Other parameters included:
field of view $=240\times 240\text{ mm}^{2}$, displacement encoding frequency
$k_{e}=0.1\text{ cycles/mm}$, flip angle $15^{\circ}$ and echo time
$=1.08\text{ ms}$. We have $282$ strain matrices in total from $57$ patients,
in which $262$ strain matrices from $53$ patients are used for training, and
$20$ matrices from $4$ other patients are used for test. All of the ground-
truth TOS curves are manually labeled by experts.

Data augmentation and pre-processing. Utilizing the circular nature of strain
matrices, the data are augmented by shifting the matrices and labels along the
sector dimension, which generates new data while keeping the meaningful
cardiac strain pattern. Besides, mixup [16] is also used in our experiments,
which generates new data by taking the weighted average of two arbitrary input
data and the weighted average of corresponding labels. $3582$ strain matrices
and corresponding labels are generated with these methods. Together with the
raw training data, in total $3844$ data are used for training. The last $5$
frames of strain matrices are usually noisy and thus trimmed. Moreover,
horizontal zero padding and vertical duplication padding are applied to unify
the size of strain matrices.

Backbone network structure. In our experiment, the joint sub-network contain 3
successive convolutional layers with the same output channel size as 16 and
kernel size as 3. Each the classification and regression sub-network contains
3 convolutional layers, which are followed by fully connected layers with
decreasing dimensions. Each of the convolutional layers are followed by a
batch normalization layer, a max-pooling layer and a ReLU layer. Each of the
fully connected layers are followed by a batch normalization layer and a ReLU
layer, except for the last layers, which are followed by the leaky ReLU layers
as mentioned in the methodology section.

Parameter tuning. Both the regression and multi-task networks are trained on
an Nvidia 2080Ti GPU with $11$ GB RAM over maximum $1000$ epochs using an Adam
optimizer [17] with early stop. The optimal learning rate, batch size,
classification loss weight $\lambda_{cls}$ and regularizer weight $r$ are
decided through parameter-tuning process using the tree-structured parzen
estimator (TPE) [18]. The empirical optimal learning rate, batch size and $r$
for the regression network are $1E-2$, $64$, $10$ and $0.5$, respectively. For
the multi-task network, the empirical optimal learning rate, batch size,
$\lambda_{cls}$ and $r$ are $1E-3$, $64$, $10$ and $0.1$, respectively. The
optimal hyper-parameters for the active contour based method is tuned using
Bayesian optimization algorithm [19] on an Intel i5-1035G4 CPU, and the
optimal $\lambda=0.60206$, $\beta=4.8778$ and $\gamma=0.43512$.

###  4.2 TOS Prediction

Fig. 4: Left: Exemplar prediction results of different methods; right: Boxplot
of average prediction error on single strain matrix.

Accuracy. Three examples of test data prediction of compared methods are shown
in Fig. 4(a). It can be observed that our proposed multi-task network offers
the most accurate late activation prediction, in terms of both location and
magnitude. The average activation time prediction error on the test dataset of
the activation contour-based method, the regression network and our proposed
network are $207.5$ms, $97.4$ms and $66.9$ms, respectively, which are also
shown in Fig. 4(b), from which we can see our proposed method offers overall
more accurate and robust prediction. Besides, in the case 1 our method
successfully recognized the scar region, which is described in the background
section, while the other two mistakenly considered the scarring region as late
activation region.

3D activation maps. We further visualization the distribution of activation
time on the myocardium surface as 3D activation maps as are shown in Fig. 5.
The regions with large TOS values are shown in red and thus are severe late
activating regions, and the regions having small TOS values are shown in blue
and are normal regions. Our proposed multi-task network successfully captures
both the two major late activation regions at apex inferior and basal
inferolateral, and gains significant higher accuracy in both late activation
region localization and TOS value prediction.

Fig. 5: 3D Activation Maps built from prediction of different TOS prediction
algorithms.

Model interpretation with Grad-CAM. The Grad-CAM results for the regression
model and our MTL model are shown in Fig. 6. While we expect the network to
recognize the pre-stretching followed by late contraction pattern to identify
late activation region, the regression network takes almost only the deep blue
late contraction region into consideration and ignores the shallow yellow pre-
stretching region when detecting LMA. This explains why it can be misled by
the myocardial scarring pattern. On the contrary, the multi-task network
focuses on almost all frames near the late activation region, which includes
both pre-stretching and late contraction regions.

Fig. 6: First row: Strain matrices and TOS curves; second row: Strain matrix
overlapped with inverse Grad-CAM, where the highlighted regions contribute
most to the models’ decision.

##  5 Discussion

Our proposed method can be further improved by utilizing additional medical
information. For example, we can introduce the scar location information and
convert the binary classification task into a multi-class classification task,
which detects LMA and scar region in the same time. This will not only further
increase the performance on LMA estimation, but also performing valuable scar
detection task in the same time.

Meanwhile, our current 3D activation maps of the LV surface is estimated by
bilinear interpolation. However, due to the large space between neighbor
scanning slices, slight prediction error can cause large different on the 3D
activation map. Therefore, we will investigate a better reconstruction to
reduce such errors.

##  6 Conclusion

This paper presented a novel multi-task learning network for an improved
cardiac LMA detection. Our newly introduced classification task of identifying
non-late-activating regions substantially improves the performance of LMA
detection. We evaluated the proposed multi-task network on the cardiac
activation time prediction task on CMR images. Experimental results
demonstrated that our proposed method outperforms the state-of-the-art [7, 9].
We further examined the LMA feature patterns captured by our model and the
baseline algorithms. Results suggested that our method focuses better on the
real late activation pattern of pre-stretching followed by contraction. We
commit to make all data and code produced in this work publicly available
online upon the acceptance of this manuscript.

## References

  * [1] William T Abraham, Westby G Fisher, Andrew L Smith, David B Delurgio, Angel R Leon, Evan Loh, Dusan Z Kocovic, Milton Packer, Alfredo L Clavell, David L Hayes, et al.,  “Cardiac resynchronization in chronic heart failure,”  New England Journal of Medicine, vol. 346, no. 24, pp. 1845–1853, 2002.
  * [2] Arthur J Moss, W Jackson Hall, David S Cannom, Helmut Klein, Mary W Brown, James P Daubert, NA Mark Estes III, Elyse Foster, Henry Greenberg, Steven L Higgins, et al.,  “Cardiac-resynchronization therapy for the prevention of heart-failure events,”  New England Journal of Medicine, vol. 361, no. 14, pp. 1329–1338, 2009.
  * [3] Eugene S Chung, Angel R Leon, Luigi Tavazzi, Jing-Ping Sun, Petros Nihoyannopoulos, John Merlino, William T Abraham, Stefano Ghio, Christophe Leclercq, Jeroen J Bax, et al.,  “Results of the predictors of response to crt (prospect) trial,”  Echocardiography, vol. 2608, pp. 2616, 2008.
  * [4] Derek V Exner, Angelo Auricchio, and Jagmeet P Singh,  “Contemporary and future trends in cardiac resynchronization therapy to enhance response,”  Heart rhythm, vol. 9, no. 8, pp. S27–S35, 2012.
  * [5] Kenneth C Bilchick, Sujith Kuruvilla, Yasmin S Hamirani, Raghav Ramachandran, Samantha A Clarke, Katherine M Parker, George J Stukenborg, Pamela Mason, John D Ferguson, J Randall Moorman, et al.,  “Impact of mechanical activation, scar, and electrical timing on cardiac resynchronization therapy response and clinical outcomes,”  Journal of the American College of Cardiology, vol. 63, no. 16, pp. 1657–1666, 2014.
  * [6] Raghav Ramachandran, Xiao Chen, Christopher M Kramer, Frederick H Epstein, and Kenneth C Bilchick,  “Singular value decomposition applied to cardiac strain from mr imaging for selection of optimal cardiac resynchronization therapy candidates,”  Radiology, vol. 275, no. 2, pp. 413–420, 2015.
  * [7] Daniel A Auger, Kenneth C Bilchick, Jorge A Gonzalez, Sophia X Cui, Jeffrey W Holmes, Christopher M Kramer, Michael Salerno, and Frederick H Epstein,  “Imaging left-ventricular mechanical activation in heart failure patients using cine dense mri: Validation and implications for cardiac resynchronization therapy,”  Journal of Magnetic Resonance Imaging, vol. 46, no. 3, pp. 887–896, 2017.
  * [8] Michael Kass, Andrew Witkin, and Demetri Terzopoulos,  “Snakes: Active contour models‖ international journals of computer vision,”  1988\.
  * [9] Jiarui Xing, Sona Ghadimi, Mohamad Abdi, Kenneth C. Bilchick, Frederick H. Epstein, and Miaomiao Zhang,  “Deep networks to automatically detect late-activating regions of the heart,”  in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), 2021, pp. 1902–1906.
  * [10] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra,  “Grad-cam: Visual explanations from deep networks via gradient-based localization,”  in Proceedings of the IEEE international conference on computer vision, 2017, pp. 618–626.
  * [11] Bradley T Wyman, William C Hunter, Frits W Prinzen, and Elliot R McVeigh,  “Mapping propagation of mechanical activation in the paced heart with mri tagging,”  American Journal of Physiology-Heart and Circulatory Physiology, vol. 276, no. 3, pp. H881–H891, 1999.
  * [12] Loren P Budge, Adam S Helms, Michael Salerno, Christopher M Kramer, Frederick H Epstein, and Kenneth C Bilchick,  “Mr cine dense dyssynchrony parameters for the evaluation of heart failure: comparison with myocardial tissue tagging,”  JACC: Cardiovascular Imaging, vol. 5, no. 8, pp. 789–797, 2012\.
  * [13] Jonathan Baxter,  “A bayesian/information theoretic model of learning to learn via multiple task sampling,”  Machine learning, vol. 28, no. 1, pp. 7–39, 1997.
  * [14] Rich Caruana,  “Multitask learning,”  Machine learning, vol. 28, no. 1, pp. 41–75, 1997.
  * [15] Andrew L Maas, Awni Y Hannun, and Andrew Y Ng,  “Rectifier nonlinearities improve neural network acoustic models,”  in Proc. icml, 2013, vol. 30, p. 3.
  * [16] Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz,  “mixup: Beyond empirical risk minimization,”  arXiv preprint arXiv:1710.09412, 2017.
  * [17] Diederik P Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.
  * [18] James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl,  “Algorithms for hyper-parameter optimization,”  Advances in neural information processing systems, vol. 24, 2011\.
  * [19] Jasper Snoek, Hugo Larochelle, and Ryan P Adams,  “Practical bayesian optimization of machine learning algorithms,”  Advances in neural information processing systems, vol. 25, 2012\.