Least-square approach for singular value decompositions of scattering problems

It was recently observed that chiral two-body interactions can be efficiently represented using matrix factorization techniques such as the singular value decomposition. However, the exploitation of these low-rank structures in a few- or many-body framework is nontrivial and requires reformulations that explicitly utilize the decomposition format. In this work, we present a general least-square approach that is applicable to different few- and many-body frameworks and allows for an efficient reduction to a low number of singular values in the least-square iteration. We verify the feasibility of the least-square approach by solving the Lippmann-Schwinger equation in factorized form. The resulting low-rank approximations of the $T$ matrix are found to fully capture scattering observables. Potential applications of the least-square approach to other frameworks with the goal of employing tensor factorization techniques are discussed.

In this context second-quantized representations of two-and many-body operators provide the fundamental input for all basis-expansion methods, e.g., manybody perturbation theory (MBPT) [18][19][20][21][22], coupled cluster (CC) theory [23][24][25], self-consistent Green's function (SCGF) theory [26][27][28], and the in-medium similarity renormalization group (IMSRG) [29][30][31].While the use of a single-particle basis is very convenient in practice, the operator representation in this basis requires an extensive number of basis functions to enable robust extractions of nuclear observables.As such, finding and employing alternative operator bases provides a promising alternative to more efficiently represent the underlying objects.Recently, the use of low-rank operator expansions obtained from a truncated singular value decomposition was shown to provide excellent approximations to chiral two-nucleon interactions [32,33].Based on such low-rank approximations it was shown that two-nucleon scattering as well as ground-state properties of mediummass nuclei and the nuclear-matter energy can be very well reproduced from low-rank approximations of chiral interactions.
One of the major advantages of such low-rank approximations is their potential to reduce the required computational resources with respect to storage and the operation cost associated to tensor contractions such as matrix multiplications.However, fully exploiting the structure of factorized many-body operators requires a reformulation of the underlying many-body approach in terms of the decomposition components themselves.While this strategy has been extensively studied in quantum chemistry (see, e.g., Refs.[34][35][36][37][38][39][40]), in nuclear physics we are just starting to explore such ideas in many-body calculations.Finally, the use of factorized tensor representations is at the heart of the density matrix renormalization group (DMRG), which has been used with great success in condensed matter physics and quantum chemistry [41][42][43].Recently, the DMRG ansatz has also been employed in nuclear physics applications [44][45][46][47].In addition, there have been various applications employing factorization techniques in fitting procedures or as diagnostic tools [48][49][50][51].
Ultimately, factorization techniques may provide a way of extending the reach of ab initio nuclear structure calculations to heavier and more exotic systems.The computational demands of such calculations are due to i) the increase in model-space dimension necessary to obtain converged calculations of heavy nuclei [16,17], ii) the need for refined truncation schemes in the many-body expansion [31,52,53], and iii) the use of symmetryunrestricted bases to account for nuclear deformation effects in open-shell nuclei [53][54][55][56].Exploiting the low-rank properties of nuclear interactions can help to push the present frontiers to access significantly larger many-body spaces and thus better capture correlations.This strategy is complementary to ongoing efforts to compress calculations using importance truncation methods [57][58][59][60] and to construct improved bases with superior conver-gence properties [53,[61][62][63].
In this work, we present a novel strategy that builds upon a least-square minimization of the decomposition error of the unknown tensor in a given framework.The equations we obtain operate exclusively on the decomposed factors without reconstructing the full operators at any point.Additionally, they are independent of the details of the few-or many-body method.For this reason, the least-square approach is a general strategy that can be used to reformulate few-and many-body methods to exploit tensor factorization techniques.As a proof of concept we apply the least-square approach to the Lippmann-Schwinger equation and the full T matrix.This paper is organized as follows.In Sec.II, the leastsquare approach is introduced.Section III provides the application to the Lippmann-Schwinger equation including numerical results for low-rank T matrices.Finally, we conclude with an outlook on future perspectives in Sec.IV.

II. LEAST-SQUARE FACTORIZATION A. General rationale
In the following, we aim at deriving a factorized form of algebraic equations of the form where T denotes the unknown tensor object, V denotes the (nuclear) interaction, and the function f (•) encodes the specifics of the underlying few-or many-body framework.The ellipsis indicates the possible presence of additional (method-specific) tensors in a given framework.We start from a factorized representation of the manybody tensor where the objects X (i) define the factors of the decomposition (see, e.g., Ref. [64] for a review on tensor decompositions).Obtaining computational benefits from such a factorization requires the reformulation of the fewor many-body formalism, fully operating on the factors themselves instead of the initial (undecomposed) tensors.
Practically, this yields a new set of equations X (m) = g m (X (1) , ..., X (m) , V, ...) , where the update functions g i depend on the chosen decomposition.In addition, we can also employ a factorized form for the interaction V using a potentially different tensor format.

B. Tensor format
The general strategy laid out here can be applied to different tensor formats.In this work, we focus on the singular value decomposition (SVD) of an N × N matrix M , which we take to be real for simplicity, where (•) † denotes the Hermitian adjoint.The diagonal matrix Σ = diag(s 1 , ..., s N ) contains the ordered set of nonnegative singular values s i , and the left and right matrices of singular vectors, L and R, are unitary.By keeping only the leading R SVD singular values s 1 , ..., s RSVD and the corresponding columns of the L and R matrices we obtain the truncated SVD of the matrix M (indicated by the tilde) which provides a rank-R SVD approximation to the initial matrix.In the following, we assume a factorized form for the two-body interaction and similarly for the unknown many-body tensor1 While many other matrix decompositions exist, e.g., eigenvalue and Cholesky decompositions, the SVD format is particularly versatile since it requires neither normality nor positive definiteness of the matrix.
In quantum chemistry, factorizations into even more tensors have been applied in the context of MBPT and CC calculations [35,36,40].The employed tensor hypercontraction (THC) formats are governed by a larger amount of decomposition factors, i.e., m = 5 in Eq. ( 2) as opposed to m = 3 for the case of the SVD.For a discussion of the THC format in the context of nuclear theory, see Refs.[58,65].

C. Minimization procedure
We follow a least-square approach that minimizes the distance of the decomposed tensor to its original counterpart [36].By introducing the error tensor ∆ T ≡ T − T we define the cost function ⌃T < l a t e x i t s h a 1 _ b a s e 6 4 = " s J e S m 7 f 5 6 ⌃T < l a t e x i t s h a 1 _ b a s e 6 4 = " s J e S m 7 f 5 6 LT < l a t e x i t s h a 1 _ b a s e 6 4 = " C k p 8 k t X 0 q Z u p w a 8 LT < l a t e x i t s h a 1 _ b a s e 6 4 = " C k p 8 k t X 0 q Z u p w a 8 T < l a t e x i t s h a 1 _ b a s e 6 4 = " l C i 1 Y K r M 9 c w Z a 9 l L d 5 N G r r 5 g J 6 0 = " > A A A B 8 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R V s M e C F 4 9 V + g V p L J v N J l 2 6 2 Q 2 7 G 6 G E / g w v H h T x 6 q / x 5 r 9 A I G A c / w C m + W t l 6 s d + t j 2 V q y i p l T + A P r 8 w c C + J E L < / l a t e x i t > T † < l a t e x i t s h a 1 _ b a s e 6 4 = " Z S 5 Q Q 6 g + + n V w 5 e y k u i Y e 4 W 1 n W i o = " > A A A B 8 H i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 9 m t g j 0 W v H i s 0 C 9 p 1 5 L N Z r e h S X Z J s k I p / R V e P C j i 1 Z / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e k H K m j e t + O 4 W N z a 3 t n e J u a W / / 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " l q V c T X e Z 7 t I Z a q g g / M f 6 T i 6 q n e V s u 1 q z y O I j p E R + g E 2 e g c 1 d A N q q M m o u g R P a N X 9 G Y 8 G S / G u / E x a y 0 Y + c w + + g P j 8 w e F w J Z J < / l a t e x i t > cost T = Figure 1.Diagrammatic representation of the cost function, Eq. ( 9). - ⌃T < l a t e x i t s h a 1 _ b a s e 6 4 = " s J e S m 7 f 5  where • Fro.denotes the Frobenius norm, This can be diagrammatically represented as shown in Fig. 1, where orange symbols indicate Hermitian adjoints of the blue symbols and connecting lines correspond to tensor contractions.The factorized working equations are obtained by optimizing cost T , i.e., setting partial derivatives with respect to the decomposition factors to zero where X ∈ {L T , Σ T , R † T }.Because the function cost T is real-valued and analytic, derivatives with respect to X and X † are linearly dependent and only one set must be taken into account.Diagrammatically, performing a derivative ∂cost T /∂X corresponds to the removal of the corresponding tensor vertex X in the tensor network (see Fig. 2 for the derivative with respect to X = L † T ).For the various decomposition factors we obtain from Eq. ( 10) the following set of derivatives for SVD-based decompositions: LT < l a t e x i t s h a 1 _ b a s e 6 4 = " C k p 8 k t X 0 q Z u p w a 8 where f (T ) encodes the (non-factorized) working equation, Eq. (13).The specific example of the Lippmann-Schwinger equation will be discussed in Sec.III.

D. Master equations
Since the derivative ∂cost T /∂X † is linear in X, all factors other than X can be contracted in a so-called environment matrix A X .Figure 3 shows the example of setting the derivative with respect to L † T to zero, where the contraction of the colored area gives the environment matrix A L associated with L T .We are left with the solution of a linear problem X • A X = B X , where B X corresponds to the first terms in Eqs.(12) and we note that the environment matrix can be a left and/or right factor.Thus, the update step can be written as The explicit expressions for the environment matrices are In contrast to the left and right matrices the derivative ∂cost T /∂Σ † T produces two environment matrices A Σ1 and A Σ2 .Note that the environment matrices are tensorformat specific and do not depend on the many-body approach itself (which is encoded in the tensors B X ).Thus, the system of Eqs. ( 12) constitutes a set of master equations governing any SVD-structured framework using algebraic equations.Finally, the update steps for the different factor matrices are given by where all evaluations in the tensor-structured framework can be performed using efficient linear algebra operations.In the following, we will refer to the update equations, Eqs.(15), as the SVD-factorized least square (SVD-LS) equations.

E. Explicit orthogonalization
The solution of Eqs. ( 15) is not constrained to conserve unitarity of the left and right matrices, i.e., L † Unitarity can be explicitly enforced through an additional orthogonalization from a QR factorization, where Q R , Q L are unitary matrices and R L , R R are upper triangular matrices.Modified factor matrices are obtained via such that the left and right matrices are manifestly unitary.To restore the diagonality of Σ T we perform an additional (untruncated) SVD where the unitary matrices LT and RT are absorbed into the left and right matrices, respectively, and we obtain the final decomposition Practically, the orthogonalization is performed in each iteration step.Due to the small size of the corresponding matrices the computational overhead is negligible.The restoration of a proper SVD format has the formal advantage of simplifying the master equations of the least-square approach.Due to the diagonality of the Σ T matrix the environment matrices [Eqs.(14)] can be analytically inverted giving rise to the update step The simplified update step in Eqs.(20) reduces the numerical sensitivity to ill-conditioned environment matrices in the formation of their inverses.In the following, we refer to the SVD-LS approach as the one with explicit orthogonalization, Eqs.(20).

F. Computational advantages
The computational benefit of factorization techniques comes from the decreased memory requirements and lower number of floating point operations in the tensor contractions, both of which depend on the dimension of the matrix objects.
The evaluation cost of a matrix-matrix product scales as C = O(N 3 ) for dense N × N matrices each with an associated memory cost M = O(N 2 ).In presence of an SVD-factorized matrix the SVD-LS update equations induce an operation count of C = O(N 2 R SVD ) with memory cost M = O(N R SVD ) for the factorization components.The full scaling is recovered in the limit of an exact decomposition, which in the case of an SVD corresponds to keeping all N singular values, lim We note that this statement is only true asymptotically, since the SVD induces a storage overhead in absence of a truncation.However, this will only affect the prefactor and not the scaling exponent.In practice, obtaining benefits through factorizations relies on the low-rank properties of the tensors such as the input interaction, i.e., how accurate low-rank approximations are at R SVD N .

III. TWO-BODY SCATTERING
A. Lippmann-Schwinger equation In the following, we systematically apply the leastsquare approach to the Lippmann-Schwinger equation where G 0 denotes the (diagonal) free Green's function and T the T matrix, which we take to be right-side halfon-shell.Equation ( 22) is of the general form of an algebraic update equation, Eq. ( 13), with For the initialization of the T matrix factors the firstorder Born approximation T (0) = V is employed leading to X (0) T = X V for X ∈ {L, Σ, R † }, thus requiring the same target rank for the T matrix and the potential.
The two-nucleon (NN) potential (and the T matrix) are represented in a partial-wave basis with the final and initial orbital angular momenta l and l , the two-body spin S, the two-body total angular momentum J, the two-body isospin T with projection M T , and the absolute values of the outgoing and incoming relative momenta k and k .In each partial-wave channel, the NN potential is represented using N = 100 momentum mesh points up to k max = k max = 6.0 fm −1 .Similarly, the Lippmann-Schwinger equation is solved in a partial-wave-decomposed form (with 2 /m = 1), kα|T where α, α , α are collective labels for the partial-wave quantum numbers.

B. A priori decomposition analysis
Before turning to the study of the SVD-LS approach, we begin with a study of the low-rank properties of the T matrix obtained from direct-inversion techniques of the Lippmann-Schwinger equation.In the following, we employ the N 3 LO NN potential from Entem, Machleidt, and Nosyk (EMN) with a cutoff Λ = 500 MeV [7].Note that the singular value spectrum is qualitatively similar for different orders, different cutoffs, and with similarity renormalization group (SRG) evolution [32,33].
Figure 4 shows a comparison of the singular spectrum of the initial NN potential and the final right-side half-onshell T matrix for different partial-wave channels.Similarly to Ref. [32], we divide the rank in the coupled channel by a factor of two to be able to compare channels with different matrix dimensions.It is evident that the initial low-rank properties directly propagate to the T matrix and the T matrix itself is dominated by very few components in the SVD expansion.For the T matrix in the 1 S 0 channel there is a strong enhancement of the leading singular value s 1 by a factor of ten going from V to the T matrix due to the large scattering length.

C. Numerical convergence of the least-square approach
The self-consistent solution of the SVD-LS equations is obtained by consecutive updates of L T , Σ T , R † T while keeping all other factors fixed.Convergence is gauged by the relative norm of the difference between consecutive factor matrices, where the superscript indicates the iteration number.
Figure 5 shows the rate of convergence for the decomposition factors as a function of iteration number for the 1 S 0 and the coupled 3 P 2 -3 F 2 partial wave.Clearly the large-scattering-length 1 S 0 channel requires a significantly larger number of iterations compared to the weaker 3 P 2 -3 F 2 channel.Moreover, the convergence for the 1 S 0 channel strongly depends on the initial rank of the potential, with a significant increase of iterations needed until convergence beyond R SVD ≈ 10.We attribute this to numerical instabilities in the inversion of the environment matrices due to small singular values at higher rank (e.g., s 20 10 −5 ).However, these high-rank components are not important for an accurate reproduction of the NN T matrix.For the 1 S 0 channel the rate of convergence of the decomposition factors Σ T and R † T is slower than for L T , in particular at low rank R SVD 3.This behavior is likely related to the use of the right-side half-on-shell T matrix, k|T (E = k 2 )|k , so that the iteration is sensitive to the energy dependence of the free Green's function G The partial-wave dependence of the rate of convergence can be understood from an analysis of the integral kernel K(E) = V G 0 (E) that enters the Lippmann-Schwinger equation.In an iterative approach, the final T matrix is obtained as the infinite Born series so that K(E) drives the numerical stability of the leastsquare approach.This can be quantified in terms of the spectral radius of the integral kernel where λ i are the eigenvalues of K(E).As the Born series constitutes a geometric series, eigenvalues |λ i | 1 will prevent an iterative approach from converging.While eigenvalues |λ i | close to unity will not necessarily prevent convergence, they induce a much lower rate of convergence in practice.For the deuteron 3 S 1 -3 D 1 channel, the presence of the bound state leads to an eigenvalue |λ i | = 1 at the deuteron binding energy E = −2.2245MeV, so that we do not present results for the deuteron channel.
This problem of non-convergence in an iterative approach can be circumvented using direct inversion techniques, which is easily possible for two-body scattering because of the small matrix dimension.However, already in the three-body sector one relies on iterative schemes and thus naturally encounters diverging Born series (see, e.g., Ref. [66]).This can be resolved by employing Padé resummation techniques on the individual terms of the Born series, which enables a robust extraction of scattering observables [67].

D. Diagnostic of the low-rank T matrix solution
We continue our analysis with the characterization of the solution of the SVD-LS approach.The exact T exact matrix is obtained from the full-rank, R SVD = N , interaction and by solving the Lippmann-Schwinger equation via direct inversion [67].We compare the low-rank T matrix from the least-square approach, T RSVD SVD-LS = L T Σ T R † T , to the exact T matrix.As error measure we study in the following absolute and relative differences of the matrix object, A matrix plot of the low-rank T matrix compared to the exact solution is provided in Fig. 6.At rank R SVD = 1 and 2 we observe a sizeable difference from the full-rank  T matrix, in particular in the low-momentum regime k, k 2 fm −1 .Once the rank is increased, deviations decrease systematically yielding only minor differences for R SVD = 3 and excellent agreement at R SVD = 5.
Using an initial NN interaction, there are two ways to obtain the low-rank T matrix: i) Decompose and reconstruct: Perform a low-rank approximation for the initial potential and use the truncated potential to obtain the low-rank T matrix from direct inversion.In this case, the T matrix factors are obtained from an explicit SVD of the resulting T matrix.
ii) Least-square approach: Perform a low-rank approximation for the initial potential and use its decomposition factors as input for the least-square approach described in Sec.II.The final low-rank T matrix is then given by the converged factors after the least-square iteration.
At fixed rank R SVD , both strategies yield an equivalent final solution, up to unitary transformations among the left and right matrices due to the re-orthogonalization.
Table I shows the quality of low-rank T matrices compared to the exact results.From rank R SVD = 1 to 2 there is a slight increase in relative error, but for larger ranks the relative error systematically decreases as the rank is increased.The leading singular values of the T matrix do not remain constant as the rank is increased due to the nonlinear dependence between the singular values of V and T .At rank R SVD = 10 the singular values stabilize and agree with the exact singular values at full rank, R SVD = 100.

