# End-to-End Learned Block-Based Image Compression with Block-Level Masked
Convolutions and Asymptotic Closed Loop Training

Fatih Kamisli<EMAIL_ADDRESS>Department of Electrical and Electronics
Engineering, Middle East Technical University, Turkey

###### Abstract

Learned image compression research has achieved state-of-the-art compression
performance with auto-encoder based neural network architectures, where the
image is mapped via convolutional neural networks (CNN) into a latent
representation that is quantized and processed again with CNN to obtain the
reconstructed image. CNN operate on entire input images. On the other hand,
traditional state-of-the-art image and video compression methods process
images with a block-by-block processing approach for various reasons. Very
recently, work on learned image compression with block based approaches have
also appeared, which use the auto-encoder architecture on large blocks of the
input image and introduce additional neural networks that perform
intra/spatial prediction and deblocking/post-processing functions. This paper
explores an alternative learned block-based image compression approach in
which neither an explicit intra prediction neural network nor an explicit
deblocking neural network is used. A single auto-encoder neural network with
block-level masked convolutions is used and the block size is much smaller
(8x8). By using block-level masked convolutions, each block is processed using
reconstructed neighboring left and upper blocks both at the encoder and
decoder. Hence, the mutual information between adjacent blocks is exploited
during compression and each block is reconstructed using neighboring blocks,
resolving the need for explicit intra prediction and deblocking neural
networks. Since the explored system is a closed loop system, a special
optimization procedure, the asymptotic closed loop design, is used with
standard stochastic gradient descent based training. The experimental results
indicate competitive image compression performance.

###### keywords:

Image compression , Deep learning

##  1 Introduction

The111Codes are shared at https://github.com/metu-kamisli/Learned-block-based-
image-compression success of deep neural networks (DNN) in computer vision
applications has attracted great attention [1] and many problems in image
processing have been revisited with the tools that DNN provide [2, 3, 4, 5,
6]. One such problem is image compression. Artificial neural network (ANN)
based image compression approaches flourished quickly and their compression
performance have caught up in recent years with that of traditional image and
video compression systems [7].

Early traditional lossy image compression methods like JPEG [8] and JPEG2000
[9] use a linear decorrelating transform to map pixel values to a transform
domain, apply scalar quantization and process the quantized transform
coefficients with the inverse transform to reconstruct the compressed image
[10]. More recent image compression systems, such as Webp, BPG and AVIF, are
based on intra-frame compression tools of video compression standards, such as
HEVC [11] or AV1 [12]. These systems operate with a block-by-block processing
approach and use spatial prediction prior to the linear decorrelating
transform and are designed as multiple-mode compression systems. In
particular, there are multiple spatial prediction modes, multiple linear
transforms and multiple block sizes along with other multiple-mode tools.
These multiple-mode tools allow the compression algorithm to adapt to varying
spatial characteristics in images and improve the compression performance.
Note that although the processing in each spatial prediction mode or transform
is linear, the collection of all multiple-mode tools and a mode decision
algorithm at the encoder constitute overall a nonlinear compression system.

ANN based compression systems, on the other hand, typically do not have
multiple modes of operation but are inherently nonlinear (or adapt to image
content) as they use ANN. They mostly use an auto-encoder based architecture
which uses convolutional neural networks (CNN) as nonlinear analysis and
synthesis transforms [13, 14]. The nonlinear analysis transform maps pixels
values to a latent domain, where scalar quantization is performed. The
quantized latent domain variables are then mapped back with the nonlinear
synthesis transform to the pixel domain to reconstruct the compressed image.
These nonlinear transforms are learned by minimizing a rate-distortion cost
over a large set of training images. To have efficient entropy coding, the
probability distribution of the quantized latent domain variables is also
learned with an ANN along with the analysis and synthesis transforms. The
entire compression system can be trained end-to-end to minimize a rate-
distortion cost function consisting of any differential distortion metric,
such mean-squared-error (MSE) or multi-scale structural similarity index
measure (MS-SSIM) [15, 16], and the estimated entropy (rate) of the quantized
latent variables via a differential relaxation method to avoid the non-
differential quantization operation [13, 14].

Most ANN based compression systems use CNN to form the nonlinear analysis and
synthesis transforms (and entropy models). CNN operate on the entire input
image [17]. On the other hand, traditional state-of-the-art image and video
compression methods use a block-by-block processing approach for various
reasons [18, 19, 20]. The primary reason is that the motion-compensated inter-
frame coding method that dominates the video compression system design is
inherently a block-based system and other aspects of video compression systems
are also designed to operate block-by-block. In addition, switching between
temporal and spatial (i.e inter and intra) prediction on a block-by-block
basis is easy to achieve and improves compression performance.

Very recently, few papers on learned image compression with block based
approaches have also appeared [21, 22, 23]. They split the input image into
large blocks (256x256 [21], 128x128 [21, 23], 64x64 [21, 22], 32x32 [22]) and
use the described auto-encoder architecture [24] on each block but introduce
additional ANN that perform spatial/intra prediction [21, 22] (to exploit
mutual information between adjacent blocks) and post-processing/deblocking
[21, 22, 23] functions (to combine the reconstructed image blocks smoothly
with each other).

This paper explores an alternative learned block-based image compression
approach in which neither an explicit intra/spatial prediction ANN nor a
deblocking ANN is used. A single auto-encoder ANN with block-level masked
convolutions is used and the block size is much smaller (8x8). By using block-
level masked convolutions in the auto-encoder, each block is processed using
neighboring left and upper blocks both at the encoder and decoder. Hence, the
mutual information between adjacent blocks is exploited during compression and
there is no need for an explicit intra/spatial prediction ANN. At the decoder,
each block is reconstructed using neighboring blocks and there is no need for
an explicit post-processing/deblocking ANN.

One advantage of the approach explored in this paper, relative to other block
based learned compression approaches in the literature, is that it has a
simpler and cleaner architecture consisting of one auto-encoder based ANN.
Another potential advantage is that it can be easily extended to video (i.e.
inter-frame) compression [18] by incorporating motion-compensated reference
blocks instead of the spatial neighbor blocks. Hence, as a byproduct,
switching between inter-frame and intra-frame compression can be achieved
simply and the small block-size of 8x8 is likely to allow for better
compression-efficiency relative to switching between the compression methods
in approaches that use larger block sizes.

The experimental results indicate that the explored learned block-based image
compression approach in this paper provides competitive compression
performance. Compared to HEVC (BPG) intra coding, the performance is slightly
inferior at low bitrates but superior at higher bitrates. Compared to state-
of-the-art learned compression approaches which do not operate on a block-by-
block basis, the performance is slightly inferior or similar at low bitrates
but catches up or is superior at higher bitrates. More details are provided in
Section 5.

The remainder of the paper is organized as follows. Section 2 reviews related
work. Section 3 presents the explored learned block-based compression
approach. The end-to-end training of the explored system, including the
Asymptotic Closed Loop design and other details are discussed in Section 4.
Section 5 presents experimental results and comparisons with other work.
Finally, Section 6 concludes the paper.

##  2 Related Work

###  2.1 End-to-end Learned Image Compression Review

While research on learned image compression can be traced back many years [25,
26] the recent surge of interest started with the works of Balle et. al [13,
14, 24], and since then, many researchers used the compression architecture in
[24] as starting point and extended it in many ways.

The baseline ANN based compression model in [13] is shown in Figure 1 a). The
image $x$ is processed with the nonlinear analysis function $g_{a}(.)$, which
consists of 3 successive layers of CNN with downsampling factors 2 (or 4) and
nonlinear generalized divisive normalization (GDN) operations, to obtain the
latent representation $y$. The latent $y$ is quantized and the quantized
latent $\hat{y}$ is lossless coded with an arithmetic coder. (The $U|Q$ symbol
denotes quantization to nearest integer in compression/inference and addition
of uniform noise during training.) The decoder decodes the bitstream to obtain
$\hat{y}$ and processes it with the nonlinear synthesis function $g_{s}(.)$,
which also consists of 3 successive layers of CNN with an upsampling factor of
2 (or) and a nonlinear inverse of generalized divisive normalization (IGDN)
operation, to reconstruct the compressed image $\hat{x}$.

(a) Baseline model [13]

(b) Hyper-prior based model [14]

(c) Joint auto-regressive and hyper-prior based model [24]

Figure 1: Operational diagrams of learned compression models

The analysis and synthesis transforms as well as the probability distribution
of the quantized latent $\hat{y}$ are learned by minimizing a rate-distortion
cost $L$ over a set of training images:

$\displaystyle L$ $\displaystyle=R(\hat{y})+\lambda D(x,\hat{x})$ (1)
$\displaystyle=E[-log_{2}p_{\hat{y}}(\hat{y})]+\lambda\cdot E[d(x,\hat{x})]$
(2)

Here, $\lambda$ is a Lagrangian multiplier that controls the trade-off between
rate and distortion and $d(.,.)$ is any differentiable distortion metric such
as mean-squared-error (MSE) or multi-scale structural similarity index measure
(MS-SSIM) [16]. The rate term $E[-log_{2}p_{\hat{y}}(\hat{y})]$ is obtained
from the probability distribution $p_{\hat{y}}(\hat{y})$ of $\hat{y}$ assuming
it is factorized over the distribution of elements $\hat{y}_{i}$ as

$\displaystyle p_{\hat{y}}(\hat{y})=\prod_{i}p_{\hat{y}_{i}}(\hat{y}_{i})$ (3)

with the elements $\hat{y}_{i}$ in the same channel being identically
distributed [13, 14].

Since quantization operation has zero derivatives almost everywhere, the rate
term is relaxed with $E[-log_{2}p_{\tilde{y}}(\tilde{y})]$ where the
probability density function $p_{\tilde{y}}(\tilde{y})$ is similarly
factorized over the densities of the elements $\tilde{y}_{i}$, which are
obtained by adding noise $u_{i}$ uniformly distributed in
$[-\frac{1}{2},\frac{1}{2}]$ to the latent variables $y_{i}$:

$\displaystyle\tilde{y}_{i}=y_{i}+u_{i}$ (4)

This relaxation allows to use $p_{\tilde{y}_{i}}(\tilde{y}_{i})$ as a proxy
for $p_{\hat{y}_{i}}(\hat{y}_{i})$ during training since for integers $k$,
$p_{\tilde{y}_{i}}(k)=p_{\hat{y}_{i}}(k)$ i.e. the probability mass and
density functions are tied together. The density values
$p_{\tilde{y}_{i}}(\tilde{y}_{i})$ are obtained from an ANN that models/learns
the cumulative distribution function (CDF) $F_{y_{i}}(.)$ of $y_{i}$ based on
Equation (4) as [13, 14]:

