# Overview of the CAIL 2023 Argument Mining Track

Jingcong Liang
Fudan University
<EMAIL_ADDRESS>
&Junlong Wang, Xinyu Zhai
Dalian University of Technology
<EMAIL_ADDRESS>
&Yungui Zhuang
Tgnet
<EMAIL_ADDRESS>
&Yiyang Zheng
Shanghai University
<EMAIL_ADDRESS>
&Xin Xu

<EMAIL_ADDRESS>
&Xiandong Ran
Xi’an Jiaotong University
<EMAIL_ADDRESS>
&Xiaozheng Dong, Honghui Rong
Lenovo KnowDee Intelligence
<EMAIL_ADDRESS>
&Yanlun Liu
Hunan Uviersity
<EMAIL_ADDRESS>
&Hao Chen
Xiamen Huatian International Vocation Institute
<EMAIL_ADDRESS>
&Yuhan Wei
Northwest Polytechnical University
<EMAIL_ADDRESS>
&Donghai Li
Beijing Thunisoft Yuandian Information Service Co., Ltd.
<EMAIL_ADDRESS>
&Jiajie Peng
Northwest Polytechnical University
<EMAIL_ADDRESS>
&Xuanjing Huang
Fudan University
<EMAIL_ADDRESS>
&Chongde Shi
Beijing Thunisoft Yuandian Information Service Co., Ltd.
<EMAIL_ADDRESS>
&Yansong Feng
Peking University
<EMAIL_ADDRESS>
&Yun Song
Northwest University of Political Science and Law
<EMAIL_ADDRESS>
&Zhongyu Wei11footnotemark: 1
Fudan University
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

We give a detailed overview of the CAIL 2023 Argument Mining Track, one of the
Chinese AI and Law Challenge (CAIL) 2023 tracks. The main goal of the track is
to identify and extract interacting argument pairs in trial dialogs. It mainly
uses summarized judgment documents but can also refer to trial recordings. The
track consists of two stages, and we introduce the tasks designed for each
stage; we also extend the data from previous events into a new dataset —
CAIL2023-ArgMine — with annotated new cases from various causes of action. We
outline several submissions that achieve the best results, including their
methods for different stages. While all submissions rely on language models,
they have incorporated strategies that may benefit future work in this field.

_K_ eywords Argument Mining $\cdot$ Judgment Document Comprehension

##  1 Introduction

One of the core tasks for judges during a trial is to understand the arguments
and testimonies from the two sides, as they need to base their verdicts on
them appropriately. Typically, the panel summarizes the trial dialog into a
judgment document, on which this judicial argument comprehension process
heavily depends [1]; however, the panel has to perform the task manually by
reading the entire document. The call for automating this process and other
judicial tasks has yielded various studies and models involving both
statistical methods [2, 3] and natural language processing (NLP) techniques
[4, 5, 6], including the latest large language models [7, 8, 9, 10].

However, judicial argument comprehension differs from other tasks in the
judicial field, as it requires extracting arguments from the trial dialog (and
its summarized form) and understanding their interaction, especially between
sides. The arguing procedure greatly resembles general argumentation
procedures like debates, where such automatic comprehension is also possible
[11, 12]. More specifically, the judicial argument comprehension task largely
overlaps the domain of argument mining [13]. The development of NLP has also
resulted in many studies that aimed to automate argument mining, covering
argument structure prediction [14, 15, 16] and interaction identification [17,
18, 19, 20, 21, 22, 23, 24].

Based on the tasks designed by [24], [25] hosted the SMP-CAIL2020-Argmine
Challenge111http://cail.cipsc.org.cn/task_summit.html?raceID=3&cail_tag=2020,
a track of the Chinese AI and Law Challenge (CAIL) 2020. The track addressed
the importance of adapting argument-mining techniques to judicial scenarios
and assisting judicial argument comprehension. Participants must develop
models to identify and extract _interacting arguments_ in trial dialogs. Here,
two arguments interact if they come from different sides (plaintiff and
defense) and make claims about the same subject, as shown in Table 1.
Interacting arguments can be consistent or partially consistent, yet in most
cases, one disputes another.

Table 1: Examples of interacting argument pairs between the plaintiff and defense in different trials. Note that in the last example, the arguments are partially consistent. No. | Side |  Argument
---|---|---
1 | Plaintiff |  同年3月14日凌晨，被告人何**回到家中向高*谎称李某某走失，并阻止高*报警但未果。
Defense |  当晚11点至12点其回家后，便与女儿高*报警。
2 | Plaintiff |  被告在家很少操持家务，也很少做饭。
Defense |  被告一直在家伺候公婆和孩子的生活起居，还要兼顾浴池的生意。
3 | Plaintiff |  2015年12月26日，被告王**、周**因果园种植需要，从原告处借款3万元，约定月利率1.5％，借款期限6个月。
Defense |  被告王**借原告款属实，但借款时原告当场扣除了3000元。

As the direct successor to the SMP-CAIL2020-Argmine Challenge, the latest CAIL
2023 Argument Mining
Track222http://cail.cipsc.org.cn/task_summit.html?raceID=5&cail_tag=2023
extended the datasets and introduced new tasks closer to real scenarios. We
divided the track into two stages, examining interacting argument
identification and extraction abilities. More specifically, the first stage
requires the model to choose the argument from five candidates by the defense
that interacts with the given plaintiff’s argument. The second stage requires
the model to extract interacting argument pairs from judgment documents
directly.

In this paper, we report and analyze the best submissions to the track,
including their methods and evaluation measures in both stages. All these
submissions used language models as their backbone, yet they have used
different approaches to improve their performance further. We also discuss the
findings and reflections from these submissions for future study in this area.

##  2 Related Work

Automatic processing of judicial tasks, such as judgment document analysis,
has been researched since the 1950s. Early studies mainly focused on
quantitative and statistical methods [26, 2, 3]; more recent works also
introduced machine learning approaches [27]. Meanwhile, the rapid development
of natural language processing (NLP) has greatly aided the current progress in
this field. For example, [4] and [6] trained $N$-gram based classifiers for
judicial verdict prediction; [5] and [28] utilized various features to predict
criminal case charges. The latest studies also utilized the power of large
language models (LLM) to tackle general judicial tasks and queries [7, 8]. In
2018, [29] held the Chinese AI Law Challenge (CAIL), initially focusing on
judgment prediction, along with a large-scale legal dataset [30]. Since then,
more judicial tasks and datasets were included, such as similar case matching
[31] and judicial document classification [32].

The specific task of argumentative analysis of judgment documents was
introduced in 2020 for its unique characteristic primarily related to argument
mining [25]. Argument mining research focuses on defining and evaluating new
algorithms for identifying and structurally analyzing human argumentation.
Typical tasks in this area include argument segmentation and extraction [33,
34, 35, 36], argument classification [15, 37, 38] and argument relation mining
[39, 40, 41]. However, the most related task to the CAIL 2023 Argument Mining
Track is interactive argument detection, a combination of argument extraction
and relation mining [42, 23]. Specifically, [24] proposed identifying
interacting argument pairs in online debate forums, which was the direct
inspiration for the tasks in our track. Since then, various methods have been
proposed to tackle this task [43, 44].

##  3 Task Details

###  3.1 Task Description

The CAIL 2023 Argument Mining Track focused on recognizing interacting
argument pairs in trial dialogs. As mentioned in Section 1, we define
interacting argument pairs as two arguments, by the plaintiff and defense,
respectively, that argue over the same subject. A typical judgment document
lists the plaintiff’s claims and evidence first, followed by the defense’s
rebuttals. Therefore, in most cases, each of the defense’s arguments is
recorded as acknowledging or denying one or more arguments by the plaintiff,
forming one or more interacting argument pairs. As a result, the track aimed
to find these interacting argument pairs in judgment documents.

We split the track into two stages. The first stage requires identifying
pairing arguments from the given candidates, and the second requires
extracting argument pairs from judgment documents. Due to the low density of
interacting argument pairs, argument pair extraction is more complicated than
identification. In this way, we can distinguish models with various levels of
ability.

More specifically, in the first stage, given an argument by the plaintiff, we
provide five candidate arguments sampled from the defense side. One of these
arguments forms an interacting argument paired with the plaintiff’s argument.
We guarantee that only one argument can form such a pair, and the model should
identify and choose it. Since this task only involves pre-extracted arguments
from the judgment document, it is relatively easy but not very practical in
real scenarios.

In the second stage, we only provide the judgment documents; the model must
extract all interacting argument pairs between the plaintiff and defense. This
task is more difficult since only a few arguments can form such pairs, causing
the labels to be highly imbalanced. We also provided trial recordings as
multimodal data for a portion of the documents as a bonus; these recordings
were transcribed, with the speaker of each speech annotated. Samples with and
without trial recordings are separated, resulting in two subtasks evaluated
separately (using the same metrics).

###  3.2 Dataset

In the SMP-CAIL2020-Argmine Challenge, [25] annotated a dataset with $4,476$
interacting argument pairs from $1,069$ judgment documents. Their annotation
included the case type (civil, criminal, etc.), the cause of action, the
plaintiff and defense, and a list of interacting argument pairs. They also
filtered cases that lacked substantial argumentation or were too short or
long. In CAIL 2021, they extended the original dataset with more cases.
However, the original dataset and the extended version lacked variety
regarding the cause of action, especially civil ones.

Following their annotation schema, we refined their dataset and extended it
with more trial cases from various crime types. Details of the data collection
and annotation procedure are available in Appendix A, including annotation
statistics. We name the new dataset CAIL2023-Argmine, and we list its data
statistics in Table 2. CAIL2023-Argmine is twice as large as the original
dataset and covers more causes of action, especially civil cases. The argument
pair density remains low, again emphasizing the challenge of the track tasks.

Table 2: Overall statistics of the original dataset and CAIL2023-Argmine. 2021: Extended version of the original dataset for CAIL 2021; 2021*: The previous dataset without cases not belonging to any specific cause of action (the last three statistics also exclude them); 2023: CAIL2023-Argmine. Statistics | 2021 | 2021* | 2023
---|---|---|---
# Documents | $1,872$ | $1,449$ | $3,620$
# Argument Pairs | $10,077$ | $7,618$ | $20,009$
Pair Density (%) | $6.60$ | $6.20$ | $8.88$
# Causes of Actions | $5$ | $5$ | $13$
# Criminal Cases | $1,304$ | $1,304$ | $2,291$
# Civil Cases | $145$ | $145$ | $1,329$

Table 3 lists data statistics per cause of action in CAIL2023-Argmine. While
the cases in the original dataset mainly came from a few causes of action like
the crime of intentional injury (“故意伤害罪”) or causing traffic casualties
(“交通肇事罪”), new cases are distributed more uniformly and cover more causes of
action. Furthermore, the original dataset only contains one civil cause of
action — maritime disputes (“海事海商纠纷”) — with only $145$ documents. In
CAIL2023-Argmine, we annotated nearly $1,200$ new documents from $4$ different
civil causes of action, greatly extending this part of the data. Another
observation is that the argument pair density in the new cases is higher than
in the old ones, especially in criminal cases.

Table 3: Data statistics per cause of action. 2021: the extended version of the dataset in the SMP-CAIL2020-Argmine Challenge; 2023*: new cases included in CAIL2023-Argmine. Causes of action that end with “罪” are criminal causes, while others are civil ones. Note that “其他” is not a valid cause of action, and we excluded corresponding cases in CAIL2023-Argmine. Dataset | Cause of Action | # Docs | # Sentences | # Pairs | Density (%)
---|---|---|---|---|---
2021 | 故意伤害罪 | $956$ | $16,601$ | $4,183$ | $5.44$
交通肇事罪 | $315$ | $5,237$ | $2,394$ | $10.52$
故意杀人罪 | $29$ | $405$ | $152$ | $10.22$
虐待罪 | $4$ | $101$ | $30$ | $6.13$
海事海商纠纷 | $145$ | $2,639$ | $859$ | $7.59$
其他* | $423$ | $7,922$ | $2,459$ | $6.20$
2023* | 故意杀人罪 | $272$ | $4,005$ | $1,366$ | $9.25$
抢劫罪 | $253$ | $3,579$ | $1,698$ | $13.41$
强奸罪 | $239$ | $3,271$ | $1,501$ | $13.22$
诈骗罪 | $168$ | $2,601$ | $1,255$ | $12.54$
盗窃罪 | $55$ | $800$ | $469$ | $16.23$
借款合同纠纷 | $299$ | $4,157$ | $1,353$ | $9.27$
离婚纠纷 | $298$ | $4,612$ | $1,498$ | $8.23$
生命权、身体权、健康权纠纷 | $294$ | $4,515$ | $1,727$ | $9.73$
买卖合同纠纷 | $293$ | $4,248$ | $1,524$ | $9.88$

In addition to judgment documents, we prepared trial audio and transcriptions
for the multimodal bonus task in the second stage. Due to limited open access
to trial recordings, we could only obtain trial videos for some cases in
CAIL2023-Argmine. We extracted the audio from the videos and transcribed them
into dialog texts using speech recognition models. Since the model output did
not contain speaker information, we further annotated the speaker of each
speech in the dialog. Table 4 shows statistics of these multimodal data.

