# Efficient algorithms for collecting the statistics of large-scale IP address
data

Jie Chen  Hua Mao<EMAIL_ADDRESS>College of Computer Science,
Sichuan University, Chengdu 610065, P.R. China  Department of Computer and
Information Sciences, Northumbria University, Newcastle, NE1 8ST, U.K.

###### Abstract

Compiling the statistics of large-scale IP address data is an essential task
in network traffic measurement. The statistical results are used to evaluate
the potential impact of user behaviors on network traffic. This requires
algorithms that are capable of storing and retrieving a high volume of IP
addresses within time and memory constraints. In this paper, we present two
efficient algorithms for collecting the statistics of large-scale IP addresses
that balance time efficiency and memory consumption. The proposed solutions
take into account the sparse nature of the statistics of IP addresses while
building the hash function and maintain a dynamic balance among layered memory
blocks. There are two layers in the first proposed method, each of which
contains a limited number of memory blocks. Each memory block contains 256
elements of size $256\times 8$ bytes for a 64-bit system. In contrast to
built-in hash mapping functions, the proposed solution completely avoids
expensive hash collisions while retaining the linear time complexity of hash-
based solutions. Moreover, the mechanism dynamically determines the hash index
length according to the range of IP addresses, and can balance the time and
memory constraints. In addition, we propose an efficient parallel scheme to
speed up the collection of statistics. The experimental results on several
synthetic datasets show that the proposed method substantially outperforms the
baselines with respect to time and memory space efficiency.

###### keywords:

large-scale IP addresses , memory blocks , hash table , sorting , network
traffic

††journal: Information Processing Letterst1t1footnotetext: This work was
partially supported by National Natural Science Foundation of China (NSFC)
under Grants 61303015 and Sichuan Science and Technology Program under Grant
2021YJ0078

##  1 Introduction

In recent years, the amount of network traffic has increased significantly
because of the rapid development of emerging network services, such as video
streaming, instant messaging, and online payment services. It is crucial to
evaluate the potential impact of the features of user behavior on network
traffic management. User behavior features are usually extracted from IP
packets, which contain IP addresses. There is a close correspondence between
informative user behavior features and the IP addresses that frequently appear
in the packets. Hence, how to effectively obtain the statistics of large-scale
IP address data in a timely manner, such as every few minutes, is a
challenging problem in network traffic measurement. Obtaining the statistics
of large-scale IP address data typically consists of two tasks: counting the
number of occurrences of each IP address and sorting the results in a specific
order.

A number of statistical algorithms for solving this problem have been studied
in the last few decades [1, 2, 3, 4, 5]. A classic divide-and-conquer strategy
has been proposed in which IP addresses are first divided into multiple
subsets. Then, each subset is individually computed using a statistics
collection method, such as a top $k$ trie algorithm [6]. Finally, to merge and
sort the results of the multiple subsets, sorting algorithms (e.g., bubble
sort, insertion sort, merge sort, selection sort, or quick sort) are used [7,
8]. However, the reduction of computational cost is an intractable problem in
the sorting procedure [9, 10, 11]. For example, the average complexity and
worst-case complexity of the bubble, insertion, and selection sort algorithms
are $O\left({{n^{2}}}\right)$ [8, 12], where $n$ represents the number of
unsorted records. This indicates that the merging and sorting step of the
multiple subsets occupies a large amount of memory and is extremely
computationally costly to perform when collecting the statistics of millions
or tens of millions of IP addresses. Therefore, statistics collection
algorithms using the divide-and-conquer strategy still face several challenges
caused by the rapid increase of the large-scale records, such as bounded
memory and computational cost restrictions.

A hash table is an effective method for collecting the statistics of IP
addresses [13]. It uses a hash function to compute a hash codes for an array
of buckets with the statistical results. The hash function assigns each key to
a unique bucket for each IP address. Unfortunately, the hash function can
generate the same hash code for more than one IP address. With the increase in
the generation of big data, millions or tens of millions of records have
become ubiquitous in network traffic. Therefore, this approach could cause
several hash collisions, especially for a large number of IP addresses.
Although many strategies can be employed to avoid collisions, such as linear
probing, quadratic probing, and double hashing, they require extra storage
space and computation.

The statistics collection algorithm should be stable, effective, and efficient
for large-scale records. To overcome the disadvantages of general statistics
collection methods, a number of parallel techniques have been developed for
large-scale records by optimizing the efficiency and complexity. For example,
these algorithms have been extended to the corresponding parallel structures
for parallel hardware architectures, such as many-core and multi-core
platforms [14, 15]. However, these parallel statistics collection algorithms
are still the parallel computing variants of the original algorithms, and they
continue to ignore the importance of the original characteristics of the
unsorted records.

In this paper, we present two efficient algorithms for collecting the
statistics of large-scale IP address data. We can obtain the frequently
occurring IP addresses from the statistics, which can be regarded as a pre-
processing step of user behavior analysis in network traffic management.
Because of the increasing volume and speed of network traffic, it has become
expensive and impractical to handle all IP addresses contained the IP packets.
By taking full advantage of the successive characteristics of memory addresses
and the fixed range of each individual part of an IP address, we design two
relationship mapping mechanisms between memory blocks and IP addresses for a
four-dimensional sparse matrix. The sparse matrix stores the number of
occurrences of the individual IP addresses, in which the positions of the rows
and columns are employed to represent the mapping relationship between the
memory blocks and IP addresses. Specifically, we construct a two-layer memory
block (TLMB) to implement the first mapping mechanism for the IP addresses. In
addition, we employ a single shared memory block (SSMB) for all IP addresses
to implement the other mapping mechanism of the IP addresses. The mechanisms
of the mapping relationship effectively remove the information about trivial
user behaviors that are irrelevant for statistical analysis. The proposed
methods can be extended to the corresponding parallel versions for specific
hardware architectures. Extensive experiments on several synthetic datasets
show the effectiveness of the proposed method.

Our contributions are summarized as follows.

  1. 1.

We present two efficient methods for collecting the statistics of large-scale
IP address data that use two different relationship mapping mechanisms, TLMB
and SSMB, between memory blocks and IP addresses for a four-dimensional sparse
matrix.

  2. 2.

The computational cost of TLMB is linearly proportional to the number of IP
addresses with a limited memory resource, and the memory use of SSMB remains
almost unchanged with a reasonable computational cost, regardless of the
number of IP addresses.

  3. 3.

A parallel computation optimization scheme for multiple computers is proposed
to effectively improve the computational efficiency and dramatically reduce
memory use.

  4. 4.

Our extensive experimental results using synthetic datasets demonstrate that
our proposed method shows clear superior performance in comparison with the
baselines, striking a balance between computational cost and memory use.

The remainder of this paper is organized as follows. We briefly describe some
of the original sorting techniques as well as various parallel sorting
algorithms in Section 2. In Section 3, we introduce the proposed method. The
experimental results are presented in Section 4. Finally, we draw the
conclusions of the study in Section 5.

##  2 Background

###  2.1 Classical sort techniques

Bubble sort is a classical sorting algorithm in which each element in a list
is compared with its neighboring elements and swapped until they are in the
desired order [16]. Bubble sort leads to $\left({n-1}\right)$ number of passes
and $\frac{{n\left({n-1}\right)}}{2}$ number of iterations if $n$ elements are
given. Insertion sort is a simple and efficient sorting algorithm that
iteratively takes one element and finds its appropriate position in the sorted
list by comparing it with neighboring elements. It becomes less efficient as
the number of records increases. Selection sort determined the smallest number
in an unsorted list and swaps it with the first number in the sorted list.
Then, it finds the next smallest element from the remaining list and swaps
with the second element in the sorted list. Consequently, the number of sorted
elements at the top of the list increases while the rest remain unsorted.
Merge sort, which is based on the divide-and-conquer principle, repeatedly
divides the array into two halves and then combines them in a sorted manner
For more details of classical sort algorithms, such as quick sort and heap
sort, we refer the reader to the comments in [17, 7]. It is well known that
computational cost and required memory are the primary concerns in sorting
algorithms [18, 19]. Table 1 shows the time complexity of the best case,
average case, and worst case of several classical sorting algorithms [7, 20,
21].

Table 1: Time complexity of various sorting algorithms. Algorithm | $Best\left(n\right)$ | $Average\left(n\right)$ | $Worst\left(n\right)$
---|---|---|---
Bubble sort | $O\left({n}\right)$ | $O\left({n^{2}}\right)$ | $O\left({n^{2}}\right)$
Insertion sort | $O\left({n}\right)$ | $O\left({n^{2}}\right)$ | $O\left({n^{2}}\right)$
Selection sort | $O\left({n^{2}}\right)$ | $O\left({n^{2}}\right)$ | $O\left({n^{2}}\right)$
Merge Sort | $O\left({n\log n}\right)$ | $O\left({n\log n}\right)$ | $O\left({n\log n}\right)$
Quick Sort | $O\left({n\log n}\right)$ | $O\left({n\log n}\right)$ | $O\left({n^{2}}\right)$
Heap Sort | $O\left({n\log n}\right)$ | $O\left({n\log n}\right)$ | $O\left({n\log n}\right)$

###  2.2 Parallel hardware architectures

The hardware architecture of modern processors usually consists of more than
two independent central processing units (CPUs) or graphics processing units
(GPUs). Parallel software platforms can be implemented using high-level
programming frameworks for specific hardware architectures [22]. The Compute
Unified Device Architecture (CUDA) is a parallel computing platform for
general computing on GPUs. Most parallel sorting algorithms are variants of
standard, well-known sorting algorithms adapted to GPU hardware architecture.
For example, Cederman designed a quick sort for the GPU platform [23], and
Peters proposed an adaptive bitonic sorting algorithm with a bitonic tree for
GPUs [24]. The parallel computation of sorting algorithms is considered to be
the most efficient way of sorting elements on parallel hardware architectures
[15].

Figure 1: Example of a memory block of size $256\times 8$ bytes containing 256
elements.

##  3 Proposed method

###  3.1 Problem Formulation