$\displaystyle
p_{\tilde{y}_{i}}(\tilde{y}_{i})=F_{y_{i}}(\tilde{y}_{i}+\frac{1}{2})-F_{y_{i}}(\tilde{y}_{i}-\frac{1}{2})$
(5)

With this relaxation, the entire compression system can be trained end-to-end
to determine the parameters of ANN forming $g_{a}(.)$, $g_{s}(.)$ and
$F_{y}(.)$ by minimizing the rate-distortion cost in Equation (2) over a set
of training images.

The baseline ANN based compression model was improved in [14] by transmitting
a learned hyper-prior $\hat{z}$ (similar to side or forward information in
compression theory) that will allow the decoder to decode the latent $\hat{y}$
with a probability model that is spatially adapted to the coded image (see
Figure 1 b). The hyper-prior $z$ is generated from the latent $y$ using a
hyper-encoder $h_{a}(.)$ ANN, then quantized and lossless coded. The decoder
decodes the bitstream of the hyper-prior and then processes it with the hyper-
decoder $h_{s}(.)$ ANN to generate per-variable mean and standard deviations
($\mu_{i},\sigma_{i}$) for the conditional probability distribution
$p_{\hat{y}|\hat{z}}(\hat{y}|\hat{z})$ of the latent $\hat{y}$ conditioned on
the hyper-prior $\hat{z}$.

The rate-distortion cost to be minimized becomes

$\displaystyle L$ $\displaystyle=R(\hat{z})+R(\hat{y})+\lambda D(x,\hat{x})$
(6)
$\displaystyle=E[-log_{2}p_{\hat{z}}(\hat{z})]+E[-log_{2}p_{\hat{y}|\hat{z}}(\hat{y}|\hat{z})]+\lambda\cdot
E[d(x,\hat{x})]$ (7)

where the rate of the hyper-latent $E[-log_{2}p_{\hat{z}}(\hat{z})]$ is now
relaxed and modeled as the latent in the base model. The rate of the latent
$E[-log_{2}p_{\hat{y}|\hat{z}}(\hat{y}|\hat{z})]$ is now relaxed via
$E[-log_{2}p_{\tilde{y}|\tilde{z}}(\tilde{y}|\tilde{z})]$ during training
where the conditional density $p_{\tilde{y}|\tilde{z}}(\tilde{y}|\tilde{z})$
is factorized over the densities
$p_{\tilde{y}_{i}|\tilde{z}}(\tilde{y}_{i}|\tilde{z})$ of the elements
$\tilde{y}_{i}$, which are obtained by adding noise uniformly distributed in
$[-\frac{1}{2},\frac{1}{2}]$ to the latent variables $y_{i}$. The conditional
density values $p_{\tilde{y}_{i}|\tilde{z}}(\tilde{y}_{i}|\tilde{z})$ are
obtained using conditional CDFs $F_{y_{i}|\tilde{z}}(.)$, similar to Equation
(5), and $F_{y_{i}|\tilde{z}}(.)$ are modeled with the Gaussian CDF whose mean
and variance is given by $h_{s}(.)$, i.e.

$\displaystyle
p_{y_{i}|\tilde{z}}(y_{i}|\tilde{z})\sim\mathcal{N}(\mu_{i},\sigma_{i}).$ (8)

The hyper-prior based compression model was also improved in [24] by forming
the conditional density of the latent variables not only from the hyper-prior
but also from the causal context $\tilde{y}_{<i}$ of each latent variable:

$\displaystyle
p_{y_{i}|\tilde{z},\tilde{y}_{<i}}(y_{i}|\tilde{z},\tilde{y}_{<i})\sim\mathcal{N}(\mu_{i},\sigma_{i})$
(9)

This joint (auto-regressive and hierarchical) priors model (see Figure 1 c)
further improves the compression performance, however, due to the spatial
dependency introduced by the causal context $\tilde{y}_{<i}$ the decoding
process can not by fully parallelized anymore [24].

Many researchers have used the joint (auto-regressive and hierarchical) priors
model compression architecture in [24] and modified some of its components in
several ways to obtain better compression performance [27, 28, 29]. For
example, Cheng et al. [27] use attention modules and residual connections in
the auto encoder architecture and replace the Gaussian distribution model in
Equation (9) with a Gaussian mixture model. Akbari et al. [28] replace the
standard CNN layers with generalized octave convolution layers to work with
multi-frequency feature maps. Yilmaz et al. [29] replace the CNN and GDN
layers in the auto encoder with self-organized operational neural network
layers.

###  2.2 Block Based Learned Image Compression

Recently, papers on learned image compression with block based approaches have
also appeared [21, 22, 23]. They split the input image into large blocks
(256x256 [21], 128x128 [21, 23], 64x64 [21, 22], 32x32 [22]) and use the above
summarized joint model architecture in [24] on each block but introduce
additional ANN that perform spatial/intra prediction and post-
processing/deblocking functions. The spatial/intra prediction ANN are used to
exploit the correlation from the previously reconstructed left and upper
blocks to improve the compression of the current block [21, 22]. The post-
processing/deblocking ANN are used to remove blocking artifacts that occur at
the block boundaries mostly due to independent processing of the blocks [21,
22, 23].

##  3 Explored Block Based Image Compression Approach

This paper explores an alternative learned block-based image compression
approach in which neither an explicit intra/spatial prediction ANN nor a
deblocking ANN is used. A single auto-encoder based ANN with block-level
masked convolutions is used and the block size is much smaller (8x8). By using
block-level masked convolutions in the auto-encoder ANN, each block is
processed using neighboring left and upper blocks both at the encoder and
decoder. Hence, the mutual information between adjacent blocks is exploited
during compression and there is no need for an explicit intra/spatial
prediction ANN. At the decoder, each block is reconstructed using neighboring
blocks and there is no need for an explicit post-processing/deblocking ANN.

The details of the system architecture and its operation during
compression/inference is discussed next. Section 4 discusses the details of
training the system.

###  3.1 ANN Architecture and Operation During Compression/Inference

The operational diagram of the system for compression/inference is shown in
Figure 3. First, the image $\mathbf{x}$ is processed with the $B2C$ (block to
channel) subsystem, which scans the pixels inside a $3$x$B$x$B$ block and
arranges them along the channel dimension, i.e. an (RGB) image of size
$3$x$H$x$W$ is converted to a tensor of size
$3B^{2}$x$\frac{H}{B}$x$\frac{W}{B}$. This conversion allows to use the CNN
routines in common deep learning frameworks for block-level convolutions.

Figure 2: Operational diagram for compression/inference. The encoder contains
all subsystems (except AD) in the figure while the decoder contains the
subsystems inside the red dashed lines. Both the encoder and decoder operate
closed loop in a block-by-block manner.

Figure 3: Fusion of original block $x_{j}$ and reconstructed blocks in the
causal neighborhood $\hat{x}_{<j}$ at the encoder with CNN in the first layer
of $t_{a}(.)$. Masked convolution [30] is used to use only the causal
neighborhood of, i.e. left and upper, reconstructed blocks.

The system uses a block-by-block processing (in raster scan order) approach as
traditional image and video compression systems. Hence, one block, shown with
$x_{j}$ in Figure 3, is processed through the entire system and its
reconstruction, denoted with $\hat{x}_{j}$, is stored in the decoded blocks
buffer to be used for the processing of the next blocks.

At the encoder, the processing is performed as follows. In the first step, the
CNN processing result of an original block $x_{j}$ (of size $3B^{2}$x$1$x$1$)
is added to the CNN processing result of its previously reconstructed causal
neighborhood of blocks $\hat{x}_{<j}$ (of size $3B^{2}$x$3$x$3$) as shown in
Figure 3. The result is then further processed with GDN (Generalized Divisive
Normalization [31, 13]) activation functions and CNN layers to produce the
latent variables block $y_{j}$. This processing at the encoder is denoted with
$t_{a}(.)$ in Figure 3 and can be seen as a nonlinear predictive analysis
transform, playing the role of intra prediction and transforming the
prediction residual in traditional image compression.

$\displaystyle y_{j}=t_{a}(x_{j},\hat{x}_{<j})$ (10)

The latent block is quantized (by rounding each element to the nearest
integer) to produce the quantized latent block $\hat{y}_{j}$, which is
lossless compressed into a bitstream by an arithmetic encoder (AE) as shown in
Figure 3. The probability distribution that the AE uses is determined by the
parameters $\mu_{j}$ and $\sigma_{j}$ that are produced by an ANN, denoted by
$n(.)$ in Figure 3, which processes the previously reconstructed causal
neighborhood of blocks $\hat{x}_{<j}$. This ANN consists of CNN with masked
convolutions (to capture only the causal context) or regular convolutions
(with kernel size 1x1) and LeakyReLU activation functions. Note that this ANN
based subsystem $n(.)$ is common to both the encoder and decoder as they need
to generate the same probability distribution
$p_{\hat{y}_{j}|\hat{x}_{<j}}(\hat{y}_{j}|\hat{x}_{<j})$ that both the AE and
the arithmetic decoder (AD) must use. These processing steps can be summarized
as follows:

$\displaystyle\hat{y}_{j}=Q(y_{j})$ (11)
$\displaystyle(\mu_{j},\sigma_{j})=n(\hat{x}_{<j})$ (12) $\displaystyle
p_{\hat{y}_{j}|\hat{x}_{<j}}(\hat{y}_{j}|\hat{x}_{<j})\sim(\mu_{j},\sigma_{j}).$
(13)

The quantized latent variable block $\hat{y}_{j}$ is processed together with
previously reconstructed causal neighborhood of blocks $\hat{x}_{<j}$ via the
reconstruction ANN, denoted with $t_{s}(.)$ in Figure 3, to obtain the
reconstruction block $\hat{x}_{j}$.

$\displaystyle\hat{x}_{j}=t_{s}(\hat{y}_{j},\hat{x}_{<j})$ (14)

In the first step of the reconstruction ANN, the CNN processing result of the
quantized latent block $\hat{y}_{j}$ (of size $M$x$1$x$1$) is added to the CNN
processing result of the previously reconstructed causal neighborhood of
blocks $\hat{x}_{<j}$ (of size $3B^{2}$x$3$x$3$), similar to the first step of
processing at the encoder. The result is then further processed with IGDN
(Inverse Generalized Divisive Normalization [31, 13]) and CNN layers to
produce the reconstructed block $\hat{x}_{j}$. This reconstruction process at
the decoder, denoted with $t_{s}(.)$ in Figure 3, can be seen as a nonlinear
predictive synthesis transform, playing the role of the following operations
in traditional image compression: inverse transforming the quantized transform
coefficients, performing intra prediction and adding their results. The
reconstructed block $\hat{x}_{j}$ is stored in reconstructed blocks buffer and
the compression/decompression of the next block can start.

