# Segment-Level Vectorized Beam Search Based on
Partially Autoregressive Inference

###### Abstract

Attention-based encoder-decoder models with autoregressive (AR) decoding have
proven to be the dominant approach for automatic speech recognition (ASR) due
to their superior accuracy. However, they often suffer from slow inference.
This is primarily attributed to the incremental calculation of the decoder.
This work proposes a partially AR framework, which employs segment-level
vectorized beam search for improving the inference speed of an ASR model based
on the hybrid connectionist temporal classification (CTC) attention-based
architecture. It first generates an initial hypothesis using greedy CTC
decoding, identifying low-confidence tokens based on their output
probabilities. We then utilize the decoder to perform segment-level vectorized
beam search on these tokens, re-predicting in parallel with minimal decoder
calculations. Experimental results show that our method is $12$ to $13$ times
faster in inference on the LibriSpeech corpus over AR decoding whilst
preserving high accuracy.

Index Terms—  Decoding algorithm, autoregressive, semi-autoregressive, hybrid
CTC/attention, beam search

##  1 Introduction

Due to recent advances in deep learning, automatic speech recognition (ASR)
has witnessed remarkable achievements [1, 2, 3]. ASR plays an essential role
in facilitating human-computer interaction by providing an interface for
converting audio to text and demonstrating substantial applicability in real-
world scenarios. In particular, the RNN-Transducer model [4], which operates
relatively fast and can be extended to streaming speech recognition, is widely
used in real-world applications. Recent research in ASR has also made
significant progress in achieving higher accuracy through Attention-based
Encoder-Decoder (AED) models [5, 6]. AED models have also been utilized in
models such as Whisper [7] and speech translation [8], and their usefulness
has been reevaluated. However, there are various trade-offs that can
potentially limit its application in certain scenarios. For example, one can
construct an ASR system with high accuracy through large and complex models,
but this comes at the expense of increased computational cost and inference
time.

There has been extensive research devoted to alleviating the trade-off between
recognition accuracy and inference speed. Inspired by the great success in
neural machine translation [9, 10], non-autoregressive (NAR) models have been
actively studied in the context of ASR, with the aim of achieving fast
inference [11, 12, 13, 14]. Compared to the conventional autoregressive (AR)
models [15, 16], which generates output at each step conditioned on the
previously generated outputs, NAR models can produce multiple outputs
simultaneously. This parallel computation accelerates the inference process of
ASR, resulting in a significant reduction in inference time compared to AR
models [17]. However, it remains a challenge to achieve the same level of
recognition accuracy as AR models. In addition, NAR models require a complex
model structure or a unique training strategy for the successful
implementation of parallel generation during inference [17].

Regarding AR decoding, the decoder is trained to learn linguistic information.
This allows us to utilize the relationship between the current token and the
previous token. Furthermore, it is common to enhance accuracy by implementing
the beam search algorithm, which is a heuristic search for the best
hypothesis. However, because AR requires the previous tokens to estimate the
next token, it is not possible to parallelize the inference of a single audio.
As a consequence, achieving NAR-level speed through inference parallelization
is challenging.

In this paper, we focus on the difference in trade-off balance between AR and
NAR, and propose a new decoding method, the partially autoregressive (PAR)
method as a fast and accurate decoding method. By utilizing the NAR approach
and segment-level vectorized beam search, we can compensate for the weaknesses
of both AR and NAR. This results in a better trade-off between accuracy and
latency. In particular, we show that by optimizing the inference operations of
the pre-trained hybrid connectionist temporal classification (CTC) and
attention model, we are able to achieve NAR-level inference speed while
maintaining its high accuracy, as well as not needing additional training of
the model.

In this paper, we make the following contributions:

  * •

We propose a new decoding framework that combines AR and NAR.

  * •

We demonstrate a better balance of accuracy-latency trade-off without
additional training.

(a) Autoregressive (AR)

(b) Non-autoregressive (NAR)

(c) Partially autoregressive (PAR)

Fig. 1: Overview and comparison of AR, NAR, and PAR decoding. ``¡sos¿— denotes
the “start-of-sequence” symbol,and the mask token is denoted by ``#— or red
characters. PAR is a hybrid of AR and NAR methods, in which the masking
process is applied first, followed by segment-level vectorized beam search.

##  2 Related Works

Numerous studies have been conducted to balance the accuracy-latency trade-
off. In general, there are two main approaches; reducing latency while
preserving good accuracy or improving accuracy while preserving fast inference
speed.

To reduce latency, there are several methods to lower the computational cost
during inference. One such method is the pruning technique [18, 19, 20], which
reduces the number of parameters in a trained model by identifying subnetworks
with better performance and fewer parameters. Another method is knowledge
distillation [21, 22], which trains a smaller model to reduce the number of
parameters required during inference. To further reduce the number of search
iterations, hypotheses can be predicted in a batch [23], or the search process
can be stopped prematurely [24]. In addition, other approaches focus on
utilizing machine resources more efficiently[25, 26, 27], rather than relying
solely on the model architecture, to achieve faster processing. In this study,
we parallelized the AR decoding of single audio to achieve high-speed
processing. The process is similar to batch processing of hypotheses, but in
this case, the audio is segmented and parallelized in addition to the
hypotheses.

On the other hand, we can also improve the accuracy of a low-latency model,
such as a non-autoregressive (NAR) model, though at a lower accuracy than
typical AR models. For example, by replacing and inserting the model output
with a mask sequence [11, 12, 13, 28], NAR can more accurately predict the
target sequence compared to a standard NAR model. We can also improve the
accuracy by utilizing an external language model[29, 30]. Additionally,
several investigations have been undertaken to improve both the inference
speed and accuracy, such as parallelizing the decoding process in the
streaming situation[31]. In our work, it is not possible to perform parallel
processing of all tokens as investigated above. However, high accuracy can
still be achieved by utilizing AR decoding in areas where parallel processing
leads to decreased accuracy.

Research has investigated the combination of AR and NAR methods [32]. The
researchers trained a dual-mode Transformer decoder that can be used for both
NAR and AR-style processes. They first applied NAR-style decoding to generate
several hypotheses, then used AR-style rescoring to select the best
hypothesis. Our approach indeed uses AR decoding after NAR decoding. However,
we utilize NAR decoding in our work to parallelize AR decoding within a single
audio, so the purpose of NAR decoding is different.

(d) AR
---
(e) NAR
Fig. 2: Average inference time and proportion of time spent on the encoder,
decoder, and CTC computation for (a) AR, as well as the encoder and decoder
computation for (b) NAR architectures.

##  3 Background

In this section, to provide a clear understanding of how the proposed approach
combines the benefits of autoregressive (AR) and non-autoregressive (NAR)
decoding, we begin with a brief overview of the conventional encoder-decoder-
based models, including hybrid CTC-attention [33] and Mask-CTC [13].
Additionally, we investigate the role of each encoder and decoder in
influencing the overall inference time for each model.