Table 4: Statistics of multimodal data in CAIL2023-Argmine. There are more recordings than documents because some cases have multiple trial recordings. Statistics | Value
---|---
# Documents | $134$
# Argument Pairs | $688$
Pair Density (%) | $8.49$
# Recordings | $146$
Average Duration | $50^{\prime}12^{\prime\prime}866$
# Speeches | $22,654$

###  3.3 Evaluation and Baseline

For the first stage, we evaluate the submissions by calculating the accuracy
of candidate prediction over the test set:

$S_{1}=\mathrm{Acc}=\frac{\sum_{i=1}^{n}I\\{y_{i}=\hat{y}_{i}\\}}{n}$ (1)

where $y_{i}$ and $\hat{y}_{i}$ are the true label and prediction of the
$i$-th sample in the test set respectively.

However, the accuracy was no longer suitable for the second stage due to
imbalanced labels (pairing or not), as seen in Table 2. Therefore, we used the
$F_{1}$ score as the evaluation metric for this stage:

$F_{1}=\frac{2\mathrm{TP}}{2\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}$ (2)

where $\mathrm{TP}$ is the number of correct predictions for _all_ cases,
$\mathrm{FP}$ is the number of incorrect predictions, and $\mathrm{FN}$ is the
number of missed pairs. Specifically, we calculated $F_{1}$ scores for
multimodal and text-only samples separately, obtaining $F_{1}^{\mathrm{m}}$
and $F_{1}^{\mathrm{t}}$ respectively. We took their weighed average as the
score of this stage:

$S_{2}=0.8F_{1}^{\mathrm{t}}+0.2F_{1}^{\mathrm{m}}$ (3)

Finally, we ranked teams that participated in both stages using the following
integrated score:

$S=0.3S_{1}+0.7S_{2}$ (4)

We provided a BERT-based baseline model [45] for both stages during the
track333https://github.com/china-ai-law-challenge/CAIL2023/tree/main/lblj.
This baseline model treats both tasks as binary classification. It pairs the
plaintiff’s argument to each of the defense’s candidate arguments and predicts
whether they form an interacting argument pair or not; Specifically, it does
not use judgment documents, trial transcripts, or multimodal data. We trained
the model with the training set of each stage to obtain the final baseline for
that stage.

In the first stage, the model achieved an accuracy of $0.60$; however, in the
second stage, it only achieved $F_{1}^{\mathrm{m}}=0.42$ and
$F_{1}^{\mathrm{t}}=0.30$. The evaluation result further proved that the task
in the second stage is more difficult than in the first stage.

##  4 Track Results

In the CAIL 2023 Argument Mining Track, $35$ teams from universities and
enterprises have submitted their models during the first stage. Among them,
$28$ teams successfully superseded our baseline model and were able to
participate in the second stage. Finally, $10$ teams have submitted valid
models for the second stage.

###  4.1 Submissions

This section provides descriptions of the top $7$ teams to whom we awarded
prizes. The rest of the teams were not required to submit technical reports
about their models, yet we will still show their results in Section 4.2.

DUT-large fine-tuned multiple pre-trained language models and employed a soft
voting approach to ensemble them. In the first stage, they concatenated the
plaintiff’s argument with each of the defense’s candidate arguments as the
model input. To enhance robustness, they introduced adversarial training
through Projected Gradient Descent (PGD) and R-Drop regularization. In the
second stage, they used a similar approach by treating the task as a binary
classification task. They involved sentence-BERT models [46] for this stage.
They employed several techniques, such as pseudo-labeling, focal loss, and the
Fast Gradient Method (FGM) for adversarial training.

X used BERT-based models for both stages but under different approaches. In
the first stage, they fine-tuned a MacBERT model [47] integrated with a multi-
choice setup to match arguments. They also leveraged the sentence-BERT model
for sentence pair binary classification in the second stage.

zyy utilized Qwen-14B [48] for both stages. The model was fine-tuned by Low-
Rank Adaptation (LoRA) [49], with NefTune applied to enhance performance. They
created prompts for each stage, respectively. In the first stage, the prompt
included the plaintiff’s argument and all five candidates and instructed the
model to choose the interacting one. In the second stage, they provided
arguments from both sides, prompting the model to find all matching pairs
between them.

xxxin fine-tuned a MacBERT model for the first stage. In addition to argument
pairs, they included designated tokens to identify the cause of action
(including its category and specific cause) and the plaintiff’s and defense’s
arguments. Each batch consisted of argument pairs with the same plaintiff’s
argument but with different defense arguments. In the second stage, they used
multiple methods, such as pseudo-labels, negative sampling, and adversarial
training, to deal with imbalanced labels.

AG adopted a list-wise learning approach, transforming both stages into a
sorting task. In the first stage, they used the provided candidates directly
as a list, and in the second stage, they constructed the list using Cartesian
products. They fine-tuned the Lawformer model [50], utilizing adversarial
learning with FGM.

福气by introduced multi-model voting for both stages. In the first stage, they
employed a combination of a multiple choice question-answering model, a text-
matching model capturing argument relation, and a re-ranking model to refine
the selection. They did not include the question-answering model and the re-
ranking model for the second stage; instead, they enhanced the text-matching
model with negative samples.

酸菜饺子&玉米饺子 fine-tuned a RoBERTa [51] model on large law corpora. They also
adopted FGM to enhance the robustness of the model.

###  4.2 Results

Table 5 demonstrates the CAIL 2023 Argument Mining Track results, including
the scores of both stages and the final score to rank teams. Although all
teams performed better than the baseline in the first stage, not all could
beat the baseline in the second stage. DUT-large significantly outperformed
all other teams in both stages, followed by X and zyy, which are very close.
Other teams mentioned in Section 4.1 were on par with X and zyy in one stage
but not in the other (“福气by” and “酸菜饺子&玉米饺子” in the first, xxxin and AG in the
second).

Table 5: CAIL 2023 Argument Mining Track final results, sorted by the final score $S$. Only teams with valid submissions for the second stage are listed. We round all scores to two digits. Team | $S_{1}$ | $S_{2}$ | $S$ | Rank
---|---|---|---|---
DUT-large | $0.72$ | $0.50$ | $0.56$ | $1$
X | $0.68$ | $0.45$ | $0.52$ | $2$
zyy | $0.69$ | $0.44$ | $0.52$ | $3$
xxxin | $0.63$ | $0.45$ | $0.51$ | $4$
AG | $0.65$ | $0.44$ | $0.50$ | $5$
福气by | $0.69$ | $0.42$ | $0.50$ | $6$
酸菜饺子&玉米饺子 | $0.68$ | $0.35$ | $0.45$ | $7$
婷之队 | $0.67$ | $0.30$ | $0.41$ | $8$
争议观点队 | $0.63$ | $0.28$ | $0.39$ | $9$
sxu-wzt | $0.64$ | $0.09$ | $0.26$ | $10$

The correlation coefficient of $S_{1}$ and $S_{2}$ is $0.5191$, and the
distribution of the two scores of different teams is illustrated in Figure 1.
Although most teams shared models across the two stages, the scores of the two
tasks are not strictly correlated. This is probably because the two tasks,
though deeply related, vary in some aspects such as the density of interactive
argument pairs.

Figure 1: Score distribution of teams with valid submissions for the second
stage.

##  5 Analysis

###  5.1 Common Techniques

Upon analyzing the models proposed by top teams, we discovered some common
techniques that many teams have adopted:

#### Task Transformation

In both stages, all teams have transformed the task into more standard ones.
In the first stage, teams treated the task as binary classification or
multiple choice question answering. In the second stage, binary classification
was more prevalent, while some teams also turned the task into a ranking
problem. Such transformation allowed teams to adopt mature solutions for
predefined standard tasks that are easier to implement and train.

#### Pretrained Language Model

All teams have utilized pre-trained language models (PLM) as their backbone,
reflecting the power of such models in NLP tasks. The selected models varied
among teams, from smaller BERT-sized models (sBERT, RoBERTa) to larger LLMs
(Qwen-14B). Some teams introduced models specifically pre-trained on law
corpora (Lawformer) or more general Chinese corpora (ERNIE) [52] to match the
tasks more closely.

#### Domain Specific Fine-tuning

Regardless of the backbone PLM selection, all teams have fine-tuned their
backbones on the training set we provided and other similar data. This domain
adaptation approach could enhance the proficiency of PLMs within the specified
scenario and task. Teams adopted different input/output formats and training
strategies depending on the model used and how they transformed the tasks.

###  5.2 Promising Tricks

Besides the above methods, teams have proposed different tricks to improve
their results. We list some of them below:

#### Data Augmentation

Some teams have involved data augmentation techniques, including pseudo-
labeling and negative sampling. These methods could alleviate the negative
effect of imbalanced labels, especially in the second stage. Meanwhile, they
could improve the generalizability and robustness of the model.

#### Adversarial Training

Besides data augmentation, several teams have used adversarial training to
fight against imbalanced labels. For instance, DUT-large, AG, and “酸菜饺子&玉米饺子”
used the Fast Gradient Method (FGM) to implement adversarial training;
Moreover, DUT-large adopted Projected Gradient Descent (PGD) for adversarial
training in the first stage.

#### Loss Function and Regularization

Since most teams treated the tasks as binary classification, they generally
chose cross entropy as the loss function. However, other loss functions were
also introduced, such as focal loss. Teams also applied regularization methods
like R-Drop to enhance model robustness.

#### Inter-task Transfer Learning

Several teams exploited their model from the first stage during the second
stage, transferring it to the new scenario. For example, AG unified the tasks
from the two stages and performed further fine-tuning to adapt to the second
stage.

#### Large Langauge Model

As mentioned, zyy utilized Qwen-14B, a large language model (LLM) for text
generation. Therefore, they have involved unique fine-tuning strategies for
LLMs such as LoRA and NefTune. They also designed prompt templates for the
task in the two stages.

#### Model Ensembling

Many teams trained multiple models and compared their performance during each
stage. Several teams have employed model ensembling to exploit all model
outputs for further improvements. They designed different methods for the
final prediction, such as direct or soft voting.

Table 6: $F_{1}$ scores of the two subtasks in the second stage. We only include teams whose model performed better (achieving a higher $S_{2}$) than our baseline in this stage. Team | $F_{1}^{\mathrm{t}}$ | $F_{1}^{\mathrm{m}}$ | $S_{2}$
---|---|---|---
DUT-large | $0.52$ | $0.44$ | $0.50$
xxxin | $0.47$ | $0.39$ | $0.45$
X | $0.47$ | $0.37$ | $0.45$
zyy | $0.45$ | $0.42$ | $0.44$
AG | $0.46$ | $0.36$ | $0.44$
福气by | $0.44$ | $0.33$ | $0.42$
Baseline | $0.42$ | $0.30$ | $0.39$

###  5.3 Multimodal Bonus Task

As stated in Section 3, we split out a multimodal bonus task in the second
stage, providing trial recordings and transcriptions for the documents in this
task. However, none of the top $7$ teams listed in Section 4.1 utilized the
multimodal data, only extracting interacting argument pairs from the given
judgment document. We conjecture that the large gap between judgment documents
and trial recordings (including transcriptions) could have made it more
difficult to exploit the multimodal resource.

Nevertheless, Table 6 shows the $F_{1}$ scores of teams listed in Table 5 that
superseded our baseline in the second stage. Most teams utilized BERT-like
PLMs and followed a similar score distribution as our baseline, where
$F_{1}^{\mathrm{m}}$ is significantly lower than $F_{1}^{\mathrm{t}}$. This
phenomenon follows our expectation because the cases from the bonus task had
different causes of action compared to other cases in this stage. In contrast,
zyy has gained an $F_{1}^{\mathrm{m}}$ much closer to $F_{1}^{\mathrm{t}}$,
using a true LLM (Qwen-14B). This phenomenon suggests that modern LLMs are
more robust in cross-domain tasks.

##  6 Conclusion

This paper presented the CAIL 2023 Argument Mining Track, a track of the
Chinese AI and Law Challenge (CAIL) 2023 that focused on identifying and
extracting interacting argument pairs in judgment documents. We divided the
track into two consecutive stages; we received $10$ valid submissions for the
second stage from $28$ teams that passed the threshold in the first stage. We
introduced the tasks and corresponding data for each stage and outlined the
best submissions. Compared to our provided baseline, teams have generally
performed better, using various techniques that are beneficial for similar
tasks in judicial scenarios.

However, there is still a long way to achieve the goal of fully automated
judicial argument comprehension. The tasks we proposed in the two stages are
simplified from real-world demands, potentially limiting their application in
actual trials. Moreover, it remains unclear how much multimodal data from
trial recordings can assist such tasks, as the top teams did not exploit the
multimodal data we provided.

Nevertheless, we thank all participating teams for taking the time to
participate in this challenging track, especially the top teams in the second
stage that provided technical reports of their submissions. We believe that
future research can use the tasks, data, and submissions from our track to
build more powerful automatic systems for judicial use.

## Acknowledgements

This work is supported by National Natural Science Foundation of China (No.
62176058) and National Key R & D Program of China (2023YFF1204800). The
project’s computational resources are supported by CFFF platform of Fudan
University.