The decoder works in a similar manner with block-by-block processing and
consists of the sub-systems in the red dashed box in Figure 3. First, the ANN
$n(.)$ processes causal context $\hat{x}_{<j}$ to obtain probability
distribution parameters $\mu_{j}$ and $\sigma_{j}$ that are used by the AD to
obtain quantized latent block $\hat{y}_{j}$. Then $\hat{y}_{j}$ and causal
context $\hat{x}_{<j}$ are processed jointly by ANN $t_{s}(.)$ to reconstruct
the block $\hat{x}_{j}$.

Finally, note that the block-by-block processing prevents CNN computations to
be performed on the entire images or feature maps, i.e. one block has to be
processed through the entire system including all CNN and then the next block
processing can start. Such computations are also present in some other
learning based approaches [24] and may not fully capitalize on the
computational resources of standard GPU hardware or software. More details are
discussed in Section 5.3.

##  4 End-to-End Training with Asymptotic Closed Loop Design

###  4.1 Asymptotic Closed Loop (ACL) Optimization

The explored block based learned image compression system is a closed loop
system. In other words, the compression result $\hat{x}_{j}$ of one block is
used for the compression of the next blocks. This poses a fundamental problem
for the end-to-end training of the system since initially the compression
system with randomly initialized ANN weights cannot produce proper
reconstructed blocks $\hat{x}_{<j}$ using which the system can be trained.

A naive approach would be to use the original (uncompressed) blocks in the
causal neighborhood $x_{<j}$ instead of the reconstructed blocks
$\hat{x}_{<j}$ during training. This is termed the open-loop optimization
procedure in [32] and produces inferior compression performance when such
trained system is used in the actual compression/inference, which operates
closed loop. In particular, the compression/inference performance is
catastrophic at low bit-rates where the reconstructed blocks and original
blocks differ significantly.

We adopt the Asymptotic Closed Loop (ACL) design procedure from [32], which
was proposed for predictive vector quantizer design in signal compression and
was also applied to predictor design in video compression [33]. The ACL
procedure provides stable open-loop training while ultimately optimizing the
compression system for closed loop operation.

The ACL design algorithm proceeds in iterations where in each iteration the
compression system is trained open loop (with standard stochastic gradient
descent) until convergence. The operational diagram of the compression system
in open loop training is shown in Figure 5. In the $k^{th}$ ACL iteration, the
reconstructions $\mathbf{\hat{x}}^{(k)}$ that the system uses in the training
are not taken from the output of the system but from a secondary training set,
which was produced with the converged system at the end of the $(k-1)^{th}$
ACL iteration. In the initial iteration ($k=0$), the system uses the original
images also for the reconstructions, i.e. $\mathbf{\hat{x}}^{(0)}=\mathbf{x}$.

Figure 4: Operational diagram during training. System operates open loop.

Figure 5: Sample validation set RD Cost curves of open loop and closed loop
operations over SGD epochs

Another difference between the open loop operation of the system in the ACL
iterations and its closed loop operation during actual compression/inference
is that the block-by-block processing is not necessary in the open loop
operation. In other words, entire images can be processed through the system
at once, which speeds up training. Hence, the original blocks buffer and
decoded blocks buffer in Figure 3 are not present in Figure 5. A summary of
the used ACL design algorithm is given in Algorithm 1.

With every ACL iteration, the system is trained with reconstructions
$\mathbf{\hat{x}}^{(k)}$ that are closer to its produced reconstructions
$\mathbf{\tilde{x}}$ and the closed loop compression performance improves. We
observe convergence after several ACL iterations. Sample curves of the
validation RD costs of the open loop and closed loop operations of the system
are shown in Figure 5. The curves are plotted for many SGD training epochs,
which cover several ACL iterations. The end of each ACL iteration is shown
with vertical dashed lines. It can be seen that in the initial ($k=0$)
iteration, there is a big gap between the open and closed loop RD costs, which
reduces with following ACL iterations.

Algorithm 1 Asymptotic Closed Loop (ACL) training

Define training and validation images sets :
$\\{\mathbf{x}\\},\\{\mathbf{x_{VAL}}\\}$

Initialize reconstructed images set :
$\\{\mathbf{\hat{x}}^{(0)}\\}=\\{\mathbf{x}\\}$

Initialize compression system parameters : $\Phi^{(0)}$ $\leftarrow$ random

for ACL iteration k $=0$ to $\infty$ do

$\mathbf{1.}$ Train compression system open loop with SGD until convergence :
$\Phi^{(k+1)}\leftarrow
train(\Phi^{(k)},\\{\mathbf{x}\\},\\{\mathbf{\hat{x}}^{(k)}\\})$

$\mathbf{2.}$ Get compression systems closed loop validation cost :
$RDCost^{(k+1)}_{VAL}\leftarrow\mathbf{\Phi_{CL}}^{(k+1)}(\\{\mathbf{x_{VAL}}\\})$

$\mathbf{3.}$ Break if $RDCost_{VAL}$ values converge

$\mathbf{4.}$ Generate reconstructed images set for next iteration :
$\\{\mathbf{\hat{x}}^{(k+1)}\\}\leftarrow\mathbf{\Phi_{OL}}^{(k+1)}(\\{\mathbf{x}\\},\\{\mathbf{\hat{x}}^{(k)}\\})$

end for

###  4.2 Rate-Distortion Cost for SGD Optimization of System Parameters

In each ACL iteration, the compression system parameters are optimized to
minimize the following rate-distortion (RD) cost $L^{(k)}$ over a set of
training and reconstruction images, $\\{\mathbf{x}\\}$ and
$\\{\mathbf{x^{(k)}}\\}$ respectively.

$\displaystyle L^{(k)}$ $\displaystyle=R(\hat{y})+\lambda D(x,\tilde{x})$ (15)
$\displaystyle=E[-log_{2}p_{\hat{y}|\hat{x}^{(k)}}(\hat{y}|\hat{x}^{(k)})]+\lambda\cdot
E[d(x,\tilde{x})]$ (16)

Here, $\lambda$ is a Lagrangian multiplier that controls the trade-off between
rate and distortion and $d(.,.)$ is a differentiable distortion metric such as
mean-squared-error (MSE) or multi-scale structural similarity index measure
(MS-SSIM). The rate $R(\hat{y})$ of the latent is obtained from the
conditional distribution $p_{\hat{y}|\hat{x}^{(k)}}(\hat{y}|\hat{x}^{(k)})$ of
the latent $\hat{y}$ conditioned on the causal neighborhood of the
reconstructions in $\hat{x}^{(k)}$ and there is no hyper-latent transmission
as in [24] or the block-based approaches using its main architecture [21, 22,
23].

Similar to [14, 24], the rate of the latent is relaxed via
$E[-log_{2}p_{\tilde{y}|\hat{x}^{(k)}}(\tilde{y}|\hat{x}^{(k)})]$ where the
conditional density $p_{\tilde{y}|\hat{x}^{(k)}}(\tilde{y}|\hat{x}^{(k)})$ is
factorized over the densities
$p_{\tilde{y}_{i}|\hat{x}^{(k)}}(\tilde{y}_{i}|\hat{x}^{(k)})$ of the elements
$\tilde{y}_{i}$, which are obtained by adding noise uniformly distributed in
$[-\frac{1}{2},\frac{1}{2}]$ to the latent variables $y_{i}$. The conditional
density values $p_{\tilde{y}_{i}|\hat{x}^{(k)}}(\tilde{y}_{i}|\hat{x}^{(k)})$
are obtained using conditional CDFs $F_{y_{i}|\hat{x}^{(k)}}(.)$, similar to
Equation (5), and $F_{y_{i}|\hat{x}^{(k)}}(.)$ are modeled with the Gaussian
CDF whose mean and variance is given by the ANN $n(.)$ that processes the
causal context of the reconstruction $\hat{x}^{(k)}$ (see Figure 5), i.e.

$\displaystyle
p_{y_{i}|\hat{x}^{(k)}}(y_{i}|\hat{x}^{(k)})\sim\mathcal{N}(\mu_{i},\sigma_{i}).$
(17)

With this relaxation, the entire compression system can be trained end-to-end
to determine the parameters of ANN forming $t_{a}(.)$, $t_{s}(.)$ and $n(.)$
by minimizing the rate-distortion cost in Equation (16) over a set of training
images.

###  4.3 Fine Tuning of System Parameters

Upon convergence of the ACL iterations, it is observed that the MSE of the
reconstructions in successive iterations are very similar ( i.e.
$d(\mathbf{x_{VAL}},\mathbf{\hat{x}_{VAL}^{(k)}})\sim
d(\mathbf{x_{VAL}},\mathbf{\hat{x}_{VAL}^{(k+1)}})$) but the closed loop
validation RD cost

$RDCost^{(k+1)}_{VAL,CL}\leftarrow\mathbf{\Phi_{CL}}^{(k+1)}(\\{\mathbf{x_{VAL}}\\})$

can still be higher than the open loop validation RD cost

$RDCost^{(k+1)}_{VAL,OL}\leftarrow\mathbf{\Phi_{OL}}^{(k+1)}(\\{\mathbf{x_{VAL}}\\},\\{\mathbf{\hat{x}_{VAL}}^{(k)}\\}).$

This small discrepancy hints that although at convergence of ACL iterations,
successive iterations provide reconstructions with similar distortions, their
statistics may still be slightly different since the closed loop operation
does give slightly inferior RD cost than open loop operation. Thus, a fine-
tuning stage was added where the reconstructions produced by the system
($\tilde{x}$ in Figure 5) were fed back to the system in place of
$\hat{x}^{(k)}$ and the cost to be minimized was formed as average of the two
RD costs from the first and second runs of the system as follows

$\displaystyle L^{(k)}=$
$\displaystyle\frac{1}{2}(E[-log_{2}p_{\tilde{y}_{0}|\hat{x}^{(k)}}(\tilde{y}_{0}|\hat{x}^{(k)})]+\lambda\cdot
E[d(x,\tilde{x}_{0})])~{}+$
$\displaystyle\frac{1}{2}(E[-log_{2}p_{\tilde{y}_{1}|\tilde{x}_{0}}(\tilde{y}_{1}|\tilde{x}_{0})]+\lambda\cdot
E[d(x,\tilde{x}_{1})]).$ (18)

Here, $\tilde{x}_{0}$ and $\tilde{x}_{1}$ are the reconstructions and
$\tilde{y}_{0}$ and $\tilde{y}_{1}$ are the latents from the first and second
runs of the system.

In summary, the system parameters were fine tuned for a few more ACL
iterations following Algorithm 1 except that in Step 1 the RD cost that was
minimized is the one in Equation (18). This improved the RD cost of the system
in closed loop operation. In Figure 5, the fine-tuning ACL iterations start
with the $5^{th}$ iteration after epoch 252.