Table 1:  Decoding example from the LibriSpeech test-clean set (1089-134686-0002). The target sequence is initially predicted by gCTC and replaced with masks (“#”). Consecutive masks are merged into one mask. Red tokens indicate the best hypotheses from the segment-level vectorized beam search at each iteration. We set max_iteration to $5$ and have five sentences for this sample. However, we only show iterations 1 and 2 since there is no difference from the third iteration.  Masked sequence | after early night# the yellow lamps would light up here and there the squalid quarter of the#el#
---|---
iteration=1 |  after early nightfall the yellow lamps would light up here and there the squalid quarter of the braels
iteration=2 |  after early nightfall the yellow lamps would light up here and there the squalid quarter of the brothels
Ground truth | after early nightfall the yellow lamps would light up here and there the squalid quarter of the brothels

###  3.1 Autoregressive ASR with Hybrid CTC/Attention

AR decoding is a recursive method for estimating target sequences. In the
example shown in Fig. 1(a), the token `s` is estimated first and then used to
estimate the next token, `e`. This simple left-to-right beam search similar to
[15] is widely used in ASR to search for the most likely transcriptions [16].

However, due to the iterative nature of beam search, this can greatly slow the
inference time of AR decoding. For example, using the hybrid CTC/attention
architecture, the beam search requires decoder computation and CTC rescoring
process for each search iteration, hence inference time increases if the
search process is iterated a large number of times. This effect is shown in
Fig. 2, which shows the proportion of time spent on the encoder and decoding
process, including decoder and CTC prefix score computation, as well as the
average time spent during inference for various lengths of audio input. This
shows that as the audio length and inference time increase, the proportion of
time spent in the decoding process also increases.

###  3.2 Non-autoregressive ASR with Mask-CTC

NAR decoding is a method that avoids recursively estimating the target
sequences to address the problem of slow inference. There are various
architectures for this method [17], but this work focuses on the Mask-CTC [13]
model. In the Mask-CTC model, we first estimate the target sequences with
greedy CTC decoding (gCTC) output, and then a mask is applied based on the CTC
probabilities for each token. As illustrated in Fig. 1(b), the token `e`, `a`,
and `m` is masked due to its lower probability. Then, the masked token `#` is
estimated using the masked language model decoder. Since the decoder is only
required for the masked token, and the number of decoder calculations is
fixed, the number of decoder computations is significantly reduced compared to
the AR method.

Fig. 2 illustrates the proportion of time spent on the encoder and decoder,
and the inference time for the Mask-CTC model. We set the number of the
decoder iterations to $10$ for measurement. Compared to Fig. 2, we can see
that the encoder’s share of NAR is larger than AR. Considering the number of
encoder computations is always $1$, a large encoder’s share means a shorter
computation time. As audio length increases, the inference time also
increases, but the impact on inference time is small. Therefore the difference
in inference time between AR and NAR can be seen as the difference in the
proportion of the encoder and decoding process.

However, there is an accuracy issue with the Mask-CTC-based NAR decoding. If
the number of masked tokens is different from the actual number of tokens, the
accuracy degrades significantly. For example, if the correct sequence is `s,
e, a` and the masked sequence is `s, #, #, a`, the result of Mask-CTC becomes
incorrect, since it tends to assign the same numbers of tokens to the masks,
e.g., `s, e, e, a` or `s, e, a, a`. Similarly, if the masked sequence is `s,
e, #, a`, it will have an insertion error. In our preliminary experiments, we
observed that almost 40% of the masked sequences do not match the proper
length. Additionally, we have observed that insertion errors caused by the
absence of a target token, occur with a probability of 2%.

##  4 Partially Autoregressive framework

###  4.1 Partially Autoregressive Inference

To address the issues inherent in AR and NAR decoding, we propose a partially
autoregressive (PAR) decoding method. Both AR and NAR models use the CTC and
decoder components, but there is a significant difference in the trade-off
balance based on usage. AR uses an iterative process to predict the target
sequence, so it does not have an accuracy issue with target sequence length.
NAR first uses gCTC results to reduce the number of tokens that need to be
predicted with the decoder, resulting in fast inference. In PAR, we combine
NAR-style CTC usage and AR-style decoder usage to utilize these two advantages
fully.

The architecture of PAR is illustrated in Fig. 1(c). PAR first generates a
sequence of tokens by gCTC approach, then apply the mask process using the
Mask-CTC method [13]. Finally, it predicts the tokens corresponding to a mask
token using beam search, similarly to the AR approach. To reduce the
computational complexity, we propose the segment-level vectorized beam search,
which significantly reduces the number of search iterations. By applying this
beam search, we can solve the NAR accuracy issue related to the incorrect
target length.

As an example, let us focus on a case where only one mask is in the sequence
for simplicity, and let the $P_{\mathrm{thres}}$ represent the threshold value
ranging from $0$ to $1$. Initially, we use the gCTC result, obtained without
the AR process. However, the resulting text can contain errors due to the
conditional independence assumption. We expect these errors can be corrected
by the AR process, where we utilize the pre-trained decoder and beam search.
To determine which tokens to update using the AR process, we use
$P_{\mathrm{thres}}$ as a filter for posterior probability. Tokens with
probability less than $P_{\mathrm{thres}}$ will be considered less confident
and replaced with the mask token. Consecutive mask tokens are merged into a
single mask because we will estimate the sequences in the AR process. Finally,
we use beam search to estimate the tokens.

We iterate the beam search for max_iteration times for any audio length, where
max_iteration denotes the maximum number of tokens for one mask. The value of
max_iteration is highly dependent on the threshold $P_{\mathrm{thres}}$, where
a higher value is required when $P_{\mathrm{thres}}$ is closer to 1 and more
tokens are replaced with mask tokens.

Algorithm 1 Segment-level Vectorized Beam Search

Beam size: $B$ and the number of masks: $S$

Encoder output: $X$

Initialize decoder cache: $C$

Initialize a $S$-length list for ended hypotheses: $F_{S}$

Initialize hypotheses $Y_{S,B}$ and end-of-sequences $E_{S}$

for $i=1$ to $max\\_iteration$ do

Calculate probability $prob=Decoder(Y_{S,B},X,C)$

Update hypotheses $Y_{S,B}$ by top-k method.

for $s=1$ to $S$ do

for $b=1$ to $B$ do

if last token of $Y_{s,b}$ is $E_{s}$ then

$F_{s}.push(Y_{s,b})$

end if

end for

end for

end for

for $s=1$ to $S$ do

Replace $s$-place mask in the masked sequence with best hypothesis in $F_{s}$

end for

###  4.2 Segment-level Vectorized Beam Search

The beam search process explained in the previous section focused on the case
where there is a single mask token. However, in practice, multiple mask tokens
may exist for a given sequence. To handle multiple masks, we parallelize the
decoding of masks by extending the vectorized beam search[23]. The vectorized
beam search is an extension of traditional beam search that allows for the
calculation of $B$ beams simultaneously. In an offline scenario, it can also
be extended to the parallel computation of $S$ utterances. In other words,
$S\times B$ hypotheses can be calculated as a single batch at each step. In
our work, we treat the value $S$ as the number of masks to enable parallel
computation of all mask tokens.

The overview of the proposed segment-level vectorized beam search is described
in Algorithm 1. First, we initialize $S\times B$ hypotheses $Y_{S,B}$ from the
gCTC result and masked sequence. Initially, there is only one hypothesis, so
$Y_{s,1}$ contains the gCTC result up to the corresponding mask as the
hypothesis. The rest of the hypotheses, from $Y_{s,2}$ to $Y_{s,B}$, are
initialized with dummy hypotheses. Additionally, we store the next token for
each mask in an end-of-sequence list, $E_{S}$. Next, we calculate the
probabilities of the next token for each hypothesis and update $Y_{S,B}$. We
then check, for each mask $s$, if $E_{s}$ is predicted as the next token for
each of the hypothesis from $Y_{s,1}$ to $Y_{s,B}$. We push the ended
hypothesis to the list of ended hypotheses, ${F_{s}}$. Finally, after
iterating $max\\_iteration$ times, we replace each mask token with the best
hypothesis in $F_{s}$. To simplify the implementation, we apply padding to the
missing hypotheses so that the number of hypotheses remains constant, similar
to [34].

Table 1 shows an example of the decoding process of a sample in the
LibriSpeech test-clean set. After calculating the gCTC result and merging the
consecutive masks, there are three masks in this example. Since the beam size
$B$ is set to $10$, there are $30$ hypotheses. In the first iteration, we
correctly estimated the first and third mask tokens. Since we used the BPE
token for this model, the red characters `fall`, `bra`, and `s` were estimated
in a single iteration. In the second iteration, we successfully predicted the
second mask, which has two corresponding tokens: `bro` and `th`. This example
demonstrates that it is possible to predict multiple tokens from a single mask
token, enabling us to handle the accuracy issue in NAR due to incorrect target
length. Moreover, the overall number of iterations in beam search is
significantly lower than in AR by utilizing the segment-level vectorized beam
search, indicating that PAR is effective in avoiding AR issues.

We determine whether to end the beam search by observing the next token of the
mask token. For example, in the example in Fig. 1(c), when the first mask
detects a space token `_` and the second mask detects `b`, the hypothesis is
terminated and pruned from the overall search.

Table 2: Dataset descriptions. The order of evaluation sets corresponds to the
error results in Table 3. †Data split based on ESPnet [35].

Dataset | Language | Hours | Token | Metric | Evaluation Sets
---|---|---|---|---|---
AISHELL-1 [36] | zh | 170 | char | CER | dev / test
JSUT [37] | ja | 10 | char | CER |  dev / test ($\dagger$)
LS-100 [38] | en | 100 | BPE | WER | dev-{clean,other} / test-{clean,other}
LS-960 [38] | en | 960 | BPE | WER | dev-{clean,other} / test-{clean,other}
TED-LIUM2 [39] | en | 210 | BPE | WER | dev / test

Table 3: Comparison of WER, CER, RTF, elapsed time, and memory usage for each
model. RTF, inference time, and memory usage compare mean and standard
deviation (in brackets) values. The symbols ($\downarrow$) and ($\uparrow$)
indicate a lower or higher number is preferable, respectively.

| AR |  | PAR |
---|---|---|---|---
Dataset |  RTF ($\downarrow$) |  Error [%] ($\downarrow$) | Memory Usage [MB] |  |  RTF ($\downarrow$) |  Error [%] ($\downarrow$) | Memory Usage [MB] |  Speedup ($\uparrow$)
AISHELL-1 | 0.027 (0.006) | 4.6 / 5.0 | 176.3 (4.2) |  | 0.010 (0.006) | 4.6 / 4.9 | 176.3 (5.9) |  2.70$\times$
JSUT | 0.129 (0.021) | 11.8 / 13.2 | 208.4 (6.2) |  | 0.018 (0.008) | 12.0 / 13.3 | 199.4 (6.4) |  7.17$\times$
LS-100 | 0.111 (0.038) | 6.2 / 16.8 / 6.4 / 17.1 | 197.5 (25.3) |  | 0.009 (0.006) | 6.5 / 17.2 / 6.7 / 17.7 | 210.6 (81.6) |  12.33$\times$
LS-960 | 0.110 (0.039) | 2.2 / 5.2 / 2.5 / 5.2 | 528.4 (44.8) |  | 0.008 (0.006) | 2.2 / 5.5 / 2.5 / 5.6 | 550.5 (132.6) |  13.75$\times$
TED-LIUM2 | 0.182 (0.060) | 7.3 / 7.1 | 281.1 (47.7) |  | 0.012 (0.008) | 7.6 / 7.3 | 474.8 (298.9) |  15.17$\times$

##  5 Experiment

###  5.1 Experimental setup

####  5.1.1 Models

To test the effectiveness of PAR, we compare both its performance and
inference speeds of various algorithms to AR and NAR. As PAR decoding can be
seen as an optimization of AR decoding, PAR can be tested using pre-trained AR
models. As a result, for both AR and PAR inference, we used E-Branchformer
models [6] that were pre-trained on the following datasets and are publicly
available on the HuggingFace hub:
AISHELL-1111https://huggingface.co/pyf98/aishell_e_branchformer,
JSUT222https://huggingface.co/pyf98/jsut_e_branchformer
LibriSpeech-100h333https://huggingface.co/pyf98/librispeech_100_e_branchformer,
LibriSpeech-960h444https://huggingface.co/asapp/e_branchformer_librispeech,
TED-LIUM2555https://huggingface.co/pyf98/tedlium2_e_branchformer. Note that we
used two types of models trained with LibriSpeech: a model trained with the
entire dataset of $960$ hours of audio (LS-960), and another model trained
with the subset of $100$ hours of audio (LS-100). The model trained with the
LS-100 dataset is considered less accurate than the one trained with LS-960,
so we used LS-100 to investigate the effect of gCTC accuracy on PAR decoding.
For comparison against NAR decoding, we trained models using the Conformer [5]
architecture on LS-100. We prepared AR and NAR models with approximately 30
million parameters for a fair evaluation. We used the Conformer model to test
if there would be any differences caused by encoder architecture. All the
models were trained and evaluated using the ESPnet toolkit [35].

####  5.1.2 Decoding setup

For AR decoding, we set the beam size to $10$, the CTC weight to $0.3$, and
employed the vectorized beam search [23]. For PAR decoding, the beam size is
set to $10$, but note that for PAR inference, we do not compute the CTC prefix
score, so the CTC weight is set to $0$. We also set the $P_{\mathrm{thres}}$
to $0.95$ and max_iteration to $5$ for PAR decoding. The decoding speed was
measured using a single RTX 2080 Ti.

####  5.1.3 Evaluation Datasets

We used several datasets of different languages for evaluation, described in
Table 2. For this study, we used LibriSpeech and TED-LIUM2 as English
datasets, JSUT as a Japanese dataset, and AISHELL-1 as a Mandarin dataset. For
LibriSpeech, we evaluated dev and test sets, each containing clean and other
sets. Each of the five evaluation datasets was evaluated using the AR/PAR
model trained on the equivalent dataset. For example, the E-branchformer model
pre-trained on the AISHELL-1 dataset was evaluated using the AISHELL-1
evaluation set.

####  5.1.4 Metrics

We evaluated accuracy with word error rate (WER) or character error rate (CER)
according to the dataset. We used the real-time factor (RTF) to measure the
inference speed. The maximum memory usage of the GPU during inference is also
measured to observe the effect of mask parallel decoding. For the RTF and
memory usage, we took the mean of all evaluation sets. The speedup column
compares the RTF values with PAR compared to AR.

Fig. 3: Average inference time and proportion of time spent on the encoder and
decoder computation during PAR decoding. The decoder’s share is greatly
reduced from AR.

###  5.2 Results

####  5.2.1 AR and PAR

The evaluation results comparing AR and PAR are shown in Table 3. Focusing on
accuracy, it is evident that compared to AR, PAR shows a similar WER or CER
for all models and evaluation datasets. The beam search process can properly
update the tokens that gCTC cannot accurately estimate. While the accuracy has
slightly degraded, this is due to the accuracy of the gCTC result. We will
explain it further in detail in the following limitations section 5.3.1. In
terms of the RTF, we have achieved approximately $10$ times faster inference
compared to AR. In particular, since the inference speed does not largely
depend on the audio length, the standard deviation is approximately $17.5\%$
of that of AR decoding. This feature is more effective for longer audio and in
the test-clean dataset, and we observed $89.7\times$ speed up for $29$ second
audio as the maximum speedup. Note that this speedup depends on several
factors, such as the number of masks or input audio length.

Fig. 3 shows the proportion of time spent on the encoder and decoder process
and the inference time for each audio length, evaluated on the LS-960 dataset.
Compared to Fig. 2, the decoder-to-encoder ratio is significantly reduced.
Furthermore, we can observe that the increase in inference time with respect
to audio length is small, similar to that of NAR. This is because, in this
work, we set the max_iteration to $5$, which means that the number of beam
search iterations is limited to a maximum of 5. Moreover, if there are no
masks, the inference is simply a gCTC result, where we do not run decoder
computation. Therefore, the AR process does not have a high computational cost
for long audio inputs. We can confirm that the inference speed of PAR depends
on the computation speed of each model, such as the encoder or decoder,
considering that the inference time slightly increases with audio length and
the number of search iterations is fixed. The computation speed of the encoder
or decoder depends on the audio length, but its impact on the overall
inference speed is small.

We investigated the correlation between the WER and RTF as shown in Fig. 4, by
changing the beam size from $1$ to $20$, and measured the WER and RTF using
the test-clean dataset from Librispeech on the E-Branchformer based pre-
trained models for LS-100 and LS-960 datasets. Comparing the AR models, the
RTF differs greatly because the model sizes are different, however, the RTF of
the PAR method remained almost the same. This is because the inference time of
PAR depends on the computation time of the models itself. Therefore, the PAR
curves have a significantly lower RTF compared to the AR curves, and the
change in RTF is negligible as the beam-size changes, hence the narrower
width. Although the RTF remains constant, we still see the differences in WER
for the PAR models, especially with the pre-trained model of the LS-100
dataset compared to the model for the LS-960 dataset. This is due to the
accuracy degradation problem we mention in Section 5.3.1.

Fig. 4: The comparison of WER and RTF measured using the AR and PAR methods.
We used the models trained with LS-100 and LS-960 datasets and measured by
changing the beam size between $1$ and $20$.

####  5.2.2 NAR and PAR

Comparing NAR and PAR, we can see that PAR is not as fast as NAR, as shown in
Table 4. The number of decoder iterations is $10$, which is larger than
max_iteration for PAR, but it performs faster than PAR. One reason for this is
that the size of the decoder input is different. In this work, we added a
dummy hypothesis as mentioned in Section 4.2. As a result, the batch size of
decoder input for PAR decoding is $S\times B$, where the batch size for NAR
decoding is $S$. Therefore, it seems that the computation time per decoder is
shorter for Mask-CTC NAR, and the resulting decoder processing time becomes
faster.

In terms of accuracy, we can see that PAR outperforms NAR. Improvement in
performance is similar to that of E-branchformer models in both RTF and WER,
and it was confirmed that there was no difference in improvement due to
architecture differences. From these results, it is evident that applying the
PAR method can solve the built-in accuracy issue for the NAR method we
mentioned in Section 3.2. In particular, we can confirm that PAR obtains a
speed between AR and NAR but achieves similar accuracy to AR. Therefore, a new
trade-off balance that is not present in AR and NAR has been realized.

Table 4: Comparison with NAR. The speedup shows the speedup from the AR decoding. | Model |  RTF ($\downarrow$) |  WER [%] ($\downarrow$) |  Speedup ($\uparrow$)
---|---|---|---|---
AR | CTC/Attention | 0.198 (0.080) | 6.7 / 18.3 / 7.0 / 18.6 |  1.00$\times$
NAR | CTC | 0.005 (0.004) | 7.7 / 21.0 / 7.9 / 21.4 |  39.60$\times$
Mask-CTC | 0.008 (0.005) | 7.1 / 20.8 / 7.5 / 21.0 |  24.75$\times$
PAR | CTC/Attention | 0.014 (0.009) | 6.2 / 18.5 / 6.6 / 18.7 |  14.14$\times$

###  5.3 Limitations

####  5.3.1 Accuracy with PAR

The accuracy of PAR can be degraded if the gCTC result is not accurate. If the
result of gCTC is incorrect with high confidence, we cannot use the AR process
to refine the gCTC result. From the comparison of the LS-100 and LS-960 models
in Table 3, we can see that more accuracy degradation can occur with the
LS-100 model.

The accuracy may also degrade at a higher $P_{\mathrm{thres}}$ because the
number of target tokens per mask may exceed max_iteration. Since we stop the
beam search after max_iteration iterations, if the number of target tokens
exceeds max_iteration, we cannot predict the entire sequence for one mask.
Therefore, the accuracy may be degraded if the $P_{\mathrm{thres}}$ is closer
to $1.0$. Fig. 5 describes the relationship between WER and
$P_{\mathrm{thres}}$ for max_iteration. Using a max_iteration of $5$, we
observe the accuracy degrades at higher $P_{\mathrm{thres}}$ levels. This
issue is caused by the lack of beam search iterations, so to solve this
problem, we need to increase the max_iteration. In Fig. 5, we increased the
max_iteration to $8$ and observed a more accurate result.

Fig. 5: The relationship between the WER and $P_{\mathrm{thres}}$. We
evaluated by changing the $P_{\mathrm{thres}}$ from $0.95$ to $0.999$. We used
the E-Branchformer-based pre-trained model for LS-960.

####  5.3.2 Memory usage

It is important to note the standard deviation of memory usage in Table 3
increases greatly compared to AR. Since the decoder process in PAR is computed
simultaneously for all masks, we need more GPU memory if there are many masks.
Therefore if the number of masks increases due to long audio inputs, high
$P_{\mathrm{thres}}$, or low accuracy of the gCTC result, we may get an out-
of-memory error as GPU memory is exceeded during inference. Considering that
the inference of all masks does not depend on each other, it is possible to
alleviate this issue by using multiple GPUs to perform inference.

####  5.3.3 Segment-level Vectorized Beam Search

If a masked sequence contains multiple masks, the predicted tokens for the
second or later masks may not be accurate due to the inaccuracy of the tokens
predicted for the first mask. For example, in Fig. 1(c), the masked sequence
is `se#_cuc#ber`, where we have two masks at the positions of tokens `e` and
`am`. Since we use the gCTC result to predict masked tokens, the decoder input
for estimating the mask for the `am` part becomes `see_cuc`, which is
incorrect, instead of the correct `sea_cuc`. This incorrect decoder input
might impact the accuracy of the prediction. However, considering that even
the AR decode may also have an incorrect decoder input, we believe PAR can
search the same hypothesis with AR by utilizing the beam search.

##  6 Conclusion

In this work, we propose a partially autoregressive framework to obtain a new
trade-off balance between accuracy and latency. With our novel architecture
design to compensate for the weakness inherent in each NAR and AR, PAR is a
decoding method that takes advantage of the strengths of these two methods. In
our experiments, we observed that the AR model can be inferred at NAR-level
speeds without sacrificing accuracy. Notably, the LS-960 pre-trained model
achieved a $13.75\times$ speedup with the same WER on the test-clean set. We
believe that using PAR can significantly improve the usability of the
traditional hybrid CTC/Attention model. Our framework has limitations, such as
memory usage issues, that restrict the scenarios where PAR can be applied. In
future work, we plan to extend our framework to edge devices with limited
computational resources.

## References

  * [1] Geoffrey Hinton, Li Deng, Dong Yu, George E. Dahl, Abdel-rahman Mohamed, Navdeep Jaitly, Andrew Senior, Vincent Vanhoucke, Patrick Nguyen, Tara N. Sainath, and Brian Kingsbury,  “Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups,”  IEEE Signal Process. Mag., vol. 29, no. 6, pp. 82–97, 2012.
  * [2] Alex Graves, Abdelrahman Mohamed, and Geoffrey Hinton,  “Speech recognition with deep recurrent neural networks,”  in Proc. ICASSP, 2013, pp. 6645–6649.
  * [3] Rohit Prabhavalkar, Takaaki Hori, Tara N Sainath, Ralf Schlüter, and Shinji Watanabe,  “End-to-end speech recognition: A survey,”  arXiv preprint arXiv:2303.03329, 2023.
  * [4] Alex Graves,  “Sequence transduction with recurrent neural networks,”  in Proc. ICML, 2012.
  * [5] Anmol Gulati, James Qin, Chung-Cheng Chiu, Niki Parmar, Yu Zhang, Jiahui Yu, Wei Han, Shibo Wang, Zhengdong Zhang, Yonghui Wu, and Ruoming Pang,  “Conformer: Convolution-augmented transformer for speech recognition,”  in Proc. Interspeech, 2020, pp. 5036–5040.
  * [6] Kwangyoun Kim, Felix Wu, Yifan Peng, Jing Pan, Prashant Sridhar, Kyu J. Han, and Shinji Watanabe,  “E-Branchformer: Branchformer with enhanced merging for speech recognition,”  in Proc. SLT, 2023, pp. 84–91.
  * [7] Alec Radford, Jong Wook Kim, Tao Xu, Greg Brockman, Christine McLeavey, and Ilya Sutskever,  “Robust speech recognition via large-scale weak supervision,” 2022.
  * [8] Hirofumi Inaguma, Shun Kiyono, Kevin Duh, Shigeki Karita, Nelson Yalta, Tomoki Hayashi, and Shinji Watanabe,  “ESPnet-ST: All-in-one speech translation toolkit,”  in Proc. ACL, Online, 06 2020, pp. 302–311, Association for Computational Linguistics.
  * [9] Jiatao Gu, James Bradbury, Caiming Xiong, Victor OK Li, and Richard Socher,  “Non-autoregressive neural machine translation,”  in Proc. ICML, 2018.
  * [10] Yisheng Xiao, Lijun Wu, Junliang Guo, Juntao Li, Min Zhang, Tao Qin, and Tie-yan Liu,  “A survey on non-autoregressive generation for neural machine translation and beyond,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023\.
  * [11] Nanxin Chen, Shinji Watanabe, Jesús Villalba, Piotr Żelasko, and Najim Dehak,  “Non-autoregressive transformer for speech recognition,”  IEEE Signal Processing Letters, vol. 28, pp. 121–125, 2021.
  * [12] William Chan, Chitwan Saharia, Geoffrey Hinton, Mohammad Norouzi, and Navdeep Jaitly,  “Imputer: Sequence modelling via imputation and dynamic programming,”  in Proc. ICML, 2020, pp. 1403–1413.
  * [13] Yosuke Higuchi, Shinji Watanabe, Nanxin Chen, Tetsuji Ogawa, and Tetsunori Kobayashi,  “Mask CTC: Non-Autoregressive End-to-End ASR with CTC and Mask Predict,”  in Proc. Interspeech, 2020, pp. 3655–3659.
  * [14] Ethan A Chi, Julian Salazar, and Katrin Kirchhoff,  “Align-Refine: Non-autoregressive speech recognition via iterative realignment,”  in Proc. NAACL-HLT, 2021, pp. 1920–1927.
  * [15] Ilya Sutskever, Oriol Vinyals, and Quoc V Le,  “Sequence to sequence learning with neural networks,”  in Proc. NeurIPS, Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, Eds. 2014, vol. 27, Curran Associates, Inc.
  * [16] William Chan, Navdeep Jaitly, Quoc Le, and Oriol Vinyals,  “Listen, attend and spell: A neural network for large vocabulary conversational speech recognition,”  in Proc. ICASSP, 2016, pp. 4960–4964.
  * [17] Yosuke Higuchi, Nanxin Chen, Yuya Fujita, Hirofumi Inaguma, Tatsuya Komatsu, Jaesong Lee, Jumon Nozaki, Tianzi Wang, and Shinji Watanabe,  “A comparative study on non-autoregressive modelings for speech-to-text generation,”  in Proc. ASRU, 2021, pp. 47–54.
  * [18] Yann LeCun, John Denker, and Sara Solla,  “Optimal brain damage,”  in Proc. NeurIPS, 01 1989, vol. 2, pp. 598–605.
  * [19] Frankle Jonathan and Carbin Michael,  “The lottery ticket hypothesis: Finding sparse, trainable neural networks,”  in Proc. ICLR, 2019.
  * [20] Cheng-I Jeff Lai, Yang Zhang, Alexander H Liu, Shiyu Chang, Yi-Lun Liao, Yung-Sung Chuang, Kaizhi Qian, Sameer Khurana, David Cox, and Jim Glass,  “PARP: Prune, adjust and re-prune for self-supervised speech recognition,”  Proc. NeurIPS, vol. 34, pp. 21256–21272, 2021.
  * [21] Geoffrey Hinton, Oriol Vinyals, and Jeffrey Dean,  “Distilling the knowledge in a neural network,”  in NIPS Deep Learning and Representation Learning Workshop, 2015\.
  * [22] Heng-Jui Chang, Shu-wen Yang, and Hung-yi Lee,  “DistilHuBERT: Speech representation learning by layer-wise distillation of hidden-unit bert,”  in Proc. ICASSP, 2022, pp. 7087–7091.
  * [23] Hiroshi Seki, Takaaki Hori, Shinji Watanabe, Niko Moritz, and Jonathan Le Roux,  “Vectorized beam search for CTC-attention-based speech recognition,”  in Proc. Interspeech, 2019, pp. 3825–3829.
  * [24] Jiahui Yu, Chung-Cheng Chiu, Bo Li, Shuo yiin Chang, Tara N Sainath, Yanzhang (Ryan) He, Arun Narayanan, Wei Han, Anmol Gulati, Yonghui Wu, and Ruoming Pang,  “FastEmit: Low-latency streaming ASR with sequence-level emission regularization,”  in Proc. ICASSP, 2021.
  * [25] Zhuoyuan Yao, Di Wu, Xiong Wang, Binbin Zhang, Fan Yu, Chao Yang, Zhendong Peng, Xiaoyu Chen, Lei Xie, and Xin Lei,  “WeNet: Production oriented streaming and non-streaming end-to-end speech recognition toolkit,”  in Proc. Interspeech, 2021, pp. 4054–4058.
  * [26] Binbin Zhang, Di Wu, Zhendong Peng, Xingchen Song, Zhuoyuan Yao, Hang Lv, Lei Xie, Chao Yang, Fuping Pan, and Jianwei Niu,  “WeNet 2.0: More productive end-to-end speech recognition toolkit,”  in Proc. Interspeech, 2022, pp. 1661–1665.
  * [27] Masao Someki, Yosuke Higuchi, Tomoki Hayashi, and Shinji Watanabe,  “Espnet-ONNX: Bridging a gap between research and production,”  in Proc. APSIPA ASC, 2022, pp. 420–427.
  * [28] Yosuke Higuchi, Hirofumi Inaguma, Shinji Watanabe, Tetsuji Ogawa, and Tetsunori Kobayashi,  “Improved Mask-CTC for non-autoregressive end-to-end ASR,”  in Proc. ICASSP, 2021, pp. 8363–8367.
  * [29] Hayato Futami, Hirofumi Inaguma, Sei Ueno, Masato Mimura, Shinsuke Sakai, and Tatsuya Kawahara,  “Non-autoregressive error correction for CTC-based ASR with phone-conditioned masked LM,”  in Proc. Interspeech, 2022, pp. 3889–3893.
  * [30] Yosuke Higuchi, Brian Yan, Siddhant Arora, Tetsuji Ogawa, Tetsunori Kobayashi, and Shinji Watanabe,  “BERT meets CTC: New formulation of end-to-end speech recognition with pre-trained masked language model,”  in Proc. Findings of EMNLP, 2022, pp. 5486–5503.
  * [31] Jay Mahadeokar, Yangyang Shi, Ke Li, Duc Le, Jiedan Zhu, Vikas Chandra, Ozlem Kalinli, and Michael L. Seltzer,  “Streaming parallel transducer beam search with fast-slow cascaded encoders,”  in Proc. Interspeech, 2022.
  * [32] Zhengkun Tian, Jiangyan Yi, Jianhua Tao, Shuai Zhang, and Zhengqi Wen,  “Hybrid autoregressive and non-autoregressive transformer models for speech recognition,”  IEEE Signal Process. Lett., vol. 29, pp. 762–766, 2022.
  * [33] Shinji Watanabe, Takaaki Hori, Suyoun Kim, John R Hershey, and Tomoki Hayashi,  “Hybrid CTC/attention architecture for end-to-end speech recognition,”  IEEE Journal of Selected Topics in Signal Processing, vol. 11, no. 8, pp. 1240–1253, 2017.
  * [34] Liu Xiaoyu, Lin Eric, and Ning Emma,  “Journey to optimize large scale transformer model inference with ONNX runtime,” https://cloudblogs.microsoft.com/opensource/2021/06/30/journey-to-optimize-large-scale-transformer-model-inference-with-onnx-runtime, 2021,  [Online; Accessed on June-27-2023].
  * [35] Shinji Watanabe, Takaaki Hori, Shigeki Karita, Tomoki Hayashi, Jiro Nishitoba, Yuya Unno, Nelson Enrique Yalta Soplin, Jahn Heymann, Matthew Wiesner, Nanxin Chen, Adithya Renduchintala, and Tsubasa Ochiai,  “ESPnet: End-to-end speech processing toolkit,”  in Proc. Interspeech, 2018, pp. 2207–2211.
  * [36] Hui Bu, Jiayu Du, Xingyu Na, Bengu Wu, and Hao Zheng,  “AISHELL-1: An open-source Mandarin speech corpus and a speech recognition baseline,”  in Proc. O-COCOSDA, 2017, pp. 1–5.
  * [37] Ryosuke Sonobe, Shinnosuke Takamichi, and Hiroshi Saruwatari,  “JSUT corpus: free large-scale Japanese speech corpus for end-to-end speech synthesis,”  arXiv preprint arXiv:1711.00354, 2017.
  * [38] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur,  “Librispeech: An ASR corpus based on public domain audio books,”  in Proc. ICASSP, 2015, pp. 5206–5210.
  * [39] Anthony Rousseau, Paul Deléglise, and Yannick Estève,  “Enhancing the TED-LIUM corpus with selected data for language modeling and more TED talks,”  in Proc. LREC, 2014, pp. 3935–3939.


# Critical correlations of Ising and Yang-Lee critical points from Tensor RG

Sankhya Basu  Physics program and Initiative for the Theoretical Sciences, The
Graduate Center, CUNY, New York, NY 10016, USA  Department of Physics and
Astronomy, College of Staten Island, CUNY, Staten Island, NY 10314, USA  Vadim
Oganesyan  Physics program and Initiative for the Theoretical Sciences, The
Graduate Center, CUNY, New York, NY 10016, USA  Department of Physics and
Astronomy, College of Staten Island, CUNY, Staten Island, NY 10314, USA
Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue,
New York, New York 10010, USA

###### Abstract

We examine feasibility of accurate estimations of universal critical data
using tensor renormalization group (TRG) algorithm introduced by Levin and
Nave LevinNaveTRG . Specifically, we compute critical exponents
$\gamma,\gamma/\nu,\delta,\eta$ and amplitude ratio $A$ for the magnetic
susceptibility from one- and two-point correlation functions for three
critical points in two dimensions – isotropic and anisotropic Ising models and
the Yang-Lee critical point at finite imaginary magnetic field. While TRG
performs quantitaviely well in all three cases already at smaller bond
dimension, $D=16$, the latter two appear to show more rapid improvement in
bond dimension, e.g. we are able to reproduce exactly known results to better
than one percent at bond dimension $D=24$. We comment on the relationship
between these results and earlier results on conformal dimensions, fixed
points of tensor RG, and also compare computational costs of tensor
renormalization vs. conventional Monte-Carlo sampling.

##  I Introduction

Formulation and subsequent development of stochastic (Monte-Carlo) approaches
to the many-body problem is arguably the most significant advancement in
computational science to-date. As envisaged in the seminal paper by Metropolis
and coauthors metropolis1 ; metropolis2 , the method became a defacto general
purpose tool, and a must-include topic in even basic coursework. Although path
integral formulation of quantum mechanics enabled significant progress in
adapting stochastic techniques to studying quantum problems, the inherent
limitations of this approach continue to stimulate research into alternatives.
Perhaps the most successful of these is Density-matrix renormalization group
(DMRG) in one-dimensional problems Schollwock2005 ; DMRG92  and related tensor
extensions to higher dimensions Verstraete2004 , whose closely related time
evolution variants (TEBD, tDMRG) can access short time dynamics – entanglement
growth severely limits the extent of excitation energy and/or time amenable to
DMRG techniques in 1D (and, more severely, 2 and 3 dimensional counterparts).
DMRG and other variational schemes built around correctly capturing patterns
of quantum entanglement remain the go-to method for low dimensional quantum
many-body problems. As it turns out some of the basic insights and tools for
distilling/managing entanglement are useful for designing or refining other
approximation schemes for many-body correlations. An example of that is the
tensor renormalization group (TRG) algorithm LevinNaveTRG  – a systematic,
quantitatively accurate way to implement the celebrated idea of real space
coarse graining transformation due to Migdal and Kadanoff Migdal1 ;
MigdalKadanoff . TRG essentially solves the bookkeeping problem inherent in
the Migdal-Kadanoff scheme, relying on singular value decompositions to
efficiently compress interactions among block-spins (see Sec. II.1 for further
details). Importantly, the TRG scheme is not restricted to purely statistical
problems (characterized by positive Gibbs weights), e.g. in our recent work
Basuetal2021  we discovered symmetry broken steady states in non-unitary
dynamics for one dimensional Ising chains, intimately connected with analytic
continuation of equilibrium 2D Ising model to complex temperature, and also
the complex field Yang-Lee problem GarciaWeiYL ; YangLee52 ; LeeYang52  we
study in this work.

While several improvements have been developed SRG_PRL ; SRG_PRB ; HOTRG ;
EvenblyTNR ; LoopTNR , the plain-vanilla TRG scheme LevinNaveTRG  may be
thought of as analogous to the basic Metropolis algorithm; both of them
establishing the essential idea of the method, demonstrating its proof-of-
principle usefulness _and_ remaining the simplest working implementation of
the method sufficient for a quick start on a given problem. TRG does not
suffer from critical slowing down that hampers conventional Monte-Carlo,
rather being limited by comparatively more subtle problems. Much of the
existing literature on improving TRG focuses on its inability to fully remove
short range correlations especially near critical points; however this
limitation appears not to hamper the ability to approximate _universal_
properties in a convergent fashion, as suggested by our results.

The main thrust of this paper is to establish this point systematically by
computing and analysing two-point correlation functions for large separations.
Our direct computations readily yield estimates of critical exponents on the
order of one percent or even better, consistent with some of the previous
works on tensor renormalization interpretting finite size effects using a
conformal field theory ansatzTEFR ; EvenblyTNR ; LoopTNR , but well outside of
what typical Monte Carlo studies are capable of. Importantly, we benchmark TRG
away from well trodden isotropic 2D Ising model and find much clearer evidence
for convergence in bond dimension. We also aim to promote a view that for many
reasons having to do with ease of implementation, power to access exceedingly
large distances _and_ potential versatility of being adaptable to study
dynamical problems, TRG should be viewed as a viable alternative to Monte-
Carlo, and introduced to students early on. With that in mind, the next
Section II summarizes the basic ingredients of the analysis in a self-
contained manner, including a detailed disucssion of errors. Section III
contains systematic analysis of critical exponents $\gamma$, $\gamma/\nu$,
$\eta$, $\delta=1/\sigma$ as well as the susceptibility amplitude ratio $A$
for isotropic and anisotropic square lattice Ising models, and critical
exponents $\sigma$ and $\eta$ at the Yang-Lee critical point. We conclude with
a summary and a few afterthoughts in Sec. IV.

##  II Methods

###  II.1 Operator Renormalization with TRG

Figure 1: Step 1 of TRG. Replacing microscopic degrees of freedom (DOF) living
on verices in favor of tensor network and edge DOF. In this implementation
each bond’s transfer matrix $M$ (panel a) is decomposed using SVD (panel b)
with components recombined into tensors $T$ (panel c).

The goal of all TRG (and its descendants) is to express partition sums as
products of factors ($S_{j}=\pm 1$)

$\displaystyle Z_{0}(\beta,J)=Tr\left[e^{-\beta H}\right]=\Pi_{j}z_{0,j}$ (1)
$\displaystyle Z_{1}(\beta,J,r)=Tr\left[S_{r}e^{-\beta
H}\right]=\Pi_{j}z_{1,j}$ (2) $\displaystyle
Z_{2}(\beta,J,r_{1},r_{2})=Tr\left[S_{r_{1}}S_{r_{2}}e^{-\beta
H}\right]=\Pi_{j}z_{2,j},\rm{etc}$ (3)

at inverse temperature $\beta$ (setting $k_{B}$ to unity), with the standard
nearest neighbor spin-spin interactions ($J_{ij}$) and field ($h$) and
Hamiltonian

$H=-\sum_{\langle i,j\rangle}J_{ij}S_{i}S_{j}-h\sum_{i}S_{i}\;.$ (4)

Free energy, magnetization and the two-point correlator can be computed
straightforwardly

$\displaystyle-\beta f$ $\displaystyle\equiv\log Z_{0}/L^{2}=\sum_{j}\log
z_{0,j}/L^{2},$ (5) $\displaystyle m$ $\displaystyle\equiv
Z_{1}/Z_{0}=\Pi_{j}z_{1,j}/z_{0,j},$ (6) $\displaystyle C(r_{1},r_{2})$
$\displaystyle\equiv Z_{2}/Z_{0}=\Pi_{j}z_{2,j}/z_{0,j},$ (7) $\displaystyle
C_{C}(r)$ $\displaystyle=C(r_{1},r_{2}+r)-m^{2},$ (8)

where we also defined the connected correlator $C_{C}$. Importantly, the
multiplicative nature of the method results in relatively compact expressions,
with number of terms scaling as logarithm of the system size, with typical
thermodynamic observables exhibiting rapid convergence in system size, as
expected.

(a) (b)

(c)

Figure 2: Steps 2 and 3 of TRG. a) tensor networks at the beginning
(horizontal lattice with solid interconnected lines) and at the end (lattice
rotated by 45∘ with dotted interconnected lines) of 1 step of renormalization,
panel (c) for the formula for renormalized tensors $\tilde{T}$; b) use of SVD
to compress initial tensors $T$ and introduce new DOF. For computing
correlation functions some of these may be spin operators instead of local
Gibbs weight $T$; c) tracing out old DOF to arrive at the new coarse grained
tensor network with renormalized tensors $\tilde{T}$ denoted as black filled
circles in panel (a).