Assume that IP flow data $FD$ is a sequence of IP records, that is,
$FD=\left\\{\left({p_{1},x_{1}}\right),\left({p_{2},x_{2}}\right),...,\left({p_{n},x_{n}}\right)\right\\}\in{\mathbb{R}^{m\times
n}}$ and ${n\geq 1{e^{6}}}$, where each pair of elements
$\left({p_{i},x_{i}}\right)$ $\left({i\in\left[{1,n}\right]}\right)$ consists
of an IP address $p_{i}$ and a set of corresponding user behavior attributes
$x_{i}$. Given a finite set of IP addresses
$P=\left\\{{{p_{1}},{p_{2}},...,{p_{n}}}\right\\}\in{\mathbb{R}^{n}}$, the
purpose of the IP address statistics task is to efficiently determine the
first $k$ IP addresses of the most frequent occurrences in $P$, where $k\ll
n$.

###  3.2 IP Address Statistics Task

A standard IP address is composed of four decimal numbers ranging from 0 to
255 which are separated by dot symbols. An individual IP address is logically
divided into four parts by splitting it with respect to each dot symbol, and
each part of an IP address has an integer value. We create a four-dimensional
array for the statistics of IP addresses, where the length of each dimension
in the array is 256. Each element of the array can be employed to store the
number of occurrences of the IP address according to the relationship mapping
between the index of each dimension of the array and the integer value of the
corresponding part of the IP address. For example, consider the individual IP
address 1.2.3.4 and the four-dimensional array $fd\\_{array}$. The number of
occurrences of this IP address is stored in ${fd}\\_{array}[1][2][3][4]$.
However, individual IP addresses in the host logs often make up a small
proportion of all IP addresses. The array can be considered to be sparse
because most of its elements are zeros. Consequently, we can carefully design
a four-dimensional sparse matrix to store the number of an individual IP
address by taking full advantage of the successive characteristics of array
addresses and the fixed range of an individual part of an IP address.

Algorithm 1 TLMB

0: A finite set of IP addresses
$X=[{x_{1}},{x_{2}},...,{x_{n}}]\in{\mathbb{R}^{m\times n}}$, number $k>1$.

1: Construct $256\times 256$ memory blocks of size 128 MB for the first layer;

2: for $i=1:n$ do

3: Assume that $a$, $b$, $c$ and $d$ each represent one of the integer values
of the four parts of IP address $x_{i}$.

4: Calculate the index of the memory block in the first layer: $p=a\times
255\times 255+b\times 255$.

5: if the value of the $j$-th element is null in the $p$-th memory block then

6: Create a memory block in the second layer, set all elements of the memory
block to zero, and store the starting address of the memory block in the
$j$-th element.

7: else

8: Obtain the starting address of the memory block $m$ in the second layer by
finding the $j$-th element of the $p$-th memory block.

9: end if

10: Add 1 to the $d$-th element of memory block $m$ in the second layer.

11: end for

12: Construct a minimum heap of size $k$ using each non-zero element of the
memory blocks in the second layer.

12:

13: Traverse the nodes of the heap to obtain the $k$ IP addresses and their
number of occurrences.

###  3.3 IP Usage Storage and Retrieval Strategies for the Four-dimensional
Sparse Matrix

Assume that each memory address of a 64-bit system can be stored in an element
8 bytes in size, and the number of occurrences of an individual IP address is
no more than $2^{64}$. There is an array of size 256 elements that consists of
$256\times 8$ bytes of memory. The array is regarded as a memory block that
contains a contiguous address space, as shown in Fig. 1. In other words, the
addresses of all bytes of the array are sequential in the memory block.
Therefore, the position of the array can be indexed by the integer value of a
particular part of the IP address. We present two efficient methods to collect
the statistics of large-scale IP address data, each of which contains a
relationship mapping mechanism between memory blocks and IP addresses for the
four-dimensional sparse matrix.

Figure 2: An example of the mapping relationships between memory blocks and an
IP address.

####  3.3.1 First method: TLMB

The first proposed mapping mechanism of IP addresses is TLMB. The four parts
of the IP address are represented in four layers, where each layer is made up
of one or more memory blocks. The first layer only contains one memory block,
whereas the second layer contains 256 memory blocks. Each memory block
contains 256 elements. Each element of the memory block in the first layer is
employed to store the starting addresses of the corresponding 256 memory
blocks in the second layer. Similarly, the third layer contains $256\times
256$ memory blocks, the size of which is 128 MB in memory. Then, the element
of each memory block in the third layer stores the starting address of the
corresponding memory block in the fourth layer. This would be 32 GB in size if
we adopted a pre-allocation strategy for all memory blocks in the four layers.
Hence, we present an alternative pre-allocation strategy for the memory
blocks. A memory block will be allocated only when the first three parts of an
initial IP address have been given. In particular, pre-allocating a big memory
block of size 128 MB containing $256\times 256$ contiguous memory blocks is
feasible in a modern computer. Consequently, the first two layers can be
removed from this architecture if the third layer has contiguous memory blocks
of 128 MB.

We formally present a storage strategy for IP addresses that consists of two
layers that consist of a limited number of memory blocks. The first layer
contains $256\times 256$ memory blocks. The first three parts of the IP
address can be mapped into the corresponding position of the element in a
particular memory block of the first layer according to the individual values
of the three parts. We allocate a memory block in the other layer for the IP
address when its first three parts are initially given. Each element of a
memory block in this layer stores the number of occurrences of the
corresponding IP address. Figure 2 shows an example of the relationship
mapping between the memory blocks of two layers and an IP addresses. The
positions of the first three dimensionalities of the sparse matrix can be
mapped to the elements of the memory blocks in the first layer. The positions
of the final dimensionality of the sparse matrix can be indexed by combining
the starting address of the memory block in the second layer with the integer
value of the four parts of the IP address.

We traverse all elements of the memory blocks of the second layer to obtain
the maximum number of occurrences of elements if $k=1$. Otherwise, we
construct a minimum heap of size $k$. The statistical results of the first $k$
IP addresses are saved in the heap, which is a special binary tree and
implemented by an array of size $k$. The construction of the heap is completed
by traversing all elements of the memory blocks of the second layer. The tree
node in the leap contains two important attributes: an IP address and its
number of occurrences. The final IP addresses and numbers of occurrences can
be obtained by a traversal of the nodes of the heap. The complete procedure
for determining the first $k$ of the most frequent IP addresses from host logs
is outlined in Algorithm 1.

####  3.3.2 Second method: SSMB

We also designed an SSMB that stores all IP addresses and their statistics. IP
addresses are logically divided into at most 256 subsets according to the
value of the first part of each individual IP address. We construct a single
memory block of $256\times 256\times 256$ elements that is $256\times
256\times 256\times 8$ bytes in size, that is, 128 MB. The memory block is
shared by all subsets. For each subset, the last three parts of the IP address
can be mapped into the corresponding position of the element in the shared
single memory block, which is always initialized at the beginning of the
relationship mapping. Then, the element of this memory block stores the
statistics of the IP addresses in this subset.

Algorithm 2 SSMB

0: A finite set of IP addresses
$X=[{x_{1}},{x_{2}},...,{x_{n}}]\in{\mathbb{R}^{m\times n}}$, number $k>1$.

1: Construct a memory block of size 128 MB saved by all IP addresses;

2: All IP addresses are logically partitioned into $q$ subsets according to
the first part of each IP address.

3: for $i=1:n$ do

4: Assume that $a$, $b$, $c$ and $d$ each represent one of the integer values
of the four parts of IP address $x_{i}$.

5: for $j=1:q$ do

6: if $j==q$ then

7: Calculate the position of the memory block in the first layer: $p=b\times
255\times 255+c\times 255+d$.

8: Add 1 to the $p$-th element of the memory block.

9: end if

10: end for

11: Construct a minimum heap of size $k$ using the each non-zero elements of
the memory block.

12: end for

12:

13: Traverse of the nodes of the heap to obtain the $k$ IP addresses and their
number of occurrences.

We further perform a round traversal of the memory block to initialize a
minimum heap of size $k$ after the relationship mapping has been completed in
the first subset. Then, we continue to perform a round traversal of the memory
block to adjust the heap after the relationship mapping has been completed for
the subsequent subsets. Finally, we obtain the first $k$ most frequently
occurring IP addresses in the heap. The complete procedure for finding the
first $k$ of the most frequency IP addresses from host logs is outlined in
Algorithm 2.

###  3.4 Memory Use and Complexity Analysis

We first evaluate the memory use and computational complexity of the first
proposed method, TLMB. The size of each memory block is $256\times 8$ bytes,
and there are $256\times 256$ memory blocks in the first layer. Hence, the
size of the memory blocks in the first layer is 128 MB. Assume that the number
of the distinct first three parts of the IP addresses is $s$. The number of
memory blocks in the second layer is linearly proportional to $s$. Moreover,
the memory size of the minimum heap is $\left(k+8\right)$ bytes, where each
tree node contains two attributes, that is, an IP address and the number of
occurrences. The total size of the memory of the proposed method is
approximately the sum of the three parts, that is, 128 MB, $s$ KB, and
$\left(k+8\right)$ bytes. The computational complexity of the two layers for
calculating the IP address statistics is ${O}\left(n\right)$ in Algorithm 1,
where $n$ is the number of IP addresses. In addition, the computational
complexity of constructing a minimum heap of size $k$ is ${O}\left({k\log
k}\right)$, where $k$ is the number of tree nodes in the heap. Consequently,
the overall computational complexity of the proposed algorithm is
${O}\left({k\log k+n}\right)$.

We next evaluate the memory use and computational complexity of the second
proposed method, SSMB. The memory size of SSMB is 128 MB. Assume that the
number of the distinct first parts of the IP addresses is $q$. The
computational complexity of the mapping mechanism of the IP addresses in
Algorithm 1 is ${O}\left(qn\right)$, where $n$ is the number of IP addresses.
Similarly, the computational complexity of constructing a minimum heap of size
$k$ is ${O}\left({k\log k}\right)$ for each subset. Consequently, the overall
computational complexity of the proposed algorithm is ${O}\left({q\left(k\log
k+n\right)}\right)$.

###  3.5 Parallel computation optimization techniques

####  3.5.1 Parallel computation mechanism of TLMB