E. Two-nucleon phase shifts
We finally turn to the description of two-nucleon phase shifts based on low-rank T matrices.The phase shifts are obtained from the least-square T matrix using the converged factor matrices T = L T Σ T R † T .Figure 7 shows the phase shifts for the 1 S 0 and the coupled 3 P 2 -3 F 2 partial waves using low-rank T matrices with R SVD = 1, . . ., 5. At very low rank, the phase shifts significantly deviate from the exact results.This is most pronounced in the 3 F 2 partial wave and for the mixing angle 2 , while the 3 P 2 channel shows very little sensitivity to the rank.The enhanced sensitivity in the mixing angle has also been observed for the deuteron channel in Ref. [32].The quality of the approximation is systematically improved in all channels and at rank R SVD = 5 the phase shifts are excellently reproduced up to laboratory energies E lab 300 MeV.Since the low-rank T matrices were already shown to reproduce the exact T matrix very well (see Fig. 6), it is clear that derived quantities yield a similarly good approximation error.

IV. CONCLUSIONS AND OUTLOOK
In this work, we have presented a new strategy to solve algebraic equations from a factorization ansatz.Following a least-square approach, the master equations are derived from a stationarity condition of the cost function based on the error tensor.We have derived a set of update equations for the individual factor matrices.This strategy is general and can be used to exploit tensor factorization techniques in few-and many-body calculations.Moreover, the SVD-LS equations can be adapted to other tensor formats as previously shown in quantum chemistry applications [36].The feasibility of the least-square approach is practically demonstrated for the Lippmann-Schwinger equation.By employing an SVD form for the potential, a factorized form for the T matrix was obtained.Using an additional explicit orthogonalization during the self-consistent iterations enables the recovery of a proper SVD format for the T matrix itself.The right-side half-on-shell T matrix and twonucleon phase-shifts are in excellent agreement with fullrank calculations, reflecting the low-rank structures of chiral interactions [32,33] and their propagation to the associated T matrices.We note, however, that the transformation to single-particle bases may lower the efficacy of the SVD approximation of nuclear potentials in the many-body calculations due to the coupling with centerof-mass degrees of freedom as discussed in Ref. [33].
While we have discussed computational benefits in terms of abstract scaling laws, their demonstration in our two-body application is difficult due to the very short runtimes of the order of hundreds of miliseconds.Future applications to larger-scale simulations will allow for more systematic studies of the underlying computational benefits.
Our work also suggests extensions to nuclear few-and many-body calculations.An interesting application will be to explore the least-square approach in nonperturbative many-body methods, such as CC or IMSRG calculations of medium-mass nuclei.These have the necessary algebraic equations for the coupled-cluster amplitudes or the evolution operator, respectively, which could benefit from factorization methods through the least-square approach.
N e P C j i 1 b / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e m H C m j e d 9 O 4 W N z a 3 t n e J u a W / / 4 P D I P T 5 p 6 z N e P C j i 1 b / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e m H C m j e d 9 O 4 W N z a 3 t n e J u a W / / 4 P D I P T 5 p 6 z N e P C j i 1 b / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e m H C m j e d 9 O 4 W N z a 3 t n e J u a W / / 4 P D I P T 5 p 6 z