There are three basic steps (see Figs. 1 and 2) in going from a microscopic
description of the problem (i.e. $H$) to expressions for $z_{n,j}$: (i)
expressing a given problem as a tensor network by tracing out microscopic
degrees of freedom (DOF) and introducing new auxiliary DOF; (ii) approximate
compression of individual tensors in the network which yields a set of DOF
akin to Hubbard-Stratonovich fields optimally chosen to encode local
correlations; (iii) tracing out of previous set of auxiliary degrees of
freedom (can be done exactly) to arrive at the next tensor network with local
geometry equivalent to the previous one, but with fewer, more optimal “coarse-
grained” DOF. Steps (ii) and (iii) are iterated until the entire network is
coarse grained to just very few tensors that can then be traced over directly.

None of the steps are unique, e.g. in step (i) the simple Ising model we may
be rewritten using plaquette variables (following the seminal paper by Levin
and Nave LevinNaveTRG  or in terms of bond variables, using a singular value
decomposition (SVD) of the local transfer matrices $M_{x}$ (in the
$x$–direction)and $M_{y}$ (in the $y$–direction) SRG_PRL ; SRG_PRB . The
former turns out inferior as it only works at zero field and can lead to slow
artifacts in performance on bipartite lattices BasuUnpub  and we use the
latter method to construct our 2D tensor network described below. The transfer
matrices $M_{x}$ and $M_{y}$ are constructed from the Hamiltonian (here we
allow for the possibility of anisotropic couplings $J_{x}\neq J_{y}$) given by
(Fig.1(a))

$\displaystyle M_{x}=\begin{pmatrix}e^{\beta(J_{x}+h/2)}&e^{-\beta J_{x}}\\\
e^{-\beta J_{x}}&e^{\beta(J_{x}-h/2)}\end{pmatrix}$ (9) $\displaystyle
M_{y}=\begin{pmatrix}e^{\beta(J_{y}+h/2)}&e^{-\beta J_{y}}\\\ e^{-\beta
J_{y}}&e^{\beta(J_{y}-\beta h/2)}\end{pmatrix}$ (10)

The matrices $U_{0x},V^{\dagger}_{0x}(U_{0y},V^{\dagger}_{0y})$ (Fig. 1(b))
are constructed using SVD of the transfer matrix $M_{x}(M_{y})$, which are
then subsequently contracted as shown in Fig.1(c) to derive the 2D tensor
network (see Fig. 1(d)).

The “plain-vanilla” original TRG algorithm LevinNaveTRG  uses a simple SVD
based “rewiring” truncation transformation for step (ii) and an exact trace in
step (iii), which is what we stick with in this work, for simplicity and is
shown in Fig. 2. One caveat we must note is that the tensors $T_{A}$ (blue)
and tensors $T_{B}$ (red) in Fig. 2(a) are identical tensors, but labeled
differently due to the distinct ways in which one can perform the SVD
decomposition on tensor $T$, as shown in Fig. 2(b) to obtain the rank–$3$
tensors $S^{1},S^{2},S^{3}$, and $S^{4}$. The ‘$S$’ tensors are then
contracted in the manner shown in Fig. 2(c) integrating out the old DOF to
obtain the renormalized tensor $\tilde{T}$. In purely statistical problems it
is conceptually advantageous to employ the so called Takagi decomposition that
produces positive definite renormalized tensors that are amenable to direct
block-spin interpretationBerkeretal ; bal2017

The overall size of the lattice does not appear explicitly in the method, but
rather in the number of iterations performed. For a fixed system size of
interest the final tensor therefore has elements scaling exponentially in
system size, which is not numerically tractable. Instead, a simple
normalization scheme is imposed to rescale the renormalized tensor $\tilde{T}$
after each iteration to line it up with the previous $T$, with rescaling
factors $z_{m,j}$ collected in this process, except for the last one, which
comes from tracing out the renormalized final tensor network (typically made
up of very few tensors).

TRG can achieve an exceptionally high degree of accuracy at small or moderate
computational cost, especially in two-dimensional problems. This is generally
understood as a consequence of individual coarsegraining steps being local in
space and well controlled through application of SVD _and_ rapid convergence
to relatively simple fixed point tensors, e.g. in the paramagnetic phase, and
so the individual errors in the method are small both at the beginning of the
flow and also late in renormalization. As the only explicitly approximate step
is the compression step (ii) using SVD we shall refer to the numerical errors
of the algorithm broadly as “finite bond dimension” errors. For some
quantities these errors are heavily dominated by initial steps in
renormalization, due to volume factors (the size of the lattice shrinks
exponentially fast under renormalization), e.g. the expression for the free
energy above can be rewritten as
$f\to\sum_{j^{\prime}=0}^{j}f_{j^{\prime}}/2^{j^{\prime}}$ with $f_{j}\lesssim
1$ to exhibit this predominance of short distance contributions. As we are
more interested in universal properties, such short distance dominated errors
are not especially worrisome or interesting. They do appear, most prominently,
as non-universal bond dimension dependent shifts of the critical temperature.
The pattern of short distance errors appears erratic leading to these errors,
while small, also being of random sign and poorly convergent.

The more delicate, and in a certain sense fundamental, issue is whether any
given RG scheme reaches the correct fixed point tensor. It is empirically
known that the original TRG scheme does so sufficiently far away from the
critical point, correctly capturing paramagnetic (PM) and ferromagnetic (FM)
phases of the Ising model corresponding to infinite and zero temperature fixed
points, respectively. There appears a range of temperatures (that presumably
shrinks with bond dimension) near the critical temperature $T_{c}$ where the
fixed point tensors of the transformation are non-critical so-called corner
double line (CDL) tensors, distinct from the PM/FM onesTEFR ; EvenblyTNR
(this is also apparently true of some of the other schemes, e.g. HOTRGlyu2021
). It is unclear, based on published previous work, how these spurious fixed
points of tensor renormalization depend on the initial choice of a mapping to
a tensor network, lattice symmetry and especially presence (and magnitude) of
small symmetry breaking fields.

There have been several improvements to TRG to address (remove) this artifact
through better local coarse graining steps and also work to implement global
optimization of the coarse graining process. The detailed nature of these
unphysical fixed points, their significance to long distance correlations
_and_ their robustness against weak symmetry breaking fields have only been
partially elucidated. Our analysis of _universal_ content of critical points
in this work helps fill that gap, especially in regard to the last aspect,
i.e. introduction of very small symmetry breaking field. We expect (and find)
that long distance properties have finite bond dimension errors that are of a
fixed sign as they originate (presumably) from how well a finite dimensional
TRG fixed point tensor approximates the true critical fixed point tensor that
is infinite dimensional. Unlike variational tensor based approaches using
matrix product states (MPS), where finite bond dimension implies finite
correlation length rounding of the phase transition, it is not clear (to us)
whether/how similar length generation takes place in the present context.
Empirically, based on the data we present below (see Figs. 4 and 5) for the
two-point correlation function, such length scale, if present, is either
astronomically large or does not appear in simple two-point correlations.

Lastly, we survey prior efforts and methodology to study critical exponents
using tensor renormalization. By far the most common approach is to assume
that tensor renormalization correctly reproduces the expected conformal field
theory operator structure and use the goarsegrained tensor of a finite lattice
to identify conformal dimensions of operators from finite size correction to
the partition function by following Cardy’s ansatz, i.e. without explicitly
computing correlation functionsTEFR ; EvenblyTNR ; LoopTNR . Alternately, in
one of the earliest attempts, Hinczewski and BerkerBerkeretal  estimated the
thermal eigenvalue of the RG flow by directly computing the tangent map of TRG
to find its Lyapunov spectrum. This beautiful application of basic
renormalization group seemed to work well for small bond dimension but then
authors encountered strange flow instabilities and had to resort to the more
brute force method of differentiating the free energy and examining the
resultant specific heat and the associated exponent $\alpha$ (see also a much
more recent effort along these lines by Lyu and collaboratorslyu2021 ). The
specific heat singularity was also used to locate $T_{c}$ in other
worksBerkeretal . Off-critical correlation functions were computed away from
$T_{c}$ by Nakamoto and TakedaTakedaTRGcorrelation  and used to estimate
exponents $\beta$ and $\nu$. Importantly, the calculation of the two-point
function simplifies to a concatenation of two one point function flows for
simple exponential families of spatial distances, e.g.
$r=(0,2^{j}),(2^{j},0),(2^{j},2^{j})$ on the square lattice with steps up to
$j$ tracking the renormalization of bulk tensors and that of the spin, then
fusing them into a new composite operator and renormalizing that for larger
$j$. Such exponentially separated set of points is also best for studying
powerlaws using log-log plots and so we exploit this double advantage in our
analysis below. In summary, while the CFT method is by far the most powerful
and elegant of these, we opted for explicit renormalization of one- and two-
point operators to examine observables directly and study criticality with
minimal set of assumptions, anchored in direct observation of scale invariant
behavior.

###  II.2 Locating the Critical Point and Estimating Exponents

The critical point is a singularity in system’s macroscopic properties. It is
a low dimensional manifold in the infinite dimensional space of all
parameters, corresponding to the so called relevant operators. For the Ising
transition there are only two – temperature and field, so we write
$m(T_{c}-T,h)\sim(T_{c}-T)^{\beta},\sim\rm{sign}(h)|h|^{\sigma}$ for $h\to
0^{+}$ and $T=T_{c}$, respectively. The susceptibility $\chi\equiv\lim_{h\to
0}\dfrac{\partial m}{\partial h}$, follows $\chi_{>}(T>T_{c})\to
A_{>}(T-T_{c})^{-\gamma}$ (and similarly for $T<T_{C}$) and
$\chi_{\text{max}}(T=T_{c})\sim L^{\gamma/\nu}$, where $L$ is the linear
system size. The two point correlator at the critical point $T=T_{c}$ follows
$C(r_{1}-r_{2})\sim|r_{1}-r_{2}|^{-\eta}$ (in two-dimensions, also in the
generic case; see the Yang-Lee case below for an exceptional case). Previous
studies have examined exponents $\beta$, $\nu$ and $\alpha$ by directly
computing the observables (as we do here)Berkeretal ; LevinNaveTRG , this
paper will tackle $\eta,\gamma,\sigma$, as well as the amplitude ratio
$A=A_{>}/A_{<}$ which is also universaldelfino1998universal . Since the
numerically found exponents will turn out close to their known values we will
not be using hyperscaling to further test the consistency of the results,
which otherwise would have been a useful check.

We shall use the divergence of the susceptibilty to locate the crtical
temperature $T_{c}(D)$, which shifts compared to its true location due to
finite bond-dimension $D$ approximation employed. Locating the critical point
is the required first step in any systematic study of the critical point. Past
work showing breakdown of TRG seems not to have focused on this point
carefullyEvenblyTNR  or used other critical “markers”, e.g. the specific heat
singularityBerkeretal  (which may be adequate for the case of the 2D Ising
model but not generally) or the divergence of correlation length in the
paramagnetic phase TakedaTRGcorrelation  . The divergence in the
susceptibility is both generic for symmetry breaking transitions and tends to
be the most prominent, and hence is the natural choice for locating the
critical point in a finite sized sample. We also note that our $T_{c}$
estimates are similar or somewhat better to those obtained by other
methodsBerkeretal ; TakedaTRGcorrelation  despite comparatively more modest
bond dimension values, which further supports our approach to locating the
critical point. Recalling that iteration index in TRG is the logarithm of the
sample size (for any bond dimension) we expect the susceptibility to diverge
exponentially with TRG iteration at $T=T_{c}$, hence, we can estimate the
relevant combination of critical exponents, $\gamma/\nu$, from the rate of
this exponential growth.We also use onset of spontaneous magnetization
(properly computed in the presence of small field – more below) as an
alternative estimate of transition point and observe that the two methods
correlate well, but not perfectly, of course. Lastly, we expect and observe
scale invariant decay in the 2-point correlation function at criticality, up
to distances that depend sensitively on residual relevant operators, e.g.
imperfect tuning of $T-T_{c}$ or small applied field. Thus, to summarize, our
convidence in locating $T_{c}(D)$ rests on consensus among three computed
properties – agreement between magnetization onset (plots not shown), peak in
$\chi(T)$ with powerlaw divergence of critical susceptibility
$\chi_{\text{max}}(T=T_{c})\sim L^{\gamma/\nu}$, _and_ powerlaw decay of the
critical two-point function with separation $C(r)\propto r^{-\eta}$, listed
here in the order of increasing sensitivity, e.g. deviations from powerlaws
appear on shorter scales in the correlation function than
susceptibility111Note, Monte-Carlo simulations employ magnetization (Binder)
cumulants to locate the critical point – these are not accessible in the
simplest TRG formulation employed here. However, typical lattice sizes
amenable to Monte-Carlo sampling are comparatively small and which makes our
direct method less reliable, presumably.  . Importantly, all this analysis is
conducted in the presence of a small symmetry breaking field, typically
$h\lesssim 10^{-12}$. We have heuristic evidence that the presence of this
field improves convergence of TRG without disrupting powerlaws out to
distances $\gtrsim 2^{12}$.

Having established the location of the critical point (for a particular bond
dimension), we can readily obtain estimates of the critical exponent $\gamma$,
the amplitude ratio $A\equiv A_{>}/A_{<}$, and the combination of critical
exponents $\gamma/\nu$, from the finite size cutoff of susceptibility
precisely at $T=T_{c}$. We also compute the two-point correlation function at
large separation, which is a particularly direct measure of criticality, from
which we extract the $\eta$ exponent. Finally, we extract $\delta$ from the
nonlinear magnetization at $T=T_{c}$, to round out the roster of standard
exponents that have not been computed previously using TRG (to the best of our
knowledge).

Table 1: Summary of universal and non-universal data on the isotropic Ising
critical point, see text for definition of absolute error $\Delta$ Table 2:
Summary of universal and non-universal data on the anisotropy 10 Ising
critical point. Figure 3: Convergence trend for exponents $\delta$, $\eta$
(closed symbols) and $\gamma$, $\gamma/\nu$ (open symbols) for the anisotropy
10 model. (Inset) Visual test of the empirical $D^{-5}$ behavior.

###  II.3 Errors

There are two sources of technical error in this analysis – (i) finite
resolution in the values of magnetic field $h$ and temperature $T$ (relevant
operators), which may be refined as needed and (ii) finite bond dimension,
which are not thoroughly understood and commonly dealt with by switching to
more elaborate (and accurate) coarsegraining schemes. As all of these are
systematic errors, one cannot rely on added runs to improve results, e.g. as
with Monte-Carlo simulations. The latter, of course suffer from their own
problems inherent in sampling large lattices, including but not restricted to
critical slowing down, all of which are entirely absent here. While finite
bond dimension can be related to finite length rounding of critical
correlations in variational MPS based approaches to classical
criticalityburgelman2022contrasting , the exisence of such finite length
rounding is far from obvious in TRG and other multiplicative renormalization
approaches. Thus it is somewhat unclear (to us) how to relate technical errors
that are clearly present in the computation to the physical errors (if any) in
conclusions drawn. We shall return to this point in Sec. IV after presenting
the results in the next Section (see also discussion above, in Sec. II.1).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: (color online) Left (right) column are data from the TRG flow with
bond dimension $D=20$ on isotropic (anisotropy 10) Ising model in
infinitesimal magnetic field $h=10^{-12}$. (a,b) Inverse temperature scan of
susceptibility $\chi(\beta)$ showing the critical divergence at $\beta_{C}$
(red dashed line) slightly displaced from the true value
$\beta_{C}^{\text{Onsager}}=1/T_{C}$ (blue dashed line). Inset: The plot of
$\log\chi$ versus the iteration index of the TRG algorithm at $\beta_{C}$
clearly shows the expected finite-size scaling of the singularity, from which
we extract $\gamma/\nu$; (c,d) Critical and slightly off-critical correlations
follow powerlaw decay out to scales determined by residual relevant
perturbations, here corresponding to deviations around separations $=2^{13}$;
(e,f) Non-linear critical magnetization displays a clear powerlaw with an
associated critical exponent $\delta$; (g,h) Off-critical susceptibility
clearly showing scale invariant behavior assoicated with exponent $\gamma$ and
also amplitude ratio $A\approx 45$ for the isotropic and $A\approx 48$ for the
anisotropic, roughly consistent with the the knownmccoyetal76amplitudes ;
delfino1998universal  exact value of $\approx 38$.

##  III results

This Section contains all the main results of the paper. For each bond
dimension $D$ we perform a systematic scan of temperature and field (for the
Ising-Onsager critical point) and complex field (for the Yang-Lee critical
point) and collect the data into log-log plots (see Figs. 4 and 5) from which
critical temperatures and critical exponents are extracted (and amplitude
ratios for susceptibility) (see Tables 1 and 2). For the Onsager-Ising
critical point we focus on contrasting the results from isotropic and
anisotropy 10 square lattice simulations (left and right columns in Fig. 4).
The Yang-Lee critical point appears considerably more challenging for TRG but
still exhibits a clear convergence trend in bond dimension. We only consider
quantities for which exact values are known, e.g. critical temperatures and
exponents, and so for any such quantity X, we define a dimensionless error
$\Delta X\equiv X^{\text{approx}}/X^{\text{exact}}-1$. For critical (inverse)
temperature $\beta_{c}$ we use a slightly more elaborate notation, with a
superscript to denote the method used to extract the estimate of critical
inverse temperature $\beta_{c}$ – $\chi$ and $m$ for susceptibility and
magnetization, respectively.

###  III.1 The Onsager-Ising critical point

Isotropic 2D Ising model was used as the default model for benchmarking most
advances in tensor renormalizationLevinNaveTRG ; SRG_PRL ; SRG_PRB ; HOTRG ;
EvenblyTNR ; LoopTNR . It has been appreciated, however, that this model is
somehow non-generic, e.g. it has artificially small corrections to scaling and
therefore may not be ideal for studying convergence properties of an
approximate method. However, much of past work focused on non-universal
observables, e.g. the free energy and critical temperature. By contrast, we
focus on universal quantitities and find (see Table 1) at best very
preliminary signs of convergence in the limited range of bond dimension
explored here.

The 2D Ising model in zero field is exactly solvable for arbitrary anisotropy,
it exhibits identical universal scaling at the Onsager critical point
regadless of anisotropy. We do expect for a strongly anisotropic Ising model
to exhibit a crossover from one to two dimensional behavior scaling which may
make it somewhat challenging to capture for an approximate computational
scheme or at the very least have a richer scaling function. Additionally, as
we recently found Basuetal2021 , while the isotropic Ising critical point is a
simple critical endpoint of a first order transition in the complex
temperature plane, the anisotropic Ising critical point is a multicritical
corner-point of a critical ordered phase, with modulated quasi-periodic order
Basuetal2021 . This does not change any universal properties along the real
temperature axis, as far as we know). Thus, the immediate proximity of a
critical phase (albeit at complex temperature) is another reason to be
interested in repeating the analysis of critical scaling to contrast the
performance of TRG between isotropic and anisotropic Ising models. The results
in Table 2 for anisotropy 10 were partly expected, i.e. larger deviations from
exact values at smaller bond dimensions. Remarkably (and surprisingly),
however, we also observe a clear convergence trend in universal exponents in
bond dimension, see Fig. 3, in some cases to an order of magnitude better than
for the isotropic model, e.g. for $\gamma$ and $\gamma/\nu$. Empirically we
observe $~{}1/D$ and $1/D^{5}$ convergence for pairs $\eta,\delta$ and
$\gamma,\gamma/\nu$ of exponents, respectively. We do not have any clear
understanding of this result nor theory for what one might expect for the
asymptotic convergence in bond dimension, as $1/D\to 0$ (see also Sec. IV).
Lastly, we note that both sets of computations reproduce one of the most
unusual aspects of the Ising-Onsager critical point – the unusually large
amplitude ratio for the susceptibility (see bottom panels of Fig. 4 and
corresponding caption).