## References

  * [1] Adrian Vermeule.  Judicial history.  The Yale Law Journal, 108(6):1311–1354, 1999.
  * [2] S. Sidney Ulmer.  Quantitative analysis of judicial processes: Some practical and theoretical applications.  Law and Contemporary Problems, 28(1):164–184, 1963.
  * [3] Stuart S. Nagel.  Applying correlation analysis to case prediction.  Texas Law Review, 42:1006, 1963.
  * [4] Octavia-Maria Sulea, Marcos Zampieri, Shervin Malmasi, Mihaela Vela, Liviu P. Dinu, and Josef van Genabith.  Exploring the use of text classification in the legal domain.  October 2017.
  * [5] Daniel Martin Katz, Michael J. Bommarito, II, and Josh Blackman.  A general approach for predicting the behavior of the supreme court of the united states.  PLOS ONE, 12(4):1–18, 04 2017.
  * [6] Yi-Hung Liu and Yen-Liang Chen.  A two-phase sentiment analysis approach for judgement prediction.  Journal of Information Science, 44(5):594–607, 2018.
  * [7] Jiaxi Cui, Zongjian Li, Yang Yan, Bohua Chen, and Li Yuan.  Chatlaw: Open-source legal large language model with integrated external knowledge bases.  June 2023.
  * [8] Shengbin Yue, Wei Chen, Siyuan Wang, Bingxuan Li, Chenchen Shen, Shujun Liu, Yuxuan Zhou, Yao Xiao, Song Yun, Xuanjing Huang, and Zhongyu Wei.  Disc-lawllm: Fine-tuning large language models for intelligent legal services.  September 2023.
  * [9] Xinyi Mou, Zhongyu Wei, and Xuanjing Huang.  Unveiling the truth and facilitating change: Towards agent-based large-scale social movement simulation.  February 2024.
  * [10] Xinnong Zhang, Haoyu Kuang, Xinyi Mou, Hanjia Lyu, Kun Wu, Siming Chen, Jiebo Luo, Xuanjing Huang, and Zhongyu Wei.  Somelvlm: A large vision language model for social media processing.  February 2024.
  * [11] Jiayu Lin, Rong Ye, Meng Han, Qi Zhang, Ruofei Lai, Xinyu Zhang, Zhao Cao, Xuanjing Huang, and Zhongyu Wei.  Argue with me tersely: Towards sentence-level counter-argument generation.  In Houda Bouamor, Juan Pino, and Kalika Bali, editors, Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing, pages 16705–16720, Singapore, December 2023. Association for Computational Linguistics.
  * [12] Jingcong Liang, Rong Ye, Meng Han, Ruofei Lai, Xinyu Zhang, Xuanjing Huang, and Zhongyu Wei.  Debatrix: Multi-dimensional debate judge with iterative chronological analysis based on llm.  March 2024.
  * [13] John Lawrence and Chris Reed.  Argument mining: A survey.  Computational Linguistics, 45(4):765–818, December 2019.
  * [14] Li Wang, Marco Lui, Su Nam Kim, Joakim Nivre, and Timothy Baldwin.  Predicting thread discourse structure over technical web forums.  In Regina Barzilay and Mark Johnson, editors, Proceedings of the 2011 Conference on Empirical Methods in Natural Language Processing, pages 13–25, Edinburgh, Scotland, UK., July 2011. Association for Computational Linguistics.
  * [15] Christian Stab and Iryna Gurevych.  Identifying argumentative discourse structures in persuasive essays.  In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 46–56, Doha, Qatar, October 2014. Association for Computational Linguistics.
  * [16] Jiangming Liu, Shay B. Cohen, and Mirella Lapata.  Discourse representation parsing for sentences and documents.  In Anna Korhonen, David Traum, and Lluís Màrquez, editors, Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 6248–6262, Florence, Italy, July 2019. Association for Computational Linguistics.
  * [17] Kaveh Taghipour and Hwee Tou Ng.  A neural approach to automated essay scoring.  In Jian Su, Kevin Duh, and Xavier Carreras, editors, Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing, pages 1882–1891, Austin, Texas, November 2016. Association for Computational Linguistics.
  * [18] Zhongyu Wei, Yang Liu, and Yi Li.  Is this post persuasive? ranking argumentative comments in online forum.  In Katrin Erk and Noah A. Smith, editors, Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers), pages 195–200, Berlin, Germany, August 2016. Association for Computational Linguistics.
  * [19] Chenhao Tan, Vlad Niculae, Cristian Danescu-Niculescu-Mizil, and Lillian Lee.  Winning arguments: Interaction dynamics and persuasion strategies in good-faith online discussions.  In Proceedings of the 25th International Conference on World Wide Web, WWW ’16, page 613–624, Republic and Canton of Geneva, CHE, 2016. International World Wide Web Conferences Steering Committee.
  * [20] Ivan Habernal and Iryna Gurevych.  Which argument is more convincing? analyzing and predicting convincingness of web arguments using bidirectional LSTM.  In Katrin Erk and Noah A. Smith, editors, Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 1589–1599, Berlin, Germany, August 2016. Association for Computational Linguistics.
  * [21] Fei Dong, Yue Zhang, and Jie Yang.  Attention-based recurrent convolutional neural network for automatic essay scoring.  In Roger Levy and Lucia Specia, editors, Proceedings of the 21st Conference on Computational Natural Language Learning (CoNLL 2017), pages 153–162, Vancouver, Canada, August 2017. Association for Computational Linguistics.
  * [22] Lu Ji, Zhongyu Wei, Xiangkun Hu, Yang Liu, Qi Zhang, and Xuanjing Huang.  Incorporating argument-level interactions for persuasion comments evaluation using co-attention model.  In Emily M. Bender, Leon Derczynski, and Pierre Isabelle, editors, Proceedings of the 27th International Conference on Computational Linguistics, pages 3703–3714, Santa Fe, New Mexico, USA, August 2018. Association for Computational Linguistics.
  * [23] Liying Cheng, Lidong Bing, Qian Yu, Wei Lu, and Luo Si.  APE: Argument pair extraction from peer review and rebuttal via multi-task learning.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 7000–7011, Online, November 2020. Association for Computational Linguistics.
  * [24] Lu Ji, Zhongyu Wei, Jing Li, Qi Zhang, and Xuanjing Huang.  Discrete argument representation learning for interactive argument pair identification.  In Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, pages 5467–5478, Online, June 2021. Association for Computational Linguistics.
  * [25] Jian Yuan, Zhongyu Wei, Yixu Gao, Wei Chen, Yun Song, Donghua Zhao, Jinglei Ma, Zhen Hu, Shaokun Zou, Donghai Li, and Xuanjing Huang.  Overview of SMP-CAIL2020-Argmine: The Interactive Argument-Pair Extraction in Judgement Document Challenge.  Data Intelligence, 3(2):287–307, June 2021.
  * [26] Fred Kort.  Predicting supreme court decisions mathematically: A quantitative analysis of the “right to counsel” cases.  American Political Science Review, 51(1):1–12, 1957.
  * [27] B. E. Lauderdale and T. S. Clark.  The supreme court’s many median justices.  American Political Science Review, 106(4):847–866, 2012.
  * [28] Bingfeng Luo, Yansong Feng, Jianbo Xu, Xiang Zhang, and Dongyan Zhao.  Learning to predict charges for criminal cases with legal basis.  In Martha Palmer, Rebecca Hwa, and Sebastian Riedel, editors, Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, pages 2727–2736, Copenhagen, Denmark, September 2017. Association for Computational Linguistics.
  * [29] Haoxi Zhong, Chaojun Xiao, Zhipeng Guo, Cunchao Tu, Zhiyuan Liu, Maosong Sun, Yansong Feng, Xianpei Han, Zhen Hu, Heng Wang, and Jianfeng Xu.  Overview of cail2018: Legal judgment prediction competition.  October 2018.
  * [30] Chaojun Xiao, Haoxi Zhong, Zhipeng Guo, Cunchao Tu, Zhiyuan Liu, Maosong Sun, Yansong Feng, Xianpei Han, Zhen Hu, Heng Wang, and Jianfeng Xu.  Cail2018: A large-scale legal dataset for judgment prediction.  July 2018.
  * [31] Chaojun Xiao, Haoxi Zhong, Zhipeng Guo, Cunchao Tu, Zhiyuan Liu, Maosong Sun, Tianyang Zhang, Xianpei Han, Zhen Hu, Heng Wang, and Jianfeng Xu.  Cail2019-scm: A dataset of similar case matching in legal domain.  November 2019.
  * [32] Chao-Lin Liu and Chwen-Dar Hsieh.  Exploring phrase-based classification of judicial documents for criminal charges in chinese.  In Proceedings of the 16th International Conference on Foundations of Intelligent Systems, ISMIS’06, page 681–690, Berlin, Heidelberg, 2006. Springer-Verlag.
  * [33] Ran Levy, Yonatan Bilu, Daniel Hershcovich, Ehud Aharoni, and Noam Slonim.  Context dependent claim detection.  In Junichi Tsujii and Jan Hajic, editors, Proceedings of COLING 2014, the 25th International Conference on Computational Linguistics: Technical Papers, pages 1489–1500, Dublin, Ireland, August 2014. Dublin City University and Association for Computational Linguistics.
  * [34] Ran Levy, Shai Gretz, Benjamin Sznajder, Shay Hummel, Ranit Aharonov, and Noam Slonim.  Unsupervised corpus–wide claim detection.  In Ivan Habernal, Iryna Gurevych, Kevin Ashley, Claire Cardie, Nancy Green, Diane Litman, Georgios Petasis, Chris Reed, Noam Slonim, and Vern Walker, editors, Proceedings of the 4th Workshop on Argument Mining, pages 79–84, Copenhagen, Denmark, September 2017. Association for Computational Linguistics.
  * [35] Yamen Ajjour, Wei-Fan Chen, Johannes Kiesel, Henning Wachsmuth, and Benno Stein.  Unit segmentation of argumentative texts.  In Ivan Habernal, Iryna Gurevych, Kevin Ashley, Claire Cardie, Nancy Green, Diane Litman, Georgios Petasis, Chris Reed, Noam Slonim, and Vern Walker, editors, Proceedings of the 4th Workshop on Argument Mining, pages 118–128, Copenhagen, Denmark, September 2017. Association for Computational Linguistics.
  * [36] Xinyi Mou, Zhongyu Wei, Changjian Jiang, and Jiajie Peng.  A two stage adaptation framework for frame detection via prompt learning.  In Nicoletta Calzolari, Chu-Ren Huang, Hansaem Kim, James Pustejovsky, Leo Wanner, Key-Sun Choi, Pum-Mo Ryu, Hsin-Hsi Chen, Lucia Donatelli, Heng Ji, Sadao Kurohashi, Patrizia Paggio, Nianwen Xue, Seokhwan Kim, Younggyun Hahm, Zhong He, Tony Kyungil Lee, Enrico Santus, Francis Bond, and Seung-Hoon Na, editors, Proceedings of the 29th International Conference on Computational Linguistics, pages 2968–2978, Gyeongju, Republic of Korea, October 2022. International Committee on Computational Linguistics.
  * [37] Vlad Niculae, Joonsuk Park, and Claire Cardie.  Argument mining with structured SVMs and RNNs.  In Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 985–995, Vancouver, Canada, July 2017. Association for Computational Linguistics.
  * [38] Orith Toledo-Ronen, Matan Orbach, Yonatan Bilu, Artem Spector, and Noam Slonim.  Multilingual argument mining: Datasets and analysis.  In Findings of the Association for Computational Linguistics: EMNLP 2020, pages 303–317, Online, November 2020. Association for Computational Linguistics.
  * [39] Andreas Peldszus and Manfred Stede.  Joint prediction in MST-style discourse parsing for argumentation mining.  In Lluís Màrquez, Chris Callison-Burch, and Jian Su, editors, Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing, pages 938–948, Lisbon, Portugal, September 2015. Association for Computational Linguistics.
  * [40] Isaac Persing and Vincent Ng.  End-to-end argumentation mining in student essays.  In Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, pages 1384–1394, San Diego, California, June 2016. Association for Computational Linguistics.
  * [41] Jianzhu Bao, Chuang Fan, Jipeng Wu, Yixue Dang, Jiachen Du, and Ruifeng Xu.  A neural transition-based model for argumentation mining.  In Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), pages 6354–6364, Online, August 2021. Association for Computational Linguistics.
  * [42] Roxanne El Baff, Henning Wachsmuth, Khalid Al Khatib, and Benno Stein.  Analyzing the Persuasive Effect of Style in News Editorial Argumentation.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 3154–3160, Online, July 2020. Association for Computational Linguistics.
  * [43] Jian Yuan, Zhongyu Wei, Donghua Zhao, Qi Zhang, and Changjian Jiang.  Leveraging argumentation knowledge graph for interactive argument pair identification.  In Chengqing Zong, Fei Xia, Wenjie Li, and Roberto Navigli, editors, Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021, pages 2310–2319, Online, August 2021. Association for Computational Linguistics.
  * [44] Jingcong Liang, Rong Ye, Meng Han, Qi Zhang, Ruofei Lai, Xinyu Zhang, Zhao Cao, Xuanjing Huang, and Zhongyu Wei.  Hi-ArG: Exploring the integration of hierarchical argumentation graphs in language pretraining.  In Houda Bouamor, Juan Pino, and Kalika Bali, editors, Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing, pages 14606–14620, Singapore, December 2023. Association for Computational Linguistics.
  * [45] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In Jill Burstein, Christy Doran, and Thamar Solorio, editors, Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota, June 2019. Association for Computational Linguistics.
  * [46] Nils Reimers and Iryna Gurevych.  Sentence-BERT: Sentence embeddings using Siamese BERT-networks.  In Kentaro Inui, Jing Jiang, Vincent Ng, and Xiaojun Wan, editors, Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP), pages 3982–3992, Hong Kong, China, November 2019. Association for Computational Linguistics.
  * [47] Yiming Cui, Wanxiang Che, Ting Liu, Bing Qin, and Ziqing Yang.  Pre-training with whole word masking for chinese bert.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 29:3504–3514, 2021.
  * [48] Jinze Bai, Shuai Bai, Yunfei Chu, Zeyu Cui, Kai Dang, Xiaodong Deng, Yang Fan, Wenbin Ge, Yu Han, Fei Huang, Binyuan Hui, Luo Ji, Mei Li, Junyang Lin, Runji Lin, Dayiheng Liu, Gao Liu, Chengqiang Lu, Keming Lu, Jianxin Ma, Rui Men, Xingzhang Ren, Xuancheng Ren, Chuanqi Tan, Sinan Tan, Jianhong Tu, Peng Wang, Shijie Wang, Wei Wang, Shengguang Wu, Benfeng Xu, Jin Xu, An Yang, Hao Yang, Jian Yang, Shusheng Yang, Yang Yao, Bowen Yu, Hongyi Yuan, Zheng Yuan, Jianwei Zhang, Xingxuan Zhang, Yichang Zhang, Zhenru Zhang, Chang Zhou, Jingren Zhou, Xiaohuan Zhou, and Tianhang Zhu.  Qwen technical report.  September 2023.
  * [49] Edward J. Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen.  LoRA: Low-rank adaptation of large language models.  In International Conference on Learning Representations, 2022.
  * [50] Chaojun Xiao, Xueyu Hu, Zhiyuan Liu, Cunchao Tu, and Maosong Sun.  Lawformer: A pre-trained language model for chinese legal long documents.  AI Open, 2:79–84, 2021.
  * [51] Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized bert pretraining approach.  July 2019.
  * [52] Yu Sun, Shuohuan Wang, Shikun Feng, Siyu Ding, Chao Pang, Junyuan Shang, Jiaxiang Liu, Xuyi Chen, Yanbin Zhao, Yuxiang Lu, Weixin Liu, Zhihua Wu, Weibao Gong, Jianzhong Liang, Zhizhou Shang, Peng Sun, Wei Liu, Xuan Ouyang, Dianhai Yu, Hao Tian, Hua Wu, and Haifeng Wang.  Ernie 3.0: Large-scale knowledge enhanced pre-training for language understanding and generation.  July 2021.