###  4.4 System Hyper Parameters

As in many learned image compression systems, a new set of system parameters
is learned for each different rate-distortion (RD) trade-off by changing the
$\lambda$ parameter in Equation (16) [13, 14, 24]. The explored system is
trained for 8 different $\lambda$ values, corresponding to those in [34], to
cover a sufficiently large range of reconstructed picture qualities.

Two sets of hyper parameters were used for the explored system, which mainly
differ in the number of channels used for the CNN in the system. One set of
hyper parameters was used for the four smaller $\lambda$ values, and another
set with larger learning capacity for the four larger $\lambda$ values where
the system needs to reconstruct higher fidelity images. The hyper parameters
are summarized in Table 1.

Table 1: Two sets of hyper parameters (HPS) to define the system | HPS1 | HPS2
---|---|---
N | 768 | 1152
M | 96 | 128
$N_{1}$ | $\frac{7}{8}N$
$N_{2}$ | $\frac{6}{8}N$
$M_{1}$ | $\frac{6}{4}N$
$M_{2}$ | $\frac{5}{4}N$
$M_{3}$ | $\frac{4}{4}N$
$K_{2}$x$K_{2}$ | 1x1 | 3x3

###  4.5 Block Size

A very important hyper parameter of the explored learned block based image
compression system is the block size. The system supports a single block size,
and a block size of 8x8 is used in the experimental results that are presented
in the next section. This block size provides competitive overall compression
results and reasonable computational/model complexity. A larger block size,
such as 16x16, requires higher number of channels for the CNN to provide
efficient compression and increases computational/model complexity
significantly. A smaller block size, such as 4x4, fails to provide competitive
compression results at lower bitrates since compression of large smooth
regions becomes quite inefficient with 4x4 block size. In summary, the fixed
block size of 8x8 is used for the main results presented in the next section,
however some results and comments for 16x16 and 4x4 block size systems are
also provided.

##  5 Experimental Results

###  5.1 Training Set and Further Details

For training, the CLIC professional and mobile datasets [35] were used with
randomly cropped patches of 256x256 and a batch size of 4 for each SGD
iteration. The Adam optimizer [36] was used with an initial learning rate of
$10^{-4}$ together with a learning rate scheduler that reduces the learning
rate by a factor of $0.8$ (up to a minimum of $2$x$10^{-5}$) when the cost
plateaus.

The entire system and associated training algorithms were implemented with
PyTorch and numerically stable Gaussian CDF implementations from [34] and the
entropy coder/decoder of [37] with [34]’s implementation were utilized. The
codes to repeat our results are shared at [38].

###  5.2 Compression Results

To test the compression performance of the explored learned block based image
compression method, the Kodak [39] and the Tecnick datasets [40] were used.
The Kodak dataset contains a set of 24 uncompressed images with 512x768
resolution (and 4:4:4 RGB color sub sampling) commonly used to evaluate image
compression methods. The Tecnick dataset contains 100 images of 1200x1200 (RGB
4:4:4) resolution.

The compression performance results are shown in Figure 6. The peak-signal-to-
noise-ratio (PSNR) and bits-per-pixel (BPP) values in the figure were
calculated by averaging the PSNR and BPP values of all images in each dataset.
The compression performance results of compared methods are taken from the
publicly available results reported at [41] which are prepared by the authors
of [14, 24].

(a) Kodak dataset [39]

(b) Tecnick dataset [40]

Figure 6: Compression performance results on Kodak and Tecnick datasets.

It can be seen from Figure 6 that the explored learned block based image
compression method provides competitive compression performance. In
particular, it provides better compression performance than traditional
compression systems JPEG and JPEG2000. Compared to BPG, the explored method
provides better compression performance at higher bitrates but slightly
inferior performance at lower bitrates. Compared to the learned compression
method of [14] (Balle et al. 2018), the explored method is on par at the low
bitrate end but provides superior results at higher bitrates. Compared to the
learned compression method of [24] (Minnen et al. 2018), the explored method
provides inferior, though competitive, results at lower bitrates but catches
up at the high bitrate end. Some visual results are provided in Figure 7. More
visual results are available at [38].

Figure 7: Reconstructed pictures at around 0.22 bpp (Zoom in the pdf file for
details.) From left to right with (PSNR, bpp): Original, JPEG2000 (26.84,
0.223), BPG (28.72, 0.222), Balle et al. [14] (28.59, 0.224), Minnen et al.
[24] (28.96, 0.208), Our method (28.59, 0.223).

There are two aspects which can hinder the performance of the explored learned
block based image compression method. The first one is that the explored
method is a block based method which processes blocks in a causal order and
therefore does not allow (for the compression of a block) fusing information
from the right and lower blocks. This is unlike many learned compression
methods [14, 24]. However, as mentioned before, block based compression can
also have desirable features, such as easy combination with video compression
systems [18, 19] (as discussed in Section 1) or internal memory requirement
advantages in hardware implementation [42, 43].

The second aspect is that the explored system operates at a single block size
of 8x8. However, in large smooth regions, the system is likely to provide
better compression with a larger block size and in detailed texture regions,
it is likely to do better with a smaller block size. To investigate this, the
explored system was trained with a (fixed) block size of 16x16 for very low
bitrates (i.e. lowest $\lambda$ value) and with a (fixed) block size of 4x4
for high bitrates (i.e. highest $\lambda$ value). The compression results are
shown in Figure 8. With a block size of 16x16, the explored system provides
slightly better compression at low bitrates. With a block size of 4x4, the
explored system provides similar (but not better) compression performance at
high bitrates. These results indicate that supporting multiple block sizes (as
traditional compression methods) in the explored learned block based system is
likely to increase compression performance, however, it requires careful joint
training of multiple block size systems with ACL design and is part of our
future research studies.

(a) 16x16 block size

(b) 4x4 block size

Figure 8: Compression performance results on Kodak dataset with different
block sizes of the explored system.

###  5.3 Computational Complexity Analysis

A computational complexity analysis is provided in Table 2 using the required
number of parameters (of the ANN) and the encoding and decoding times. The
encoding and decoding times of our method and those of Balle et al. [14] and
Minnen et al. [24] (with the implementations in [34]), all learning based
methods, were measured with our computer and run on the GPU222Nvidia RTX 3060
Ti, while those of BPG, JPEG and JPEG2000, all traditional methods, were taken
from [34] and run on the CPU. The results in Table 2 should taken as a rough
estimate of computational complexity. A thorough analysis depends on many
factors and is out of the scope of this paper.

Table 2: Number of parameters and encoding/decoding times | Number of | Run-time
---|---|---
| Parameters | (Enc., Dec.)
Our method (HPS1) | 8.6 M | 27.4 , 37.2 s
Our method (HPS2) | 24.7 M | 29.2 , 41.4 s
Balle et al. 2018 (HPS1) | 4.9 M | 0.06 , 0.03 s
Balle at al. 2018 (HPS2) | 11.3 M | 0.09 , 0.04 s
Minnen et al. 2018 (HPS1) | 12.6 M | 4.6 , 8.8 s
Minnen at al. 2018 (HPS2) | 21.8 M | 4.0 , 8.1 s
BPG | NA | 3.7 , 0.07 s
JPEG2000 | NA | 0.5 , 0.5 s

The number of parameters of all learning based methods in Table 2 are
comparable, however, the encoding and decoding times can differ by several
orders of magnitude. Traditional methods do not use CNN or other
computationally heavy algorithms and thus are typically fast on and optimized
for the CPU. Learning based methods require typically orders of magnitude more
computations and thus are preferred to run on the GPU. Balle et al.’s method
[14] has very low encoding and decoding times since it can fully exploit the
GPU’s parallel computation resources. On the other hand, Minnen et al.’s
method [24] can not fully exploit the GPU’s parallel computation resources
since it requires some serial computations (for the auto-regressive parts of
their model) which can not be efficiently computed with standard GPU hardware
and software. Similarly, our method requires a block-by-block processing
approach, which means CNN computations can not be performed on the entire
image at once, which again can not be efficiently computed with standard GPU
hardware and software.

We note that the total number of parameters and therefore the total number
computations of Minnen et al.’s method and our method are comparable to those
of Balle’s method. This indicates that other parallelization methods such as
wavefront parallel processing, well-known in traditional video compression
[44], or better optimized GPU hardware/software may enable also faster
encoding and decoding times for these methods.

##  6 Conclusions

This paper explored an alternative learned block-based image compression
approach in which, unlike other block based learned compression approaches in
the literature, neither an explicit intra prediction nor a deblocking network
neural network is used. A single auto-encoder neural network with block-level
masked convolutions is used and the block size is much smaller (8x8). Since
the explored system is a closed loop system, a special optimization procedure,
the asymptotic closed loop design, was used with standard stochastic gradient
descent based training. Experimental results indicate that the explored
learned block based approach provides competitive image compression
performance. Extending the approach to support multiple block sizes is likely
to further increase the compression performance, which is part of our future
research directions together with extension to inter-frame compression.