(a)

(b)

(c)

(d)

Figure 5: (color online) The top panel (row) shows the plots of $\log\chi$ as
a function of $\log H$ exhibiting powerlaw behavior, where $H=h-ih_{C}$ for
(a) $D=16$; and (b) $D=20$ from where we extract the critical exponent
$\sigma$. In the bottom panel, we show the plots of $\log G(r,H)$ as a
function of $\log r$ for (c) $D=16$; and (d) $D=20$. The decay for the
connected correlator $G(r,H)$ is slow with critical exponent $\eta$. The
critical exponents thus obtained are in good agreement with the analytical
values $\sigma^{\text{Ana}}=-1/6$, and $\eta^{\text{Ana}}=-4/5$. See Sec.
III.2 for details.

###  III.2 The Yang-Lee critical point

First order transitions may be thought of as branch cuts in the (possibly
enlarged) space of couplings. The seminal work by Yang and Lee details the
trajectory of one such branch cut, corresponding to a ferromagnet in a field,
$h$. Using the fundamental theorem of algebra they were able to establish
explicitly that the familiar transition that takes place at $h=0$ in the
ferromagnetic phase $(T<T_{c})$ recedes to finite purely complex fields for
$(T>T_{c}$), with corresponding branch points (critical endpoints) located at
$h_{c}=\pm i|h_{c}|$.222The formalism of analytic continuation of couplings
proved useful, e.g. in establishing the existence of essential (Griffiths)
singularites in disorded magnets. The Yang-Lee endpoints are novel, non-
Onsager critical points. They are sometimes referred toFisherPRL78  as “proto-
critical” points for having a single relevant perturbation, the complex
magnetic field. They were identified with a complex scalar $\varphi^{3}$ field
theory by FisherFisherPRL78  and, in two dimensions, the minimal model
$\mathcal{M}(5,2)$ by CardyCardyPRL85 . Hence, the exact critical exponents
$\sigma=-1/6$ and $\eta=-4/5$ are known. Importantly for what follows, the
standard two-point scaling ansatz in the vicinity of the Yang-Lee critical
point reads

