Bypassing the computational bottleneck of quantum-embedding theories for strong electron correlations with machine learning

A cardinal obstacle to performing quantum-mechanical simulations of strongly-correlated matter is that, with the theoretical tools presently available, sufficiently-accurate computations are often too expensive to be ever feasible. Here we design a computational framework combining quantum-embedding (QE) methods with machine learning. This allows us to bypass altogether the most computationally-expensive components of QE algorithms, making their overall cost comparable to bare Density Functional Theory (DFT). We perform benchmark calculations of a series of actinide systems, where our method describes accurately the correlation effects, reducing by orders of magnitude the computational cost. We argue that, by producing a larger-scale set of training data, it will be possible to apply our method to systems with arbitrary stoichiometries and crystal structures, paving the way to virtually infinite applications in condensed matter physics, chemistry and materials science.


I. INTRODUCTION
The atomic energy scales emerging in strongly correlated systems [1][2][3] can induce a broad spectrum of spectacular effects, ranging from arresting the electronic motion [1] to causing high-temperature superconductivity [4], unlocking access to new topological phases and influencing dramatically the potential-energy surfaces (PES) of molecules and solids [5][6][7][8][9][10]. Therefore, the need and the potential effects for science and society of extending to strongly correlated systems the computational materials-by-design paradigm can hardly be overstated [2]. The substantial progress achieved in the past decade in calculating the electronic structure of strongly correlated materials largely owes to the idea of combining mean-field (MF) theories, such as approximations to DFT [11][12][13][14][15][16][17] with QE [2,18,19] theoretical frameworks. Well-known examples are Dynamical Mean Field Theory (DMFT) [20][21][22][23][24][25][26][27] and Density Matrix Embedding Theory (DMET) [28,29]. As shown in Ref. [30], also the multi-orbital Gutzwiller approximation (GA) [30][31][32][33][34], which is a variational framework (equivalent to the Rotationally Invariant Slave Boson (RISB) [35][36][37] at the MF level [38,39]), can be formulated as a QE scheme featuring recursive ground-state calculations of impurity models with a finite bath called "Embedding Hamiltonians" (EH). Therefore, even if the principles underlying The fundamental idea underlying all QE theoretical frameworks consists in replacing the original (typically unfeasible) problem of directly simulating these systems with the more manageable task of solving equations for a series of EHs, composed by fractions of the material (impurities) and effective-medium degrees of freedom (self-consistently determined for describing the interaction of the impurities with their environment). The current state-of-the-art approach to tackle QE simulations is based on solving the EH equations recursively utilizing many-body techniques [42][43][44][45]. On the other hand, due to the quantum-mechanical nature of the EH, its solution ultimately has a computational cost that grows exponentially with the number of impurity degrees of freedom. Because of this reason, the practical application of these tools to complex materials is often too computationally demanding to be ever feasible.
Here we show that this problem can be efficiently tackled from a completely different perspective: capitalizing on the fact that the form of the EH is universal (i.e., it does not depend on the specific stoichiometry and crystal structure of the material considered), we bypass altogether the computationally-expensive recursive solution of the EH by "training a machine" to solve this problem once and for all. To accomplish this goal, we develop a computational framework combining machine-learning (ML) techniques, such as Kernel Ridge Regression (KRR), with a mathematical method named Standard approach (left): each iteration requires to solve numerically the EH and calculate the observables F (X) for different descriptors X. Proposed approach (right): a ML algorithm, previously trained to learn the universal function F (X), allows us to bypass the computation of the EH. "n-mode representation" [46,47] -previously used for effectively reducing the dimensionality of large-scale regression problems (e.g., for reducing the number of points required for constructing high-dimensional PESs in quantum chemistry [48,49]).
Note the fundamentally-different nature of our method -which employs ML inside the solution of the full quantum problem,-with respect to the many current uses of ML for learning pre-existing solutions [50] (e.g., for applications to different materials and structures).
We illustrate the power of our method by performing benchmark calculations of a series of actinides. Utilizing our method, we were able to calculate -at a computational cost comparable to bare DFT-the discontinuous behavior of the equilibrium volumes of the actindes as a function of their atomic number Z (actinide transition) [51], which is a phenomenon originated by a complex interplay between structural degrees of freedom, relativistic effects, atom-and orbital-selective electron correlations [30,[52][53][54][55].

II. QE ALGORITHMIC STRUCTURE
The problem of applying QE methodologies (such as the DFT+E methods) to realistic solids and molecules ultimately reduces to solving recursively multi-orbital Hamiltonians represented as follows: where k is the momentum conjugate to the unit-cell label R, the electronic shells of the atoms within the unit cell are labeled by i, j = 1, .., η, and the corresponding spinorbitals are labeled by α = 1, .., M i , β = 1, .., M j . For later convenience, with no loss of generality, we assume that the first term is non-local (i.e, that k k,ii = 0 ∀ i) and thatĤ loc includes both the one-body and the two-body local parts ofĤ: where E i describe the on-site energies (such as the crystal-field energies and the spin-orbit coupling (SOC)) andĤ int Ri depends on the Slater-Condon parameters [56], i.e., the Hubbard interaction strength U i and the Hund's coupling constant J i .
The basic algorithmic structure of all QE methods to solve the Hamiltonian [Eq. (1)] is schematically illustrated in the left side of Fig. 1. A series of EH, represented as: are constructed for describing the coupling of the impurity with its environment in a MF fashion. HereB Ri (x i ) represents an effective medium coupled with the subsystem (impurity) [Eq. (2)], which is encoded in a series of parameters x i . Determining the self-consistent parameters x i requires to calculate multiple times a series of quantities (varying for different QE methodologies) for the Hamiltonian in Eq. (3), that we schematically represent as F i .
For concreteness, here we focus on the GA. As shown in Refs. [30,37] (see also the supplemental material), this method can be regarded as a QE framework where: matrices, the latin labels a, b correspond to the bath degrees of freedom f and the output function F i is the single-particle density matrix: where the labels A, B = 1, .., 2M i run over both the impurity and the bath degrees of freedom. Therefore, consistently with the general algorithmic structure schematically represented in the left side of Fig. 1, solving the GA equations requires to evaluate recursively Eq. (5) as a function of the EH descriptors: To simplify the notation, from now on we will omit the electronic-shell label i.

Computational complexity of the EH problem
In GA ab-initio calculations it is typically necessary to deal with impurities consisting of M = 10 degrees of freedom (for d-electron shells) or M = 14 degrees of freedom (for f-electron shells). Since the bath of the EH has the same number of degrees of freedom of the impurity, see Eq. (4), the dimension of the EH space is D = 2 2M .
Note that the dimension of the EH system scales as D = 2 2M also in DMET. In fact, the differences between these 2 methods stem exclusively from their different selfconsistency conditions [30,40,41]. Within the ghost GA framework (g-GA), which is a more accurate extension of the GA [34], the number of effective-medium degrees of freedom is still finite, but larger than bare GA. Therefore, the dimension D of the EH system is exponentially higher. Finally, in DMFT [20] the number of effectivemedium degrees of freedom (and, therefore, the EH dimension D) is infinite.
In all of the theoretical tools mentioned above, the computational bottleneck is solving recursively the EH equations. In fact, this is the only reason why the cost of QE methods generally exceeds by orders of magnitude the cost of mean-field approaches, such as classic approximations to DFT. The computational framework described in the next section will allow us to bypass altogether this problem.

III. COMBINING ML WITH THE n-MODE REPRESENTATION
Rather than trying to develop more efficient computational tools for solving the EH equations, in this work we will capitalize on the universality of the function [Eq. (5)], utilizing KRR and the n-mode expansion for learning it once and for all, see the right side of Fig. 1.
The strategy of utilizing ML for bypassing expensive calculations of universal maps, often referred to as "surrogate modeling," is widely used in physics, chemistry and materials science [57][58][59][60][61][62][63]. The main obstacle to applying classic ML algorithms (such as KRR) for learning multi-variable functions is that it requires a number of training data points that scales as: where d is the number of input variables and m is the number of mesh subdivisions for each dimension. This problem is often referred to as the "exponential curse". Without use of symmetry, the number of input variables in F (X) in Eq. (5) is d = 1 + 4M 2 , which is 401 for d-electron shells and 785 for f-electron shells. Therefore, direct applications of ML methods would be extremely costly for learning this function.
To overcome this problem, here we combine KRR with the "n-mode representation," which is a technique previously explored in different contexts (and under different names), e.g., for reducing the number of points required to construct high-dimensional PESs [46,48,49,64] and for facilitating the solution of the Schrödinger equation in quantum chemical methods [65,66]. The basic idea underlying the n-mode representation is to construct approximations to the high-dimensional function F (X) in Figure 2. Representation of the proposed nKRR approach. The universal EH functions F (X) is approximated with the n-mode representation up to the desired order n (e.g., n = 2 in the picture). The KRR method is used to fit the corresponding lower-dimensional cut-functions, which are recombined into an approximation to F (X).
terms of the so-called "cut-functions," such as: which are the restrictions of F (X) to hyperplanes where subsets of the components of X are set to 0. At a given order n of the expansion, F (X) is approximated utilizing only cut functions of up to n variables. The exact function is recovered when n equals the total number of variables d, and the series converges very rapidly as a function of n in many relevant cases [46,48,49,[64][65][66]. Within our context of application, the main consequence of the n-mode expansion is that, since the effective dimensionality is limited to that of the needed order n, the input-output mapping can be determined with a number of data points that scales only as: i.e., it scales polynomially as a function of d, rather than exponentially (Eq. (9)). Specifically, the hereby proposed nKRR methodology consists of the following steps: • Breaking down the universal functions F (X) in lower-dimensional cut-functions, using the n-mode expansion.
• Learning the corresponding lower-dimensional cut functions, up to the desired order, using KRR.
• Combining the cut-functions into the desired nmode approximation.
A schematic representation of this approach is shown in Fig. 2. For completeness, a short introduction to the KRR method and the n-mode representation is provided in the supplemental material.
We point out that, as opposed to other dimensionalityreduction techniques (where the number of input variables is decreased), the nKRR method allows us to take into account from the outset all descriptors of the EH system -in a manner such that the effective dimensionality is substantially reduced. In the next section we will also capitalize on general physical arguments inherent in the specific structure of Eq. (4). This will allow us to derive a convenient parametrization of the EH in DFT+GA calculations, speeding up dramatically the convergence of the n-mode representation.
We want to point out that, besides the GA, the nKRR methodology described above could as well be implemented in combination with DMET or more accurate QE methods, such as the g-GA [34] and DMFT.

IV. APPLICATION TO ACTINIDE SYSTEMS
Here we describe in detail our implementation of the nKRR method for actinide systems. For simplicity, we will focus on the case of a generic EH consisting of felectron shells in an isotropic medium -which is typically a good approximation for actinide systems, where the dominant role of the SOC allows us to average over the crystal-field splittings. Under these assumptions, using group-theoretical considerations, it can be shown [37,67] that the 14 × 14 matrices E, D, λ c are diagonal and fully determined by their respective j = 5/2 and j = 7/2 components E j , D j , λ c j , where j is the label of the total angular momentum for an f-electron shell. Furthermore, as discussed in the supplemental material, the conservation of the total number of electrons implies that the trace of the single-particle density matrix of the EH is M . Therefore, the only independent descriptors of the EH are the interaction parameters U, J and the following variables: In fact, with no loss of generality, we can set j∈{5/7,7/2} (E j − λ c j ) = 0, as changing this variable corresponds to applying a chemical-potential shift in Eq. (4) (which would be redundant, as the number of electrons M in the EH is fixed).
Furthermore, the behavior of F (X) (Eq. (5)) is fully determined by the following functions: where j z is the quantum label of the third component of the total angular momentum for each j.
Consistently with Ref. [30], here we set the screened Hubbard interaction U = 4.5 eV and the Hund's coupling constant J = 0.36 eV. Therefore, the only free embedding descriptors are X 1 , .., X 5 .

Parametrization of F (X) and training data set
In this section we describe in detail our procedure for setting up the nKRR method, that we are going to utilize for performing DFT+GA calculations of systems involving Pa, U, Np, Pu and Am.
As pointed out in the supplemental material, the speed of convergence of the n-mode representation can be improved by a suitable change of variables, as (by construction) the accuracy of the approximation tends to be higher in the proximity of the domain of the cut functions utilized at the chosen order of truncation. Furthermore, it is pivotal to ensure that the grid of training data points utilized for learning the cut functions [Eq. (10)] is sufficiently large, as ML methods can be predictive only within the training-data range.
It is particularly convenient for our applications to exploit the possibility of expressing the function of Eq. (13) in terms of shifted variables: For setting up the components ofX and the training data set we utilized the following procedure, which is based on physical considerations inherent in the properties of the atomic impurities of interest.
To determineX 1 , we have pre-calculated the behavior of F 1 as a function of X 1 , at fixed X j = 0 ∀ j = 2, .., 5. Note that, in this limit, the impurityĤ loc is isolated from the bath. Therefore, the values of F 1 are quantized and correspond to the nominal (integer) f-electron occupations along the actinide series. For each actinide, we have setX 1 as the middle point of the interval of X 1 values such that F 1 (X 1 , 0, .., 0) equals the corresponding nominal occupation, see Fig. 3. We note that X 2 describes the impurity SOC, which is essentially an atomic property, i.e., it is typically almost independent of the environment in DFT and DFT+GA calculations. Therefore, for each actinide we have setX 2 based on the to the nominal atomic values, which we pre-calculated using LDA. Since X 3 , X 4 , X 5 describe the EH bath and its coupling with the impurity, their range is generally system dependent. Therefore, the choice ofX 3 ,X 4 ,X 5 is essentially arbitrary. In our calculations we have set them based on a single DFT+GA calculation of δ-Pu at its experimental equilibrium volume.
The range of the training data set was estimated by performing LDA+GA calculations of δ-Pu at ±35% of its experimental equilibrium volume. This choice proved to be sufficient for performing all calculations performed in this work. Note that our implementation interactively queries the user if EH parameters beyond the training range are explored in a calculation. Whenever this happens, new training data can be generated and stored in a database. This type of iterative supervised learning -which is often called "active learning procedure,"-allows one to assess the validity of the simulations and to extend systematically the range of applicability of the nKRR algorithm.
The numerical values of the components ofX and the data mesh of the components of Y -obtained with the procedure outlined above-are reported in the supplemental material.

V. BENCHMARK CALCULATIONS
To assess the power of our method, we performed LDA+GA benchmark calculations of different actinide solids, utilizing the nKRR method described above to solve the EH Hamiltonian. We will refer to this framework as the GA+nKRR, while we will call GA+ED the standard DFT+GA approach resulting from using ED as an EH solver.
A particularly interesting property of the actinide series is the anomalous dependence of their equilibrium volumes as a function of the atomic number. In fact, while the equilibrium volume of the lighter actinides varies con-tinuously as a function Z (from Pa to Pu), it displays a pronounced discontinuity between Pu and Am. This volume anomaly is often called actinides transition [51,68], and it is originated by a complex interplay between structural degrees of freedom with SOC, atom-and orbital-selective electron correlations [30,[52][53][54][55]. Therefore, capturing this behavior constitutes a very strict benchmark of our method.
In panel Fig. 4(a) we show the equilibrium volumes of the low-temperature allotropes of Pa, U, Np, Pu and Am. The GA+nKRR (n = 3) and GA+ED calculations (performed at U = 4.5 eV and J = 0.36 eV, as in Ref. [30],) are shown in comparison with LDA and the experiments [52,69]. In Fig. 4(b1-b5) we also show the corresponding energy-volume curves -whose minima are the points shown in Fig. 4(a). Remarkably, for all systems considered, the GA+nKRR method is substantially more accurate than DFT already for n = 2, while it becomes essentially as accurate as GA+ED for n = 3.
We point out that each ED solution of the EH Hamiltonian takes about 10 minutes and requires 10GB of RAM on average, while it takes only about 0.1 seconds and 50MB of RAM within the nKRR framework. Because of this reason, our GA+nKRR method is essentially as expensive as bare DFT. In particular, the cost in terms of computational time and RAM is mainly determined by the EH solver within the GA+ED framework. Instead, within GA+nKRR the computational bottleneck is determined by the DFT operations, such as constructing the Kohn-Sham Hamiltonian and calculating the electron density at each iteration (see the supplemental material). For the calculations performed in this work, where the DFT part was performed using the all-electron scheme implemented in WIEN2k [70], utilizing the nKRR method for solving the EH equations reduced the computational time by a factor of 10 to 100 (depending on the system). The relative computational gain of applying our methodology would presumably be even higher by utilizing less computationally-demanding implementations of DFT.

VI. CONCLUSIONS
In summary, in this work we proposed a new computational framework for simulating strongly-correlated electron systems, which offers the unique possibility of stepping up substantially the accuracy with respect to meanfield theories (such as classic approximations to DFT and DFT+U [56]), at a comparable computational cost. This was accomplished combining QE theoretical frameworks with a fit-for-purpose ML technique (the nKRR), where the learning problem is facilitated by the n-mode expansion. The fact that our method reduces the complexity of the learning problem from exponential to polynomial makes it realistically possible to extend it to systems with arbitrary stoichiometry and crystal structures. Realizing this program will pave the way to virtually infinite ap-   For completeness, here we briefly summarize the equations underlying the formulation of the GA as a QE scheme, which was previously derived in Refs. [1,2].
The GA solution is obtained by calculating the saddle-points of the following Lagrange function [2]: where:Ĥ R, D, λ, λ c , ∆ are complex block-matrices whose respective µ is the chemical potential, N is the total number of electrons in the system (normalized to the number of k-points N ), k,ij are matrices constituted by M i × M j blocks labeled by i, j with entries αβ k,ij and |Ψ 0 is the most general single-particle wavefunction within the space ofĤ qp , see Eq. (2). By construction [1], the "embedding states" |Φ i are assumed to lie within the M i -particle subspace ofĤ emb i , see Eq. (3), i.e., they satisfy the following equation: Physical observables can be calculated from the parameters of the theory realizing the saddle-point of Eq. (1). In particular, the total energy of the system equals the saddle-point value of L N . The expectation values of any local operatorÔ {c † Riα , c Riα } can be calculated as follows [1]: where |Φ i is the ground state of the self-consistent EHĤ emb The local self energy, instead, is expressed in terms of the self-consistent parameters R i and λ i as follows [3]: The most expensive operation necessary for obtaining the saddle-point of Eq. (1) is to compute recursively the ground-state of the EH, see Eq. (3). The nKRR methodology proposed in this work allows us to bypass altogether this operation.

GA+DFT
Within DFT+GA, the Kohn-Sham parameters k,ij and E are updated at each charge iteration and determined selfconsistently. Evaluating these parameters and calculating the electron density at each charge iteration is essentially as expensive as in all classic DFT implementations.
Our LDA and LDA+GA calculations were performed utilizing the DFT code WIEN2k [4]. The LDA+GA solver was implemented following Ref. [2]. The LAPW interface between WIEN2k and the RISB was implemented as described in Ref. [5], utilizing the fully-localized limit (FFL) double-counting functional [6]. All calculations were performed setting RKmax = 9. The convergence with respect to the number of k-points was verified for all systems considered.

II. THE n-MODE REPRESENTATION
As mentioned in in the main text, applying directly KRR for learning a multivariate function F (X 1 , X 2 , . . . , X d ) from data built on a mesh with m data points per axis requires N ∼ m d function evaluations. Therefore, the complexity of the learning problem grows exponentially as a function of the number of variables d.
The basic idea underlying the n-mode representation is to represent a high-dimensional function F in terms of the following lower-dimensional functions: etc.., where: F 0 = F (0, 0, . . . , 0, 0, 0, . . . , 0, 0, 0, . . . , 0) etc.., are the so-called "cut-functions," which are restrictions of F (X) to hyperplanes where subsets of the components of X are set to 0. Specifically, the n-mode representation of F (X) is given by the following equation: It can be readily verified that, when all terms are retained, Eq. (9) is an exact identity (which is a major advantage compared to other approximations such as the Taylor expansion). Truncating this series up to a given order n < d provides us with an approximation that is exact only over the domains of the order-n cut-functions, while elsewhere it is an approximation that tends to be more accurate in the proximity of the domains of the order-n cut functions.
By inspecting Eq. (9) we note that, building a mesh of training data with m points per axis, the nKRR method up to order n requires only the following number of data points: which scales polynomially (as d n ) as a function of d, rather than exponentially. Note that, as opposed to other dimensionality-reduction methods -where the number of effective input variables is decreased,-the n-mode expansion includes from the outset all variables.

Variable shifts
Whenever it is possible to estimate the range of input values where F has to be evaluated for a particular application (e.g., based on physical arguments inherent in the particular context of application), performing a suitable change of variables X = v(Y ) can reduce significantly the necessary truncation order n. In fact, a change of variables can be designed in such a way that the relevant range of input values is as close as possible to the domains of the cut functions of G(Y ) = F (v(Y )), where the n-mode expansion is more accurate (by construction).
In particular, in the main text we have exploited this freedom to shift the origin with a change of variables of the form: which improved considerably the speed of convergence of the n-mode expansion, facilitating the learning problem.
The numerical values of the components ofX 1 andX 2 are reported in Table I, whileX 3 ,X 4 ,X 5 have been all set to −0.694 eV (based on a single DFT+GA calculation of δ-Pu at its experimental equilibrium volume). While in this work we have applied the n-mode representation within the context of quantum embedding methods, this method has been originally designed and explored in different contexts (and under different names). In particular, as mentioned in the main text, the n-mode representation has recently gained significant attention for representing and computing potential energy surfaces (PES) of molecules [7][8][9][10].
Also the so-called "incremental method" [11] and "many-body expansion" for computing electronic energies of molecules can be considered as a variations to the n-mode expansions [10]. These tools have both been exploited for obtaining and representing accurate electronic energies of large molecules and molecular clusters [12]. Furthermore, the incremental method has also been applied previously with a focus on strong electron correlation [13,14]. Specifically, these works applied the incremental expansion directly to the solution of the Schrödinger equation for solving specific chemical problems. The incremental methods for single point electronic energies has been also combined with the PES n-mode expansion to obtain a double incremental expansion of the PES paving the way for obtaining linear scaling construction of PESs -an otherwise hard-to-imagine result [10]. Finally, the n-mode representation can be seen as one variant of high-dimensional model representation (HDMR) and is sometimes denoted cut-HDMR. In turn, HDMR is closely related to the ANOVA method of statistics [15], and it has been applied, from this side, for sparse-grid methods in high-dimensional problems [16]. Recently cut-HDMR is receiving significant attention in other fields, for example machine learning in engineering [17].

III. IMPLEMENTATION OF KRR METHOD
Kernel ridge regression (KRR) is a non-parametric form of regression. Here we describe the specific procedure utilized in the calculations presented in this work.
Given a continuous function F (X), the KRR method provides us with an approximation represented as follows: where X ∈ R d , {X l ∈ R d | l = 1, . . . , N } is a set of points belonging to the domain of F (known as "feature vectors"), and k is the so-called "kernel" function. Specifically, in this this work we utilized the so-called the "radial basis function (RBF)" kernel (also known as the Gaussian kernel), which is defined as follows: where X 2 = N m=1 |X m | 2 is the standard Euclidean norm. The procedure for determining the coefficients α l and the kernel width parameter σ is the following: • A "training data" set {Y l = F (X l ) ∈ R d | l = 1, . . . , N } is constructed by evaluating F on the feature vectors {X l ∈ R d | l = 1, . . . , N }.
• The following minimization is performed: where [K σ ] lm = k σ (X l , X m ) is the kernel matrix, [Y ] m = Y m is the training-data-set vector and the right member of Eq. (14) is called the "cost function".
• The coefficients α l of Eq. (12) are: where σ and λ, typically named "hyperparameters," are determined using the so-called "cross-validation" (CV) method, that is an empirical protocol designed to optimize the predictive power of the KRR model. In particular, λ is a regularization parameters that can be used to remove singularities in Eq. (15) and avoiding over-fitting.
The standard k-fold procedure consists in dividing the data set into a finite number b of batches (folds) with similar size. The ML solver is trained using b − 1 of these batches at a given pair of λ and σ values. The (trained) ML solver is subsequently used for predicting the data points in the excluded fold, and the errors in the predictions are measured. This procedure is repeated, excluding each fold once and keeping a running total of the error for the given pair of hyperparameters. A mesh of possible combinations of λ and σ are tested as described above, and the pair that yields the smallest error is chosen.
In our work we employed a variation to the classic k-fold CV outlined above. Specifically, we have restricted the optimization of the hyperparameters λ and σ to values such that λ 10 −6 . The purpose of this cutoff was to reduce spurious high-frequency components in the KRR function approximation [Eq. (12)]. This improved the overall stability of the GA+nKRR algorithm, without compromising the accuracy of the nKRR solver.

Data-set mesh
By construction, the meshes of the EH input variables Y i -utilized in our calculations for training the KRR solver-are symmetric around 0. The respective number of data points per axis m i , as well as the minimum Y min i and maximum Y max i values of the corresponding intervals, are reported in Table II.