## References

  * [1] Y. Guo, Y. Liu, A. Oerlemans, S. Lao, S. Wu, M. S. Lew, Deep learning for visual understanding: A review, Neurocomputing 187 (2016) 27–48.
  * [2] C. Tian, Y. Xu, L. Fei, K. Yan, Deep learning for image denoising: A survey, in: International Conference on Genetic and Evolutionary Computing, Springer, 2018, pp. 563–572.
  * [3] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, C. I. Sánchez, A survey on deep learning in medical image analysis, Medical image analysis 42 (2017) 60–88.
  * [4] Z. Wang, J. Chen, S. C. Hoi, Deep learning for image super-resolution: A survey, IEEE transactions on pattern analysis and machine intelligence 43 (10) (2020) 3365–3387.
  * [5] S. Minaee, Y. Y. Boykov, F. Porikli, A. J. Plaza, N. Kehtarnavaz, D. Terzopoulos, Image segmentation using deep learning: A survey, IEEE transactions on pattern analysis and machine intelligence.
  * [6] S. Ma, X. Zhang, C. Jia, Z. Zhao, S. Wang, S. Wang, Image and video compression with neural networks: A review, IEEE Transactions on Circuits and Systems for Video Technology 30 (6) (2019) 1683–1698.
  * [7] Y. Hu, W. Yang, Z. Ma, J. Liu, Learning end-to-end lossy image compression: A benchmark, IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * [8] G. K. Wallace, The jpeg still picture compression standard, IEEE transactions on consumer electronics 38 (1) (1992) xviii–xxxiv.
  * [9] M. Rabbani, R. Joshi, An overview of the jpeg 2000 still image compression standard, Signal processing: Image communication 17 (1) (2002) 3–48.
  * [10] V. K. Goyal, Theoretical foundations of transform coding, IEEE Signal Processing Magazine 18 (5) (2001) 9–21.
  * [11] J. Lainema, F. Bossen, W.-J. Han, J. Min, K. Ugur, Intra coding of the hevc standard, IEEE transactions on circuits and systems for video technology 22 (12) (2012) 1792–1801.
  * [12] Y. Chen, D. Murherjee, J. Han, A. Grange, Y. Xu, Z. Liu, S. Parker, C. Chen, H. Su, U. Joshi, et al., An overview of core coding tools in the av1 video codec, in: 2018 Picture Coding Symposium (PCS), IEEE, 2018, pp. 41–45.
  * [13] J. Ballé, V. Laparra, E. P. Simoncelli, End-to-end optimized image compression, arXiv preprint arXiv:1611.01704.
  * [14] J. Ballé, D. Minnen, S. Singh, S. J. Hwang, N. Johnston, Variational image compression with a scale hyperprior, arXiv preprint arXiv:1802.01436.
  * [15] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE transactions on image processing 13 (4) (2004) 600–612.
  * [16] Z. Wang, E. P. Simoncelli, A. C. Bovik, Multiscale structural similarity for image quality assessment, in: The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003, Vol. 2, Ieee, 2003, pp. 1398–1402.
  * [17] J. Gu, Z. Wang, J. Kuen, L. Ma, A. Shahroudy, B. Shuai, T. Liu, X. Wang, G. Wang, J. Cai, T. Chen, Recent advances in convolutional neural networks, Pattern Recognition 77 (2018) 354–377.  doi:https://doi.org/10.1016/j.patcog.2017.10.013.
URL https://www.sciencedirect.com/science/article/pii/S0031320317304120

  * [18] D. J. Le Gall, The mpeg video compression algorithm, Signal Processing: Image Communication 4 (2) (1992) 129–140.  doi:https://doi.org/10.1016/0923-5965(92)90019-C.
URL https://www.sciencedirect.com/science/article/pii/092359659290019C

  * [19] T. Wiegand, G. Sullivan, G. Bjontegaard, A. Luthra, Overview of the H.264/AVC video coding standard, Circuits and Systems for Video Technology, IEEE Transactions on 13 (7) (2003) 560–576.  doi:10.1109/TCSVT.2003.815165.
  * [20] G. Sullivan, J. Ohm, W.-J. Han, T. Wiegand, Overview of the High Efficiency Video coding (HEVC) standard, Circuits and Systems for Video Technology, IEEE Transactions on 22 (12) (2012) 1649–1668.  doi:10.1109/TCSVT.2012.2221191.
  * [21] Y. Wu, X. Li, Z. Zhang, X. Jin, Z. Chen, Learned block-based hybrid image compression, IEEE Transactions on Circuits and Systems for Video Technology (2021) 1–1doi:10.1109/TCSVT.2021.3119660.
  * [22] Z. Yuan, H. Liu, D. Mukherjee, B. Adsumilli, Y. Wang, Block-based learned image coding with convolutional autoencoder and intra-prediction aided entropy coding, in: 2021 Picture Coding Symposium (PCS), IEEE, 2021, pp. 1–5.
  * [23] Z. Zhao, C. Jia, S. Wang, S. Ma, J. Yang, Learned image compression using adaptive block-wise encoding and reconstruction network, in: 2021 IEEE International Symposium on Circuits and Systems (ISCAS), IEEE, 2021, pp. 1–5.
  * [24] D. Minnen, J. Ballé, G. D. Toderici, Joint autoregressive and hierarchical priors for learned image compression, Advances in neural information processing systems 31.
  * [25] J. Jiang, Image compression with neural networks–a survey, Signal processing: image Communication 14 (9) (1999) 737–760.
  * [26] G. Toderici, S. M. O’Malley, S. J. Hwang, D. Vincent, D. Minnen, S. Baluja, M. Covell, R. Sukthankar, Variable rate image compression with recurrent neural networks, arXiv preprint arXiv:1511.06085.
  * [27] Z. Cheng, H. Sun, M. Takeuchi, J. Katto, Learned image compression with discretized gaussian mixture likelihoods and attention modules, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 7939–7948.
  * [28] M. Akbari, J. Liang, J. Han, C. Tu, Generalized octave convolutions for learned multi-frequency image compression, arXiv preprint arXiv:2002.10032.
  * [29] M. A. Yílmaz, O. Kelesş, H. Güven, A. M. Tekalp, J. Malik, S. Kíranyaz, Self-organized variational autoencoders (self-vae) for learned image compression, in: 2021 IEEE International Conference on Image Processing (ICIP), IEEE, 2021, pp. 3732–3736.
  * [30] A. Van den Oord, N. Kalchbrenner, L. Espeholt, O. Vinyals, A. Graves, et al., Conditional image generation with pixelcnn decoders, Advances in neural information processing systems 29.
  * [31] J. Ballé, V. Laparra, E. P. Simoncelli, Density modeling of images using a generalized normalization transformation, arXiv preprint arXiv:1511.06281.
  * [32] H. Khalil, K. Rose, S. L. Regunathan, The asymptotic closed-loop approach to predictive vector quantizer design with application in video coding, IEEE transactions on image processing 10 (1) (2001) 15–23.
  * [33] B. Vishwanath, T. Nanjundaswamy, K. Rose, Effective prediction modes design for adaptive compression with application in video coding, IEEE Transactions on Image Processing 31 (2021) 636–647.
  * [34] J. Bégaint, F. Racapé, S. Feltman, A. Pushparaja, Compressai: a pytorch library and evaluation platform for end-to-end compression research, arXiv preprint arXiv:2011.03029.
  * [35] R. T. L. T. J. B. E. A. N. J. F. M. George Toderici, Wenzhe Shi, Workshop and challenge on learned image compression (clic2020) (2020).
URL http://www.compression.cc

  * [36] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, in: Y. Bengio, Y. LeCun (Eds.), 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
URL http://arxiv.org/abs/1412.6980

  * [37] J. Duda, Asymmetric numeral systems: entropy coding combining speed of huffman coding with compression rate of arithmetic coding (2014).  arXiv:1311.2540.
  * [38] F. Kamisli, Learned-block-based-image-compression.
URL https://github.com/metu-kamisli/Learned-block-based-image-compression

  * [39] E. Kodak, Kodak lossless true color image suite (PhotoCD PCD0992).
URL http://r0k.us/graphics/kodak

  * [40] N. Asuni, A. Giachetti, TESTIMAGES: A large-scale archive for testing visual devices and basic image processing algorithms (SAMPLING 1200 RGB set), in: STAG: Smart Tools and Apps for Graphics, 2014.
URL
https://sourceforge.net/projects/testimages/files/OLD/OLD_SAMPLING/testimages.zip

  * [41] Tensorflow compression.
URL https://github.com/tensorflow/compression

  * [42] F. Bossen, B. Bross, K. Suhring, D. Flynn, Hevc complexity and implementation analysis, IEEE Transactions on Circuits and Systems for Video Technology 22 (12) (2012) 1685–1696.  doi:10.1109/TCSVT.2012.2221255.
  * [43] B.-F. Wu, C.-F. Lin, Memory-efficient architecture for jpeg 2000 coprocessor with large tile image, IEEE Transactions on Circuits and Systems II: Express Briefs 53 (4) (2006) 304–308.  doi:10.1109/TCSII.2005.862042.
  * [44] C. C. Chi, M. Alvarez-Mesa, B. Juurlink, G. Clare, F. Henry, S. Pateux, T. Schierl, Parallel scalability and efficiency of hevc parallelization approaches, IEEE Transactions on circuits and systems for video technology 22 (12) (2012) 1827–1838.


# Light-induced ideal Weyl semimetal in HgTe via nonlinear phononics

Dongbin Shin<EMAIL_ADDRESS>Department of Physics and Photon Science,
Gwangju Institute of Science and Technology (GIST), Gwangju 61005, Republic of
Korea  Max Planck Institute for the Structure and Dynamics of Matter and
Center for Free Electron Laser Science, 22761 Hamburg, Germany  Angel Rubio
<EMAIL_ADDRESS>Max Planck Institute for the Structure and Dynamics of
Matter and Center for Free Electron Laser Science, 22761 Hamburg, Germany
Nano-Bio Spectroscopy Group, Departamento de Fisica de Materiales, Universidad
del País Vasco UPV/EHU- 20018 San Sebastián, Spain  Center for Computational
Quantum Physics (CCQ), The Flatiron Institute, 162 Fifth avenue, New York NY
10010.  Peizhe Tang<EMAIL_ADDRESS>School of Materials Science and
Engineering, Beihang University, Beijing 100191, P. R. China  Max Planck
Institute for the Structure and Dynamics of Matter and Center for Free
Electron Laser Science, 22761 Hamburg, Germany

###### Abstract

Interactions between light and matter allow the realization of out-of-
equilibrium states in quantum solids. In particular, nonlinear phononics is
one of the efficient approaches to realizing the stationary electronic state
in non-equilibrium. Herein, by using extended $ab~{}initio$ molecular
dynamics, we identify that long-lived light-driven quasi-stationary geometry
could stabilize the topological nature in the material family of HgTe
compounds. We show that coherent excitation of the infrared-active phonon mode
results in a distortion of the atomic geometry with a lifetime of several
picoseconds. We show that four Weyl points are located exactly at the Fermi
level in this non-equilibrium geometry, making it an ideal long-lived
metastable Weyl semimetal. We propose that such a metastable topological phase
can be identified by photoelectron spectroscopy of the Fermi arc surface
states or ultrafast pump-probe transport measurements of the nonlinear Hall
effect.