##  Appendix A Annotation Details

The dataset for SMP-CAIL2020-Argmine (including its extension in 2021) only
covered $5$ causes of action. Meanwhile, many cases did not have a specific
cause of action (tagged as “其他”). First, we excluded the above vague cases in
CAIL2023-Argmine to ensure the quality of crime features. Next, to increase
the diversity of the dataset, we selected $9$ causes of action that had no or
few cases in the original dataset, including criminal and civil ones, as
listed in Table 7. We collected public judgment documents and filtered them
according to the rules proposed by [25]; for each cause of action, if there
are more than $300$ valid cases, we sampled $300$ of them.

Table 7: Statistics of the annotation for new cases in CAIL2023-Argmine. Causes of action that end with “罪” are criminal causes, while others are civil ones. “Agreement” and “Disagreement” denote agreement and disagreement argument pairs, respectively. We define agreement and disagreement pairs as in [25]: we only count argument pairs where the defense fully acknowledges the plaintiff as agreement ones. Cause of Action | # Documents | # Agreement | # Disagreement | Cohen’s $\kappa$
---|---|---|---|---
故意杀人罪 | $272$ | $604$ | $762$ | $0.18$
抢劫罪 | $253$ | $794$ | $904$ | $0.26$
强奸罪 | $239$ | $414$ | $1,087$ | $0.28$
诈骗罪 | $168$ | $602$ | $653$ | $0.25$
盗窃罪 | $55$ | $319$ | $150$ | $0.21$
借款合同纠纷 | $299$ | $246$ | $1,107$ | $0.15$
离婚纠纷 | $298$ | $282$ | $1,216$ | $0.24$
生命权、身体权、健康权纠纷 | $294$ | $82$ | $1,645$ | $0.28$
买卖合同纠纷 | $293$ | $209$ | $1,315$ | $0.22$
Summary | $2,171$ | $3,552$ | $8,839$ | $0.24$

Next, we arranged for a group of undergraduate and graduate students majoring
in law to annotate these documents. We only kept the part containing
testimonies from both sides and extracted the sentences from the document
beforehand to unify the argument span. The annotation schema was the same as
in [25], where annotators needed to choose the case type, the cause of action,
the plaintiff and defense, and a list of interacting argument pairs, except
that they only needed to select pre-extracted sentences as argument
candidates. Two annotators annotated each document. If the list of selected
argument pairs differed between the annotators, an adjudicator would decide on
the valid pairs among all candidates. We invited a law professor as the
adjudicator to ensure the annotation quality.

After annotating all documents, we further filtered documents whose annotated
cause of action (by both annotators) did not match its presumed one because it
indicates that the cause of action could be vague in this case. Table 7 shows
annotation statistics of the final documents per cause of action. In most
cases, there are more disagree pairs than agree pairs, especially in civil
cases; the only exception is lancery (“盗窃罪”), where the evidence is usually
undoubtedly clear. The Cohen’s $\kappa$ is not very high, making adjudication
critical during annotation.


# Optimized Implementation for Calculation and Fast-Update of Pfaffians
Installed to the Open-Source Fermionic Variational Solver mVMC

RuQing G. Xu<EMAIL_ADDRESS>Tsuyoshi Okubo  Synge Todo
Masatoshi Imada  Department of Physics, University of Tokyo, Hongo 7-3-1,
Bunkyo, Tokyo, 113-0033, Japan  Institute for Physics of Intelligence,
University of Tokyo, Hongo 7-3-1, Bunkyo, Tokyo, 113-0033, Japan  JST, PRESTO,
Honcho 4-1-8, Kawaguchi, Saitama, 332-0012, Japan  Institute for Solid State
Physics, University of Tokyo, Kashiwanoha 5-1-5, Kashiwa, Chiba,277-8581,
Japan  Deptartment of Applied Physics, Waseda University, Okubo 3-4-1,
Shinjuku, Tokyo, 169-8555, Japan  Toyota RIKEN, Yokomichi 41-1, Nagakute,
Aichi, 480-1118, Japan

###### Abstract

In this article, we present a high performance, portable and well templated
implementation for computing and fast-updating Pfaffian and inverse of an
even-ranked skew-symmetric (antisymmetric) matrix. It is achieved with a skew-
symmetric, blocked variant of the Parlett-Reid algorithm and a blocked update
scheme based on the Woodbury matrix identity. Installation of this framework
into the geminal-wavefunction-based _many-variable Variational Monte Carlo_
(mVMC) code boosts sampling performance to up to more than $6$ times without
changing Markov chain’s behavior. The implementation is based on an extension
of the _BLAS-like instantiation software_ (BLIS) framework which has optimized
kernel for many state-of-the-art processors including Intel Skylake-X, AMD
EPYC Rome and Fujitsu A64FX.

###### keywords:

Variational Monte Carlo, Ground-state Method, Quantum Lattice Model, Skew-
symmetric Matrix, Pfaffian, Blocked Algorithm, LAPACK, BLAS

###### MSC:

[2020] 15A90, 65Z05, 81-08

††journal: Computater Physics
Communicationsmytitlenotemytitlenotefootnotetext: Intro for software package
Pfaffine and its installation to mVMC with additional fast-update techniques.

##  1 Program Summary

Program title: |  Pfaffine and PfUpdates library for mVMC[1]
---|---
CPC Library link to program: |  TBA
Developer’s repository link: |  https://github.com/issp-center-dev/mVMC/tree/master/src/pfupdates
|  https://github.com/xrq-phys/Pfaffine
Code Ocean capsule: |  TBA
Licensing provisions: |  MPL-2.0 (for new Library part)
Programming language: |  C++14 (for new Library part)
Nature of problem: |  Finding a method for computing and updating Pfaffian and inverse of a skew-symmetric matrix that yields a high performance on modern processor architectures.
Solution method: |  Deploying a blocked version of the Parlett-Reid algorithm with BLIS serving as assembly-level backend. Updating is approached using a modified Woodbury matrix identity.

##  2 Introduction

Numerical variational methods for correlated many-body systems have played an
important role in studies of novel condensed-matter phenomena such as high-
temperature superconductivity and spin liquids[2, 3, 4, 5, 6, 7, 8, 9].
Compared to exponentially expensive methods like exact diagonalization[10],
variational methods retains advantage that it demands polynomial computational
complexity, yet is able to provide results of decent accuracy. For Fermionic
systems, there are several popular formulations of variational approaches.
Density matrix renormalization group (DMRG)[2], which later evolved into the
matrix product state (MPS)[11] ansatz, has achieved tremendous success in one-
dimensional systems but suffers from lack of expressiveness in higher
dimensions, while its extension for two-dimensional systems called projected
entangled pair state (PEPS)[12, 13, 14] proved to be very challenging to
optimize due to local minima problems[15, 16, 17]. On the contrary, the
variational Monte Carlo method is easy to evaluate and optimize, while still
capable to express those novel quantum states well[18, 19, 20, 21, 22].

During the past few years, several packages for variational Monte Carlo have
been developed and opened to the public[23, 24], among them is mVMC[1] which
focuses on the study of lattice-based systems. In mVMC, the variational
wavefunction with $N_{\mathrm{pair}}$ spin-up Fermions and the same number of
spin-down Fermions is taken to be the form:

$\left|{\Psi}\right\rangle=\mathcal{L}_{s}\mathcal{L}_{p}\mathcal{P}_{J}\mathcal{P}_{G}\left|{\phi}\right\rangle.$
(1)

In eq. 1, the core part $\left|{\phi}\right\rangle$ is given by:

$\left|{\phi}\right\rangle=\left(\sum_{i<j}f_{ij}c_{i}^{\dagger}c_{j}^{\dagger}\right)^{N_{\mathrm{pair}}}\left|0\right\rangle,$
(2)

usually referred as a “pair product” (PP) or “geminal wavefunction” with
$c_{i}^{\dagger}$ as the Fermion creation operator on the $i$th spin-orbital,
while $\mathcal{P}_{J},\mathcal{P}_{G}$ are real-space correlation operators
with Jastrow and Gutzwiller factors, respectively and
$\mathcal{L}_{s},\mathcal{L}_{p}$ are quantum-number projection operators for
spin and momentum, respectively, to restore the original symmetry of the
Hamiltonian. Readers interested in
$\mathcal{P}_{J},\mathcal{P}_{G},\mathcal{L}_{s}$ and $\mathcal{L}_{p}$ are
advised to refer to [18] and the rest sections of this article only focus on
PP. With mean-field based $\left|{\phi}\right\rangle$[18] providing a very
good initial state and $\mathcal{P}_{J},\mathcal{P}_{G}$ further considering
effects introduced by interactions like Coulomb repulsion, mVMC is able to
grab important features of many-body quantum mechanical systems[19, 20],
providing excellent approximations to their ground states as well as low-
energy excited states [25, 26].

Within the process of calculation, mVMC samples real space configurations
$\left\\{\left|{x^{(k)}}\right\rangle\right\\}$ using a Markov chain (MC)
process. The probability for some configuration $\left|{x}\right\rangle$ to be
sampled is proportional to its overlap with the wavefunction ansatz, i.e.:

$\mathrm{Prob}\left(\left|{x}\right\rangle\textrm{ is
sampled}\right)=\frac{{\left\langle{x}\middle|{\Psi}\right\rangle}^{2}}{\left\langle{\Psi}\middle|{\Psi}\right\rangle}.$
(3)

According to the detailed balance condition of the Markov chain, eq. 3 yields
a stochastic process dominated by:

$\mathrm{Prob}\left(\left|{x}\right\rangle\longrightarrow\left|{x^{\prime}}\right\rangle\right)=\min\left(1,\frac{{\left\langle{x^{\prime}}\middle|{\Psi}\right\rangle}^{2}}{{\left\langle{x}\middle|{\Psi}\right\rangle^{2}}}\right).$
(4)

As $\mathcal{P}_{J}$ and $\mathcal{P}_{G}$ are diagonal in configuration
representation, while $\mathcal{L}_{p}$ and $\mathcal{L}_{s}$ expands to
summed spatial operators, $\left\langle{x}\middle|{\Psi}\right\rangle$ expands
to:

