Compressing the two-particle Green’s function using wavelets Theory and application to the Hubbard atom

Precise algorithms capable of providing controlled solutions in the presence of strong interactions are transforming the landscape of quantum many-body physics. Particularly exciting breakthroughs are enabling the computation of non-zero temperature correlation functions. However, computational challenges arise due to constraints in resources and memory limitations, especially in scenarios involving complex Green’s functions and lattice effects. Leveraging the principles of signal processing and data compression, this paper explores the wavelet decomposition as a versatile and efficient method for obtaining compact and resource-efficient representations of the many-body theory of interacting systems. The effectiveness of the wavelet decomposition is illustrated through its application to the representation of generalized sus-ceptibilities and self-energies in a prototypical interacting fermionic system, namely the Hubbard model at half-filling in its atomic limit. These results are the first proof-of-principle application of the wavelet compression within the realm of many-body physics and demonstrate the potential of this wavelet-based compression scheme for understanding the physics of correlated electron systems.


Introduction
The emergence of highly precise algorithms capable of providing controlled solutions even within the realm of strong interactions is transforming the field of quantum many-body physics [1][2][3][4][5][6][7].In various practical scenarios though, memory limitations and constraints on computational resources make achieving the necessary precision challenging.Particular challenges include the storage problems associated with two particle Green functions [8], which depend on three independent frequencies and three independent momenta, and even for single-particle Green's functions when many degrees of freedom are involved, such as in the presence of spin-orbit coupling or in multi-orbitals applications, as often encountered, for instance, in ab-initio many-body calculations [9].The large size of these quantities also imposes constraints on the ability to manipulate and operate on these quantities.The Matsubara formalism commonly used for thermodynamic properties of systems at non-zero temperature [10] leads to a polynomial growth in size with respect to inverse temperature β, rapidly surpassing the capabilities of even state-of-the-art high-performance computing architectures.These challenges underscore the need for developing efficient compression schemes that allow, on the one hand, to conveniently store the many-body quantities, while still permitting, on the other hand, to perform operations typical to many-body physics.
By taking advantage of the specific mathematical structure of the imaginary time Green's function G(τ ) given by the spectral Lehmann representation, two novel representation strategies leverage the fact that the analytic continuation kernel or Lehmann kernel can be approximated to high accuracy by a low rank decomposition to derive compressed representations of G(τ ).The first approach, known as intermediate representation (IR) [11,12], uses the singular value decomposition (SVD) of the Lehmann kernel to construct orthogonal bases for the imaginary time Green's function.
a These authors contributed equally to this work.b e-mail: emin.moghadas@tuwien.ac.at (corresponding author) arXiv:2402.13030v2[cond-mat.str-el]4 Sep 2024 The IR has proven effective in numerous applications, demonstrating success also in scenarios that involve two-particle quantities [8,[13][14][15][16][17][18].The second compression method, which is called the discrete Lehmann representation (DLR), expands the imaginary time Green's function onto a basis set composed of exponentials, such that the corresponding expansion can be conceptualized as a discrete version of the Lehmann representation, utilizing an effective spectral density composed of a sum of δ functions [19,20].Since the IR and DLR schemes are based on a low rank representation of the spectral Lehmann kernel, they are generally restricted to the imaginary time or Matsubara domains [19,20].
Here we aim at a compression method that is agnostic to the underlying representation of the data, but is still able to yield an effective and compact representation.In this context, for the treatment of quantum-field theory correlation functions, the quantics tensor train (QTT) ansatz [21] has been recently proposed.At the same time, the discrete wavelet transform (DWT) emerges as a complementary candidate, as its versatility allows to find compact representations of not only imaginary time or Matsubara data but also quantities defined for real times or frequencies as well as momentum dependent objects.
The wavelet decomposition is a powerful mathematical technique that has gained significant attention within the framework of signal processing and data compression [22][23][24][25][26][27][28][29].Unlike conventional methods such as the Fourier transform, the wavelet decomposition is able to resolve both frequency and location information in a signal, making it particularly suitable at representing complex and dynamic data, while efficiently concentrating information in a few significant coefficients, thus allowing for the removal of less essential details.This adaptability enables wavelets to compress data with minimal loss of important features, and has made them one of the dominant methods for various applications, from image and audio compression to the analysis of time-series data [29][30][31][32][33][34][35][36][37][38].In this paper, we show the versatility and efficiency of the wavelet compression to address the ever-growing demand for more compact and resource-efficient data representations in the many-body theory of interacting systems.
To prove the efficiency of the wavelet decomposition in producing compact representations of many-body quantities, such as generalized susceptibilities and the self-energy of an interacting fermionic system, we focus our attention on a prototypical correlated platform: the Hubbard model at half-filling in its atomic limit (referred to as Hubbard atom in the following) [39].Despite its very simple form, that is reflected in the possibility of generating data from closed analytical expressions at negligible computational cost [40], this model exhibits complex dynamical two-particle correlations, making it the perfect playground to test the validity and effectiveness of our wavelet-based compression scheme.
The paper is structured as follows: in Sec. 2, we first provide an introduction about the most important features of the Hubbard atom and subsequently in Sec. 3, a primer on the mathematical basics of the wavelet transform, detailing how it can be utilized to implement a compression method through the pruning of a fraction of the smallest coefficients in terms of magnitude.We also present and explore two metrics, the compression ratio in memory and the Structural Similarity Index Measure (SSIM), demonstrating their utility in measuring the strength of the compression and giving an introduction into the mathematical foundation behind the SSIM.Finally, in Sec. 4, we will apply the method to the generalized spin and charge susceptibilities as well as the self-energy of the Hubbard atom.In particular, we show that even when pruning a very large fraction of the detail coefficients, the original data can be reconstructed very accurately.This result suggests that only a few features in the underlying data may be critical for describing the physical phenomena of correlated Hubbard atoms.Although specific to this toy model, our work acts as a foundation for further research in this direction.While in fact the focus of this study has been predominantly the methodological development, future work may concentrate on the application of our method for a more in-depth analysis of the underlying physics of correlated electron systems beyond the atomic limit, where the competition of kinetic and potential energies generates a whole richer set of physical phenomena.

Model and formalism
The model we consider in this work, i.e., the Hubbard atom, represents a basic building block of several fundamental Hamiltonians used for the description of many-electron physics.In particular, the Hubbard atom can be viewed as the limit of vanishing electronic hopping of the celebrated Hubbard model [39,41] or, equivalently, of vanishing hybridization in the not less famous Anderson impurity model [42].The atomic limit (AL) of the Hubbard model or Anderson Impurity model describes a single site with an instantaneous electrostatic repulsion U in the presence of an external chemical potential µ.Its explicit expression in second quantization formalism reads: where the density operators of electrons with spin σ =↑, ↓ can be expressed as nσ = ĉ † σ ĉσ in terms of the creation/annihilation operators ĉ † σ /ĉ σ of an electron with spin σ =↑, ↓ on the atomic site considered.Throughout this work, we will set µ = U/2, which corresponds to fixing the average occupation of the system to 1 electron (half-filling), i.e., to a regime where correlation effects are expected to be the strongest.
Due to the Pauli principle, the Hilbert space of the model consists of four states, corresponding to the state without electrons, i.e. |0⟩, to the two states with one electron in one of its two spin configurations | ↑⟩,| ↓⟩, and to a double occupancy of electrons with opposite spins | ↑↓⟩.Evidently, this reduced Hilbert-space allows for an analytical solution and makes it possible to exactly compute several many-electron quantities typically used in the quantum many-body description of correlated electrons.At the same time, the presence of the interaction term U allows to capture, albeit in a simplified fashion, nontrivial aspects of the many-electron physics in strong-coupling regimes, whose theoretical description is often quite hard.In particular, among several studies, where the Hubbard atom has been used as a starting point for more complex analyses [43,44], we should also recall that its solution [40,43,45] turned out to be crucial to investigate the different manifestations of the breakdown of the self-consistent perturbation theory for the many-electron problem, such as the divergences of the two-particle irreducible vertex functions [46][47][48][49][50][51][52][53][54][55] as well as the misleading convergence of bold perturbative series [55][56][57][58][59], and to clarify their intrinsic interconnection [60].
Here we focus on the generalized susceptibilities, which are directly related to the one-(G (1) ) and two-particle (G (2) ) Green's functions defined by G (1)  σ (τ and where the external brackets define the grand-canonical thermal averaging The specific quantities of interest for our work are Matsubara frequency dependent quantities defined by the Fourier transform as follows: where the fermionic Matsubara frequencies ν are defined as (2n + 1) π β , ∀n ∈ Z.We also introduce the electronic self-energy defined as where G (1) 0,σ (ν) = 1 iν is the Green's function of the corresponding noninteracting case (U = µ = 0).As no magnetic field is applied to our system (h = 0), Ĥ is perfectly SU (2)-symmetric in the spin-space, and all one-particle quantities, like Σ σ (ν), are, thus, spin independent.For the half-filled Hubbard atom our convention for the chemical potential given below Eq. 1 implies that the exact expression for the self-energy is We are also interested in the susceptibilities defined from the two-particle correlator as: where ν, ν ′ as well ω are respectively fermionic and bosonic Matsubara frequencies, the latter being defined as 2nπ β , ∀n ∈ Z, and The above definition of the generalized susceptibility in terms of this specific combination of one-and two-particle Green's function is particularly useful because after performing the summation on both fermionic Matsubara frequencies (ν, ν ′ ), it directly yields the physical response of the system to an external perturbation in the framework of the linear response theory.
In particular, the following spin combinations define the generalized charge (c)/spin (s) susceptibilities of our system: where the +(−) sign should be used for the charge (spin) susceptibility, respectively.Hence, the corresponding physical charge and spin susceptibilities can be directly computed as We recall that, in the (particularly relevant) static case (ω = 0), the two response functions yield respectively the isothermal compressibility (χ c ∝ ∂n ∂µ ) and the magnetic susceptibility (χ s ∝ ∂m ∂h ) of the system considered.In the following, we will mostly focus on the static case, which amounts to set ω = 0 in Eq. ( 6), and will henceforth omit the bosonic frequency index.In this way the generalized susceptibilities of the Hubbard atom can be regarded as two-dimensional matrices in the fermionic Matsubara frequency space.

Wavelet transform and methodological details
The wavelet transform is a powerful tool for the analysis of frequency components in signals.Similar to the Fourier transform, the signal will be decomposed in a different basis set in frequency space where it is represented by coefficients obtained from the transform.Unlike the Fourier transform, the wavelet transform makes use of temporally enclosed, wave-like functions as its basis.These wavelet packets can be translated to gain information on the location in time of certain frequencies while scaling (stretching) wavelets allows a resolution in the frequency domain [26,29,30].Since its discovery, the advantages of the wavelet transform over other methods in frequency analysis have made it a popular tool in a wide range of areas such as image processing, data compression, or speech recognition [29][30][31][32][33][34][35][36][37].
The transform computes coefficients by convolving sets of wavelet functions over a signal.Wavelets are wave-like functions ψ ∈ L 2 (R) that share the following properties: We can build a wavelet basis D by scaling and translating such a mother wavelet ψ: where u is used for translating the wavelet and s determines the scale.The transformation of a function f ∈ L 2 (R) is then defined as: with where f ⋆ ψs (u) is the convolution between f and ψs (u).The Fourier transform, which henceforth will be denoted by a circumflex symbol over the quantity of interest, can be used to obtain ψs in the frequency domain From Eq. 10b one can see that ψ(0) = +∞ −∞ ψ(t)dt = 0. Hence the wavelet function ψ can be interpreted as a high pass filter [29].
The inverse wavelet transform of a function f (t) can be defined as where C ψ is given by The above equation is also known as weak admissability condition [29].
We can introduce a scaling function ϕ, also known as father wavelet, which captures the low-frequency components of the signal representing slower changes in the data.In other words this corresponds to an accumulation of wavelet functions for scales s larger than 1.The scaling function is essential for the reconstruction of the original signal.It provides the basis for the low-frequency components or approximation at each level of resolution, while the wavelet functions capture the detail components.The modulus of the Fourier transform of ϕ can be computed using the Fourier transform of the mother wavelet as: from which the scaling function can then be derived via an anti Fourier transform [29].Using Eq. ( 16) and Eq.(10b) we see that showing that the scaling function can be interpreted as a low pass filter.Furthermore, defining the scaling function as allows us, in analogy to Eq. ( 12), obtain the low frequency approximation of f at scale s via a convolution with the scaling function f ⋆ φ * s (u) [29].In order to practically execute computations the discrete wavelet transform (DWT) [26,29,30], is used.Using the DWT, a discrete signal f can be rewritten as follows [61]: Observe, that we have two sets of basis functions {ϕ j0,k [n]} k∈Z and {ψ j,k [n]} (j,k)∈Z 2 where j and k denote the indices associated to scaling and translating the wavelets.The functions ϕ[n], ψ[n] and f [n] are defined on a discrete interval [0, N − 1].As before, the function ϕ corresponds to the scaling function, while ψ refers to the wavelet function.Since these two functions are orthogonal to each other, we can compute the coefficients W i as discrete convolutions of the functions ψ and ϕ over the discretized signal f [61]: As mentioned in Sec. 2, since we are mainly interested in the static (ω = 0) case, we can represent the two-particle quantities as matrices depending on the two fermionic frequencies ν, ν ′ .Consistently, a two-dimensional DWT can be defined by applying the one-dimensional DWT first on the rows of our signal and subsequently on the columns [26,29,30].In this way we replace the original one-dimensional scaling and wavelet functions by their two-dimensional counterparts {ϕ, ψ h , ψ v , ψ d }, where the former is now the two-dimensional scaling function responsible for capturing the low frequency components of the data and the latter three, which denote the two-dimensional wavelet functions, are responsible for the details in the horizontal, vertical and diagonal directions of the signal, respectively.Hence, the corresponding two-dimensional signal f , indexed by two indices m and n, can now be represented by [61] where now an additional index r has been added for translations along the second dimension.As before, the coefficients W can be obtained by convolving the signal f with the scaling and wavelet functions [61]: where i ∈ {h, v, d}.If we consider the specific case of our interest, i.e.. the two-dimensional DWT on the generalized susceptibilities, we can obtain the corresponding explicit expressions simply by making the substitutions f → χ and m, n → ν, ν ′ in Eqs. ( 21) and (22).Evidently, the DWT outlined above can be easily generalized to arbitrary dimensions (e.g. to 3D, if one needs to compress bosonic-frequency dependent susceptibilities) by simply applying the transform consecutively along each dimension of the signal [29].
The wavelet transform can be applied to a signal multiple times recursively.The depth of this recursion is denoted by the level of the decomposition.Each computation of the decomposition at level l + 1 takes the approximation coefficients of the current level l and applies the decomposition to this approximation [24,26,29,30,61].In Fig. 2 the wavelet transform applied to an image for different levels is displayed.Fig. 2b shows a wavelet decomposition of level 1.The top left quarter is the approximation, while the top right, bottom left and bottom right quarters correspond to the horizontal, vertical and diagonal detail coefficients, respectively.Fig. 2c and Fig. 2d illustrate the recursion where the approximations that can be seen in the top left corners of Fig. 2b and Fig. 2c are replaced by, again, another set of approximation and detail coefficients.
In the following, we introduce the key methodological details required for obtaining the results presented in the next section.First, we describe how the wavelet transform is used to compress the data we are working with.Afterwards, we will highlight which metrics can be exploited to quantify the compression.

Compression using wavelet transform
We use the discrete wavelet transform as a means of compression as implemented in the open source package PyWavelets [62].In contrast to other popular methods such as the Principal Component Analysis (PCA) [63] or the Singular Value Decomposition (SVD) [64] we rely on a wavelet basis that is agnostic to the input data.
The compression method we chose is rather simple and brute force.For a decomposition of level l, we leave the highest level approximation as it is.All the detail coefficients will be truncated based on their position in the distribution of detail coefficients.More specifically, we compute the q-quantile d q over the distribution of absolute values of the detail coefficients and set all detail coefficients d where |d| ≤ d q to 0. The reconstruction of the sparsified coefficients then yields the compressed reconstruction of the original data.The proposed procedure is illustrated in Fig. 3.

Quantification of compression
In order to compare the effect of the compression on the physical quantities as well as to assess the compression strength with varying degrees of decomposition levels and threshold quantiles, we first need to introduce a reliable measure.For our purposes, we choose two different metrics for the quantification of the wavelet-based compression.
The first method we use, permits us to measure the efficiency of our compression algorithm in terms of memory The proposed compression algorithm using the wavelet transform.After the wavelet decomposition of the original data, the absolute high-frequency detail coefficients are thresholded by a given q-quantile: all coefficients whose magnitude is smaller than the quantile will be set to zero.The wavelet reconstruction yields then the compressed data.Figure taken from [65].
requirements.The second metric, which is based on established methods in image processing, is used to quantify to which degree the compressed image deviates from the original one.

Memory usage
The original vertices can be represented as dense matrices.Almost each element of the matrix is nonzero.In memory, this dense matrix is stored as a contiguous block of numeric values, each corresponding to one entry in the matrix.If we have an n × m matrix, we occupy nm blocks of the size of the datatype that stores the value of each entry.In our case, the values are represented in double-precision floating point format occupying 8 bytes.When compressing the matrix using the wavelet transform, we expect the majority of entries to become zero transforming the dense matrix into a rather sparse one, depending on the strength of the compression.Since most of the elements then share the same value of 0, the matrix can be stored in a more memory efficient format.One of the most simple formats to achieve this efficiency is the COO format as described in [66], where instead of storing a numeric value for each entry, we store tuples only for those values that are nonzero.To characterize a nonzero entry and reconstruct a dense matrix representation, we require three values for each nonzero entry: the row and the column in the original representation and the value itself.
To quantify the level of compression, we simply compute a compression factor from the space taken up in memory.The size in memory of the original matrix in dense matrix representation m o and memory size of the reconstructed matrix in sparse matrix format m r form a compression ratio mr mo .In the remainder of this work, however, we are interested in the percentage of compression, that the algorithm achieves.This relative compression ratio can be expressed via Note, that the representation of rather dense matrices in sparse matrix format might be larger than in the dense matrix format due to overhead in the data structure.

Structural similarity index measure (SSIM)
The sheer memory size of the vertices is essentially a measure of sparsity in the matrix.To get better insight into the compression, we adopt the Structural Similarity Index Measure (SSIM) into our analysis.The SSIM, which was first introduced by Wang et al. in Ref. [67], is a measure for the similarity of two images and is widely used in image processing.[68,69] The SSIM is defined as where x and y refer to the coordinates in the input images A and B. The quantities l, c and s denote the so-called luminance, contrast and structure, respectively.They are computed using an 11 × 11 pixel Gaussian filter around x and y with a standard deviation of σ = 1.5 [67].
The SSIM essentially splits the task of comparing the similarity of two images into three subtasks: (i) a comparison of the luminance, (ii) a comparison of the contrast, and (iii) a comparison of the structure of the images before combining their results to obtain the final index [67].The luminance is defined as where µ A and µ B are the mean values computed on the input images and are functions of x and y.Similarly, we define σ A and σ B as the variances on the input images to compute the contrast c as Finally, the structure s is defined in a very similar way as where σ AB now represents the covariance.For simplicity we set α = β = γ = 1 and C 3 = C2 2 to obtain the SSIM [67,69,70]: Here, L corresponds to the dynamic range values for each pixel.For an 8-bit image we have L = 255.C 1 and C 2 are then defined as [67] K 1 and K 2 were set to K 1 = 0.01 and K 2 = 0.03.

Quantification of the compression
We begin our discussion with the analysis of the compression strength of the DWT applied to the generalized spin and charge susceptibilities of the Hubbard atom.To this end, we will compress the atomic limit vertices using various decomposition levels and threshold quantiles and afterwards compare the reconstructed quantities to the original ones.Using the two metrics discussed in Sec.3.2, namely the relative compression ratio and the SSIM, we analyse the performance of the DWT in terms of sheer memory requirements as well as reconstruction quality.The two-particle Figure 6: Illustration of the unfolding procedure for constructing a perfectly symmetric reconstructed vertex.We start by extracting the green rectangle from the reconstructed matrix and unfold it using the relations from Eq. ( 29) quantities of the atomic limit were evaluated for β = 1000 and U = 1 where a frequency grid size of 512 × 512 was chosen.The results are presented in Fig. 4. The left panel shows the relative compression ratio, or 1 − m r /m o , as a function of decomposition level and threshold quantile.Evidently, with increasing threshold quantile and decomposition level the memory compression increases.We note, that since we truncate a fixed percentage of detail coefficients during our compression procedure, the rate of compression is independent of temperature as well as of which channel, i.e. spin or charge, we consider.The central and right panels show the SSIM scores for the spin and charge channel respectively again as a function of decomposition level and threshold quantile.We observe a clear-cut difference in performance between the two quantities.In particular, whereas the charge susceptibility displays a nearly perfect score across the board, the reconstruction of the spin susceptibility visibly deteriorates for higher decomposition levels and threshold quantiles.The above calculations were also conducted for lower β values, namely for β = {10, 100}.In these cases, however, both for the spin and charge susceptibilities the wavelet compression achieved an extremely high SSIM score, which was almost indistinguishable from 1.In order to provide a potentially more qualitative insight into the results presented in Fig. 4, we show explicit examples of the wavelet compression for the specific case of the spin channel in Fig. 5.In the left panel we plot the spin susceptibility at U = 1, β = 1000 which are compared to two reconstructed susceptibilities at decomposition levels l = 8 and threshold quantiles q = 0.95 and q = 0.99 which are respectively shown in the upper and lower row of the central panel.In the graphs of the right panel the absolute error between the original and reconstructed quantities is presented.The difference of the scale for the two cases indicates that an increasing threshold quantile or rather a stronger compression will result in a larger reconstruction error, consistent with the SSIM scores displayed in Fig. 4.
We also notice here that differently from other schemes, our DWT method does not suffer from the shift of sharp features in the ν − ν ′ plane (see diagonal structure in Fig. 2a) occurring for a finite bosonic frequency ω.

Preservation of symmetries
In addition to the overall compression performance of the DWT as discussed above, we would also like to address its ability to correctly capture the symmetries of the generalized susceptibilities.In the particular case of half filling considered here, the generalized susceptibilities can be shown to be real, symmetric and centrosymmetric matrices (for further details we refer to Refs.[5,49,51,71,72]), which can be formulated via the following equation In order to analyze whether our wavelet-based compression method is capable of preserving these symmetries, we designed the following experiment.First the generalized susceptibility is compressed and reconstructed.Subsequently, we cut out a piece and use this to construct a perfectly symmetrized susceptibility consistent with Eq. ( 29) as visualized in Fig. 6.This newly obtained perfectly symmetric susceptibility χ rec sym is then compared to the original, reconstructed χ rec to quantify the degree of symmetry violation without considering other reconstruction errors.In particular, we calculate the root mean squared error (RMSE) between these two quantities for the spin and charge susceptibilities with a frequency grid size corresponding to 512 × 512 at a decomposition level of l = 8 and a threshold quantile of q = 0.99, which for these parameters corresponds to the maximal possible compression.In addition, as a direct comparison we also perform a similar symmetry analysis with two-particle quantities obtained via state-of-the-art quantum Monte Carlo (QMC) computations [73], where we compare the original Monte Carlo vertex with a symmetrized one.The motivation behind this comparison lies in the fact, that QMC provides a numerically exact solution of the AIM.Hence, the degree of symmetry breaking can be controlled by choosing sufficiently high sampling statistics.Therefore, we can utilize high-quality QMC data as a benchmark to estimate to what degree symmetry violations occur in our wavelet-based compression scheme.To this end, we utilize data sets of the static, spin and charge susceptibilities of the Anderson impurity model (AIM) at different temperatures, namely β = {10, 60, 100, 300}, corresponding to the quantities from Fig. 1 in Ref. [52] or equivalently from Ref. [74].The resulting RMSE values are presented in Fig. 7.We immediately observe a similar if not smaller RMSE for the wavelet-compressed susceptibilities compared to QMC, where interestingly for the latter the symmetries of the spin susceptibility seem to be less violated.In the case of the DWT on the other hand, consistent with our findings in Fig. 4, the symmetry preservation of the charge as compared to the spin channel appears to be slightly better.In order to obtain a qualitative understanding of this symmetry analysis, Fig. 8 provides an illustrative example of the proposed procedure.The left and central panels respectively show the reconstructed and perfectly symmetric susceptibilities for β = 1000 and U = 1 at the same compression parameters as used for the results in Fig. 7, whereas in the right panel the absolute difference between these two is plotted.One can indeed observe that apart from the triangular piece that was cut out to symmetrize the susceptibility the error is merely numerical noise.Lastly, it is noteworthy that other bases than the Haar wavelet prove to be significantly less suitable for correctly capturing the symmetries of the half-filled atomic limit susceptibilities (see Appendix B).

Compression of the physical susceptibilities
In the current section, we focus on characteristic physical properties or observables associated to the atomic limit susceptibilities and how these are affected by our proposed wavelet-based compression algorithm.As a first step, we calculate the physical response functions χ s/c from the generalized spin and charge susceptibilities.These are then compared to the physical susceptibilities obtained from the respective compressed quantities.We are particularly interested in the temperature behaviour of T •χ s/c .Therefore, we perform the wavelet compression for various threshold quantiles to the generalized spin and charge susceptibilities at different temperature values.In this case we chose the maximal possible decomposition level of l = 10 for a frequency grid size of 2048 × 2048.The corresponding results are shown in Fig. 9.The upper panel shows the comparison between the true, uncompressed and the compressed quantities, whereas in the lower two panels the absolute errors for the spin and charge channel are displayed.We immediately observe a perfect reconstruction of the physical response functions even for the maximal threshold quantile of q = 1, where essentially only the final, single approximation coefficient is kept.The fact, that even after such a drastic truncation a perfect reconstruction can be achieved, can be explained in the following way.The scaling function of the Haar wavelet, responsible for capturing the low frequency dynamics of the signal, computes the average between two neighbouring pixels or data points.Recursively applying this scaling function to the signal at every level and subsequently reconstructing it after truncating all detail coefficients will hence result in a new signal where each point is given by the mean value of the original signal.Thus, performing the Matsubara summation to obtain the physical susceptibilities will then, even for a maximal threshold of q = 1, lead to the same result if one would perform the summation over the original uncompressed data.We also want to stress, on an unrelated note, that the failure of the spin susceptibility to converge to a value of 0.5 in the zero temperature limit is due to the fact, that we perform the Matsubara summation over a finite frequency box and do not consider the appropriate asymptotics of the generalized susceptibilities.

Eigenvalues and eigenvectors of generalized susceptibilities
As mentioned in the Introduction, in spite of its simple nature, the atomic limit of the Hubbard model is capable of capturing non trivial features of strongly correlated electron systems in the strong coupling limit.An arguably important aspect is the breakdown of the self-consistent perturbation theory, as well as its formal [46,56] and physical [50,75] manifestations.We recall that such a breakdown is signalled by the divergences of the two-particle irreducible vertex functions, which are related to the generalized susceptibilities via the Bethe-Salpeter equation (BSE) [5,40,71,72].Specifically, for certain parameter sets, some eigenvalues of the generalized susceptibility in the charge channel vanish making its matrix representation singular and the corresponding BSE non invertible.The first vanishing eigenvalue, associated to the first divergence of the irreducible vertex, marks thus the onset of the non-perturbative regime.Evidently, the eigenvalues as well as the corresponding eigenvectors' structure encode [50,60,[76][77][78][79] significant information for the physics of strongly correlated many-electron systems.For this reason, in the following numerical experiments, we examine the effect of the DWT on the eigenvalue and eigenvector structure of the generalized susceptibilities.To begin with, we will analyze the capabilities of the proposed wavelet-based compression to locate the exact point in the Hubbard atom's U − T phase-diagram where the irreducible vertex function in the charge channel diverges or, equivalently, where the eigenvalues of the generalized charge susceptibility become zero [40,48].Fig. 10 plots the number of negative eigenvalues of the charge susceptibility as a function of U and T , where for the sake of simplicity and illustrative purposes we chose a very coarse grid.The emerging transitions, hence mark the divergence lines where the eigenvalues of the charge susceptibility vanish.In order to investigate whether the DWT is able to determine the location of the vertex divergence, we fix the temperature to β = 1000 and find via binary search the interaction value where the lowest eigenvalue becomes zero.The eigenvalues are computed for the true, uncompressed as well as the compressed charge susceptibility at various threshold quantiles and a decomposition level of l = 6.In Tab. 1 we report the resulting values of U for which the first eigenvalue vanishes, where the error indicates the precision we set for the binary search algorithm.Additionally, for the sake of completeness, we also provide the relative compression as well as the SSIM score.Unsurprisingly, in this case, where we consider the charge susceptibility which showed very high and robust compression metrics in Sec.4.1.1,we observe no deviation from the true critical U -value up to our set precision.
Since, the compression in the charge channel works almost perfectly, it is perhaps more instructive to consider the spin susceptibility as well which for sufficiently low temperatures leads to a significant deterioration of the SSIM scores.Therefore, in a similar fashion as before, we will now investigate the eigenvalues of the generalized spin susceptibility at β = 1000 and U = 1.Specifically, we will monitor how the maximal positive eigenvalue, which is, to a first approximation, responsible for the Curie behaviour of the static spin response of the system [76], behaves upon compression with varying threshold quantiles.Here again, we have chosen a decomposition level of l = 6 and a frequency grid with the size 512 × 512.The results are presented in Tab. 2, where in addition to the maximal eigenvalue λ max , we also report the eigenvalue scaled by its corresponding spectral weight w max , which is defined as As before, we also show the relative compression ratio as well as the SSIM score.In the case of the spin channel, we now observe a slight deviation of the quantities under consideration from the true, uncompressed values.Nonetheless, up to a relative error in the order of magnitude of O(10 −4 ), λ max and w max are hardly affected by the compression, despite the relatively poor SSIM score.
Furthermore, to also obtain an insight into the structural behaviour of the generalized susceptibilities with varying compression strengths, we have also investigated in detail how the associated eigenvectors for the parameter sets in the two cases discussed above are affected by the DWT.In particular, we looked (i) at the eigenvector corresponding to the lowest eigenvalue of the charge susceptibility at U = 0.003628, where said eigenvalue becomes zero, as well as (ii) the eigenvector of the largest eigenvalue of the spin susceptibility at U = 1, because the first is associated to the breakdown of the perturbation expansion and the second controls the Curie magnetic response of the system.In both cases we fixed the temperature to β = 1000.We compare the eigenvectors which we obtained from the compressed susceptibilities to the original ones for various threshold quantiles and a decomposition level of l = 6, whereas the frequency grid size was again set to 512 × 512.The results are shown in Fig. 11, where the original and compressed eigenvectors of the spin and charge susceptibilities are displayed in the plots of the upper and lower row respectively.In contrast to the charge channel, whose lowest eigenvector, which is fully localized in frequency, is always perfectly reconstructed, we can observe an increasingly coarser reconstruction of the broad structure of the largest eigenvector of the spin susceptibility.This different behavior can be qualitatively explained upon inspecting the generalized susceptibilities in Matsubara frequency space more closely.In Fig. 12 we show, for the same parameter sets as before, the original spin and charge susceptibilities and their compressed counterparts, where the threshold quantile is now set to q = 0.99.One can easily note, then, that the generalized spin susceptibility is much broader distributed in frequency space than its charge counterpart, which exhibits a frequency behaviour almost reminiscent of a delta-function.These characteristics are also very distinctively manifested in the frequency structure of the respective eigenvectors we have considered.As a consequence, the wavelet transform is not able to capture all the information contained in the spin susceptibility with fewer wavelet coefficients as opposed to the charge channel, which due to its sharper features requires far less coefficients.This leads to the coarse-graining effect observed in the reconstructed eigenvectors (upper row in Fig. 11) as well as the compressed spin susceptibility (upper right panel in Fig. 12) when too many detail coefficients are truncated.We also want to highlight, that this aspect gives a reasonable explanation as to why in Fig. 4 a significantly smaller SSIM score was observed for the spin susceptibility at higher β values.Additionally, the compression metric for the memory and the SSIM are reported.The divergence in the original, uncompressed vertex is also found to be at U = 0.003628 ± 10 −7 .

4.2.3
Schwinger-Dyson-Equation with the compressed G (2)   The final set of numerical experiments will be devoted to the Schwinger-Dyson equation (SDE), also referred to as the equation of motion for the one-particle Green's function.Our goal is to analyze the effect of the wavelet-based compression of the two-particle Green's function, the central component of the SDE, on the resulting self-energy  obtained from the equation of motion.We note, that in this case the Green's function depends on three frequencies, which will require a three-dimensional DWT.The SDE for the local self-energy is given by where F ↑↓ denotes the corresponding spin components of the full vertex function [71].The above equation can be equivalently formulated in terms of the two-particle Green's function [71] Σ In order to assess the effect of the DWT on the self-energy, we will extract it from the SDE using the analytic two-particle Green's function of the Hubbard atom and compare the results with the self-energy obtained from a   11, where the threshold quantile is set to q = 0.99.wavelet-compressed G (2) .The compression will be done at various threshold quantiles and a maximal decomposition level for a frequency grid of 64 × 64 × 64.The original and reconstructed self-energies for various temperature values are reported in the lower panel of Fig. 13.In this case, we employed a three-dimensional DWT to compress the entire two-particle Green's function before plugging it back into the SDE.We generally observe larger reconstruction errors for smaller values of β, this is especially noticeable for larger threshold quantiles such as q = 1.Upon lowering the temperature, these errors are also further reduced.This, at first glance, anomalous behaviour can be explained by the analytic properties of G (2), νν ′ ω ↑↓ .The two-particle Green's function under consideration follows a Curie-Weiss law and shows a 1/T divergence along the main diagonal ν = ν ′ for β → ∞ (for the analytic expression of F νν ′ ω ↑↓ we refer to Eq. (2.223) in Ref. [71]).This leads to very sharp features along the main diagonal compared to the background ν ̸ = ν ′ , which, as we discussed before in Sec.4.2.2 when we considered the eigenvectors of the generalized susceptibilities, are significantly easier for the DWT to capture accurately with fewer coefficients.
Evidently, exploiting the knowledge we gained from compressing the physical susceptibilities in Sec.4.2.1, it is also possible to perform the compression for two-dimensional slices along the ν-frequency axis.This way, as subsequently in the SDE a Matsubara summation in the ν ′ and ω frequency direction will be performed, we can rest assured that the wavelet-based compression will yield exact results.In the upper panels of Fig. 13 the reconstructed self-energies are compressed in this manner.And indeed, we observe a perfect reconstruction even for a threshold quantile of q = 1, consistent with our findings from Fig. 9. Needless to say, this way of performing the wavelet compression, might pose Original Σ(ν) q = 0.95 q = 0.99 q = 1 Figure 13: The self-energy computed from the SDE with the original, uncompressed two-particle Green's function G (2) (ν, ν ′ , ω) compared to the one obtained from a compressed two-particle Green's function at various threshold quantiles, for different temperatures.Upper panel: The compression is performed on (ν ′ , ω)-slices along the ν-frequency range, where no subsequent Matsubara summation is performed in the SDE.Lower panel: Here the compression is performed on the entire two-particle Green's function as a function of three frequencies.
higher memory requirements than simply performing a three-dimensional DWT for the entire two-particle Green's function.

Conclusions and outlook
This study has shown that wavelet compression methods can efficiently represent generalized susceptibilities and self-energies within the Hubbard model at half-filling in its atomic limit.We have shown that the wavelet method can efficiently address computational challenges associated with storing and manipulating the non-zero temperature correlation functions in systems characterized by strong interactions, and have discussed how to characterize the efficiency of the compression.The results demonstrate that the wavelet method is versatile and resource-efficient.As a concluding remark, we can give a crude estimation on the storage capacities of the wavelet transform when it comes to two-particle quantities often encountered in many-body theory.In the general case, a two-particle correlation function will depend on three frequency and three momentum indices in addition to possibly multiple orbital degrees of freedom.Depending on the dimensionality of the lattice considered for the Hamiltonian the overall memory scaling will be proportional to O(N 4 orb N 3 ν N 3d k ) with N orb orbitals, N ν number of Matsubara frequencies for each of the three frequency indices and N k datapoints for the momentum grids in each lattice dimension.One can immediately see that this scaling will quickly become immensely challenging to handle.For instance, neglecting, for the sake of simplicity, the orbital and momentum dependence and considering only N ν = 1000 for all three Matsubara frequencies, the memory requirements to store N 3 ν datapoints as double precision floating points will already be of the order of 8Gbytes.In Sec. 4 we observed that the compression ratios where the reconstruction quality was still considerably good were between 90 % and 98 %, which would compress the considered two-particle object down to 160 -800Mbytes yielding a significant reduction of memory requirements.
The successful representation of the atomic limit opens the door to the broader usage of signal processing and data compression principles.Future directions include the adaptation of wavelet-based compression techniques to explore non-zero doping conditions, moving away from the symmetry constraints of half filling.Additionally, future studies will aim to extend its applicability to more realistic scenarios, such as the investigation of two-particle Green's functions within the Anderson impurity model or the Hubbard model at intermediate or weak coupling.A distinctive feature of the wavelet-based approach, distinguishing it from the previously mentioned IR and DLR compression methods, is its capacity for a data-agnostic application.This characteristic allows for its use in both momentum-dependent and real-frequencies-dependent Green's functions.In this respect, particularly intriguing will be its potential application in addressing interacting systems under nonequilibrium conditions, thereby broadening the scope of research possibilities in the field.
We observe that a seamless integration of wavelet methods into the existing workflow of numerical methods for the many-body problem requires the application of fundamental mathematical operations inherent to many-body physics -specifically, convolutions and inversions -within the wavelet basis [80][81][82].Achieving this understanding represents a critical step to establish the wavelet-based framework as an integral component of numerical methodologies for tackling complex many-body systems.
From a pure data science perspective, research directions may extend toward broadening the applicability of our proposed compression scheme to higher-dimensional wavelet functions, ones that go beyond mere products of one-dimensional functions and instead possess intrinsic higher-dimensional characteristics.Pursuing this research trajectory holds the potential to provide even more effective compression techniques, enhancing the overall efficiency.An alternative approach involves leveraging machine learning techniques to obtain optimal performance in the compression process.This can be achieved by incorporating learnable thresholding functions and employing dictionary learning to determine the set of basis functions that most effectively represent the underlying data [83][84][85].
In conclusion, this work serves as a pioneering proof-of-principle application of wavelet compression in the field of many-body physics, opening new possibilities for future research efforts and providing a valuable contribution to the ongoing evolution of quantum many-body physics methodologies.
being the corresponding partition function and β the inverse temperature of the system, ĉ( †)σ (τ ) = e τ Ĥc ( †)σ e −τ Ĥ, T τ denotes Wick's time-ordering operator in imaginary times τ and we have exploited the time-translational invariance of the problem considered, by setting the time argument of the last operator in the Green's function definitions equal to zero.

Fig. 1
examples for prominent wavelet functions ψ and the corresponding scaling function ϕ of different wavelets are shown.

Figure 2 :
Figure 2: Decomposition of the generalized spin susceptibility vertex at β = 100, using the Haar wavelet (the coefficients of the wavelet decomposition are normalized for better visibility).

Figure 3 :
Figure3: The proposed compression algorithm using the wavelet transform.After the wavelet decomposition of the original data, the absolute high-frequency detail coefficients are thresholded by a given q-quantile: all coefficients whose magnitude is smaller than the quantile will be set to zero.The wavelet reconstruction yields then the compressed data.Figure taken from[65].

Figure 4 :Figure 5 :
Figure 4: Left panel: Relative compression ratio 1 − m r /m o for the generalized susceptibilities χ ν,ν ′ as a function of decomposition level and threshold quantile.Note, that due to our chosen compression procedure, the compression rate is independent of temperature or which channel, i.e. spin or charge, is chosen.Central panel: SSIM as a function of the decomposition level and threshold quantile for the generalized spin susceptibility at U = 1 and an inverse temperature of β = 1000.Right panel: SSIM for the generalized charge susceptibility at U = 1 and β = 1000.For the calculations here a frequency grid size of 512 × 512 was chosen.For lower values of β the SSIM scores for the spin and charge channel show close to no deviation from 1 independently of threshold quantile or decomposition level.Original s

Figure 7 :
Figure 7: Comparison of the symmetry conservation in the wavelet decomposition with QMC.The RMSE between the reconstructed χ rec and perfectly symmetric, reconstructed χ rec sym is shown for the atomic limit at U = 1 as a function of temperature at a decomposition level of l = 8 and a threshold quantile of q = 0.99.Here we chose a grid size of 512 × 512.The results are compared to the RMSE between the charge susceptibility of the AIM at U = 5.75 and its symmetrized counterpart.The data for the AIM is obtained from Ref. [52] Fig.1 therein.

Figure 8 :
Figure 8: Example of the reconstructed χ rec (left panel) and perfectly symmetric, reconstructed χ rec sym (central panel) for the spin channel at U = 1 and β = 1000 alongside the absolute difference.For the wavelet decomposition a level of l = 8 and a threshold of q = 0.99 is chosen, whereas the frequency grid of the susceptibility has the size 512 × 512.

Figure 9 :Figure 10 :
Figure9: Physical spin and charge susceptibilities χ c/s (T ) as computed from the respective generalized quantities for various threshold quantiles.Here, a decomposition level of l = 10 and a frequency grid with a size of 2048 × 2048 is used.In the lower panels, the absolute errors with respect to the original uncompressed physical susceptibilities are displayed.

Figure 11 :
Figure 11: Upper panels: Eigenvector corresponding to the maximal eigenvalue of the reconstructed generalized spin susceptibility for U = 1 and β = 1000 at different threshold quantiles compared to the eigenvector of the original, uncompressed susceptibility.panels: Eigenvector corresponding to the lowest eigenvalue for the reconstructed charge susceptibility at various threshold quantiles compared to the eigenvector calculated from the uncompressed quantity.For the decomposition level l = 6 and a frequency grid size of 512 × 512 were chosen.Here, the eigenvectors were computed at U = 0.00362759873 and β = 1000, i.e.where the lowest eigenvalue vanishes.

Figure 12 :
Figure12: Zoom on the first 20 Matsubara frequencies of the charge and spin susceptibilities alongside their respective reconstructions from the parameter sets of Fig.11, where the threshold quantile is set to q = 0.99.

Table 1 :
Value of the Hubbard U where the lowest eigenvalue of the reconstructed generalized charge susceptibility χ c at β = 1000 vanishes, for a decomposition level of l = 6 and a frequency size of 512 × 512 for various threshold quantiles.

Table 2 :
Value of the maximal eigenvalue of the generalized charge susceptibility at U = 1 and β = 1000 vanishes, for a decomposition level of l = 6 and a frequency size of 512 × 512 for various threshold quantiles.Additionally, the compression metric for the memory and the SSIM are reported.The magnitude of the maximal eigenvalue of the original, uncompressed susceptibility is λ max = 987605.434and factoring the spectral weight gives λ max w max = 326.211.