Non-equilibrium states driven by strong light-matter interactions have
attracted a lot of attention in recent years, due to intriguing physical
phenomena without equilibrium counterparts and interesting ultrafast
applications [1, 2, 3, 4, 5, 6]. Unlike modulating ground states, the coupling
between the laser field and matter often produces non-thermal
electronic/phononic changes in their electronic structures and charge
populations. For example, Floquet engineering is expected to be realized by
applying a continuous light [7, 3, 8, 9, 4, 6, 5, 10, 11]. Recent theoretical
works have suggested that coupling between a time-periodic oscillating field
and the Dirac bands could result in a topologically non-trivial insulating
state among the Floquet spectra [12]. However, the electronic Floquet states
are difficult to be achieved and only reported on a few materials [13, 10, 14]
because of dissipation, decoherence, and heating effects [15, 14]. In
contrast, nonlinear coupling between materials and phonon via terahertz (THz)
pulses can lead to robust quasi-stationary states [4, 16, 17, 18, 19, 20, 21],
such an approach is named the nonlinear phononics. We consider an infrared
(IR) active phonon mode resonantly excited by the THz pulse, its oscillation
amplitude can be largely enhanced to an anharmonic range. Via the nonlinear
phonon interaction coupled with another phonon (e.g. $R_{c}$ mode), a time-
averaged atomic displacement can be induced along $R_{c}$’s oscillating
direction, resulting in an effective lattice distortion in ultrafast time
scale. Consequently, electronic properties can be efficiently manipulated,
realizing long-lived non-equilibrium states [22, 23]. In this context, light-
induced superconductivity has been experimentally shown in cuprate, alkali-
doped C60, and organic conductors [24, 18, 19]. Light-induced magnetic
orderings in perovskite ErFeO3 and piezomagnetic crystal CoF2 [25, 26], as
well as a light-induced ferroelectric transition in SrTiO3 [20, 27, 28, 29,
30] have been assessed from the same perspective. On the theoretical side, the
light-induced non-equilibrium states induced by nonlinear phononics have
usually been investigated using a model Hamiltonian that includes the relevant
coupling to specific modes and some phenomenological dissipative terms [20,
29, 21, 31, 32]. In this regard, first-principles calculations considering
non-equilibrium dynamics, including both couplings to the external field and
realistic dissipation channels, are needed.

As a notable exotic quantum state of matter, topological Dirac and Weyl
semimetals host gapless excitations with topological charges in the bulk and
Fermi arc states on the surface [33, 34], which have distinctive electronic
and optical responses. The electronic structures of these semimetals could be
affected by light-matter interactions. For example, the light-induced
anomalous Hall effect has been observed in the graphene [35], where the
asymmetric photo-carrier population effect induced by the circularly polarized
light contributes to an anomalous Hall signal [36]. Furthermore, the phase
transition from a Dirac semimetal to a Floquet Weyl semimetal in NaBi3 was
signified as an intentional breaking of the time-reversal symmetry by the
driven light field [37]. On the other hand, a long-lived meta-stable lattice
distortion could also result in the stationary states in non-equilibrium [38,
39]. For example, the driven phase transition is observed in bulk WTe2 [38,
40] and ZrTe5 [41]. Experimental evidence for photo-carrier-induced non-
equilibrium dynamics in Weyl materials are reported for $\beta$-WP2 [42] and
elemental tellurium [43, 44] recently, in which the possible topological phase
transition and topological band modulations have been observed. However,
because of various excrescent contributions from non-topological states near
the Fermi level in these materials mentioned above [45, 33], the ultrafast
transport and optical signals will be strongly influenced by these trivial
states and cannot be entirely ascribed to light-induced topological
characteristics.

Herein, we demonstrate a robust metastable light-induced topological phase
transition to an ideal Weyl semimetal from a trivial semimetallic state in the
family of bulk HgTe, as shown in the schema in Figs. 1(a) and 1(b). To achieve
this, we have employed a new implementation of $ab~{}initio$ molecular
dynamics (AIMD), which additionally introduces the direct interaction between
the classical ions and the oscillating external electric field, including
realistic nonlinear phonon interactions and dissipation via phonon-phonon
scatterings [20, 29, 21, 46, 47]. In this way, we identify a non-equilibrium
geometry with distortion along the combined direction of IR modes (named as
IRy and IRz shown Figs. 1(d) and 1(e)) when the IRx mode in HgTe is resonantly
excited (see Fig. 1(c)). The new steady structure of HgTe hosts an ideal Weyl
semimetal phase with a lifetime of several picoseconds. When we increase the
laser intensity, the amplitude of the atomic distortion is enhanced, further
enlarging the distance between Weyl points in bulk and the length of the Fermi
arc on the surface. In addition, similar THz field-induced dynamics are
observed in bulk HgSe and $\beta$-HgS. Because the light-driven HgTe is an
ideal Weyl semimetal without the trivial band contributions around the Fermi
level, the signatures measured by time- and angle-resolved photoemission
spectroscopy (TrARPES) and nonlinear Hall conductivity [48, 49, 50] have a
clear topological origin.

Figure 1:  Light-induced Weyl points in bulk HgTe through the nonlinear phonon
interaction between IR modes. (a) Schematic phase diagram of THz field-induced
Weyl points from trivial semi-metal HgTe and the modification of potential
energy surface for induced phonon mode under the excited IRx mode. (b)
Schematic image of THz field-induced atomic geometry distortion of HgTe
through the nonlinear phonon interaction. Three degenerated IR phonon modes
along (c) x-direction, (d) y-direction, and (e) z-direction.  Figure 2:
Nonlinear phonon interaction in HgTe and light-induced atomic geometry
distortion. (a) Nonlinear phonon dynamics of IR phonon modes simulated by AIMD
simulation. (b) Modified potential energy surface of $Q_{ind}$ by shifted
atomic geometry toward $Q_{x}$. Time profile of (c) THz field and (d) induced
phonon modes evaluated from AIMD simulations. (e) E-field strength dependence
of geometry distortion along $Q_{ind}$ in HgTe family.

Before discussing the light-induced dynamical simulations, we describe the
vibrational properties and nonlinear phonon interactions in HgTe. In ambient
conditions, bulk HgTe with the Zinc Blende structure hosts three degenerate IR
phonon modes at the $\Gamma$ point (see Figs. 1(c)-1(e)). To investigate their
nonlinear phonon interactions, we perform conventional AIMD simulations in
which the initial atomic structure is distorted along the eigenvector
direction of IRx mode with the value of
$\bf{Q}$${}_{x}(t=0)=1.2$Å$\sqrt{\mu}$. As shown in Fig. 2(a), induced
oscillations of IRy and IRz modes can be observed, indicating the anharmonic
phonon couplings between the IRx and a linearly combined oscillation of IRy
and IRz. Notably, the nonlinear phonon interaction shifts the origin of the
oscillation, corresponding to a lattice distortion along
$\hat{\bf{Q}}$${}_{ind}=0.81\hat{\bf{Q}}$${}_{y}+0.58\hat{\bf{Q}}$z where
$\hat{\bf{Q}}$y,z denotes the unit vector along the direction of equilibrium
IRy,z eigenmodes. Their nonlinear dynamics can also be revealed from the
calculated potential energy surface as shown in Fig. 2(b). At equilibrium
($\bf{Q}$${}_{x}=0$), the potential surface shows a parabolic curve with a
minimum energy point at $\bf{Q}$${}_{ind}^{min}=0$. For both positive and
negative phonon displacements of $\bf{Q}$${}_{x}=\pm 1.2$ Å$\sqrt{\mu}$, the
minimum points shift to a negative value ($\bf{Q}$${}_{ind}^{min}<0$)
regardless of the sign of $\bf{Q}$x. These results indicate that a resonant
pumping of the IRx mode could strongly modify the potential surface. Similar
nonlinear phonon interactions are obtained in HgSe and $\beta$-HgS under
resonant THz pumping (see Supplementary Materials (SM) [51]).

Anharmonicity and dissipation effects of phonons should be considered to
describe the non-equilibrium structural dynamics induced by a THz pulse. To
achieve this goal, we proceed with extended AIMD simulations that include
direct interaction between oscillating electric field and charged ions under
the Born-Oppenheimer approximation and in the periodic boundary condition with
a $3\times 3\times 3$ supercell geometry of HgTe (see SM for computational
details, including discussion on the role of screening effect) [51]. In
practice, the THz field excites the IRx phonon mode via an electric field
oscillating:
$\bf{E}$$(t)=E_{0}\sin(\omega(t-t_{0}))e^{-(t-t_{0})^{2}/\sigma^{2}}\bf{\hat{x}}$
in which $\sigma=0.8~{}ps$, $t_{0}=2.5~{}ps$ and $\omega=3.5$ THz with a
maximum amplitude of $1.25$ MV/cm as shown in Fig. 2(c). Under this resonant
laser pumping, the oscillation amplitude of IRx is increased up to $1.8$
Å$\sqrt{\mu}$ until $4$ ps. Such excitation leads to the emergence of the non-
zero amplitude of induced mode ($\bf{Q}$ind) through the nonlinear phononics,
as shown in Fig. 2(d). Notably, the origin of $\bf{Q}$${}_{ind}(t)$ is shifted
away from the equilibrium position in the steady state with a lifetime of a
few picoseconds, which is much longer than the broadening of the THz pumping
pulse ($\sigma=0.8$ ps). Therefore, a meta-stable structural distortion is
achieved through nonlinear phonon interaction between the highly excited IRx
and IRind modes. After a few picoseconds (>8 ps), excited IRx and IRind
phonons dissipate via anharmonic couplings with other phonon modes, acting as
a heating source and increasing the total temperature of this system. In the
SM, we provide detailed discussions on the role of dissipation and thermostat
on the light-induced meta-stable atomic geometry [51].

The time-averaged lattice distortion $\bf{Q}$${}_{ind}^{avg}$ depends on the
intensity of the applied THz pulse as shown in Fig. 2(e). For instance, the
stronger THz pulse ($\sim 1.75$ MV/cm) leads to a more considerable lattice
distortion ($\bf{Q}$${}^{avg}_{ind}\sim 0.84$ Å$\sqrt{\mu}$) in HgTe, compared
with the case with weaker pumping strength ($\bf{Q}$${}^{avg}_{ind}\sim 0.08$
Å$\sqrt{\mu}$ at E$=0.5$ MV/cm). The same tendency is found in HgSe and
$\beta$-HgS, which have similar nonlinear phonon interactions with that in
HgTe. These results show that the amplitude of lattice distortion in the meta-
stable state can be controlled by varying the intensity of the THz field pulse
in the HgTe family.

To better understand the role of anharmonicity and dissipation, we compare our
extended AIMD simulation with the results from the simplified lattice model
fitted to $ab~{}initio$ calculations [20, 29]. This model considers harmonic
IRx and IRind modes and their minimal nonlinear coupling and the equations of
motion are given as follows:
$\ddot{Q}_{x}+\gamma\dot{Q}_{x}+\Omega^{2}Q_{x}=-2k_{nl}Q_{x}Q_{ind}+Z^{*}E(t)$,
and $\ddot{Q}_{ind}+\gamma\dot{Q}_{ind}+\Omega^{2}Q_{ind}=-k_{nl}Q_{x}^{2}$,
where $Z^{*}$ is the mode effective charges, $k_{nl}$ is the coupling
coefficient for nonlinear phonon interaction, $\gamma$ is the dissipation
rate, and $\Omega$ is the frequency of these phonon modes (see SM for details)
[51]. We find that this model can reproduce the $ab~{}initio$ results about
the structural evolution of HgTe very well in the low range of field strength
($E$ < $1$ MV/cm) (see Fig. 2(e)). However, it fails to reproduce the
displacements in the high field strength range ($E$ > $1$ MV/cm), due to the
restrictive description of dissipation channels and anharmonicity. These
results demonstrate that our extended AIMD calculations can provide more
realistic dynamics in HgTe to investigate light-induced metastable phases than
the model simulations.