In the worst case, the distinct first three parts of the IP addresses cover
all the binary combinations. Hence, the size of the memory blocks in the
second layer is 32 GB. We present a parallel computation scheme on multiple
computers for improving the computational efficiency and reducing memory use.
Assume that there are $2^{r}$ computers available for parallel computation,
where $r$ represents a positive integer. The task of collecting IP address
statistics is then divided into multiple subtasks, which are performed by
$2^{r}$ computers, respectively, according to the first $r$ bits of the first
part of the IP addresses. The number of memory blocks reduces to
$\left(256\times 256\right)/2^{r}$ for $2^{r}$ computers. Simultaneously, the
number of memory blocks will decrease to $32/2^{r}$ GB in the second layer.
For example, the number of memory blocks in the first layer is $256\times 64$
for each computer when $r=2$, and the size of memory blocks in the second
layer is 8 GB in the worst case. In addition, if four computers perform the
task of computing IP address statistics by partitioning the first two bits of
the first part of the IP addresses, then the second layers of multiple
computers are merged into a complete second layer, where is employed to
construct a minimum heap of size $k$. Finally, the parallel computation
results can be obtained in a manner similar to the last step of Algorithm 1.

####  3.5.2 Parallel computation mechanism of SSMB

Assume that all IP addresses are logically divided into $q$ subsets according
to the value of the first part of an individual IP address. Further assume are
$q$ computers for parallel computation, where the statistics collection task
of each subset can be performed by an individual computer. A minimum heap of
size $k$ is shared among these computers. Hence, this greatly increases the
computational efficiency of the task by $q$ times.

Table 2: Statistics of the three synthetic datasets. Data sets | IP Records | Individual IP Addresses | Size
---|---|---|---
1 | 5, 000, 000 | 50, 000 | 77.5 MB
2 | 10, 000, 000 | 100, 000 | 155 MB
3 | 50, 000, 000 | 500, 000 | 775 MB
Table 3: Computational cost (s) of different methods on the three synthetic datasets. Data | $k$ | Hash Mapping | IP Mapping | Ours (TLMB) | Ours (SSMB)
---|---|---|---|---|---
1 | 10 | 1458.14 (2.32) |  15.18 (0.07) |  2.51 (0.01) | 18.43 (0.05)
100 | 1458.56 (2.76) |  15.21 (0.04) |  2.52 (0.01) | 18.43 (0.06)
2 | 10 | 2927.23 (11.64) |  17.41 (0.11) |  4.92 (0.01) | 27.59 (0.05)
100 | 2934.65 (28.93) |  17.45 (0.05) |  4.95 (0.01) | 27.65 (0.06)
3 | 10 | 14547.09 (13.95) |  35.33 (0.09) |  24.22 (0.07) | 101.01 (0.17)
100 | 14566.36 (26.78) |  35.39 (0.05) |  24.24 (0.04) | 101.20 (0.2)
Table 4: Memory use (MB) of different methods on the three synthetic datasets. Data | $k$ | Hash Mapping | IP Mapping | Ours (TLMB) | Ours (SSMB)
---|---|---|---|---|---
1 | 10 |  43.63 (0.12) | 16332.81 (0.07) | 189.61 (0.16) |  139.77 (0.13)
100 |  43.73 (0.06) | 16332.77 (0.05) | 189.64 (0.07) |  139.79 (0.07)
2 | 10 |  74.8 (0.12) | 16332.75 (0.12) | 239.11 (0.1) |  139.74 (0.13)
100 |  74.83 (0.05) | 16332.77 (0.07) | 239.06 (0.15) |  139.7 (0.12)
3 | 10 |  324.45 (0.99) | 16332.59 (0.06) | 628.33 (0.07) |  139.6 (0.15)
100 |  324.4 (1.02) | 16332.6 (0.16) | 628.27 (0.01) |  139.71 (0.16)

##  4 Experiments

###  4.1 Experimental Settings

In this section, we evaluate the performance of the proposed methods
111https://github.com/chenjie20/IPStatistics on three synthetic datasets that
contain 5 million, 10 million, and 50 million randomly generated IP records.
Each individual IP address contains one or more of IP records. The average
number of IP records is 100 for each individual IP address. Parameter $k$
represents the number of frequently occurring IP addresses. The statistics of
these datasets are summarized in Table 2.

We compared the proposed method with the following baselines:

  1. $\bullet$

Hash Mapping. Each IP record is mapped into an entry with a statistical result
using a hash table 222https://github.com/activesys/libcstl. Next, the
statistical results are used to construct a minimum heap of size $k$.

  2. $\bullet$

IP Mapping. All IP records are partitioned into $q$ subsets according to the
first part of each IP address. The statistics of the IP records in each subset
are mapped into an array, whose memory is pre-allocated on a computer
according to the last three parts of each IP address. The first $k$ most
frequent IP addresses are chosen from each subset. Next, a minimum heap of
size $k$ is constructed using the $q\times k$ IP addresses.

Two metrics were employed to evaluate the sorting performance, that is,
computational cost and memory use. All experiments were implemented using the
C language on a Windows platform with an Intel i7-9700k CPU and 32 GB RAM.

###  4.2 Experiment Results

Parameter $k$ was set to 10 or 100, and we repeated each experiment 10 times.
The average computational costs and standard deviations are reported in Table
4, and the mean memory use and standard deviations are given in Table 4. The
results show that TLMB consistently outperformed all the other methods in
terms of computational cost. For example, TLMB achieves computational costs of
2.51 s and 2.52 s when $k=10$ and $k=100$, respectively. When the number of IP
addresses increases from 5 million to 50 million with $k=10$, the gap in
computational costs of TLMB and IP Mapping are 12.6 s and 11.11 s,
respectively. We also observed the same advantages when $k=100$. In addition,
SSMB shows competitive results when compared with the comparison methods in
terms of memory use. The memory use of SSMB is always approximately 139 MB,
even when the number of IP addresses changes. In contrast, the computational
cost of SSMB substantially outperforms that of Hash Mapping under different
numbers of IP addresses. IP Mapping obtained the lowest computation cost for
all numbers of IP addresses. However, the highest memory use results of IP
Mapping are consistent with expectations.

###  4.3 Discussion

The gap in computational cost between TLMB and Hash Mapping dramatically
increases when the number of IP records increases from 5 million to 50
million. This is because TLMB avoids hash collisions when an IP address is
mapped to the corresponding memory block. The computational cost of TLMB is
linearly proportional to the number of IP addresses. Moreover, the memory use
of SSMB remains almost unchanged regardless of the number of IP addresses.
This is consistent with the theory underlying the second proposed mapping
mechanism. There is a negligible effect on the computational costs of TLMB and
SSMB when $k$ increases from 10 to 100. Moreover, the changes in the memory
use of TLMB and computational cost of SSMB are tolerable in practical
applications as the number of IP addresses increases. Consequently, TLMB and
SSMB reach a reasonable balance between computational cost and memory use when
compared with Hash Mapping and IP Mapping. This reveals that the two
relationship mapping mechanisms for memory blocks and IP addresses are
effective approaches for the design of the four-dimensional sparse matrix.

##  5 Conclusion

The collection of the statistics of large-scale IP address data is one of the
most fundamental problems in network traffic measurement. In this paper, we
addressed this problem. Specifically, the two proposed methods present two
different relationship mapping mechanisms between memory blocks and IP
addresses to strike a balance between computational cost and memory use. They
can be employed to search for frequently occurring IP addresses in practical
applications. The extensive experimental results demonstrate the effectiveness
of the proposed methods.

## References

  * [1] Y. N. Jing, P. Tu, X. P. Wang, G. D. Zhang, Distributed-log-based scheme for ip traceback, in: The Fifth International Conference on Computer and Information Technology, Shanghai, China, 2005, pp. 711–715.
  * [2] E. Kapur, P. Kumar, S. Gupta, Proposal of a two way sorting algorithm and performance comparison with existing algorithms, International Journal of Computer Science, Engineering and Applications 2 (3) (2012) 61–78.
  * [3] S. T. Klein, On the connection between hamming codes, heapsort and other methods, Information Processing Letters 113 (17) (2017) 617–620.
  * [4] M. L. Fredman, An intuitive and simple bounding argument for quicksort, Information Processing Letters 3 (114) (2014) 137–139.
  * [5] A. Agapitos, S. M. Lucas., Evolving efficient recursive sorting algorithms, in: 2006 IEEE International Conference on Evolutionary Computation, Vancouver, BC, Canada, 2016, pp. 2677–2684.
  * [6] Y. N. Jing, P. Tu, X. P. Wang, G. D. Zhang, Space-efficient data structures for top-k completion, in: Proceedings of the 22nd international conference on World Wide Web, Rio de Janeiro, Brazil, 2013, pp. 583–594.
  * [7] T. Cormen, C.Leiserson, R. Rivest, C.Stein, Introduction to Algorithms, MIT press, 2009.
  * [8] J. Hammad, A comparative study between various sorting algorithms, International Journal of Computer Science and Network Security 15 (3) (2015) 358–367.
  * [9] S. J. Puglisi, W. F. Smyth, A. Turpin, The performance of linear time suffix sorting algorithms, in: Data Compression Conference, Snowbird, UT, USA, USA, 2005, pp. 358–367.
  * [10] F. A. Louza, S. Gog, G. P. Telles., Optimal suffix sorting and lcp array construction for constant alphabets, Information Processing Letters 118 (2017) 30–34.
  * [11] S. Abdel-Hafeez, A. Gordon-Ross, An efficient o($n$) comparison-free sorting algorithm, IEEE Transactions on Very Large Scale Integration (VLSI) Systems 25 (6) (2017) 1930–1942.
  * [12] M. Shabaz, A. Kumar, Sa sorting: a novel sorting technique for large-scale data, Journal of Computer Networks and Communications (2019) 1–7.
  * [13] I. Múller, P. Sanders, A. Lacurie, W. Lehner, F. Fárber, Cache-efficient aggregation: Hashing is sorting, in: Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data, Melbourne, Bictoria, Australia, 2015, p. 1123–1136.
  * [14] N. Satish, M. Harris, M. Garland, Designing efficient sorting algorithms for manycore gpus, in: 2009 IEEE International Symposium on Parallel & Distributed Processing, Rome, Italy, 2009, pp. 1–10.
  * [15] D. P. Singh, I. Joshi, J. Choudhary, Survey of gpu based sorting algorithms, International Journal of Parallel Programming 46 (6) (2018) 1017–1034.
  * [16] P. M. E. Shutler, S. W. Sim, W. Y. S. Lim, Analysis of linear time sorting algorithms, The Computer Journal 51 (4) (2008) 451–469.
  * [17] F. Idrizi, A. Rustemi, F. Dalipi, A new modified sorting algorithm: A comparison with state of the art, in: 2017 6th Mediterranean Conference on Embedded Computing, Bar, Montenegro, 2017, pp. 1–6.
  * [18] I. Rusu, Sorting signed permutations by reversals using link-cut trees, Information Processing Letters 132 (2018) 44–48.
  * [19] S. Jukna, H. Seiwert, Sorting can exponentially speed up pure dynamic programming, Information Processing Letters 159 (2020) 451–469.
  * [20] G. Kocher, N. Agrawal, Analysis and review of sorting algorithms, International Journal of Scientific Engineering and Research 2 (3) (2014) 81–84.
  * [21] N. Faujdar, S. P. Ghrera, Analysis and testing of sorting algorithms on a standard dataset, in: 2015 Fifth International Conference on Communication Systems and Network Technologies, Gwalior, India, 2015, pp. 1–10.
  * [22] S. Chen, J. Qin, Y. Xie, J. Zhao, P. A. Heng, A fast and flexible sorting algorithm with cuda, in: International Conference on Algorithms and Architectures for Parallel Processing, Taipei, Taiwan, China, 2009, pp. 281–290.
  * [23] D. Cederman, P. Tsigas, A practical quicksort algorithm for graphics processors, in: European Symposium on Algorithms, 2008, pp. 246–258.
  * [24] H. Peters, O. Schulz‐Hildebrandt, N. Luttenberger, Fast in‐place, comparison‐based sorting with cuda: a study with bitonic sort, Concurrency and Computation: Practice and Experience 23 (7) (2011) 681–693.