Figure 2 .
Figure 2. Derivative of the cost function with respect to L † T .

T < l a t e x i t s h a 1 _
b a s e 6 4 = " U m x d N E o 8 K z c J + l T k q H Y 1 T e s H E 9 8 = " > A A A B 6 H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x G w R w D X j w m k B c k S 5 i d d J I x s 7 P L z K w Q l n y B F w + K e P W T v P k 3 T p I 9 a G J B Q 1 H V T X d X E A u u j e t + O 7 m t 7 Z 3 d v f x + 4 e D w 6 P i k e H r W 1 l I 8  u I E a 3 E E d m s B g B M / w C m + O c F 6 c d + d j 2 Z p z s p l T + A P n 8 w c F J I 2 X < / l a t e x i t > ⌃T < l a t e x i t s h a 1 _ b a s e 6 4 = " s J e S m 7 f 5 6x R 5 t U 4 C d x 5 D Q B 5 K 9 U o = " > A A A B 7 3 i c b V B N S w M x E J 2 t X 7 V + r X r 0 E i y C p 7 J b B X s s e P F Y s V / Q L i W b Z t v Q J L s m W a E s / RN e P C j i 1 b / j z X 9 j 2 u 5 B W x 8 M P N 6 b Y W Z e m H C m j e d 9 O 4 W N z a 3 t n e J u a W / / 4 P D I P T 5 p 6 z