Figure 3:  Formation of Weyl points in the shifted equilibrium position by the
nonlinear phonon interaction and possible experimental observation in TrARPES.
(a) Displacement direction of the induced mode ($Q_{ind}$) in cubic geometry.
(b) The band structure of HgTe with the equilibrium geometry on $k_{z}=0$
plane. (c) The band structure of HgTe in the THz-induced non-equilibrium
geometry on the $k_{z}=35.3$ mÅ-1 plane. (d) Distance between Weyl points
depending on the displacement of the induced mode $Q_{ind}$ in HgTe family.
Fermi arc surface states in the (100) surface with (e) $Q_{ind}=0.5$
Å$\sqrt{\mu}$ and (f) $Q_{ind}=0.8$ Å$\sqrt{\mu}$. In (a), pink plane
indicates the (110) mirror plane. In (b) and (c), the color map indicates the
energy variation corresponding to its z-axis.

The THz pulse-induced structural distortion drives bulk HgTe into an ideal
Weyl semimetal with four Weyl points at the Fermi level in non-equilibrium. At
the equilibrium geometry, bulk HgTe is a gapless trivial semimetal with
inverted bands, in which the degeneracy at the $\Gamma$ point (see Fig. 3(b))
is protected by the crystal symmetry ($F\bar{4}3m$ space group). Under atomic
distortion along $\bf{Q}$ind, only one mirror symmetry $m_{1,1,0}$ (see Fig.
3(a)) is preserved in new structure ($Cm(8)$ space group). With a displacement
of $\bf{Q}$${}_{ind}=0.5$ Å$\sqrt{\mu}$, fourfold degenerated points at the
Fermi level split into four Weyl points away from the $\Gamma$ point, their
positions in the momentum space are labeled as $W_{1}^{\pm}$ and $W_{2}^{\pm}$
as shown in Fig. 3(c), in which $\pm$ stands for the sign of the topological
charges. Under mirror symmetry and time-reversal symmetry, these Weyl points
can be converted from one to others (see SM [51]), sharing the same energy
eigenvalue. Such a phase is named the ideal Weyl semimetal [52] since Weyl
points are precisely at the Fermi level and low-energy electronic states are
always topologically non-trivial. In Fig. 3(e), we plot the (100) surface
state, where four Weyl points project to different positions on the side
surface and Fermi arc surface states connect each pair with opposite chiral
charges, which can be detected in TrARPES [48]. Our result shows that the THz-
driven lattice distortion can induce a topological phase transition in bulk
HgTe.

The positions of the light-induced Weyl points depend on the strength of the
lattice distortion. In Fig 3(d), we show that the change of distance ($d_{k}$)
between two Weyl points $W_{1(2)}^{+}$ and $W_{1(2)}^{-}$ with respect to
lattice distortions. Before the critical value of $\bf{Q}$${}_{ind}\sim 0.4$
Å$\sqrt{\mu}$, the nodal lines instead of Weyl points (see Fig. 3(d)) can be
observed in the bulk and the similar result has been reported in strained
cases [52]. Above this value, for example, $\bf{Q}$${}_{ind}=1.0$
Å$\sqrt{\mu}$ corresponding to $1.75$ MV/cm maximum field strength gives rise
to a Weyl point separation of $d_{k}=52.8$ mÅ-1. The same conclusion is valid
for bulk HgSe and $\beta$-HgS (see Fig. 3(d)). The evolution of the light-
induced Weyl points with increasing pumping strength can be captured from the
evolution of the Fermi arc surface states. Compared with the case for
$\bf{Q}$${}_{ind}=0.5$ Å$\sqrt{\mu}$ shown in Fig. 3(e), the Fermi arc surface
state (see Fig. 3(f)) becomes larger under stronger THz laser pumping
($\bf{Q}$${}_{ind}=0.8$ Å$\sqrt{\mu}$ correspondingly). We anticipate that
such features will be captured in TrARPES with THz laser pumping, a powerful
tool to detect the Fermi arc surface state and its related dynamic evolution
[48, 49, 50, 10].

Figure 4:  Possible experimental observable of the light-induced Weyl points;
nonlinear Hall effect. (a) Berry curvature $\Omega_{x}$ of distorted HgTe with
$Q_{ind}=0.5$ Å$\sqrt{\mu}$ at $k_{z}=35.3$ mÅ-1 plane. (b) Berry curvature
dipole of distorted HgTe with $Q_{ind}=0.5$ Å$\sqrt{\mu}$ in 3D plot. Diagonal
terms of Berry curvature dipole (c) $D_{xx}$ and (d) $D_{yy}$ with distorted
geometry along $Q_{ind}$ at $100$ K. (e) Temperature and Fermi level
dependency on the Berry curvature dipole $D_{xx}$ with $Q_{ind}=0.5$
Å$\sqrt{\mu}$. (f) The absolute maximum value of Berry curvature dipole
$D_{xx}$ near the Fermi level at $100$ K.

Besides TrARPES, THz field-induced Weyl points can also be experimentally
observed in quantum transport measurements, such as the nonlinear Hall effect
[53, 54, 55]. The nonlinear Hall current defined as
$\bm{J}_{a}=\chi_{abc}\bm{E}_{b}\bm{E}^{*}_{c}$, with
$\chi_{abc}=\epsilon_{adc}\frac{e^{3}\tau}{2(1+i\omega\tau)}D_{bd}$ depends on
the Berry curvature dipole $D_{bd}$, in which $e$ is the electron charge,
$\tau$ is the scattering time, and $\omega$ is the frequency of the probe
laser field [53]. The total Berry curvature is zero in HgTe due to the time-
reversal symmetry. However, the inversion symmetry-breaking condition induced
by the strong light-matter interaction in HgTe leads to a non-zero Berry
curvature dipole that provides the nonlinear Hall effect (see Fig. 4) [56].

The distorted geometry of HgTe induced by nonlinear phonon interaction shows a
non-zero Berry curvature dipole. Based on the symmetry analysis, the bulk HgTe
in equilibrium shows zero Berry curvature dipole due to inversion symmetry. In
contrast, the distorted HgTe has a finite value with elements $D_{xx}=-D_{yy}$
and $D_{zz}=0$ constrained by the mirror symmetry (see SM [51]). Furthermore,
we use the density functional theory to calculate Berry curvature and Berry
curvature dipole for distorted HgTe. In Fig. 4(a), we plot Berry curvature
$\Omega$ at $k_{z}=35.3$ mÅ-1 plane of Brillouin zone, which asymmetrically
distributes around Weyl points and could result in a non-zero Berry curvature
dipole distribution $\bf{D}$$(\textbf{k})$ [55] (see Fig. 4(b)) because
$\bf{D}$$(\textbf{k})$ $\sim
f_{\textbf{k}}\nabla_{\textbf{k}}\Omega(\textbf{k})$, where $f_{\textbf{k}}$
is the occupation of Bloch state. In Fig. 4(c) and 4(d), we show the
calculated Berry curvature dipole for structures with different $\bf{Q}$ind
values, corresponding to different THz pumping intensities. Furthermore, we
find that $D_{xx}$ and $D_{yy}$ have two peaks with opposite signs around the
Fermi level, corresponding to contributions from electron and hole parts of
Weyl fermions.

We also investigated thermal effects on the Berry curvature dipole to verify
whether the dynamic nonlinear Hall effect can be observed at finite
temperatures. The temperature-dependent Berry curvature dipole is evaluated
with an electronic Fermi-Dirac distribution function at a given temperature
$T$. The results are shown in Fig. 4(e), whose temperature dependency
indicates a measurable nonlinear Hall current can be experimentally achieved
at high-temperature regions (even at room temperature). Once we fix the
temperature (see Fig. 4(f)), peak values of Berry curvature dipole increase
with the enhancement of pumping strength (larger $\bf{Q}$ind value) because of
the increasing distance between Weyl points with opposite chiral charges.
These results indicate the possibility of identifying light-induced Weyl
points in bulk HgTe via pump-probe ultrafast experiments.

In summary, we theoretically investigate a light-induced topological phase
transition in the family of HgTe materials through nonlinear phonon
interaction. We use an extended AIMD framework to show that THz laser coherent
pumping could drive a new non-equilibrium steady-state with a distorted atomic
structure via nonlinear phononics. The long-lived lattice distortion is
enhanced by increasing the pump laser strength. Bulk HgTe with the distorted
geometry hosts four Weyl points at the Fermi level, an ideal Weyl semimetal in
non-equilibrium. The same scenario has been observed in HgSe and $\beta$-HgS.
We anticipate that ultrafast pump-probe experiments can be used to validate
the proposed topological phase transitions in bulk HgTe family via detecting
surface Fermi arc states and nonlinear Hall effects. Our work provides an
opportunity to optically manipulate topological properties of quantum
materials via nonlinear phononics, which is important for future ultrafast
all-optical topological device switching. Furthermore, the topics discussed
herein can be extended to the new field of cavity materials engineering or
cavity-induced phase transition [5, 6, 30, 8, 2].

###### Acknowledgements.

We acknowledge financial support from the Grupos Consolidados (IT1453-22), the
Cluster of Excellence ‘CUI: Advanced Imaging of Matter’ of the Deutsche
Forschungsgemeinschaft (DFG) - EXC 2056 - project ID 390715994. D.S. was
supported by the National Research Foundation of Korea (NRF) grant funded by
the Korea government (MSIT) (No. RS-2023-00253716). The computational work was
supported by the National Supercomputing Center with supercomputing resources
including technical support (KSC-2023-CRE-0073). P.Z.T. was supported by the
National Natural Science Foundation of China (Grants No. 12234011 and No.
12374053). We also acknowledge support from the Max Planck–New York Center for
Non-Equilibrium Quantum Phenomena. The Flatiron Institute is a division of the
Simons Foundation. The authors thank Dante M. Kennes and Noejung Park for the
fruitful discussion.