$\left\langle{x}\middle|{\Psi}\right\rangle=\sum_{\theta}\beta_{\theta}\sum_{\mathbf{r}}\alpha_{\mathbf{r}}\left\langle{T_{\mathbf{r}}(x)}\middle|{R_{\theta}(\phi)}\right\rangle.$
(5)

Here $\left|R_{\theta}(\phi)\right\rangle$ is $\left|\phi\right\rangle$
rotated by angle $\theta$ and $\left|T_{\mathbf{r}}(x)\right\rangle$ is
$\left|x\right\rangle$ spatially translated by vector $\mathbf{r}$. Meshes for
$(\theta,\mathbf{r})$ and their corresponding weights are defined so as the
original symmetry of Hamiltonian can be correctly restored. We recommended
readers to refer to [18, 1] for further details of those quantities. Equation
5 tells us that to generate samples satisfying eq. 3 we have to compute:

  * 1.

$\left\\{\left\langle{T_{r}(x^{(0)})}\middle|{R_{\theta}(\phi)}\right\rangle\middle|r,\theta\right\\}$
for the first sample;

  * 2.

$\left\\{\left\langle{T_{r}(x^{(k)})}\middle|{R_{\theta}(\phi)}\right\rangle\middle|r,\theta\right\\}$
with
$\left\\{\left\langle{T_{r}(x^{(k-1)})}\middle|{R_{\theta}(\phi)}\right\rangle\middle|r,\theta\right\\}$
known;

which are exactly the hotspots of the whole mVMC code.

Because $\left|{R_{\theta}(\phi)}\right\rangle$ generates a constant set of
pair products
$\left|{R_{\theta}(\phi)}\right\rangle=\left(\sum_{i<j}f_{ij}^{\theta}c_{i}^{\dagger}c_{j}^{\dagger}\right)^{N_{\mathrm{pair}}}\left|{0}\right\rangle$
that does not change throughout the MC sampling process, and
$\left|{T_{r}(x)}\right\rangle$ can be viewed as from another chain of samples
with a constant displacement, we omit these transformation operators in later
discussion and only use $\left\langle{x}\middle|{\phi}\right\rangle$ which is
evaluated as:

$\left\langle{x}\middle|{\phi}\right\rangle=\mathrm{Pf}\,\mathbf{X}\textrm{
where }X_{\alpha\beta}=\left\\{\begin{array}[]{ll}0&x_{\alpha}=x_{\beta}\\\
f_{x_{\alpha}x_{\beta}}&x_{\alpha}<x_{\beta}\\\
-f_{x_{\beta}x_{\alpha}}&\text{otherwise}\end{array}\right.,$ (6)

where $x_{\alpha}$ denotes spin-orbital index of the $\alpha$-th Fermion under
real-space configuration $\left|{x}\right\rangle$.

Current implementation of mVMC[1] already features a nice threading behavior,
making it able to run on large-scale computer clusters. However, this
implementation bears defective microarchitecture usage on modern processor
hardware due to the following two reasons:

  * 1.

Heaviest calculation of the mVMC method involves operations on skew-symmetric
(antisymmetric) matrices, which is not available via the standard BLAS
interface. The alternative library, known as Pfapack77[27], is programmed
without considering pipeline-level hardware efficiency hence not optimal for
those state-of-the-art processors.

  * 2.

To avoid recalculating Pfaffian every time for each Markov chain update, a
fast-update scheme is adopted in current mVMC implementation. However, only
rank-1 update scheme is implemented. The rank-1 update strategy is, by its
name, mainly invoking rank-1 updates of matrices
$\left(\mathbf{X}=\mathbf{X}+\mathbf{uv}^{T}\right)$, which is doomed to hit
the memory bound and impact performance.

To get around these two severe hurdles upon mVMC’s performance, a new fast-
update scheme as well as a new implementation of the skew-symmetric matrix
library is necessary so that:

  * 1.

All computational hotspots are effectively handled with assemblies fully
optimized for the hosting hardware.

  * 2.

Fast-update of Pfaffian invokes matrix-matrix operations instead of rank-1
updates.

Such are exactly what our optimization does.

This article is organized in the following fashion: Section 3 reviews and
extends Wimmer’s approach in Pfapack77[27] on tridiagonalizing skew-symmetric
matrices and computing Pfaffians. Section 4 presents the blocked fast-update
scheme for further speeding up Pfaffian updates. Section 5 specifies
programming interfaces of this acceleration framework as well as our assembly-
level tuned implementation. Section 6 shows benchmarking results indicating
superiority of this new implementation and finally, Sec. 7 concludes our
development work.

##  3 Skew-symmetric Matrix and Pfaffian

Pfaffian of a size $2n\times 2n$ skew-symmetric matrix $\mathbf{X}$ is defined
as:

$\mathrm{Pf}\,\mathbf{X}=\frac{1}{2^{n}n!}\sum_{\sigma\in
S_{2n}}\mathrm{sgn}(\sigma)\prod_{i}^{n}X_{\sigma(2i-1),\sigma(2i)},$ (7)

and satisfies:

$\left(\mathrm{Pf}\,\mathbf{X}\right)^{2}=\det\mathbf{X},$

where $S_{2n}$ is the permutation group of sets with $2n$ elements. Directly
evaluating eq. 7 is of exponential complexity due to the number of elements in
$S_{2n}$. However, there are linear transformations that can be applied to
simplify computation just like the case of calculating determinants.

Before introducing tridiagonalization computation of Pfaffian, let us first go
over definition of the Gaussian elimination matrix:

$\mathbf{M}_{k}=\mathbf{I}_{2n}-\boldsymbol{\alpha}_{k}\mathbf{e}_{k}^{T},$
(8)

where $\mathbf{I}_{2n}$ is the $2n\times 2n$ identity, $\mathbf{e}_{k}$ is the
$k$-th elementary vector of Euclidean space (i.e.
$\mathbf{e}_{k}=(\underbrace{0,\cdots,0}_{k-1\text{
zeros}},1,\underbrace{0,\cdots,0}_{2n-k\text{ zeros}})^{T}$) and
$\boldsymbol{\alpha}_{k}$ is constructed with entries:

$\left(\alpha_{k}\right)_{i}=\left\\{\begin{array}[]{ll}0&i\leq k\\\
\left.X_{i,k-1}\middle/X_{k,k-1}\right.&i>k\end{array}\right.$ (9)

so that left-applying $\mathbf{M}_{k}$ to $\mathbf{X}$ eliminates all
$\left.X_{i,k-1}\right|_{i>k}$ with $\left.X_{ij}\right|_{i\leq k,\forall j}$
untouched. Applying $\mathbf{M}_{k}$ to both sides of a skew-symmetric matrix
$\mathbf{X}$ eliminates non-tridiagonal part of the $\left(k-1\right)$-th
column and row without changing its Pfaffian[27]. This will be the main weapon
for efficient Pfaffian calculation in the formulation of this section.

### Blocked Tridiagonal Decomposition

Blocked tridiagonal decomposition of skew-symmetric matrices (sktdf) was first
formulated by Wimmer in Pfapack77[27]. In their work a blocked and skew-
symmetric variant of the Parlett-Reid algorithm[28] is adopted so that:

$\displaystyle\mathbf{P}\mathbf{X}\mathbf{P}^{T}$
$\displaystyle=\mathbf{L}\mathbf{T}\mathbf{L}^{T},$ (10)
$\displaystyle\mathbf{P}$
$\displaystyle=\mathbf{P}_{n-1}\cdots\mathbf{P}_{2},$
$\displaystyle\mathbf{T}$
$\displaystyle=\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}\cdots\mathbf{M}_{\left[2::2+r^{(1)}\right]}\mathbf{X}\mathbf{M}_{\left[2::2+r^{(1)}\right]}^{T}\cdots\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}^{T}$
$\displaystyle\mathbf{L}$
$\displaystyle=\left(\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}\cdots\mathbf{M}_{\left[2::2+r^{(1)}\right]}\mathbf{P}^{T}\right)^{-1}$
$\displaystyle\mathbf{M}_{\left[i::i+r^{(t)}\right]}$
$\displaystyle=\mathbf{M}_{i+r^{(t)}}\mathbf{P}_{i+r^{(t)}}\cdots\mathbf{M}_{i+1}\mathbf{P}_{i+1}\mathbf{M}_{i}\mathbf{P}_{i}.$

Here $\mathbf{P}_{k}$ are permutation matrices that swap the $k$-th and
$\left(p(k)>k\right)$-th columns (rows) of $\mathbf{X}$ when multiplied on the
right (left) side. We call linear transformation
$\mathbf{M}_{\left[i::i+r^{(t)}\right]}$ a rank-$r^{(t)}$ pivoted Gaussian
elimination where $r^{(t)}\in\mathbb{Z}^{+}$ is a parameter for performance
tuning. It is constructed by merging multiple Gaussian elimination matrices:

$\displaystyle\mathbf{M}_{\left[k::k+1\right]}$
$\displaystyle=\mathbf{M}_{k+1}\mathbf{P}_{k+1}\mathbf{M}_{k}$
$\displaystyle=\left(\mathbf{I}_{2n}-\boldsymbol{\alpha}_{k+1}\mathbf{e}_{k+1}^{T}\right)\mathbf{P}_{k+1}\left(\mathbf{I}_{2n}-\boldsymbol{\alpha}_{k}\mathbf{e}_{k}^{T}\right)$
$\displaystyle=\left(\mathbf{I}_{2n}-\underbrace{\begin{bmatrix}\boldsymbol{\alpha}_{k+1}&\boldsymbol{\alpha}^{(p)}_{k}\end{bmatrix}}_{\mathbf{V}_{\left[k::k+1\right]}}\underbrace{\begin{bmatrix}\mathbf{e}_{k+1}^{T}\\\
\mathbf{e}_{k}^{T}\end{bmatrix}}_{\mathbf{E}_{\left[k::k+1\right]}}\right)\mathbf{P}_{k+1}$
(11) $\displaystyle\mathbf{M}_{\left[k-r::k+1\right]}$
$\displaystyle=\mathbf{M}_{k+1}\mathbf{P}_{k+1}\mathbf{M}_{\left[k-r::k\right]}$
$\displaystyle=\left(\mathbf{I}_{2n}-\boldsymbol{\alpha}_{k+1}\mathbf{e}_{k+1}^{T}\right)\mathbf{P}_{k+1}\left(\mathbf{I}_{2n}-\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}^{T}\right)\mathbf{P}_{k}\cdots\mathbf{P}_{k-r}$
$\displaystyle=\left(\mathbf{I}_{2n}-\underbrace{\begin{bmatrix}\boldsymbol{\alpha}_{k+1}&\mathbf{V}_{\left[k-r::k\right]}^{(p)}\end{bmatrix}}_{\mathbf{V}_{\left[k-r::k+1\right]}}\underbrace{\begin{bmatrix}\mathbf{e}_{k+1}^{T}\\\
\mathbf{E}_{\left[k-r::k\right]}^{T}\end{bmatrix}}_{\mathbf{E}_{\left[k-r::k+1\right]}}\right)\mathbf{P}_{k+1}\mathbf{P}_{k}\cdots\mathbf{P}_{k-r},$
(12)

where $\mathbf{V}_{\left[...k\right]}^{(p)}$ and
$\boldsymbol{\alpha}^{(p)}_{k}$ denote a matrix and column vector with $k$-th
and $p(k)$-th rows swapped. Note also that in derivation of eqs. 11 and 12 we
used $\mathbf{e}_{k+1}^{T}\boldsymbol{\alpha}_{k}=0$ and
$\mathbf{e}_{k}=\mathbf{P}_{k+r}\mathbf{e}_{k}$ implied by eq. 8. We refer to
[27] for further details of this algorithm.

Once factorization (10) is done, with the fact that neither $\mathbf{P}_{k}$
or $\mathbf{M}_{k}$ affects the matrix’s Pfaffian, $\mathrm{Pf}\,\mathbf{X}$
can be computed as:

$\mathrm{Pf}\,\mathbf{X}=\mathrm{Pf}\,\mathbf{T}=\prod_{i=1}^{n}T_{2i-1,2i}.$
(13)

### Inversion from Tridiagonalization

With eq. 13 completes Wimmer’s derivation in [27], but mVMC requires
computation of $\mathbf{X}^{-1}$ in addition to $\mathrm{Pf}\,\mathbf{X}$.
Current version of the implementation is initializing another GETRF/GETRI call
to invert $\mathbf{X}$ as a non-symmetric matrix. This treatment requires
$\sim(2n)^{3}$ multiply-accumulate operations (MACs). However, one might
notice that by inverting eq. 10, $\mathbf{X}^{-1}$ can be expressed as:

$\displaystyle\mathbf{P}\mathbf{X}^{-1}\mathbf{P}^{T}$
$\displaystyle=\left(\mathbf{L}^{-1}\right)^{T}\mathbf{T}^{-1}\mathbf{L}^{-1}$
(14)
$\displaystyle=\mathbf{P}\mathbf{M}_{\left[2::2+r^{(1)}\right]}^{T}\cdots\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}^{T}\mathbf{T}^{-1}\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}\cdots\mathbf{M}_{\left[2::2+r^{(1)}\right]}P^{T}$
$\displaystyle\Rightarrow\mathbf{X}^{-1}$
$\displaystyle=\mathbf{M}_{\left[2::2+r^{(1)}\right]}^{T}\cdots\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}^{T}\mathbf{T}^{-1}\mathbf{M}_{\left[n-1-r^{(s)}::n-1\right]}\cdots\mathbf{M}_{\left[2::2+r^{(1)}\right]},$
(15)