# Stabilizing the calculation of the self-energy in dynamical mean-field
theory using constrained residual minimization

Harrison LaBollita  Department of Physics, Arizona State University, Tempe, AZ
85287, USA  Jason Kaye  Center for Computational Quantum Physics, Flatiron
Institute, 162 5th Avenue, New York, NY 10010, USA  Center for Computational
Mathematics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Alexander Hampel  Center for Computational Quantum Physics, Flatiron
Institute, 162 5th Avenue, New York, NY 10010, USA

###### Abstract

We propose a simple and efficient method to calculate the electronic self-
energy in dynamical mean-field theory (DMFT), addressing a numerical
instability often encountered when solving the Dyson equation. Our approach
formulates the Dyson equation as a constrained optimization problem with a
simple quadratic objective. The constraints on the self-energy are obtained
via direct measurement of the leading order terms of its asymptotic expansion
within a continuous time quantum Monte Carlo framework, and the use of the
compact discrete Lehmann representation of the self-energy yields an
optimization problem in a modest number of unknowns. We benchmark our method
for the non-interacting Bethe lattice, as well as DMFT calculations for both
model systems and ab-initio applications.

##  I Introduction

Quantum impurity models, which generically describe a quantum impurity
interacting with a non-interacting fermionic bath Anderson (1961), play an
important role in the study of strongly correlated materials Gull  _et al._
(2011). In particular, they serve as auxiliary models in quantum embedding
techniques such as dynamical mean-field theory (DMFT) used to map an
interacting lattice problem to a quantum impurity problem subject to a self-
consistency condition Georges  _et al._ (1996). A common choice of quantum
impurity solver is continuous-time quantum Monte Carlo (CT-QMC) using the
hybridization expansion algorithm (CT-HYB) Werner and Millis (2006); Gull  _et
al._ (2011); Nolting (1972), in which observables of the quantum impurity,
such as the impurity Green’s function, the many-body impurity density matrix,
and higher-order correlation functions, are directly measured via sampling.

In DMFT and its ab-initio extension density functional theory (DFT) + DMFT,
the impurity Green’s function $G$ and its corresponding electronic self-energy
$\Sigma$, matrices with the dimension of the impurity space, are the primary
quantities of interest. Within the DMFT formalism, these quantities are
related via the Dyson equation, formulated in the Matsubara frequency domain
as

$G(i\nu_{n})=G_{0}(i\nu_{n})+G_{0}(i\nu_{n})\Sigma(i\nu_{n})G(i\nu_{n}),$ (1)

or equivalently

$\Sigma(i\nu_{n})=G_{0}^{-1}(i\nu_{n})-G^{-1}(i\nu_{n}).$ (2)

Here, $\nu_{n}\equiv(2n+1)\pi/\beta$ is a fermionic Matsubara frequency at
inverse temperature $\beta$, and $G_{0}$ is a bare Green’s function, which
contains the dynamical Weiss field describing the impurity bath coupling and
the non-interacting part of the impurity Hamiltonian. It is well-known that
computing the self-energy using the direct formula (2) leads to an
amplification of noise, typically arising from Monte Carlo-based impurity
solvers, which grows with increasing $\nu_{n}$ Gull  _et al._ (2011). The
mechanism of this phenomenon will be discussed in this article. Since high-
frequency information contributes significantly to quantities, such as the
Weiss field, required to satisfy the DMFT self-consistency condition,
excessive noise ultimately prevents the convergence of the DMFT loop and
charge self-consistency within DFT+DMFT. Several methods are used to address
this issue, the most common of which we briefly summarize.

An approach referred to as “tail-fitting” attempts to address the problem
directly Haule  _et al._ (2010); Seth  _et al._ (2016a). The coefficients of a
high-frequency asymptotic expansion
$\Sigma(i\nu_{n})\approx\sum_{j=0}^{N}\Sigma_{j}/(i\nu_{n})^{j}$ of the self-
energy are fit from data obtained using (2) in some fitting window
($\nu_{n}^{\mathrm{min}}$, $\nu_{n}^{\mathrm{max}}$). The window must be
chosen to cover sufficiently high frequencies so that the asymptotic expansion
is valid, but not so high that the data for the fit is overwhelmed by noise.
This asymptotic expansion is then used in place of (2) at high frequencies.
Determining the location and width of the fitting window, the number of terms
to include in the asymptotic expansion, and the point at which to “glue” the
asymptotic expansion to the representation (2) requires substantial tuning on
a case-by-case basis. This often makes it difficult to obtain well-converged
DMFT results, and represents a significant obstacle to the construction of
reliable, black-box DMFT codes.

Another approach is to expand the imaginary time Green’s function in a basis
of orthogonal polynomials, such as Legendre polynomials Boehnke  _et al._
(2011), constraining the expansion coefficients using the known high-frequency
asymptotics of the Green’s function (in particular, $G(i\nu_{n})\sim
1/(i\nu_{n})$ as $n\to\infty$). The expansion coefficients can be directly
measured, at some computational expense, within the CT-HYB algorithm.
Enforcing the correct leading-order asymptotic behavior of the Green’s
function reduces, but does not eliminate, the amplification of high-frequency
noise in the self-energy. The performance of this approach is also sensitive
to the size of the basis, which is challenging to select in a manner that fits
the data accurately without overfitting to noise.

A third approach, based on the equations of motion, is the method of improved
estimators (IE) Hafermann  _et al._ (2012); Hafermann (2014); Gunacker  _et
al._ (2016), and more recently symmetric improved estimators (sIE) Kaufmann
_et al._ (2019); Kugler (2022). These involve the derivation of higher-order
correlation functions related to the self-energy (and vertex) corrections,
which can be directly measured within CT-HYB Hafermann  _et al._ (2012). While
the sIE approach is effective in stabilizing the self-energy calculation, it
requires some implementation effort and a significant increase in
computational cost, since unlike the impurity Green’s function in CT-HYB, the
sIEs cannot be measured efficiently by partition sum sampling Kaufmann  _et
al._ (2019).

Our approach is based on the observation that the problem of high frequency
noise amplification in the self-energy is simply the result of a numerical
instability in (2): since both $G_{0}^{-1}(i\nu_{n})$ and $G^{-1}(i\nu_{n})$
grow with $\nu_{n}$, but $\Sigma(i\nu_{n})$ does not, the subtraction causes a
loss of accuracy due to catastrophic cancellation. To address this, we
therefore avoid (2) and instead obtain $\Sigma$ by using (1) to set up a
constrained minimization problem: minimize the residual $\norm{G-G_{0}\Sigma
G-G_{0}}$, in a suitable norm $\norm{\cdot}$, subject to a constraint on the
high-frequency behavior of $\Sigma$ which can be determined directly within
the CT-HYB algorithm. Discretizing all quantities using the compact discrete
Lehmann representation (DLR) Kaye  _et al._ (2022a) leads to a quadratic
minimization problem with a scalar constraint in a small number of degrees of
freedom, even at very low temperatures and large spectral widths. The
resulting constrained residual minimization (CRM) scheme is simple to
implement using existing codes, with a negligible computational cost. In our
numerical experiments, we find robust convergence of the self-energy with no
tuning of parameters, leading to accuracy at the level of the Monte Carlo
noise with the correct high-frequency behavior.

This paper is organized as follows. After introducing the required numerical
tools in Sec. II and describing our method in Secs. III.1 and III.2, we
validate our method for the exactly solvable non-interacting single-band Bethe
lattice in Sec. III.3. In Sec. IV, we then demonstrate the full procedure for
two representative examples: the interacting Bethe lattice, and an ab-initio
DMFT calculation for Sr2RuO4.

##  II Preliminaries

Our method relies on two primary tools: the high-frequency expansion of the
self-energy, and a compact representation of imaginary time and imaginary
frequency quantities called the DLR.

###  II.1 High-frequency expansion of the self-energy

A generic quantum impurity model Hamiltonian is given by

$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\mathrm{int}},$ (3)

with

$\mathcal{H}_{0}=\sum_{k\alpha}\varepsilon_{k\alpha}c^{\dagger}_{k\alpha}c_{k\alpha}\\\
+\sum_{k\alpha b}\big{(}V_{k}^{\alpha
b}c^{\dagger}_{k\alpha}d_{b}+\mathrm{h.c.}\big{)}+\sum_{ab}E_{ab}d^{\dagger}_{a}d_{b}$

and

$\mathcal{H}_{\mathrm{int}}=\frac{1}{2}\sum_{abcd}U_{abcd}d^{\dagger}_{a}d^{\dagger}_{b}d_{d}d_{c}.$