## References

  * de la Torre  _et al._ [2021] A. de la Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Rev. Mod. Phys. 93, 041002 (2021).
  * Kennes and Rubio [2022] D. M. Kennes and A. Rubio, arXiv (2022), 2204.11928 .
  * Bao  _et al._ [2022] C. Bao, P. Tang, D. Sun, and S. Zhou, Nat. Rev. Phys. 4, 33 (2022).
  * Disa  _et al._ [2021] A. S. Disa, T. F. Nova, and A. Cavalleri, Nat. Phys. 17, 1087 (2021).
  * Basov  _et al._ [2021] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, Nanophotonics 10, 549 (2021).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Nature 606, 41 (2022).
  * Oka and Kitamura [2018] T. Oka and S. Kitamura, Annu. Rev. Condens. Matter Phys. 10, 1 (2018).
  * Hubener  _et al._ [2021] H. Hubener, U. D. Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Nat. Mater. 20, 438 (2021).
  * Ruggenthaler  _et al._ [2018] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, Nat. Rev. Chem. 2, 0118 (2018).
  * Zhou  _et al._ [2023a] S. Zhou, C. Bao, B. Fan, H. Zhou, Q. Gao, H. Zhong, T. Lin, H. Liu, P. Yu, P. Tang, S. Meng, W. Duan, and S. Zhou, Nature 614, 75 (2023a).
  * Zhou  _et al._ [2023b] S. Zhou, C. Bao, B. Fan, F. Wang, H. Zhong, H. Zhang, P. Tang, W. Duan, and S. Zhou, Phys. Rev. Lett. 131, 116401 (2023b).
  * Rudner and Lindner [2020] M. S. Rudner and N. H. Lindner, Nat. Rev. Phys. 2, 229 (2020).
  * Wang  _et al._ [2013] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
  * Aeschlimann  _et al._ [2021] S. Aeschlimann, S. A. Sato, R. Krause, M. Chávez-Cervantes, U. De Giovannini, H. Hübener, S. Forti, C. Coletti, K. Hanff, K. Rossnagel, A. Rubio, and I. Gierz, Nano Lett. 21, 5028 (2021).
  * Iadecola  _et al._ [2015] T. Iadecola, T. Neupert, and C. Chamon, Phys. Rev. B 91, 235133 (2015).
  * Shin  _et al._ [2020] D. Shin, S. A. Sato, H. Hübener, U. D. Giovannini, N. Park, and A. Rubio, Npj Comput. Mater. 6, 182 (2020).
  * Först  _et al._ [2011] M. Först, C. Manzoni, S. Kaiser, Y. Tomioka, Y.-n. Tokura, R. Merlin, and A. Cavalleri, Nat. Phys. 7, 854 (2011).
  * Mitrano  _et al._ [2016] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. D. Pietro, D. Pontiroli, M. Riccò, A. Subedi, S. R. Clark, D. Jaksch, and A. Cavalleri, Nature 530, 461 (2016).
  * Buzzi  _et al._ [2020] M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. A. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, T. Siegrist, J. A. Schlueter, K. Miyagawa, K. Kanoda, M.-S. Nam, A. Ardavan, J. Coulthard, J. Tindall, F. Schlawin, D. Jaksch, and A. Cavalleri, Phys. Rev. X 10, 031028 (2020).
  * Kozina  _et al._ [2019] M. Kozina, M. Fechner, P. Marsik, T. v. Driel, J. M. Glownia, C. Bernhard, M. Radovic, D. Zhu, S. Bonetti, U. Staub, and M. C. Hoffmann, Nat. Phys. 15, 387 (2019).
  * Subedi [2015] A. Subedi, Phys. Rev. B 92, 214303 (2015).
  * Mankowsky  _et al._ [2016] R. Mankowsky, M. Först, and A. Cavalleri, Rep. Prog. Phys. 79, 064503 (2016).
  * Basov  _et al._ [2017] D. N. Basov, R. D. Averitt, and D. Hsieh, Nat. Mater. 16, 1077 (2017).
  * Fausti  _et al._ [2011] D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * Nova  _et al._ [2017] T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, Nat. Phys. 13, 132 (2017).
  * Disa  _et al._ [2020] A. S. Disa, M. Fechner, T. F. Nova, B. Liu, M. Först, D. Prabhakaran, P. G. Radaelli, and A. Cavalleri, Nat. Phys. 16, 937 (2020).
  * Li  _et al._ [2019] X. Li, T. Qiu, J. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, Science 364, 1079 (2019).
  * Nova  _et al._ [2019] T. F. Nova, A. S. Disa, M. Fechner, and A. Cavalleri, Science 364, 1075 (2019).
  * Shin  _et al._ [2022] D. Shin, S. Latini, C. Schafer, S. A. Sato, E. Baldini, U. D. Giovannini, H. Hubener, and A. Rubio, Phys. Rev. Lett. 129, 167401 (2022).
  * Latini  _et al._ [2021] S. Latini, D. Shin, S. A. Sato, C. Schäfer, U. D. Giovannini, H. Hübener, and A. Rubio, Proc. Natl. Acad. Sci. U.S.A. 118, e2105618118 (2021).
  * Pueyo and Subedi [2022] A. G. Pueyo and A. Subedi, Phys. Rev. B 106, 214305 (2022), 2208.05293 .
  * Subedi  _et al._ [2014] A. Subedi, A. Cavalleri, and A. Georges, Phys. Rev. B 89, 220301 (2014).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Liu  _et al._ [2020] J. Liu, F. Xia, D. Xiao, F. J. Garcia de Abajo, and D. Sun, Nat. Mater. 19, 830 (2020).
  * McIver  _et al._ [2020] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri, Nat. Phys. 16, 38 (2020).
  * Sato  _et al._ [2019] S. A. Sato, J. W. McIver, M. Nuske, P. Tang, G. Jotzu, B. Schulte, H. Hübener, U. D. Giovannini, L. Mathey, M. A. Sentef, A. Cavalleri, and A. Rubio, Phys. Rev. B 99, 214302 (2019).
  * Hübener  _et al._ [2017] H. Hübener, M. A. Sentef, U. De Giovannini, A. F. Kemper, and A. Rubio, Nat. Commun. 8, 1 (2017).
  * Sie  _et al._ [2019] E. J. Sie, C. M. Nyby, C. D. Pemmaraju, S. J. Park, X. Shen, J. Yang, M. C. Hoffmann, B. K. Ofori-Okai, R. Li, A. H. Reid, S. Weathersby, E. Mannebach, N. Finney, D. Rhodes, D. Chenet, A. Antony, L. Balicas, J. Hone, T. P. Devereaux, T. F. Heinz, X. Wang, and A. M. Lindenberg, Nature 565, 61 (2019).
  * Vaswani  _et al._ [2020] C. Vaswani, L.-L. Wang, D. H. Mudiyanselage, Q. Li, P. M. Lozano, G. D. Gu, D. Cheng, B. Song, L. Luo, R. H. J. Kim, C. Huang, Z. Liu, M. Mootz, I. E. Perakis, Y. Yao, K. M. Ho, and J. Wang, Phys. Rev. X 10, 021013 (2020).
  * Guan  _et al._ [2021] M.-X. Guan, E. Wang, P.-W. You, J.-T. Sun, and S. Meng, Nat. Commun. 12, 1 (2021).
  * Luo  _et al._ [2021] L. Luo, D. Cheng, B. Song, L.-L. Wang, C. Vaswani, P. M. Lozano, G. Gu, C. Huang, R. H. J. Kim, Z. Liu, J.-M. Park, Y. Yao, K. Ho, I. E. Perakis, Q. Li, and J. Wang, Nat. Mater. 20, 329 (2021).
  * Hu  _et al._ [2023] T. Hu, B. Su, L. Shi, Z. Wang, L. Yue, S. Xu, S. Zhang, Q. Liu, Q. Wu, R. Li, X. Zhou, J. Yuan, D. Wu, Z. Chen, T. Dong, and N. Wang, Adv. Opt. Mater. 11, 2202639 (2023).
  * Ning  _et al._ [2022] H. Ning, O. Mehio, C. Lian, X. Li, E. Zoghlin, P. Zhou, B. Cheng, S. D. Wilson, B. M. Wong, and D. Hsieh, Phys. Rev. B 106, 205118 (2022).
  * Jnawali  _et al._ [2020] G. Jnawali, Y. Xiang, S. M. Linser, I. A. Shojaei, R. Wang, G. Qiu, C. Lian, B. M. Wong, W. Wu, P. D. Ye, Y. Leng, H. E. Jackson, and L. M. Smith, Nat. Commun. 11, 3991 (2020).
  * Deng  _et al._ [2016] K. Deng, G. Wan, P. Deng, K. Zhang, S. Ding, E. Wang, M. Yan, H. Huang, H. Zhang, Z. Xu, J. Denlinger, A. Fedorov, H. Yang, W. Duan, H. Yao, Y. Wu, S. Fan, H. Zhang, X. Chen, and S. Zhou, Nat. Phys. 12, 1105 (2016).
  * Lu  _et al._ [2022] I.-T. Lu, J.-J. Zhou, J. Park, and M. Bernardi, Phys. Rev. Mater. 6, L010801 (2022).
  * Malic  _et al._ [2011] E. Malic, T. Winzer, E. Bobkin, and A. Knorr, Phys. Rev. B 84, 205406 (2011).
  * Xu  _et al._ [2015] S.-Y. Xu, C. Liu, S. K. Kushwaha, R. Sankar, J. W. Krizan, I. Belopolski, M. Neupane, G. Bian, N. Alidoust, T.-R. Chang, H.-T. Jeng, C.-Y. Huang, W.-F. Tsai, H. Lin, P. P. Shibayev, F.-C. Chou, R. J. Cava, and M. Z. Hasan, Science 347, 294 (2015).
  * Yang  _et al._ [2015] L. X. Yang, Z. K. Liu, Y. Sun, H. Peng, H. F. Yang, T. Zhang, B. Zhou, Y. Zhang, Y. F. Guo, M. Rahn, D. Prabhakaran, Z. Hussain, S.-K. Mo, C. Felser, B. Yan, and Y. L. Chen, Nat. Phys. 11, 728 (2015).
  * Gierz  _et al._ [2013] I. Gierz, J. C. Petersen, M. Mitrano, C. Cacho, I. C. E. Turcu, E. Springate, A. Stöhr, A. Köhler, U. Starke, and A. Cavalleri, Nat. Mater. 12, 1119 (2013).
  * [51] See supplemental material at url.
  * Ruan  _et al._ [2016] J. Ruan, S.-K. Jian, H. Yao, H. Zhang, S.-C. Zhang, and D. Xing, Nat. Commun. 7, 11136 (2016).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Ma  _et al._ [2019] Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A. M. M. Valdivia, S. Wu, Z. Du, C.-H. Hsu, S. Fang, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, E. Kaxiras, H.-Z. Lu, H. Lin, L. Fu, N. Gedik, and P. Jarillo-Herrero, Nature 565, 337 (2019).
  * Zhang  _et al._ [2018] Y. Zhang, Y. Sun, and B. Yan, Phys. Rev. B 97, 041101 (2018).
  * Gao  _et al._ [2020] Y. Gao, F. Zhang, and W. Zhang, Phys. Rev. B 102, 245116 (2020).