so that one can either invert $\mathbf{X}$ by expanding eq. 15 from the center
or by directly evaluating eq. 14 since $\mathbf{L}$ is obtainable by properly
merging and swapping all $\mathbf{V}$ entries of $\mathbf{M}$[27].

### Quantitative Cost of Computation

Evaluation of
$\mathbf{M}_{\left[k-r::k\right]}\mathbf{X}^{(k-r-1)}\mathbf{M}_{\left[k-r::k\right]}^{T}$
can be done as:

$\displaystyle\mathbf{M}_{\left[k-r::k\right]}\mathbf{X}^{(k-r-1)}\mathbf{M}_{\left[k-r::k\right]}^{T}$
$\displaystyle=\left(\mathbf{I}-\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}^{T}\right)\mathbf{X}^{(k-r-1,p)}\left(\mathbf{I}-\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}^{T}\right)^{T}$
$\displaystyle=\mathbf{X}^{(k-r-1,p)}+\mathbf{V}_{\left[k-r::k\right]}\mathbf{X}^{(k-r-1,p)T}_{\left[k-r::k\right],[::]}-\mathbf{X}^{(k-r-1,p)}_{\left[k-r::k\right],[::]}\mathbf{V}_{\left[k-r::k\right]}^{T},$
(16)

where we have used during the derivation that $\mathbf{X}$ is skew-symmetric
and $\mathbf{V}$’s columns are orthogonal to each other.
$\mathbf{X}^{(k-r-1)}$ denotes $\mathbf{X}$ after $k-r-1$ steps of Gaussian
elimination and $\mathbf{X}^{(k-r-1,p)}$ its permuted version, i.e.111 Here
for ease notation, an unblock expression is used. In real computation we
always adopt blocked manipulations. :

$\displaystyle\mathbf{X}^{(k-r-1)}$
$\displaystyle=\mathbf{M}_{k-r-1}\mathbf{P}_{k-r-1}\cdots\mathbf{M}_{2}\mathbf{P}_{2}\mathbf{X}^{(k-r-1)}\mathbf{P}_{2}^{T}\mathbf{M}_{2}^{T}\cdots\mathbf{P}_{k-r-1}^{T}\mathbf{M}_{k-r-1}^{T},$
$\displaystyle\mathbf{X}^{(k-r-1,p)}$
$\displaystyle=\mathbf{P}_{k}\cdots\mathbf{P}_{k-r}\mathbf{X}^{(k-r-1)}\mathbf{P}_{k-r}^{T}\cdots\mathbf{P}_{k}^{T}.$

As $\mathbf{X}^{(k-r-1,p)}$ always has the first $k-r-1$ columns and rows
tridiagonal, evaluation of eq. 16 only needs to consider the lower-right
portion of size $\left(2n-(k-r-1)\right)\times\left(2n-(k-r-1)\right)$. This
fact together with $\mathbf{X}$’s skew-symmetric property yields an evaluation
cost of $(2n-(k-r-1))^{2}r$. Hence the total evaluation cost to transform
$\mathbf{X}$ to $\mathbf{T}$ in eq. 10 is around $\frac{1}{3}(2n)^{3}$ MACs.

For extracting $\mathbf{X}^{-1}$, instantiation of eq. 15 requires iteratively
evaluating
$\mathbf{M}_{\left[k-r::k\right]}^{T}\mathbf{T}^{-1}_{(k-r-1)}\mathbf{M}_{\left[k-r::k\right]}$
which expands into:

$\displaystyle\mathbf{M}_{\left[k-r::k\right]}^{T}\mathbf{T}^{-1}_{(k-r-1)}\mathbf{M}_{\left[k-r::k\right]}$
$\displaystyle=\mathbf{P}_{k-r}\cdots\mathbf{P}_{k}\left(\mathbf{I}_{2n}-\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}^{T}\right)^{T}\mathbf{T}_{(k-r-1)}^{-1}\left(\mathbf{I}_{2n}-\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}^{T}\right)\mathbf{P}_{k}\cdots\mathbf{P}_{k-r}$
$\displaystyle=\mathbf{P}_{k-r}\cdots\mathbf{P}_{k}\Big{\\{}\mathbf{T}_{(k-r-1)}^{-1}-\mathbf{E}_{\left[k-r::k\right]}\left(\mathbf{T}_{(k-r-1)}^{-1}\mathbf{V}_{\left[k-r::k\right]}\right)^{T}$
$\displaystyle\quad+\left(\mathbf{T}_{(k-r-1)}^{-1}\mathbf{V}_{\left[k-r::k\right]}\right)\mathbf{E}_{\left[k-r::k\right]}+\mathbf{E}_{\left[k-r::k\right]}\mathbf{V}_{\left[k-r::k\right]}^{T}\mathbf{T}_{(k-r-1)}^{-1}\mathbf{V}_{\left[k-r::k\right]}\mathbf{E}_{\left[k-r::k\right]}\Big{\\}}\mathbf{P}_{k}\cdots\mathbf{P}_{k-r}.$
(17)

Computation of matrix-matrix products in eq. 17 requires around
$(2n)\left(2n-\left(k-r-1\right)\right)r+(2n)r^{2}$ MACs. Ignoring $(2n)r^{2}$
term which is much less than $(2n)\left(2n-\left(k-r-1\right)\right)r$, total
cost of eq. 15 becomes $\approx\frac{1}{2}(2n)^{3}$.

The alternative approach for $\mathbf{X}^{-1}$, namely eq. 14 requires
inverting $\mathbf{L}^{-1}$ first, calling a subprocedure of GETRI named TRTRI
which makes $\frac{1}{6}(2n)^{3}$ MACs. Since $\mathbf{T}$ is tridiagonal,
computation of $\mathbf{T}^{-1}\mathbf{L}^{-1}$ can be done by using the
$O\left(n^{2}\right)$-cost procedure GTSV whose cost can be neglected. The
final
$\mathbf{X}^{-1}=\left(\mathbf{L}^{-1}\right)^{T}\left(\mathbf{T}^{-1}\mathbf{L}^{-1}\right)$
step has triangular $\mathbf{L}^{-1}$ on the left hand side and only needs to
write half of $\mathbf{X}^{-1}$ since we already know that it is skew-
symmetric, resulting in $\frac{1}{6}(2n)^{3}$ MACs. Hence the whole evaluation
of eq. 14 totals up to $\frac{1}{3}(2n)^{3}$ MACs, also lower than GETRF/GETRI
deployed by the original algorithm. Both approaches for $\mathbf{X}^{-1}$ are
available in our implementation as different subversions.

##  4 Blocked Pfaffian Update from Modified Woodbury Matrix Identity

### Theory

mVMC takes configuration samples
$\left\\{\left|{x^{(k)}}\right\rangle\middle|k\in\mathbb{N}\right\\}$ with a
Markov chain (MC) process, as is shown in eq. 4. As each step
$\left|{x^{(k)}}\right\rangle\rightarrow\left|{x^{(k+1)}}\right\rangle$ within
this MC only changes the spin-orbitals of one or two electrons of
$\left|{x^{(k)}}\right\rangle$, we do not need to redo Sec. 3 calculations
each time a new sample is generated. Instead, from the Woodbury matrix
identity[29] we have:

$\displaystyle\mathrm{Pf}\left[\mathbf{X}+\mathbf{B}\mathbf{C}\mathbf{B}^{T}\right]$
$\displaystyle=\mathrm{Pf}\,\mathbf{X}\times\frac{\mathrm{Pf}\left[\mathbf{C}^{-1}+\mathbf{B}^{T}\mathbf{X}^{-1}\mathbf{B}\right]}{\mathrm{Pf}\left[\mathbf{C}^{-1}\right]}$
(18)
$\displaystyle\left(\mathbf{X}+\mathbf{B}\mathbf{C}\mathbf{B}^{T}\right)^{-1}$
$\displaystyle=\mathbf{X}^{-1}-\mathbf{X}^{-1}\mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{B}^{T}\mathbf{X}^{-1}\mathbf{B}\right)^{-1}\mathbf{B}^{T}\mathbf{X}^{-1},$
(19)

which allows Pfaffian in eq. 6 to be updated as (To give a more intuitive
image of what is happening during the MC fast-update, we also plug in here a
$4\times 4$ exemplar matrix for $\mathbf{X}$.):