Figure 3 .
Figure 3. Diagrammatic representation of the derivative ∂costT /∂L † T = 0.The colored area corresponds to the environment matrix (see text for details). ∂cost

Figure 4 .
Figure 4. Singular values (in fm) for the VNN and T matrices in different partial-wave channels for the N 3 LO EMN 500 potential.For the coupled 3 P2-3 F2 channel the rank is divided by two.

3 P 2 - 3 F 2 Figure 5 .
Figure 5. Convergence of the decomposition factors of the T matrix in the SVD-LS approach as function of iteration number in the 1 S0 channel (upper panels) and the coupled 3 P2-3 F2 partial wave (lower panels) for the N 3 LO EMN 500 potential.Results are shown for different rank approximations.

Figure 6 .
Figure 6.Absolute difference of low-rank T matrices obtained from the least-square approach compared to the exact solution of the Lippmann-Schwinger equation via direct inversion.Results are shown for the N 3 LO EMN 500 potential in the 1 S0 channel.

2 Figure 7 .
Figure 7. Two-nucleon phase shifts as a function of laboratory energy in the 1 S0 and 3 P2-3 F2 partial waves based on the low-rank T matrices obtained from the least-square approach and in comparison to the exact T matrix.Results are shown for the N 3 LO EMN 500 potential.