$C_{C}(r,\xi)=D(r/\xi)r^{d-2+\eta}$ (11)

with correlation length $\xi\sim|h-h_{c}|^{-\nu}$ and a _singular_ scaling
function $D(x)=A\exp(-Bx)/x^{2}$, where $A,B$ are constants. Importantly, the
additional singularity in $D(x)$ implies a divergent _amplitude_ of the
correlator $\sim\xi^{2}$.

Monte-Carlo techniques are powerless to access the Yang-Lee and other
similarly exotic critical points, so the analysis historically relied on a
combination of methods, including high temperature series griffiths69 ;
Kurtze79 , analysis of partition function zeroes in finite lattices,
culminating in matching to known conformal field theoriesCardyPRL85 ;
CardyMussardo . Tensor network formulation was applied to this
problemGarciaWeiYL  but did not consider two-point correlation function
$C_{C}(r)$ nor susceptibility. Our results for these two quantities are
presented in Fig. 5. In spite of the noisy data, there is clear evidence of
powerlaw dependences from which we extract $\sigma$ and $\eta$. The errors
made for bond dimension $D=16$ and $20$ are, respectively,
$\Delta_{\sigma}=-0.02,\Delta_{\eta}=0.18$ and
$\Delta_{\sigma}=0.01,\Delta_{\eta}=0.06$. These preliminary results appear
encouraging to expect a convergence trend as with the anisotropic Ising model.