$\displaystyle\left|{x}\right\rangle$
$\displaystyle=\left|{{j_{1}}{j_{2}}{j_{3}}{j_{4}}}\right\rangle$ (20)
$\displaystyle\mathbf{X}$
$\displaystyle=\left[\begin{array}[]{cccc}0&F_{j_{1}}^{j_{2}}&F_{j_{1}}^{j_{3}}&F_{j_{1}}^{j_{4}}\\\
-F_{j_{1}}^{j_{2}}&0&F_{j_{2}}^{j_{3}}&F_{j_{2}}^{j_{4}}\\\
-F_{j_{1}}^{j_{3}}&-F_{j_{2}}^{j_{3}}&0&F_{j_{3}}^{j_{4}}\\\
-F_{j_{1}}^{j_{4}}&-F_{j_{2}}^{j_{4}}&-F_{j_{3}}^{j_{4}}&0\end{array}\right]$
(25) $\displaystyle\left\langle{x}\middle|{\phi}\right\rangle$
$\displaystyle=\mathrm{Pf}\,\mathbf{X}$ (26) $\displaystyle\downarrow$
$\displaystyle\left|{x^{(1)}}\right\rangle$
$\displaystyle=\left|{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}{j_{2}}{j_{3}}{j_{4}}}\right\rangle$
(27) $\displaystyle\mathbf{X}^{(1)}$
$\displaystyle=\left[\begin{array}[]{cccc}0&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{2}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{4}}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{2}}}&0&F_{j_{2}}^{j_{3}}&F_{j_{2}}^{j_{4}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{3}}}&-F_{j_{2}}^{j_{3}}&0&F_{j_{3}}^{j_{4}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{4}}}&-F_{j_{2}}^{j_{4}}&-F_{j_{3}}^{j_{4}}&0\end{array}\right]=\mathbf{X}+\mathbf{B}^{(1)}\underbrace{\left[\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right]}_{\mathbf{C}^{(1)}}\mathbf{B^{(1)}}^{T}$ (34)
$\displaystyle B^{(1)}_{lm}$
$\displaystyle=\left\\{\begin{array}[]{ll}F_{i_{1}}^{j_{l}}-F_{j_{1}}^{j_{l}}&l\neq
1,m=1\\\ 1&l=1,m=1\\\ 0&\textrm{otherwise}\end{array}\right.$ (38)
$\displaystyle\left\langle{x^{(1)}}\middle|{\phi}\right\rangle$
$\displaystyle=\mathrm{Pf}\,\mathbf{X}^{(1)}=\mathrm{Pf}\,\mathbf{X}\times\left.\mathrm{Pf}\left(\mathbf{C}^{(1)}+\mathbf{B^{(1)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(1)}\right)\right/\mathrm{Pf}\,\mathbf{C}^{(1)}=\mathrm{Pf}\,\mathbf{X}\times\left(1-\sum_{j}\left(X^{-1}\right)_{1}^{j}u_{j}\right)$
(39) $\displaystyle\mathbf{X^{(1)}}^{-1}$
$\displaystyle=\mathbf{X}^{-1}-\mathbf{X}^{-1}\mathbf{B}^{(1)}\left(\mathbf{C}^{(1)}+\mathbf{B^{(1)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(1)}\right)^{-1}\mathbf{B^{(1)}}^{T}\mathbf{X}^{-1}$
(40) $\displaystyle\downarrow$

$\displaystyle\left|{x^{(2)}}\right\rangle$
$\displaystyle=\left|{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}{j_{3}}{j_{4}}}\right\rangle$
(41) $\displaystyle\mathbf{X}^{(2)}$
$\displaystyle=\left[\begin{array}[]{cccc}0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}F_{i_{1}}^{i_{2}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{4}}}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-F_{i_{1}}^{i_{2}}}&0&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{2}}^{j_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{2}}^{j_{4}}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{2}}^{j_{3}}}&0&F_{j_{3}}^{j_{4}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{4}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{2}}^{j_{4}}}&-F_{j_{3}}^{j_{4}}&0\end{array}\right]=\mathbf{X}+\mathbf{B}^{(2)}\mathbf{C}^{(2)}\mathbf{B^{(2)}}^{T}$
(46) $\displaystyle\left\langle{x^{(2)}}\middle|{\phi}\right\rangle$
$\displaystyle=\mathrm{Pf}\,\mathbf{X}^{(2)}=\mathrm{Pf}\,\mathbf{X}\times\left.\mathrm{Pf}\left(\mathbf{C}^{(2)}+\mathbf{B^{(2)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(2)}\right)\right/\mathrm{Pf}\,\mathbf{C}^{(2)}$
(47) $\displaystyle\mathbf{X^{(2)}}^{-1}$
$\displaystyle=\mathbf{X}^{-1}-\mathbf{X}^{-1}\mathbf{B}^{(2)}\left(\mathbf{C}^{(2)}+\mathbf{B^{(2)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(2)}\right)^{-1}\mathbf{B^{(2)}}^{T}\mathbf{X}^{-1}$
(48) $\displaystyle\downarrow$

$\displaystyle\left|{x^{(3)}}\right\rangle$
$\displaystyle=\left|{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{3}}{j_{4}}}\right\rangle$
(49) $\displaystyle\mathbf{X}^{(3)}$
$\displaystyle=\left[\begin{array}[]{cccc}0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}F_{i_{1}}^{i_{2}}}&{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}F_{i_{1}}^{i_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{1}}^{j_{4}}}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-F_{i_{1}}^{i_{2}}}&0&{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}F_{i_{2}}^{i_{3}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{2}}^{j_{4}}}\\\
{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}-F_{i_{1}}^{i_{3}}}&{\color[rgb]{0,0.44,0.57}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.44,0.57}\pgfsys@color@cmyk@stroke{0.98}{0.13}{0}{0.43}\pgfsys@color@cmyk@fill{0.98}{0.13}{0}{0.43}-F_{i_{2}}^{i_{3}}}&0&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}F_{i_{3}}^{j_{4}}}\\\
{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{1}}^{j_{4}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{2}}^{j_{4}}}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-F_{i_{3}}^{j_{4}}}&0\end{array}\right]=\mathbf{X}+\mathbf{B}^{(3)}\mathbf{C}^{(3)}\mathbf{B^{(3)}}^{T}$
(54) $\displaystyle\left\langle{x^{(3)}}\middle|{\phi}\right\rangle$
$\displaystyle=\mathrm{Pf}\,\mathbf{X}^{(3)}=\mathrm{Pf}\,\mathbf{X}\times\left.\mathrm{Pf}\left(\mathbf{C}^{(3)}+\mathbf{B^{(3)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(3)}\right)\right/\mathrm{Pf}\,\mathbf{C}^{(3)}$
(55) $\displaystyle\mathbf{X^{(3)}}^{-1}$
$\displaystyle=\mathbf{X}^{-1}-\mathbf{X}^{-1}\mathbf{B}^{(3)}\left(\mathbf{C}^{(3)}+\mathbf{B^{(3)}}^{T}\mathbf{X}^{-1}\mathbf{B}^{(3)}\right)^{-1}\mathbf{B^{(2)}}^{T}\mathbf{X}^{-1}$
(56) $\displaystyle\vdots$

Matrices $\mathbf{B}^{(r)}$ and $\mathbf{C}^{(r)}$ appearing in eqs. 39, 47,
55, 40, 48 and 56 also have a common expression for their entries:

$\displaystyle B^{(r)}_{lm}$
$\displaystyle=\left\\{\begin{array}[]{ll}F_{i_{r}}^{j_{l}}-F_{j_{r}}^{j_{l}}&l>m,m\leq
r\\\ F_{i_{r}}^{i_{l}}-F_{j_{r}}^{i_{l}}&l<m,m\leq r\\\ 1&l=m,m>r\\\
0&\textrm{otherwise},\end{array}\right.$ (61)
$\displaystyle\mathbf{C}^{(r)}_{lm}$
$\displaystyle=\begin{bmatrix}\mathbf{0}&-\mathbf{I}_{r}\\\
\mathbf{I}_{r}&\mathbf{0}\end{bmatrix}.$ (62)

Note also that we have assumed in these equations that electrons in
$\left|{x}\right\rangle$ change their positions by order from the first one.
In real implementation of mVMC where the hopping electron is selected
randomly, one needs to modify $\mathbf{B}$ accordingly to select a different
set of unit vectors corresponding to that random selection.

Table 1: Different ways for evaluating Markov Chain on $\left|{x}\right\rangle$ with $r_{\mathrm{max}}=3$ Configuration $\left|{x^{(k)}}\right\rangle$ |  . $j_{3}$$j_{2}$$j_{1}$$j_{0}$ |  . $j_{3}$$j_{2}$$j_{1}$$i_{0}$$\rightarrow$ |  . $j_{3}$$j_{2}$$i_{1}$$i_{0}$$\rightarrow$ |  . $j_{3}$$i_{2}$$i_{1}$$i_{0}$$\rightarrow$
---|---|---|---|---
Old method evaluates: | – |  eqs. 39 and 40 |  eqs. 39 and 40 |  eqs. 39 and 40
Old method modifies: | – | $\mathrm{Pf}\,\mathbf{X},\mathbf{X}^{-1}$ | $\mathrm{Pf}\,\mathbf{X},\mathbf{X}^{-1}$ | $\mathrm{Pf}\,\mathbf{X},\mathbf{X}^{-1}$
Old method calls: | SKPFA/SKTDI | SKMV, SKR2 | SKMV, SKR2 | SKMV, SKR2
New method evaluates: | – | eq. 39 | eq. 47 |  eqs. 55 and 56
New method modifies: | – | $\mathrm{Pf}\,\mathbf{X}$ | $\mathrm{Pf}\,\mathbf{X}$ | $\mathrm{Pf}\,\mathbf{X},\mathbf{X}^{-1}$
New method calls: | SKPFA/SKTDI | SKMV | SKMV | SKMM, GEMMT
before: |  after:
---|---
Start w/ $\left|{x^{(0)}}\right\rangle,k\leftarrow 0$Compute $\mathrm{Pf}\,\mathbf{X}$ and $\mathbf{X}^{-1}$Randomly select a new configuration $\left|{x^{\prime}}\right\rangle$For $\left|{x^{\prime}}\right\rangle$ compute $\mathrm{Pf}\,\mathbf{X}^{\prime}$ w/ rank-1 fast-updateCompute weight $\left\langle{x^{\prime}}\middle|{\Psi}\right\rangle$ from $\mathrm{Pf}\,\mathbf{X}^{\prime}$ yieldedAccept?$k\leftarrow k+1$$\left|{x^{(k)}}\right\rangle=\left|{x^{\prime}}\right\rangle$Rank-1 update $\left(\mathbf{X}^{(k)}\right)^{-1}$$k\geq N_{\mathrm{tot}}$?Stop$k\leftarrow k+1$$\left|{x^{(k)}}\right\rangle=\left|{x}\right\rangle$acc.rej.yesno | Start w/ $\left|{x^{(0)}}\right\rangle,k\leftarrow 0$Compute $\mathrm{Pf}\,\mathbf{X}$ and $\mathbf{X}^{-1}$$r\leftarrow 1$Randomly select a new configuration $\left|{x^{\prime}}\right\rangle$For $\left|{x^{\prime}}\right\rangle$ compute $\mathrm{Pf}\,\mathbf{X}^{\prime}$ w/ rank-$r$ fast-updateCompute weight $\left\langle{x^{\prime}}\middle|{\Psi}\right\rangle$ from $\mathrm{Pf}\,\mathbf{X}^{\prime}$ yieldedAccept?$k\leftarrow k+1,\;r\leftarrow r+1$$\left|{x^{(k)}}\right\rangle=\left|{x^{\prime}}\right\rangle$$k\geq N_{\mathrm{tot}}$?StopRank-$r$ update $\left(\mathbf{X}^{(k)}\right)^{-1}$$r<r_{\mathrm{max}}$?$k\leftarrow k+1$$\left|{x^{(k)}}\right\rangle=\left|{x}\right\rangle$acc.rej.yesnoyesno
Figure 1: Flowcharts for evaluating Markov chain on $\left|{x}\right\rangle$
before and after block-update optimization.

### Application

With eqs. 39, 47, 55, 40, 48 and 56, mVMC is able to proceed throughout the MC
without recalculating $\mathrm{Pf}\,\mathbf{X}$ upon each sample. However, in
the current version of mVMC, only rank-1 update equations (i.e. eqs. 39 and
40) are used which eventually made the program memory-bound. To get around
this problem to make full use of modern CPU architectures, we extend the way
of evaluating MC on $\left|{x^{(k)}}\right\rangle$ employed in mVMC. In our
new approach:

  * 1.

Rank-$r|_{r\in\mathbb{Z}^{+}}$ (blocked) updates are used instead of just
rank-1 updates;

  * 2.

A tuning parameter $r_{\mathrm{max}}$ is set to determine at which step should
$\mathbf{X}^{-1}$ be updated. Note that rank-$r$ update does _not_ require
$\left(\mathbf{X}^{(l)}\right)^{-1}$ from previous steps.

Table 1 gives an example of how this change looks like when operating on a
4-electron configuration with $r_{\mathrm{max}}=3$ and Fig. 1 compares the
flowcharts before and after our modification. It is important to note that our
optimization does not alter sampling process or acceptance ratio of the
original Markov chain. What we have done is to “pack up” the heaviest
computation, i.e. update of $\mathbf{X}^{-1}$, grouping each
$r_{\mathrm{max}}$ rank-1 updates into one rank-$r_{\mathrm{max}}$ update so
that the new program becomes much faster as we will see in Sec. 6.

##  5 Library Interface and Accelerated Implementation

This work is basically an optimization of the open-source program mVMC. We did
not make any changes on the user side interface defined in [1]. Instead, this
section gives a brief specification of library routines we have made to
implement procedures presented in Secs. 3 and 4.

PfaffineBLISskpfasktdfsktdiskr2kskmmPfUpdatespush()pop()merge()LAPACK Lv.BLAS
Lv.Skew-symmetricAssembly: A64FXAssembly: Skylake-XAssembly: Zen, etc. Figure
2: Two-layer structured solution for computing Pfaffian and inverse of skew-
symmetric matrix $\mathbf{X}$.

### Layered Implementation

As was stated in Sec. 2, Wimmer’s algorithm is good but their implementation
Pfapack77 is not tuned to use assembly-level optimized BLAS routines as skew-
symmetric operation is not part of the standard. To complete a high-
performance implementation for evaluating equations appeared in Secs. 3 and 4,
an extended BLAS-level implementation is needed. BLIS[30] is a fast and highly
portable framework that implements BLAS-like operations with a fresh approach
based on the GotoBLAS[31, 32] algorithm. It is able to instantiate multiple
level-3 BLAS operations with a minimal number of assembly-based kernels[33]
tuning based on an analytical model[34]. By adding a minimum amount of
interface code to BLIS, skew-symmetric BLAS operations are imported with
optimized kernels for AMD Zen, Intel Skylake-X and many other architectures,
among them a set of GEMM kernels designed for A64FX222 A64FX is the processor
name of current world top supercomputer “Fugaku”.  made by ourselves in
another joint work with Forschungszentrum Jülich. Integrating these modules we
have now built a layered approach shown in Fig. 2 for computing and block-
updating Pfaffian which is fast, well-organized and runs across different CPU
architectures.

Table 2: List of API handles of the Pfaffine library for Pfaffian and inverse computation. Here $T$ refers to different floating-point data types and $t$ their LAPACK-style prefix characters. We use “sk. matrix” as the shorthand of a skew-symmetric matrix. C++ | C (Simp.) | C/Fort. | Description
---|---|---|---
sktdf<$T$> | N/A | N/A | Tridiagonal factorize $\mathbf{X}$ (eq. 10).
sktdi<$T$> | N/A | N/A | Inverse matrix from Tridiagonal factorization (eq. 15).
skpfa<$T$> |  $t$skpfa_ | cal_$t$skpfa_ | Tridiagonal factorize and compute $\mathrm{Pf}\,\mathbf{X}$ (eq. 13),
|  |  |  $\mathbf{X}^{-1}$ computed conditionally depending on input.
skslc<$T$> | N/A | N/A | Get a column from corner-stored sk. matrix.
skswp<$T$> | N/A | N/A | Swap columns and rows of cornet-stored sk. matrix.
skr2k<$T$> | N/A | N/A |  $\mathbf{X}+\mathbf{AB}^{T}-\mathbf{BA}^{T}$ where $\mathbf{X}$ is corner-stored sk. matrix.
skmm<$T$> | N/A | N/A |  $\mathrm{XB}$ where $\mathbf{X}$ is corner-stored sk. matrix.
Table 3: Method API of C++ object updated_tdi for rank-$r$ update strategy. C++ Method | Description
---|---
updated_tdi<$T$>::updated_tdi() | Constructor. Given $f_{ij}$ and $\left|{x}\right\rangle$, compute $\mathrm{Pf}\,\mathbf{X}$ and $\mathbf{X}^{-1}$.
updated_tdi<$T$>::get_Pfa() | Get Pfaffian of current configuration.
updated_tdi<$T$>::push_update() | Push a configuration update: $\left|{x^{(r)}}\right\rangle\rightarrow\left|{x^{(r+1)}}\right\rangle$
| with one electron changed its spin-orbital position.
updated_tdi<$T$>::pop_update() | Revert the last push_update().
updated_tdi<$T$>::merge_updates() | Modify $\mathbf{X}^{-1}$ based on the last accepted configuration.

### Interface

C++ and C interfaces of our Pfaffian-inverse framework described in Fig. 2 are
listed in Tab. 2, while object-based routines for Sec. 4’s blocked update are
listed in Tab. 3. We recommend referring to our program distribution in [35]
for more API-related information.

##  6 Benchmarking Results

$5$$10$$15$$20$$25$$30$$0$$2$$4$$6$GFLOPS/sec @Intel Skylake-XReal Slater
Determinants$5$$10$$15$$20$$25$$30$$0$$5$$10$Complex Slater
Determinants$5$$10$$15$$20$$25$$30$$0$$2$$4$GFLOPS/sec @AMD
Zen2$5$$10$$15$$20$$25$$30$$0$$2$$4$$6$$8$Original mVMCBlocked
mVMC$5$$10$$15$$20$$25$$30$$0$$1$$2$Square System Length $L$
($N_{\text{Spin}}=L^{2}$)GFLOPS/sec @Fujitsu
A64FX$5$$10$$15$$20$$25$$30$$0$$2$$4$$6$ Figure 3:  Benchmark conducted on 3
different CPUs against mVMC’s Markov chain sampling process comparing original
version and our block-update optimized version. FLOPs numbers here, given by
eq. 63, show how much our method has speed up the mVMC program. Performance
draw-back of Intel Skylake-X at $L\geq 20$ is believed to be caused by the
block-update’s panel size exceeding processor’s L2 cache[36]. Future work on
deploying a hierarchy storage (e.g. one implemented in Low and van de Geijn
[37]) to $\mathbf{X}$ might address this problem and allow further speed-up.

We run a Markov chain sampling a random many-body spin wavefunction on a
$L\times L$ square lattice. Our experiment only studies the sampling process
itself against this random wavefunction so that the Hamiltonian is omitted.
The total number of samples generated by the MC is proportional to the number
of spins:

$N_{\mathrm{tot}}\propto N_{\mathrm{spin}}=L\times L,$

because MC needs more steps to decorrelate as system size grows. This
experimental setup yields:

$r_{\mathrm{acc}}N_{\mathrm{proj}}\left(N_{\mathrm{qty}}+1\right)\overbrace{N_{\mathrm{smp}}N_{\mathrm{spin}}}^{N_{\mathrm{tot}}}\left(2N_{\mathrm{spin}}^{2}+\frac{3}{4}\frac{1}{N_{\mathrm{spin}}}N_{\mathrm{spin}}^{3}\right)\equiv\omega
N_{\mathrm{spin}}^{3}\text{ MACs}$

for each system size. Here $r_{\mathrm{acc}}$ is the acceptance ratio since
rank-1 (rank-$k$) update is requested every time (every $k$ times) a
configuration update is _accepted_. $N_{\mathrm{qty}}$, $N_{\mathrm{smp}}$ and
$N_{\mathrm{proj}}$ correspond to NDataQtySmp, NVMCSample and NQ parameters in
the mVMC application. The extra
$\frac{3}{4}\frac{1}{N_{\mathrm{spin}}}N_{\mathrm{spin}}^{3}$ term inside the
bracket accounts for the recalculation of full
$\mathbf{X}\rightarrow\mathbf{X}^{-1}$ every $N_{\mathrm{spin}}$ acceptances
to remove the accumulated round-up error by fast updates. We set
$N_{\mathrm{qty}}=10$, $N_{\mathrm{smp}}=10$ and $N_{\mathrm{proj}}=16$ in our
experiment. Acceptance ratios turned out to be
$r_{\mathrm{acc}}^{\mathrm{real}}\approx 31\%$ for real Slater determinants
and $r_{\mathrm{acc}}^{\mathrm{comp}}\approx 47\%$ for complex ones regardless
of $N_{\mathrm{spin}}$. Hence the coefficients $\omega$ become
$\omega_{\text{real}}\approx 2050$ and $\omega_{\text{complex}}\approx 3000$
for real and complex Slater determinants, respectively. Since different
hardware vendors sometimes count performance data differently, here we use the
expression:

$\text{FLOPs/sec}=\left.\omega_{\text{data}}\times\text{OpCost}_{\text{data}}\times
L^{6}\middle/T_{\text{elapsed}}\right.$ (63)

to convert time measurements to floating-point performance data, where
$\text{OpCost}_{\text{data}}$ is the number of floating-point operations
(FLOPs) yielded by each MAC on corresponding data types, hence $2$ for real
Slater determinants and $8$ for complex Slater determinants.

Three CPU variants are selected to be included in this test, each having its
own significance in the HPC context:

  * 1.

Intel® Xeon® Platinum 8260 under Skylake-X microarchitecture with AVX512
vector instructions support;

  * 2.

AMD EPYC™ 7702 under Zen 2 microarchitecture with AVX2 vector instructions
support;

  * 3.

Fujitsu A64FX based on Armv8 ISA with 512-bit Scalable Vector Extension[38]
support.

Results of our blocked-update mVMC compared to the original implementation is
shown in Fig. 3. Note that all tests here are conducted single-threaded
because mVMC has parallel capability at a higher level. These results show
that performance of our blocked algorithm can ramp up quickly with the growth
of system size, yielding up to more than $6\times$ performance compared to the
original code 333The original code is compiled with corresponding vendor
compilers under default max-optimization options. on all our processors
tested, while the latter two tested architectures seem to bear potential for
higher speed-up rates on even larger system sizes. Among all test cases here,
AMD EPYC™ 7702 yields the most significant speed-up (up to $\sim 17.5\times$)
thanks to a new skinny-matrix handling strategy introduced to BLIS[39].

##  7 Conclusion

We have shown in this article a framework for computing and updating Pfaffians
that is fast, flexible and runs across different microarchitectures. mVMC
optimized with this framework has achieved significant speed up compared to
its original release. Theory as well as library implementation in this work
can also be applied to other kinds of variational methods that deals with
skew-symmetric wavefunctions, i.e. Fermions, like TurboRVB[24]. On the
implementation side, this blocking would also allow potential portings of the
related methods onto heterogeneous platforms like a general-purpose GPU
machine, although the sequential nature of Markov chains might require
additional parallelizations to unleash GPUs’ full power. Further, our work
shows a new possible pattern of open-source communities collaborating with
hardware people to deliver an assembly-level optimized code for HPC
applications.

## Acklowdgements

A64FX assembly kernels for BLIS used in this work are developed in a joint
work with Stepan Nassyr from Forschungszentrum Jülich and the Science of High-
Performance Computing (SHPC) group from the University of Texas at Austin.
A64FX test results were collected under “Program for Promoting Researches on
the Supercomputer Fugaku” by MEXT with the project title “Basic Science for
Emergence and Functionality in Quantum Matter: Innovative Strongly-Correlated
Electron Science by Integration of “Fugaku” and Frontier Experiments”
(hp200132). Our development work used also the Isambard 2 UK National Tier-2
HPC Service (http://gw4.ac.uk/isambard/) operated by GW4 and the UK Met
Office, and funded by EPSRC (EP/T022078/1). RGX acknowledge the Global Science
Graduate Course (GSGC) program of the University of Tokyo. TO and ST
acknowledge support by the Endowed Project for Quantum Software Research and
Education, the University of Tokyo (https://qsw.phys.s.u-tokyo.ac.jp/).

## References

  * Misawa  _et al._ [2019] T. Misawa, S. Morita, K. Yoshimi, M. Kawamura, Y. Motoyama, K. Ido, T. Ohgoe, M. Imada, and T. Kato, Comp. Phys. Comm. 235, 447 (2019).
  * White [1993] S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Sorella  _et al._ [2002] S. Sorella, G. B. Martins, F. Becca, C. Gazza, L. Capriotti, A. Parola, and E. Dagotto, Phys. Rev. Lett. 88, 117002 (2002).
  * Clark  _et al._ [2011] B. K. Clark, D. A. Abanin, and S. L. Sondhi, Phys. Rev. Lett. 107, 087204 (2011).
  * Iqbal  _et al._ [2013] Y. Iqbal, F. Becca, S. Sorella, and D. Poilblanc, Phys. Rev. B 87, 060405 (2013).
  * Misawa and Imada [2014] T. Misawa and M. Imada, Nat. Commun. 5, 1 (2014).
  * Hu  _et al._ [2015] W.-J. Hu, W. Zhu, Y. Zhang, S. Gong, F. Becca, and D. N. Sheng, Phys. Rev. B 91, 041124 (2015).
  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).
  * Ohgoe  _et al._ [2020] T. Ohgoe, M. Hirayama, T. Misawa, K. Ido, Y. Yamaji, and M. Imada, Phys. Rev. B 101, 045124 (2020).
  * Kawamura  _et al._ [2017] M. Kawamura, K. Yoshimi, T. Misawa, Y. Yamaji, S. Todo, and N. Kawashima, Comp. Phys. Comm. 217, 180 (2017).
  * Schollwöck [2011] U. Schollwöck, Ann. Phys. (NY) 326, 96 (2011), january 2011 Special Issue.
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, arXiv preprint cond-mat/0407066 (2004).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 143 (2008).
  * Maeshima  _et al._ [2001] N. Maeshima, Y. Hieida, Y. Akutsu, T. Nishino, and K. Okunishi, Phys. Rev. E 64, 016705 (2001).
  * Corboz [2016] P. Corboz, Phys. Rev. B 94, 035133 (2016).
  * Vanderstraeten  _et al._ [2016] L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Phys. Rev. B 94, 155123 (2016).
  * Liao  _et al._ [2019] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Phys. Rev. X 9, 031041 (2019).
  * Tahara and Imada [2008] D. Tahara and M. Imada, J. Phys. Soc. Japan 77, 114701 (2008).
  * Morita  _et al._ [2015] S. Morita, R. Kaneko, and M. Imada, J. Phys. Soc. Japan 84, 024720 (2015).
  * Ido  _et al._ [2018] K. Ido, T. Ohgoe, and M. Imada, Phys. Rev. B 97, 045138 (2018).
  * Nomura  _et al._ [2017] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Phys. Rev. B 96, 205152 (2017).
  * Nomura and Imada [2020] Y. Nomura and M. Imada, arXiv preprint arXiv:2005.14142 (2020).
  * Needs  _et al._ [2000] R. J. Needs, G. Rajagopal, M. D. Towler, P. R. C. Kent, and A. J. Williamson,  _CASINO version 1.0 User’s Manual_ (University of Cambridge, Cambridge, UK, 2000).
  * Nakano  _et al._ [2020] K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti, M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo, G. Mazzola, _et al._ , J. Chem. Phys. 152, 204121 (2020).
  * Ido  _et al._ [2020] K. Ido, M. Imada, and T. Misawa, Phys. Rev. B 101, 075124 (2020).
  * Charlebois and Imada [2020] M. Charlebois and M. Imada, Physi. Rev. X 10, 041023 (2020).
  * Wimmer [2012] M. Wimmer, ACM Trans. Math. Software 38, 1 (2012).
  * Parlett and Reid [1970] B. Parlett and J. Reid, BIT 10, 386 (1970).
  * Woodbury [1950] M. A. Woodbury, Statistical Research Group, Princeton University (1950).
  * Van Zee and van de Geijn [2015] F. G. Van Zee and R. A. van de Geijn, ACM Trans. Math. Software 41, 14:1 (2015).
  * Goto and van de Geijn [008b] K. Goto and R. van de Geijn, ACM Trans. Math. Software 35, 1 (2008b).
  * Goto and van de Geijn [008a] K. Goto and R. van de Geijn, ACM Trans. Math. Software 34 (2008a), 10.1145/1356052.1356053.
  * Van Zee  _et al._ [2016] F. G. Van Zee, T. Smith, F. D. Igual, M. Smelyanskiy, X. Zhang, M. Kistler, V. Austel, J. Gunnels, T. M. Low, B. Marker, L. Killough, and R. A. van de Geijn, ACM Trans. Math. Software 42, 12:1 (2016).
  * Low  _et al._ [2016] T. M. Low, F. D. Igual, T. M. Smith, and E. S. Quintana-Ortí, ACM Trans. Math. Software 43, 12:1 (2016).
  * Xu [2021] R. G. Xu, “Release v1.1.1 + blocked pfaffian update optimization | issp-center-dev/mVMC,”  (2021).
  * Goto and van de Geijn [2002] K. Goto and R. A. van de Geijn,  _On Reducing TLB Misses in Matrix Multiplication_ , Tech. Rep. CS-TR-02-?? (Department of Computer Sciences, The University of Texas at Austin, 2002).
  * Low and van de Geijn [2004] T. M. Low and R. van de Geijn,  _An API for Manipulating Matrices Stored by Blocks_ , Tech. Rep. TR-2004-15 (Department of Computer Sciences, The University of Texas at Austin, 2004).
  * Stephens  _et al._ [2017] N. Stephens, S. Biles, M. Boettcher, J. Eapen, M. Eyole, G. Gabrielli, M. Horsnell, G. Magklis, A. Martinez, N. Premillieu, A. Reid, A. Rico, and P. Walker, IEEE Micro 37, 26 (2017).
  * Van Zee and van de Geijn [2019] F. G. Van Zee and R. A. van de Geijn, “Level-3 GEMM performance on small matries | BLAS-like library instantiation software framework,”  (2019).