The minimum norm method for the determination of the charge density from elastic electron scattering data

Unphysical behavior in the QR algorithm based least squares determination of the expansion coefficients of the charge density obtained from limited information about the charge form factor occurs when the spread of the singular values in the matrix relating these quantities becomes too large. Setting the smallest singular values equal to zero in the singular value decomposition used in the minimum norm method yields a much more reasonable determination of the charge density. Increasing the size of the basis without increasing the range of the prior information about the charge form factor leads to ambiguities in the determination of the charge density. Numerical results in an analytic model are presented.

electron-nucleus interaction is relatively weak, multiple scattering effects can be neglected and the scattering process can be described in first order perturbation theory.The connection between the charge density and the cross section is well understood and in plane wave Born approximation Fc(q) is just the Fourier transform of p(r) which for the case of even-even nuclei, which we shall consider, is simply given by oo Fc(q) = 4n ~ dr r2jo(qr) p(r) (1) 0 wherejo is the spherical Bessel function of zero order and q is the absolute value of the three momentum transfer.Given that the experimental measurements are performed over a limited range at a finite number of values of the momentum transfer, q, a unique determination of p(r) is not possible since the resulting inverse problem is ill posed.The generally accepted procedure for determining p(r) is to expand it in a Fourier Bessel (FB) basis [1][2][3] and then determine the expansion coefficients from a least squares fit to the experimentally determined values of Fc(q).One of the problems in the aforementioned procedure is that increasing the number of terms in the expansion generally leads to unphysical oscillations in the charge density in spite of the fact that the charge form factor is well reproduced at the experimentally determined values of q.These unphysical oscillations arise, as we shall show, if one or more of the singular values of the non-square matrix to be inverted becomes abnormally small.In this case the least squares problem does not have a unique solution [41.Since the inversion in the least squares problem is generally accomplished by means of the QR algorithm [4], this fact is overlooked.If the minimum norm method, which makes use of a singular value decomposition, is used, a unique solution of the least squares problem always exists [4] and the unphysical oscillations can be eliminated.
In order to demonstrate this we make use of the following analytical model.For a charge density given by a symmetrized Fermi distribution [5] an analytical expression for the corresponding charge form factor can easily be obtained [6,7]: Only two of the parameters c~, R and d are independent since the charge density must fulfill the normalization condition 4n S drr 2 p(r) = Z. ( Since the charge density is a single-valued function defined in a finite domain it can be expanded in a basis of orthogonal functions.In the FB expansion [1][2][3] use is made of the following orthogonality relation between spherical Bessel functions in a finite domain R~ R?
(8) R~ If the form factor in Eq. ( 1) is known at q,, the coefficients a, can easily be obtained and are given by: V~(q,) a,, = 2nR~(jl(q,,Rc))2 (9) In general, however, the cross section is measured at momentum transfers different from q,.Using the expansion of the charge density given in Eq. ( 7) the charge form factor is then given by Fc(q) 4nR2~a, (-1)" ~ sin(qRc).
(10) = q .(qR-) g-_2(nn)-Now, given that n measurements are made at momentum transfers q'= (q'b q~ ..... q',) , we wish to determine the m expansion coefficients a = (al, a2, ..., a,,).In this case Eq. ( 10) leads to the following system of equations F¢ = Aa (11) where F~ ~ R" is a column vector with elements F~(~) with c~ ~ q' and the corresponding matrix A ~ R m+".If n > m and rank(A)= m, a unique solution a for which the L 2 norm tlF~ -Aal]2 is minimized exists [4].This is just the least squares problem and it is usually solved by inverting the square matrix ATA since where A r is the transpose of A. Numerically the most efficient way to accomplish this is by means of the QR algorithm which transforms ATA into a right upper triangular matrix by means of a set of orthogonal transformations.The inverse of this transformed matrix is easily obtained by simple back substitutions.However, if the rank(A) < m, there are many solutions a ~ R" for which liFe -Aalt2 is minimized.A unique solution to this least squares problem may be obtained if the additional requirement that Ilall2 be minimized is also included [4].This is nothing more than the minimum norm solution of the least squares problem and is obtained from a singular value decomposition of A. In this case A is decomposed in the following manner A = UTZV (13) where "al 0 0" 2 0 (7 r 6R m+" (14 where U ~R "+" and V~ R m+" are orthogonal and the square of the singular values a 2 are the eigenvalues of ArA with at > 0.2---Or > 0. Since this decomposition is numerically expensive, the QR algorithm is usually used to solve the least squares problem.It should also be noted that the matrix 71 = (ATA) -IAT needed to determine a in Eq. (1 1)is nothing more than the Moore-Penrose inverse of A [8,9] which often is denoted as the pseudoinverse.In general, for any n, m, the pseudoinverse of A is given in terms of the aforementioned singular value decomposition as where the non-zero diagonal elements of the matrix £-1 are a7 1.For m = n it corresponds to the inverse of A. Furthermore, a unique solution of Eq. ( 11) can be obtained with the minimum norm method for n > m when the system of equations is overdetermined [4] as well as for undetermined inverse problems with n < m.Lastly, it is interesting to note that under certain conditions Maximum Entropy methods for solving undetermined inverse problems also lead to the minimum norm method [10][11][12].
If rank(A)= m, there is no benefit in obtaining the minimum norm solution of the least squares problem as it is the same as that obtained with the QR algorithm.However, generally one does not know the rank of A. The simplest way to determine the rank of A is to determine the singular values of A. The number of non-zero singular values determines the rank of the matrix.In practice, however, if one or more of the singular values become small, the rank of the matrix A is usually less than the dimension of the square matrix ATA and as we shall show unphysical behavior will occur in the least squares fit obtained by means of the QR algorithm.In this case the problem can be avoided by using the minimum norm method and setting the small singular values to zero in the singular value decomposition of A. Unfortunately there are no hard and fast rules to determine how small a singular value should be before it is set equal to zero.Furthermore, large gaps in the values of the singular values do not always occur.
We have performed a numerical calculation for 12C in the aforementioned analytical model with Rc = 8 fm, d = 0.626 fm and R = 1.1A 1/3.As input we have taken the values of Fc(q) determined at a set of nine momentum Fig. 2. Charge density as a function ofr.The solid curve is the exact density, while the dashed curve represents the result of the minimum norm method for 30 basis states transfers (0.001 fm-t, 0.5 fro-l, 1.0 fm-1, 2.0 fm-l, 3.0 fm-i, 4.0 fm-l, 5.0 fm-i, 6.0 fro-t, 7.0 fm-1).If we require that the spread in the singular values of A be less than 105 the QR algorithm may be used to determine the set of a,, for m < 6 and the fits to the charge density look quite reasonable.For m = 7 (see Fig. 1) the smallest singular value (.00307165) is less than 10-5 the value of the largest singular value (701.129) and we have set it equal to zero in the singular value decomposition of A. In this case the use of the QR algorithm results in large unphysical oscillations in the charge density while the minimum norm method provides a reasonable fit to the charge density.In both cases we reproduce the values of Fc(q) at the nine given values of the momentum transfer.Furthermore, with the same requirement on the spread of the singular values and for the same prior information about the charge form factor, a well-behaved charge density with m = 30 may also be obtained with the minimum norm method (see Fig. 2).Here again we reproduce the values of Fc(q) at the nine given values of the momentum transfer.It is interesting to note that although the tail of the charge density agrees well with that obtained with m = 7 it does not agree well at smaller values 143 Fig. 3.The FB expansion coefficients as a function of n.The circles connected by a solid curve are the exact values of the expansion coefficients while the triangles connected by the short dashed curve designates the values obtained from the least squares method with the QR algorithm and the squares connected by the tonq dashed curve represents the result obtained from the minimum norm method for 7 basis states of r.This ambiguity arises from the fact that not enough information about the charge form factor has been provided to uniquely determine the charge density.The fact that no information has been given about the charge form factor for q _> 7 fro-i nicely demonstrates the fact that a unique determination of the charge density is not possible at smaller values of r.
In order to ascertain in a more quantitative manner the underlying reasons for the large unphysical oscillations in charge density obtained we compare the values of the expansion coefficients obtained from the least squares and minimum norm method with the exact values given by Eq. ( 9).The-results for m = 7 are given in Fig. 3.Note that the coefficients for small values of n are not well determined because of the sparse amount of information initially provided about the form factor.The exact coefficients become very small very quickly and oscillate slowly about zero.The least squares method produces much larger oscillations which persist as m is increased while the minimum norm method provides a much smoother representation of the behavior of the coefficients.The deviations in the calculated expansion coefficients, which may be substantial, lead to deviations in p as can be seen in Fig. 2. Fortunately, and this is the nice thing about the choice of the FB basis, the coefficients'with larger values of n are small and large deviations in these coefficients are not that important except at small values of r.If this were not the case the cut-off procedure used in the Fourier Bessel method would certainly not work as well as it does.Note also that setting the smallest singular value of A equal to zero is not equivalent to setting any of the values of the expansion coefficients equal to zero but rather a simple way of guaranteeing that no excessive structure is introduced into these coefficients.For this reason for larger values of n they are in better agreement In Fig. 4 we show the results for which there is a large amount of prior information as is usually the case when the charge form factor has been determined experimentally.The values of F¢(q) are given at 140 equally spaced values for 0.001 < q < 7.001 fro-1.If no assumptions are made about Fch(q) for q > 7 fm-1 reasonable and identical results are obtained for p(r) for m _< 19 from both the least squares and the minimum norm method.Note that m may take on much larger values than the suggested cut-off for rn ( <_ qma~Rc/rc) of Eq. ( 8) [1].For m _> 20 the spread in the singular becomes too large and oscillations occur in the least squares fit.In spite of the fact that no prior information is available for q > 7 fm-1 the minimum norm method provides a means of extrapolating to higher momentum which is consistent with the prior information and provides a reasonable determination of p(r).The results are, of course, model dependent as the short ranged behavior of p(r) depends on the choice of m as has been previously demonstrated.
If, as is usually the case, the experimental data, Fc, is assumed to have Gaussian errors then the uncertainties in the determination p(r) may be determined in the manner suggested in [2].In this case the variance in p(r) due to the errors in the charge form factor, AFc are given by and derivatives can be determined from Eqs. ( 7) and ( 12) (in terms of pseudoinverse ,4 used to determine the expansion coefficients).The square root of A2p(r) corresponds to the width of a Gaussian distribution which describes the uncertainty in the determination of p(r).Within the framework of Maximum Entropy methods a highly nonlinear method for treating inverse problems involving data with Gaussian errors has been given by Gull and Daniell [13].
In the present note we have shown that unphysical oscillations of the FB expanded charge density determined from information of the charge form factor occur when the spread in the singular values of the matrix relating these quantities becomes too large.These oscillations may occur in any basis and are not the result of using the FB basis [-7].When they occur, the minimum norm method with the small singular values set equal to zero should be used to obtain the expansion coefficients.Increasing the size of the basis without increasing the range of the prior information about the charge form factor leads to ambiguities in the determination of the charge density.
HGM acknowledges the financial support of the ROC-RSA Scientific Exchange Program and the Foundation for Research Development, Pretoria, as well as the hospitality of the Institute of Physics, Academia Sinica, Taiwan.YT and GDY acknowledge the support of Taiwan National Science Council under grant numbers NSC 84-2112-M-001-007 and 84o2112-M-001-025, respectively.NC is a fellow of CONICET and RR of CIC(PBA) of Argentina.RR acknowledges also a grant from Fundacion Antorchas.

Fig. 1 .
Fig.1.Charge density as a function of r.The solid curve is the exact density, while the dot-dashed curve designates the least squares fit with the QR algorithm, and the dashed curve represents the result of the minimum norm method for 7 basis states

Fig 4 .
Fig 4. Charge density as a function of r.The solid curve is the exact density, while the dot-dashed curve designates the least squares fit with the QR algorithm, and the dashed curve represents the result of the minimum norm method for 20 basis states.The minimum norm results are distinguishable from the exact calculations only at small values of r