##  IV Summary and outlook

To summarize, we computed several critical exponents of Ising and Yang-Lee
critical points using TRG method to obtain and analyse correlation functions
directly. We were especially interested to study convergence trends with bond
dimensions and documented clear convergence in anisotropic Ising model and for
the Yang-Lee critical point but not for the isotropic Ising model. At the
present we do not have a clear understanding of this observation and plan to
explore it further. One candidate culprit is the onset of spurious fixed
points of TRG near the critical point as has been previously documented TEFR ;
EvenblyTNR  for the isotropic model in zero field. It would be interesting to
explore the convergence trends of other more sophisticated methods of tensor
renormalization and also compare to the values of exponents obtained by other
methods, e.g. from scaling dimensions using CFT ansatz for the finite size
free energy. Tensor renormalization techniques can also be fruitfully applied
to statistical mechanics problems with non-positive weights as demonstrated
here on the Yang-Lee problem. There are several such interesting problems in
the literature, e.g. models with complex fugacity (similar to Yang-Lee), non-
linear sigma models with topological terms, problems with gauge fields (and
long-range interactions), quantum circuits with measurementBasuetal2021  and
we expect to see numerical progress on these in the near future. Additionally,
we expect tensor methods to make impact on problems with competing phases,
possibly separated by weak first order transitions which often present severe
difficulty to Monte-Carlo techniques. While generic three dimensional problems
are still likely too difficult for an effective application of TRG, strongly
anisotropic lattices or lattices with low local coordination may turn out more
amenable.