Here, $c^{\dagger}_{\alpha}$ and $d^{\dagger}_{a}$ are the fermion creation
operators in the bath and impurity, respectively. $E_{ab}$ captures the bare
energy structure for the quantum impurity, $U_{abcd}$ parameterizes the
electron-electron interactions, $V_{k}^{\alpha b}$ describes the hybridization
between the impurity and the bath, and $\varepsilon_{k\alpha}$ is the energy-
momentum dispersion of the bath.

The electronic self-energy of the quantum impurity problem can be expanded in
the $\nu_{n}\to\infty$ limit as

$\Sigma(i\nu_{n})=\Sigma_{0}+\frac{\Sigma_{1}}{i\nu_{n}}+\mathcal{O}\big{(}\nu_{n}^{-2}\big{)},$
(4)

where the matrices $\Sigma_{0}$ and $\Sigma_{1}$ can be obtained from the
high-frequency expansions of the full and non-interacting Green’s functions.
They are formally defined in terms of expectation values of commutators and
anti-commutators,

$\displaystyle\Sigma_{0}^{ab}$
$\displaystyle=-\langle\\{[\mathcal{H}_{\mathrm{int}},d_{a}],d_{b}^{\dagger}\\}\rangle$
(5) $\displaystyle\Sigma_{1}^{ab}$
$\displaystyle=\langle\\{[\mathcal{H}_{\mathrm{int}},[\mathcal{H}_{\mathrm{int}},d_{a}]],d_{b}^{\dagger}\rangle-(\Sigma_{0}^{ab})^{2},$

which can be computed to high accuracy within CT-HYB from the many-body
impurity density matrix at no additional computational cost Wang  _et al._
(2011); Gull  _et al._ (2011); Haule (2007). Formally, the full Hamiltonian
(3) should be used in the commutators, but all non-trivial terms from its non-
interacting part cancel in $\Sigma_{0}$ and $\Sigma_{1}$.

The high-level abstractions used in the TRIQS software library Parcollet  _et
al._ (2015) allow for the direct implementation of the commutators in (5).
Such an implementation is general, and in particular agnostic to the exact
form of $\mathcal{H}_{\mathrm{int}}$. Further details on the derivation and
implementation of (5) are given in Appendix A.

###  II.2 Discrete Lehmann representation

We give a brief overview of the DLR, and refer the reader to Refs. Kaye  _et
al._ , 2022a, b for further details.

The DLR of an imaginary time Green’s function is given by an expansion of the
form

$G(\tau)\approx\sum_{l=1}^{r}K(\tau,\omega_{l})\,\widehat{g}_{l},$ (6)

where $K(\tau,\omega)=-\frac{e^{-\tau\omega}}{1+e^{-\beta\omega}}$ is the
analytic continuation kernel appearing in the spectral Lehmann representation

$G(\tau)=\int_{-\infty}^{\infty}d\omega\,K(\tau,\omega)\rho(\omega)$ (7)

relating $G$ to its spectral density $\rho$. The DLR frequencies $\omega_{l}$
are obtained by applying a rank-revealing pivoted Gram-Schmidt algorithm to
the analytic continuation kernel, and are therefore independent of the
specific Green’s function $G$. They depend only on the desired accuracy
$\epsilon$ of the DLR expansion (6), and a dimensionless cutoff parameter
$\Lambda=\beta\omega_{\max}$, where $\omega_{\max}$ is such that
$\rho(\omega)=0$ outside $[-\omega_{\max},\omega_{\max}]$. Both $\epsilon$ and
$\Lambda$ are user-specified parameters. To determine $\Lambda$, typically
$\beta$ is given, and $\omega_{\max}$ can be estimated, or results can be
converged with respect to $\Lambda$. For any imaginary time Green’s function
$G$ obeying the cutoff $\Lambda$, there are coefficients $\widehat{g}_{l}$
such that (6) holds to accuracy $\epsilon$; in other words, the functions
$K(\tau,\omega_{l})$ span the space of all such Green’s functions. The
existence of the DLR expansion is derived from a low-rank approximation of the
analytic continuation operator (7), and crucially, it is observed that
$r=\mathcal{O}(\log(\Lambda)\log(\epsilon^{-1}))$, yielding an exceptionally
compact representation.

The DLR coefficients $\widehat{g}_{l}$ can be determined by fitting the
expansion (6) to data, or by interpolation at a collection of $r$ DLR nodes
$\tau_{k}$ also obtained using the pivoted Gram-Schmidt algorithm. The latter
requires solving an $r\times r$ linear system obtained by evaluating the
expression (6) at $\tau=\tau_{k}$. Fourier transform of (6) yields a DLR
expansion

$G(i\nu_{n})\approx\sum_{l=1}^{r}K(i\nu_{n},\omega_{l})\,\widehat{g}_{l}$ (8)

in the Matsubara frequency domain, with
$K(i\nu_{n},\omega)=(i\nu_{n}-\omega)^{-1}$ for fermionic Green’s functions.
As in imaginary time, the coefficients $\widehat{g}_{l}$ can be determined by
interpolation at a collection of $r$ DLR nodes $i\nu_{n_{k}}$. Once the DLR
coefficients $\widehat{g}_{l}$ have been obtained, the DLR expansion can be
evaluated both in the imaginary time and frequency domains, so the Fourier
transform is in effect performed analytically.

We note that the DLR is a close cousin of the intermediate representation (IR)
Shinaoka  _et al._ (2017); Chikano  _et al._ (2018), which uses a singular
value decomposition of the analytic continuation kernel to obtain an
orthogonal expansion in terms of numerically-represented basis functions,
rather than the DLR’s non-orthogonal expansion in terms of analytically-known
basis functions. The number of terms in the DLR and IR expansions at given
$\Lambda$ and $\epsilon$ is similar, with the IR containing slightly fewer. An
IR expansion can also be obtained using fitting or a similar interpolation
procedure, with interpolation nodes determined either in the manner described
above Kaye  _et al._ (2022a), or using the sparse sampling method Li  _et al._
(2020); Shinaoka  _et al._ (2022). Both the DLR and IR have been effectively
used to represent imaginary time and frequency quantities other than Green’s
functions, such as self-energies and hybridization functions Li  _et al._
(2020); Kaye and U. R. Strand (2023); Yeh  _et al._ (2022); Cai  _et al._
(2022); Shinaoka  _et al._ (2022); Sheng  _et al._ (2023); Kaye  _et al._
(2023).

##  III Constrained residual minimization algorithm

###  III.1 Instability in the direct solution of the Dyson equation

The formula (2) for $\Sigma$ is numerically unstable, causing an amplification
of noise in $G_{0}$ and $G$ at high frequencies. To demonstrate this, suppose
we perturb $G(i\nu_{n})$ by a noise function $\eta(i\nu_{n})$, and let
$\widetilde{\Sigma}$ be the approximation of $\Sigma$ obtained from (2). We
find

$\displaystyle\widetilde{\Sigma}(i\nu_{n})$
$\displaystyle=G_{0}(i\nu_{n})^{-1}-(G(i\nu_{n})+\eta(i\nu_{n}))^{-1}$
$\displaystyle=G_{0}(i\nu_{n})^{-1}-G(i\nu_{n})^{-1}(I+\eta(i\nu_{n})G(i\nu_{n})^{-1})^{-1}$
$\displaystyle=G_{0}(i\nu_{n})^{-1}-G(i\nu_{n})^{-1}\sum_{k=0}^{\infty}(\eta(i\nu_{n})G(i\nu_{n})^{-1})^{k}$
$\displaystyle=\Sigma(i\nu_{n})-G(i\nu_{n})^{-1}\sum_{k=1}^{\infty}(\eta(i\nu_{n})G(i\nu_{n})^{-1})^{k}$

as long as $\eta$ is sufficiently small. We focus on the case of scalar-valued
Green’s function and self-energy for simplicity. Since
$G(i\nu_{n})\sim(i\nu_{n})^{-1}$ as $n\to\infty$, we have
$G(i\nu_{n})^{-1}\sim i\nu_{n}$, yielding an expansion of the error of the
form

$\eta\nu_{n}^{2}+\eta^{2}\nu_{n}^{3}+\cdots.$

Here we have abused notation and ignored the frequency dependence of $\eta$.
This error growth can be viewed as a consequence of catastrophic cancellation:
the growing quantities $G_{0}^{-1}$ and $G^{-1}$ are subtracted to obtain a
bounded quantity $\Sigma$, leading to an amplification of noise. We note that
while this explanation assumes roughly uniform noise in the Matsubara
frequency domain, we do not claim this is necessarily the correct error model
for any given Monte Carlo-based impurity solver. However, it is a simple model
of the mechanism of high frequency noise amplification, which reproduces the
rate of noise growth often observed in practice.

###  III.2 Constrained residual minimization

To avoid the numerically unstable formula (2), we instead begin with (1),
defining its residual as

$R=G-G_{0}\Sigma G-G_{0}.$ (9)

Here, all quantities are considered to be numerical approximations, so that
$R\neq 0$ as it would be if (1) were satisfied exactly. Given the leading
coefficients $\Sigma_{0}$ and $\Sigma_{1}$ of the asymptotic expansion (4),
which can be computed directly as described in Section II.1, $\Sigma$ can be
obtained from $G$ and $G_{0}$ by minimizing $\norm{R}$ in a suitable norm
subject to the constraint that (4) holds. We refer to this approach as
constrained residual minimization (CRM).

We discretize (9) by evaluating it at the DLR Matsubara frequency nodes
$i\nu_{n_{k}}$. Thus, we define $r_{k}\equiv R(i\nu_{n_{k}})$, $g_{k}\equiv
G(i\nu_{n_{k}})$, $g_{0k}\equiv G_{0}(i\nu_{n_{k}})$, and
$\sigma_{k}\equiv\Sigma(i\nu_{n_{k}})$. As described in Section II.2, the DLR
expansions of $R$, $G$, and $G_{0}$ can be recovered from the values $g_{k}$
and $g_{0k}$. For $\Sigma$, we must first subtract its known constant part
$\Sigma_{0}$ for the DLR expansion to be valid. Thus, defining
$\overline{\Sigma}\equiv\Sigma-\Sigma_{0}$, we can recover the DLR expansion
of $\overline{\Sigma}$ from the values $\sigma_{k}-\Sigma_{0}$, which gives a
method of recovering $\Sigma$, and also implements the constraint on the
constant term given in (4). The discretization of (9) is therefore given by

$r_{k}\equiv g_{k}-g_{0_{k}}\sigma_{k}g_{k}-g_{0_{k}}.$ (10)

Figure 1: (a) $G(\tau)$ for the non-interacting Bethe lattice with $\beta=100$
eV-1. (b) $\imaginary G(i\nu_{n})$; note $\text{Re}\,G(i\nu_{n})=0$ by
symmetry. (c) $\imaginary\Sigma(i\nu_{n})$ computed via the direct formula (2)
and the CRM method, using data obtained in imaginary time with noise level
$\eta=10^{-4}$, along with the exact solution $\Sigma=G/4$. (d) Pointwise
error for both methods, for several choices of $\eta$, along with the
predicted error for the direct method.

We measure the residual (9) in the $L^{2}(\tau)$ norm, given by

$\norm{R}_{L^{2}(\tau)}^{2}\equiv\frac{1}{\beta}\int_{0}^{\beta}d\tau\,\norm{R(\tau)}_{F}^{2},$
(11)

where $\norm{\cdot}_{F}$ is the Frobenius norm. We show in Appendix B that
(11) can be computed from a simple quadratic form evaluated at the DLR
coefficients $\widehat{r}_{l}$ of $R$, which can be obtained from the values
$r_{k}$ by interpolation.

To implement the constraint on the first decaying term in (4), we equate the
DLR expansion of $\overline{\Sigma}$ and the desired asymptotic expansion:

$\begin{multlined}\sum_{l=1}^{r}K(i\nu_{n},\omega_{l})\widehat{\sigma}_{l}=\overline{\Sigma}(i\nu_{n})\\\
=\Sigma(i\nu_{n})-\Sigma_{0}=\frac{\Sigma_{1}}{i\nu_{n}}+\mathcal{O}(\nu_{n}^{-2}).\end{multlined}\sum_{l=1}^{r}K(i\nu_{n},\omega_{l})\widehat{\sigma}_{l}=\overline{\Sigma}(i\nu_{n})\\\
=\Sigma(i\nu_{n})-\Sigma_{0}=\frac{\Sigma_{1}}{i\nu_{n}}+\mathcal{O}(\nu_{n}^{-2}).$

Multiplying both sides by $i\nu_{n}$ and taking the limit $n\to\infty$ yields
the desired constraint

$\sum_{l=1}^{r}\widehat{\sigma}_{l}=\Sigma_{1},$

since $K(i\nu_{n},\omega)=(i\nu_{n}-\omega)^{-1}$.

The discretized CRM problem is therefore given as follows:

$\operatorname*{arg\,min}_{\sigma_{k}}\norm{R}_{L^{2}(\tau)}\quad\text{subject
to}\quad\sum_{l=1}^{r}\widehat{\sigma}_{l}=\Sigma_{1}.$ (12)

Once it has been solved, $\Sigma$ is given by the DLR expansion

$\Sigma(i\nu_{n})=\Sigma_{0}+\sum_{l=1}^{r}K(i\nu_{n},\omega_{l})\widehat{\sigma}_{l}.$

The objective function $\norm{R}_{L^{2}(\tau)}$ can be computed from
$\sigma_{k}$ using (10) and the procedure described in Appendix B, since
$g_{k}$ and $g_{0k}$ are given. The DLR coefficients $\widehat{\sigma}_{l}$
are recovered from the $\sigma_{k}$ by solving the $r\times r$ linear system

$\sum_{l=1}^{r}K(i\nu_{n_{k}},\omega_{l})\widehat{\sigma}_{l}=\sigma_{k}.$

In our numerical examples, we use the SciPy Virtanen  _et al._ (2020)
implementation of the trust region optimization method Conn  _et al._ (2000)
to solve (12). We note that the CRM problem is simply a quadratic minimization
problem with a scalar constraint, and its solution is therefore given semi-
analytically in terms of simple matrix operations. However, we find that the
trust region optimization tends to perform slightly better, due to the
appearance of ill-conditioned operations in the semi-analytic formula.

We note that since the CRM scheme seeks a balance between correctly matching
the low-frequency data via the Dyson equation, and high-frequency asymptotics
of relatively low-order, one might obtain a larger error for intermediate
frequencies in which neither of these limiting representations is sufficiently
accurate. In our numerical examples, which include realistic DMFT
calculations, we have not observed a significant contribution of this effect
to the error. A topic of our future work is the possibility of directly
measuring higher-order asymptotics, which would push the region of validity of
the asymptotic expansion to lower frequencies. In particular, we note that the
next term of the asymptotic expansion can be obtained given an approximation
of the bath hybridization by a finite number of levels, for which there exist
a variety of algorithms Georges  _et al._ (1996); Mejuto-Zaera  _et al._
(2020); Shinaoka and Nagai (2021); Huang  _et al._ (2023).

###  III.3 Benchmark: non-interacting Bethe lattice

We first test our algorithm on the non-interacting Bethe lattice, given by (1)
with $\Sigma=G/4$ and $G_{0}(i\nu_{n})=(i\nu_{n})^{-1}$. This model can be
solved analytically in real frequency Georges  _et al._ (1996), yielding the
semi-circular spectral function
$\rho(\omega)=\frac{2}{\pi}\sqrt{1-\omega^{2}}$ on $[-1,1]$. $G(\tau)$ and
$G(i\nu_{n})$ can then be computed from the Lehmann representation (7) and its
Fourier transform (obtained by replacing $K(\tau,\omega)$ with
$K(i\nu_{n},\omega)$) by numerical integration. These are shown in Figs. 1(a)
and 1(b). This directly yields a high-accuracy reference solution for
$\Sigma=G/4$.

To simulate the output of a Monte Carlo-based impurity solver in imaginary
time, we evaluate $G(\tau)$ at a large collection of random points $\tau$, and
perturb the result by random noise of magnitude $\eta$. We then obtain a DLR
expansion of $G$ by least squares fitting this data. This expansion can be
evaluated in the Matsubara frequency domain using (8). In this case, we have
$\omega_{\max}=1$, so we use the DLR cutoff parameter $\Lambda=\beta$, and set
the DLR accuracy parameter to $\epsilon=\eta$. We calculate the leading
coefficients of the asymptotic expansion of $\Sigma$ analytically:
$\Sigma_{0}=0$, and $\Sigma_{1}=1/4$. We then obtain $\Sigma$ using both the
direct formula (2) and the CRM procedure (12), as shown in Fig. 1(c) for
$\beta=100$ eV-1 and $\eta=10^{-4}$. The pointwise error of
$\Sigma(i\nu_{n})$, computed using the known result, is shown for both methods
at $\eta=10^{-4}$, $10^{-6}$, and $10^{-8}$ in Fig. 1(d). We observe the
$\mathcal{O}(\eta\,\nu_{n}^{2})$ leading order error growth derived in Section
III.1 for the self-energy obtained from (2). The CRM approach yields uniform
accuracy at the noise level $\eta$ until the asymptotic region, at which point
the correctly-captured leading-order asymptotics yields an error decreasing at
the expected $\mathcal{O}(\nu_{n}^{-2})$ rate.

We note that even if $G$ were obtained with very little noise, the high
frequency noise amplification eventually leads to incorrect tail behavior
using (2), whereas CRM yields a DLR expansion of $\Sigma$ with the correct
asymptotic behavior on the full Matsubara frequency domain.

##  IV Numerical examples

We demonstrate the CRM method for two prototypical DMFT calculations: the
interacting Bethe lattice and a three-band Hubbard model of Sr2RuO4.

Figure 2: Interacting Bethe lattice for $U/t=4$ and $\beta t=5$. (a)
$\imaginary\Sigma(i\nu_{n})$ obtained using CRM for several choices of the DLR
cutoff parameter $\Lambda$, with DLR tolerance $\epsilon=10^{-6}$ and
$N=10^{9}$ QMC samples, along with the reference $\Sigma_{\text{leg}}$
obtained using a Legendre polynomial expansion of the Green’s function. (b)
Pointwise absolute difference of $\Sigma$ and $\Sigma_{\text{leg}}$. (c)
Maximum absolute difference of $\Sigma$ and $\Sigma_{\text{leg}}$, along with
self-convergence error of $\Sigma$ versus $\Lambda$, for $N=10^{6}$ and
$10^{9}$ QMC samples. (d) Spectral function $-\imaginary G(\omega)/\pi$
approximated using Padé analytic continuation.

###  IV.1 Interacting Bethe lattice

We begin with the half-filled, single-orbital Bethe lattice, for which the
non-interacting Green’s function is determined by a semi-circular density of
states (DOS). The interaction Hamiltonian is of the simple form
$\mathcal{H}=Un_{\uparrow}n_{\downarrow}$, where we have set (in units of the
hopping $t$) $U/t=4$, and the system inverse temperature to $\beta t=5$. We
solve the corresponding quantum impurity problem using the CT-HYB impurity
solver as implemented in TRIQS/cthyb Seth  _et al._ (2016b). In order to
simplify the experiment, we perform only one DMFT cycle, which is sufficient
to test the CRM method for various choices of calculation parameters. The
density matrix $\rho_{\text{imp}}$ is measured within CT-HYB, which allows us
to compute the expectation values for the high-frequency expansion
coefficients of the self-energy given in (5).

In TRIQS/cthyb, samples of the imaginary time Green’s function are averaged
within bins of a fixed width, leading to a bin width parameter and a
corresponding binning error. In our experiments, we use $10^{4}$ equispaced
bins, which is sufficient to make the binning error negligible and to ensure
that each bin contains a sufficient number of samples. We note that properly
adjusting the number of bins requires some care, and an important future
problem is to design an algorithm which maintains the efficiency of binning
with substantially reduced errors, e.g., via projection onto a higher-order
basis. For all experiments, we set the DLR tolerance to $\epsilon=10^{-6}$,
which is below the Monte Carlo error, and we use CRM to obtain the impurity
self-energy. We compare our results to a reference self-energy
$\Sigma_{\text{leg}}$, obtained using the well-established Legendre polynomial
expansion of the Green’s function Boehnke  _et al._ (2011), with Legendre
coefficients measured directly within CT-HYB using a total of $10^{11}$ QMC
samples, a roughly $6000$ core-hour calculation.

Figure 3: Comparison between constrained residual minimization (CRM), improved
estimators (IE), and symmetric improved estimators (sIE) for the interacting
Bethe lattice example. The CT-HYB solvers were run using $10^{8}$ QMC samples.
For CRM, we used $\Lambda=20$ and $\varepsilon=10^{-6}$.

Using a CRM calculation with $N=10^{9}$ QMC samples, in Fig. 2(a) we plot
$\Im{\Sigma(i\nu_{n})}$ for several choices of the DLR cutoff parameter
$\Lambda$, along with the reference, and we plot their corresponding absolute
pointwise difference in Fig. 2(b). We observe convergence to the reference as
$\Lambda$ is increased. We note that although $\Sigma_{\text{leg}}$ was
calculated with a larger number of QMC samples, its accuracy is limited by the
presence of Monte Carlo noise which is somewhat amplified at high frequencies.
It is therefore possible that the CRM result is more accurate at high
frequencies. Fig. 2(c) shows the maximum absolute pointwise difference against
$\Lambda$ (blue curves) for $N=10^{6}$ and $N=10^{9}$ QMC samples. We also
show a “self-convergence” error (red curves), given by the maximum of
$\absolutevalue{\Sigma(i\nu_{n})-\Sigma_{\text{last}}(i\nu_{n})}$, where
$\Sigma$ is the self-energy computed with a given choice of $\Lambda$ and
$\Sigma_{\text{last}}$ is the self-energy computed with the largest value of
$\Lambda$ used, which is $\Lambda=50$. Since in practice one does not have a
known reference solution, the choice of a converged $\Lambda$ parameter can be
determined by starting with a physically-motivated guess (based on an estimate
of $\omega_{\max}$), and plotting the self-convergence error against $\Lambda$
in this manner. Both error metrics consistently suggest convergence to the QMC
noise level by $\Lambda=20$, and this noise level decreases at the expected
rate as the number of QMC samples is increased (by roughly a factor of $30$
for a $1000$-fold increase in the number of samples). Fig. 2(d) shows the
spectral function (interacting DOS) obtained via Padé analytic continuation
Ferris-Prabhu and Withers (1973) for our single-orbital model. We find a
spectral width of roughly $\omega_{\mathrm{max}}=4$ eV, consistent with the
observed convergence near $\Lambda=\beta\omega_{\max}=20$.

We lastly compare the CRM algorithm to the methods of improved estimators (IE)
and symmetric improved estimators (sIE), in which the self-energy is obtained
by measuring high-order correlation functions Hafermann  _et al._ (2012);
Hafermann (2014); Gunacker  _et al._ (2016); Kaufmann  _et al._ (2019); Kugler
(2022). We use the implementations of IE and sIE provided by the w2dynamics
code Wallerberger  _et al._ (2019), and used $N=10^{8}$ QMC samples in all
cases. Fig. 3 shows the pointwise absolute difference of $\Sigma(i\nu_{n})$
for each method and the reference self-energy $\Sigma_{\text{leg}}$, along
with that using the direct formula (2). All methods agree in the low-frequency
regime, and use of the direct formula produces the expected loss of accuracy
at high frequencies. All three alternative methods outperform the direct
formula (2), but the error of IE is still too large to be useful in practical
computations. sIE and CRM outperform the other approaches, and yield
comparable results. However, the sIE results come at a significant additional
computational cost: for the simple one-band model considered here, the cost of
measuring the sIE was roughly $10$ times larger than CRM to achieve similar
accuracy. In particular, the cost of CRM does not scale with the size of the
impurity problem, in contrast to the IE and sIE methods.

Figure 4: Comparison between the TRIQS tail-fitting algorithm and CRM for the
Sr2RuO4 example. (a,b) Real and imaginary parts of $\Sigma^{11}(i\nu_{n})$ in
the low to mid-frequency regime. The window sizes correspond to
$(\nu_{n}^{\mathrm{min}},\nu_{n}^{\mathrm{max}})=$ (2,4) (window 1), (2,6)
(window 2), and (2, 10) (window 3). The dashed black line corresponds to the
directly measured high-frequency expansion of the self-energy. (c,d) Same as
(a,b), at higher frequency. (e,f) Convergence of the chemical potential $\mu$
and local Green’s function $G$ as measured by their differences $\delta\mu$
and $\delta G$ between subsequent iterations.

###  IV.2 Three-band Hubbard model of Sr2RuO4

We next consider an ab-initio-derived three-band Hubbard model, a common
DFT+DMFT workflow, at full DMFT self-consistency. We model the correlated
Hund’s metal Sr2RuO4, where the non-interacting low-energy physics consists of
three partially filled $t_{2g}\equiv\\{d_{xy},d_{xz},d_{yz}\\}$ bands from the
Ru($4d$) site. We downfold the Kohn-Sham Hamiltonian onto three maximally
localized Wannier functions, which serve as our correlated impurity subspace.
Our low-energy model is identical to that discussed in Ref. Sheng  _et al._ ,
2023, which contains further details.

Interactions are governed by a local rotationally-invariant Hubbard-Kanamori
Hamiltonian Mizokawa and Fujimori (1995) with $U=2.3$ eV and
$J_{\mathrm{H}}=0.4$ eV. We solve the impurity problem using the CT-HYB
algorithm implemented in TRIQS/cthyb Seth  _et al._ (2016b). We use the
inverse system temperature $\beta=40$ eV-1 ($\sim$290 K). In each step of the
self-consistent DMFT calculation, the CRM method is used to obtain the self-
energy from the impurity solver. We also perform the same calculation using
the standard TRIQS tail-fitting protocol, which operates according to the
description in Sec. I: an asymptotic expansion of a few terms is obtained by
fitting the data produced by (2) on a user-specified frequency window, and
glued to the result produced by (2). We produce results for three different
choices of the fitting window.

The real and imaginary parts of the first diagonal element of the self-energy
are shown in Figs. 4(a)-(d), respectively, using CRM and the tail-fitting
algorithm implemented in TRIQS, along with the direct formula (2). The results
for the other diagonal elements are similar, and the off-diagonal elements are
zero. The self-energy produced by the CRM algorithm correctly reproduces the
low-frequency result given by (2), and has the correct high-frequency
asymptotics, $\real\Sigma(i\nu_{n})\rightarrow\Sigma_{0}$ and
$\imaginary\Sigma(i\nu_{n})\sim\Sigma_{1}/(i\nu_{n})$, as imposed by CRM.

The solutions produced using the tail-fitting method are highly dependent on
the choice of the fitting window, which must be specified by the user. For too
large of a window (window 3: $\nu_{n}^{\min}=2$, $\nu_{n}^{\max}=10$), too
much high-frequency noise is incorporated, leading to large high-frequency
errors. For too small of a window (window 1: $\nu_{n}^{\min}=2$,
$\nu_{n}^{\max}=4$), the high-frequency asymptotics are not correctly
captured. In a “Goldilocks zone” fitting window, it is possible to balance
these competing sources of inaccuracy reasonably well, but in this experiment
even a well-chosen window (window 2: $\nu_{n}^{\min}=2$, $\nu_{n}^{\max}=6$)
is outperformed by CRM in reproducing the correct high-frequency asymptotics.
Evidently, the necessity of finely hand-tuning the tail-fitting
parameters—which must ideally be done for each non-equivalent orbital
separately and, in principle, at each DMFT iteration— represents a significant
barrier to producing black-box DMFT codes, and CRM eliminates the need for
such a procedure.

We also observe that the use of CRM helps to stabilize the DMFT loop by
eliminating fluctuations in the high-frequency behavior of the self-energy. In
Figs. 4(e,f), we show the convergence during the DMFT iteration of two
important quantities, the chemical potential $\mu$ and the local Green’s
function $G$, by plotting their differences $\delta\mu$ and $\delta G$ between
consecutive iterations. The CRM approach appears to stabilize the DMFT self-
consistency, leading to smoother and more rapid convergence to a higher-
accuracy result.

##  V Conclusion

We have proposed the CRM algorithm to compute the electronic self-energy in
DMFT calculations using CT-QMC-based quantum impurity solvers in a stable
manner. CRM is agnostic to the quantum impurity solver and can be used
whenever the high-frequency behavior of the electronic self-energy is known,
which includes applications outside of DMFT. The algorithm is simple to
implement, and can therefore replace otherwise cumbersome tail-fitting
procedures. Furthermore, our results suggest that CRM yields improved
stability and convergence of physical observables, and thereby improved
accuracy of calculations within DMFT and DFT+DMFT. Combining the CRM method
with the approach of Ref. Sheng  _et al._ , 2023, which demonstrated the
replacement of dense Matsubara frequency grids with compact grids of DLR
nodes, we envision a streamlined and automated implementation of the DMFT loop
with highly efficient discretizations of all quantities and minimal parameter-
tuning.

###### Acknowledgements.

We are grateful to H. Strand for helpful discussions on alternative
algorithms. We also thank F. Kugler, O. Parcollet, A. J. Millis, and N.
Wentzell for their input. HL acknowledges support from grant NSF-DMR 2045826.
The Flatiron Institute is a division of the Simons Foundation.

##  Appendix A  High-frequency asymptotics of the self-energy

The high-frequency expansion of the interacting Green’s function is given by

$G(i\nu_{n})=\frac{G_{1}}{i\nu_{n}}+\frac{G_{2}}{(i\nu_{n})^{2}}+\frac{G_{3}}{(i\nu_{n})^{3}}+\mathcal{O}\big{(}\nu_{n}^{-4}\big{)},$
(13)

with

$G_{n}=\partial_{\tau}^{n-1}G(0^{+})-\partial_{\tau}^{n-1}G(0^{-}).$

The first few coefficients are given by Gull  _et al._ (2011)

$\displaystyle G_{1}^{ab}$
$\displaystyle=\langle\\{d_{a},d_{b}^{\dagger}\\}\rangle=\delta_{ab},$ (14)
$\displaystyle G_{2}^{ab}$
$\displaystyle=-\langle\\{[\mathcal{H},d_{a}],d_{b}^{\dagger}\\}\rangle,$
$\displaystyle G_{3}^{ab}$
$\displaystyle=\langle\\{[\mathcal{H},d_{a}],[d_{b}^{\dagger},\mathcal{H}]\\}\rangle,$

where $\mathcal{H}$ is defined in (3). The analogous expansion for the non-
interacting Green’s function is obtained using the Hamiltonian (3) without the
interaction term.

To obtain the high-frequency expansion of the self-energy, we substitute the
expansions for the full and non-interacting Green’s function into (1).
Matching terms of like powers yields

$\displaystyle\Sigma(i\nu_{n})$
$\displaystyle=G_{2}-G_{02}+\frac{G_{3}-G_{03}+G_{02}^{2}-G_{2}^{2}}{i\nu_{n}}+\mathcal{O}\big{(}\nu_{n}^{-2}\big{)}$
$\displaystyle\equiv\Sigma_{0}+\frac{\Sigma_{1}}{i\nu_{n}}+\mathcal{O}\big{(}\nu_{n}^{-2}\big{)},$

where the coefficients of the asymptotic expansion of $G_{0}$ are given, in
analogy with (13), by $G_{01}$, $G_{02}$, etc. The expressions for
$\Sigma_{0}$ and $\Sigma_{1}$ in (5) then follow from (14).

11footnotetext: The dimer example is part of the TRIQS benchmarks repository,
which can be found at
https://github.com/TRIQS/benchmarks/blob/master/Dimer/notebook.ipynb. The
Jupyter notebook provides a description of the model and its implementation
using the TRIQS software.

We test our implementation on a simple model: a dimer coupled to a bath, which
can be solved by exact diagonalization. A description of the model,
calculation parameters, and implementation in TRIQS are provided in Ref. Note,
1. We obtain a reference self-energy using the exact diagonalization library
PyED Note (2). We use the TRIQS/cthyb Seth  _et al._ (2016b) solver to sample
the many-body density matrix $\rho_{\text{imp}}$, from which we compute the
first two coefficients of the asymptotic expansion of the self-energy using
(5). Fig. 5 compares an expansion $\Sigma_{0}+\Sigma_{1}/(i\nu_{n})$ to the
reference self-energy, showing close agreement in the high-frequency regime.
This implementation is available in latest release of TRIQS/cthyb Note (3).

22footnotetext: PyED: Exact diagonalization for finite quantum systems.
https://github.com/HugoStrand/pyed.33footnotetext: TRIQS/cthyb (v3.2.0)
https://github.com/TRIQS/cthyb/tree/3.2.x Figure 5: Asymptotic expansion
$\Sigma_{0}+\Sigma_{1}/(i\nu_{n})$ of the self-energy for the dimer example,
with coefficients computed from (5) using the measured many-body density
matrix within CT-HYB, compared with a reference computed using exact
diagonalization.

##  Appendix B $L^{2}(\tau)$ norm of a DLR expansion

We show in this Appendix how to compute $||G||_{L^{2}(\tau)}$ for a DLR
expansion $G$. In the CRM method, the objective function is given by
$||R||_{L^{2}(\tau)}$, with $R$ represented by its values $r_{k}$ at the DLR
Matsubara frequency nodes; the DLR coefficients $\widehat{r}_{l}$ can be
recovered via interpolation in the DLR basis. We focus on the case of scalar-
valued $G$, since the matrix-valued case follows directly.

Let

$G(\tau)=\sum_{l=1}^{r}K(\tau,\omega_{l})\widehat{g}_{l}$ (15)

be a DLR expansion. We have

$\displaystyle||G||_{L^{2}(\tau)}^{2}$
$\displaystyle\equiv\frac{1}{\beta}\int_{0}^{\beta}d\tau\,G^{2}(\tau)$
$\displaystyle=\frac{1}{\beta}\sum_{k,l=1}^{r}\hat{g}_{k}\hat{g}_{l}\int_{0}^{\beta}d\tau\,K(\tau,\omega_{k})K(\tau,\omega_{l})$
$\displaystyle\equiv\sum_{k,l=1}^{r}\widehat{g}_{k}M_{kl}\widehat{g}_{l},$

with the last equivalence defining the matrix $M$. Thus the squared norm is
given by a quadratic form with matrix $M$, evaluated at the vector of DLR
coefficients. $M$ is straightforwardly computed, and we write two equivalent
formulas:

$M_{kl}=K(0,\omega_{k})K(0,\omega_{l})\frac{1-e^{-\beta(\omega_{k}+\omega_{l})}}{\beta(\omega_{k}+\omega_{l})}$
(16)

and

$M_{kl}=\frac{K(0,\omega_{k})K(0,\omega_{l})-K(\beta,\omega_{k})K(\beta,\omega_{l})}{\beta(\omega_{k}+\omega_{l})}.$
(17)

The $\omega_{k}=-\omega_{l}$ case is given, via the limit, as
$K(0,\omega_{k})K(0,\omega_{l})$. When $\omega_{k}\approx-\omega_{l}$, both
formulas suffer from catastrophic cancellation, but this is straightforwardly
remedied in (16), e.g., using the expm1 function exp  for the stable
calculation of $e^{x}-1$. On the other hand, when
$\beta(\omega_{k}+\omega_{l})\ll 0$, (16) is vulnerable to numerical overflow,
whereas (17) is numerically stable. We therefore use (16), evaluated using the
expm1 function, when $0<\absolutevalue{\beta(\omega_{k}+\omega_{l})}<1$, (17)
when $\absolutevalue{\beta(\omega_{k}+\omega_{l})}\geq 1$, and the exact
formula when $\omega_{k}=-\omega_{l}$.

## References

  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
  * Nolting (1972) W. Nolting, Zeitschrift für Physik A Hadrons and nuclei 255, 25 (1972).
  * Haule  _et al._ (2010) K. Haule, C.-H. Yee, and K. Kim, Phys. Rev. B 81, 195107 (2010).
  * Seth  _et al._ (2016a) P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, Comput. Phys. Commun. 200, 274 (2016a).
  * Boehnke  _et al._ (2011) L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, and O. Parcollet, Phys. Rev. B 84, 075145 (2011).
  * Hafermann  _et al._ (2012) H. Hafermann, K. R. Patton, and P. Werner, Phys. Rev. B 85, 205106 (2012).
  * Hafermann (2014) H. Hafermann, Phys. Rev. B 89, 235128 (2014).
  * Gunacker  _et al._ (2016) P. Gunacker, M. Wallerberger, T. Ribic, A. Hausoel, G. Sangiovanni, and K. Held, Phys. Rev. B 94, 125153 (2016).
  * Kaufmann  _et al._ (2019) J. Kaufmann, P. Gunacker, A. Kowalski, G. Sangiovanni, and K. Held, Phys. Rev. B 100, 075119 (2019).
  * Kugler (2022) F. B. Kugler, Phys. Rev. B 105, 245132 (2022).
  * Kaye  _et al._ (2022a) J. Kaye, K. Chen, and O. Parcollet, Phys. Rev. B 105, 235115 (2022a).
  * Wang  _et al._ (2011) X. Wang, H. T. Dang, and A. J. Millis, Phys. Rev. B 84, 073104 (2011).
  * Haule (2007) K. Haule, Phys. Rev. B 75, 155113 (2007).
  * Parcollet  _et al._ (2015) O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, Comput. Phys. Commun. 196, 398 (2015).
  * Kaye  _et al._ (2022b) J. Kaye, K. Chen, and H. U. R. Strand, Comput. Phys. Commun. 280, 108458 (2022b).
  * Shinaoka  _et al._ (2017) H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, Phys. Rev. B 96, 035147 (2017).
  * Chikano  _et al._ (2018) N. Chikano, J. Otsuki, and H. Shinaoka, Phys. Rev. B 98, 035104 (2018).
  * Li  _et al._ (2020) J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, E. Gull, and H. Shinaoka, Phys. Rev. B 101, 035144 (2020).
  * Shinaoka  _et al._ (2022) H. Shinaoka, N. Chikano, E. Gull, J. Li, T. Nomoto, J. Otsuki, M. Wallerberger, T. Wang, and K. Yoshimi, SciPost Phys. Lect. Notes , 63 (2022).
  * Kaye and U. R. Strand (2023) J. Kaye and H. U. R. Strand, Adv. Comput. Math. 49 (2023), 10.1007/s10444-023-10067-7.
  * Yeh  _et al._ (2022) C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, Phys. Rev. B 106, 235104 (2022).
  * Cai  _et al._ (2022) X. Cai, T. Wang, N. V. Prokof’ev, B. V. Svistunov, and K. Chen, Phys. Rev. B 106, L220502 (2022).
  * Sheng  _et al._ (2023) N. Sheng, A. Hampel, S. Beck, O. Parcollet, N. Wentzell, J. Kaye, and K. Chen, Phys. Rev. B 107, 245123 (2023).
  * Kaye  _et al._ (2023) J. Kaye, H. U. R. Strand, and D. Golež, (2023), arXiv:2307.08566 .
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nat. Methods 17, 261 (2020).
  * Conn  _et al._ (2000) A. R. Conn, N. I. M. Gould, and P. L. Toint,  _Trust Region Methods_ (SIAM, 2000).
  * Mejuto-Zaera  _et al._ (2020) C. Mejuto-Zaera, L. Zepeda-Núñez, M. Lindsey, N. Tubman, B. Whaley, and L. Lin, Phys. Rev. B 101, 035143 (2020).
  * Shinaoka and Nagai (2021) H. Shinaoka and Y. Nagai, Phys. Rev. B 103, 045120 (2021).
  * Huang  _et al._ (2023) Z. Huang, E. Gull, and L. Lin, Phys. Rev. B 107, 075151 (2023).
  * Seth  _et al._ (2016b) P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, Comput. Phys. Commun. 200, 274 (2016b).
  * Ferris-Prabhu and Withers (1973) A. Ferris-Prabhu and D. Withers, J. Comput. Phys. 13, 94 (1973).
  * Wallerberger  _et al._ (2019) M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held, and G. Sangiovanni, Comput. Phys. Commun. 235, 388 (2019).
  * Mizokawa and Fujimori (1995) T. Mizokawa and A. Fujimori, Phys. Rev. B 51, 12880 (1995).
  * Note (1) The dimer example is part of the TRIQS benchmarks repository, which can be found at https://github.com/TRIQS/benchmarks/blob/master/Dimer/notebook.ipynb. The Jupyter notebook provides a description of the model and its implementation using the TRIQS software.
  * Note (2) PyED: Exact diagonalization for finite quantum systems. https://github.com/HugoStrand/pyed.
  * Note (3) TRIQS/cthyb (v3.2.0) https://github.com/TRIQS/cthyb/tree/3.2.x.
  * (40) https://numpy.org/doc/stable/reference/generated/numpy.expm1.html.