###### Acknowledgements.

We would like to thank Miles Stoudenmire, Tzu-Chieh Wei, Nikko Pomata, Frank
Pollman, Chris Laumann, Aleix Bou Comas for many illuminating conversations.
The Flatiron Institute is a division of the Simons Foundation.

## References

  * (1) Michael Levin and Cody P. Nave.  Tensor renormalization group approach to two-dimensional classical lattice models.  Phys. Rev. Lett., 99:120601, Sep 2007.
  * (2) Nicholas Metropolis and S. Ulam.  The monte carlo method.  Journal of the American Statistical Association, 44(247):335–341, 1949.  PMID: 18139350.
  * (3) Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.
  * (4) U. Schollwöck.  The density-matrix renormalization group.  Rev. Mod. Phys., 77:259–315, Apr 2005.
  * (5) Steven R. White.  Density matrix formulation for quantum renormalization groups.  Phys. Rev. Lett., 69:2863–2866, Nov 1992.
  * (6) F. Verstraete and J. I. Cirac.  Renormalization algorithms for quantum-many body systems in two and higher dimensions.  7 2004.
  * (7) Alexander A. Migdal.  Recursion Equations in Gauge Theories.  Sov. Phys. JETP, 42:413, 1975.
  * (8) Leo P Kadanoff.  Notes on migdal’s recursion formulas.  Annals of Physics, 100(1):359–394, 1976.
  * (9) Sankhya Basu, Daniel P. Arovas, Sarang Gopalakrishnan, Chris A. Hooley, and Vadim Oganesyan.  Fisher zeros and persistent temporal oscillations in nonunitary quantum circuits.  Phys. Rev. Research, 4:013018, Jan 2022.
  * (10) Artur García-Saez and Tzu-Chieh Wei.  Density of yang-lee zeros in the thermodynamic limit from tensor network methods.  Phys. Rev. B, 92:125132, Sep 2015.
  * (11) C. N. Yang and T. D. Lee.  Statistical theory of equations of state and phase transitions. i. theory of condensation.  Phys. Rev., 87:404–409, Aug 1952.
  * (12) T. D. Lee and C. N. Yang.  Statistical theory of equations of state and phase transitions. ii. lattice gas and ising model.  Phys. Rev., 87:410–419, Aug 1952.
  * (13) Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y. Weng, and T. Xiang.  Second renormalization of tensor-network states.  Phys. Rev. Lett., 103:160601, Oct 2009.
  * (14) H. H. Zhao, Z. Y. Xie, Q. N. Chen, Z. C. Wei, J. W. Cai, and T. Xiang.  Renormalization of tensor-network states.  Phys. Rev. B, 81:174411, May 2010.
  * (15) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang.  Coarse-graining renormalization by higher-order singular value decomposition.  Phys. Rev. B, 86:045139, Jul 2012.
  * (16) G. Evenbly and G. Vidal.  Tensor network renormalization.  Phys. Rev. Lett., 115:180405, Oct 2015.
  * (17) Shuo Yang, Zheng-Cheng Gu, and Xiao-Gang Wen.  Loop optimization for tensor network renormalization.  Phys. Rev. Lett., 118:110504, Mar 2017.
  * (18) Zheng-Cheng Gu and Xiao-Gang Wen.  Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order.  Phys. Rev. B, 80:155131, Oct 2009.
  * (19) Sankhya Basu.  in preparation.
  * (20) Michael Hinczewski and A. Nihat Berker.  High-precision thermodynamic and critical properties from tensor renormalization-group flows.  Phys. Rev. E, 77:011104, Jan 2008.
  * (21) M. Bal, M. Mariën, J. Haegeman, and F. Verstraete.  Renormalization group flows of hamiltonians using tensor networks.  Phys. Rev. Lett., 118:250602, Jun 2017.
  * (22) Xinliang Lyu, RuQing G. Xu, and Naoki Kawashima.  Scaling dimensions from linearized tensor renormalization group transformations.  Phys. Rev. Research, 3:023048, Apr 2021.
  * (23) N. Nakamoto and S. Takeda.  Computation of correlation functions by tensor renormalization group method.  Sci. Rep. of Kanazawa Univ., 60:11, Jan 2016.
  * (24) G Delfino.  Universal amplitude ratios in the two-dimensional ising model1work supported by the european union under contract fmrx-ct96-0012.1.  Physics Letters B, 419(1):291–295, 1998.
  * (25) Note, Monte-Carlo simulations employ magnetization (Binder) cumulants to locate the critical point – these are not accessible in the simplest TRG formulation employed here. However, typical lattice sizes amenable to Monte-Carlo sampling are comparatively small and which makes our direct method less reliable, presumably.
  * (26) Lander Burgelman, Lukas Devos, Bram Vanhecke, Frank Verstraete, and Laurens Vanderstraeten.  Contrasting pseudo-criticality in the classical two-dimensional heisenberg and rp2 models: zero-temperature phase transition versus finite-temperature crossover.  arXiv preprint arXiv:2202.07597, 2022.
  * (27) Tai Tsun Wu, Barry M. McCoy, Craig A. Tracy, and Eytan Barouch.  Spin-spin correlation functions for the two-dimensional ising model: Exact theory in the scaling region.  Phys. Rev. B, 13:316–374, Jan 1976.
  * (28) The formalism of analytic continuation of couplings proved useful, e.g. in establishing the existence of essential (Griffiths) singularites in disorded magnets.
  * (29) Michael E. Fisher.  Yang-lee edge singularity and ${\phi}^{3}$ field theory.  Phys. Rev. Lett., 40:1610–1613, Jun 1978.
  * (30) John L. Cardy.  Conformal invariance and the yang-lee edge singularity in two dimensions.  Phys. Rev. Lett., 54:1354–1356, Apr 1985.
  * (31) Robert B. Griffiths.  Nonanalytic behavior above the critical point in a random ising ferromagnet.  Phys. Rev. Lett., 23:17–19, Jul 1969.
  * (32) Douglas A. Kurtze and Michael E. Fisher.  Yang-lee edge singularities at high temperatures.  Phys. Rev. B, 20:2785–2796, Oct 1979.
  * (33) John L. Cardy and G. Mussardo.  S-matrix of the yang-lee edge singularity in two dimensions.  Physics Letters B, 225(3):275–278, 1989.