Least-squares fitting with non-local model functions

and Related Phenomena, 40 (1986) 301-306 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netlierlands

Journal of Electron Spectroscopy

LEAST-SQUARES FI’ITING WITH NON-LOCAL MODEL FUNCTIONS

G.K. WERTHEIM and S.B. DICENZO

AT&T Bell Laboratories, Murray Hill, NJ 07974 (U.S.A.) (Received 7 January 1986)

ABSTRACT The implementation of an algorithm for least squares analysis with non-local model functions which minimizes the computational effort is described. A 20-parameter FORTRAN version has been run routinely on a VAX-11/785; a lo-parameter version written in C runs effectively on an AT&T PC6300.

INTRODUCTION

In a recent publication [l] we discussed the advantages of using least squares fitting to extract information from experimental photoemission data. This procedure will yield physically meaningful results provided the mathematical model function correctly embodies the physical processes which underlie the data. The use of such models is not without cost, because the required functions are not available in closed form. In general they contain one or more convolution integrals to introduce phonon and instrumental broadening. To evaluate such a function at a point requires integrations spanning the width of the broadening functions, i.e., the functions are non-local. Since least squares adjustment requires repeated evaluation of the model function, the computation tends to become burdensome. In fact, least squares fitting with such functions has the reputation of being very demanding of computer time. We will show here that this is in part due to the way in which the least squares computer program calls the model function, and describe a modification which allows least squares fitting with nonlocal functions even on a personal computer. We begin by sketching the mathematics underlying the least squares adjustment [2-41. We consider a set of N data points, yl, and a model function, f&), which is a function of n variables pi. The objective is to find the set of values for these variables which minimize R, the sum of the squares of the deviations between the data points and the model function

302

R =

2 (Yi -fi(Pj))2

(1)

I=1

At the minimum of R the partial derivatives of R with respect to the n parameters are equal to zero. Thus we have n equations of the form f(Yi-fi)

g

=

O

l
where the subscript i in the partial derivative indicates that it has been evaluated at the location of the ith data point. For a non-local function, the partial derivatives are not analytic and must be evaluated numerically. The solution to the problem, that is, the best fit to the data, consists of the n values of p which satisfy these II simultaneous equations. If f is a linear function of the parameters pi, these are simply n simultaneous, linear equations in the variablespl, which can be solved by standard methods, i.e., it is a problem in linear regression. In the cases of interest here, the dependence on pj is non-linear, and there is no general solution in closed form. In order to solve this problem by numerical means, we assume that we have an approximate solution (an initial set of values for p), and expand f in a Taylor series about this point

With this substitution

eqn. (2) takes the form

The right hand side of the n linear equations in (4) may be written as a product A Ap of an n x n matrix A and an n-dimensional vector Ap, where the elements a]* of the matrix are defined by

The solution to eqn. (4) is obtained by left-hand multiplication with the inverse of the matrix A. This yields a set of corrections, Apk, to the initial values of p. The new values of p produce a smaller value of R. Because f is not actually a linear function of pi, this new value of R is not the minimum. However, the minimum may be approached as closely as desired by an iteration of this process. The use of a Taylor series expansion requires a good set of initial parameters, if the least-squares adjustment is to avoid spurious solutions. We have found it advantageous to use an interactive graphics program that super-

303

imposes the model function on the data points and allows the user to change the input parameters until a reasonable, preliminary fit is obtained.

IMPLEMENTATION

To implement this procedure on the computer,. the left hand side of eqn. (4) as well as the matrix A are evaluated numerically. This requires the values of the function and its partial derivatives at the location of the data points. For the model functions appropriate for photoemission, neither the function nor its derivatives are available in closed form. The partial derivatives of the function are obtained by evaluating the function with the initial parameters, and with the parameters incremented, one at a time, by a small amount. To set up the A matrix the products of these derivatives evaluated at the location of each data point are summed. The matrix is then inverted by standard techniques, and the correction vector, Ap, is obtained by multiplying it with the left hand side of eqn. (4). The most time consuming part of this calculation is the setting up of the A matrix. The function must be evaluated n + 1 times at the location of the data points to obtain the partial derivatives. Since the product of each partial derivative with all the others is needed for each data point, the n values of the partial derivatives must be stored during the calculation. To avoid more extensive storage most versions of the least squares algorithm call the model function with the n + 1 sets of parameter values, for one data point at a time, requiring (n + l)N calls in all. With a non-local function this is highly inefficient, because the entire function is calculated on each call. This accounts for the excessive time required for least-squares fitting of nonlocal functions. The computational effort is drastically reduced by recasting the least squares algorithm so that the model function is called only once for a given set of parameter values, and then gives the values for all data points. The number of function or subroutine calls is thereby reduced by a factor N, and the number of calculations of the underlying function by a factor equal to the number of data points that are encompassed in the convolution integral. The only cost incurred by this change arises from the necessity to store a complete set of n x N partial derivatives while the matrix is being set up. When the original non-linear least squares programs were written, core storage was at a premium and considerable effort was made to avoid the use of large arrays; in present-day computers little is gained by adhering to this practice. Our versions of the least squares algorithm are based on the work of Marquardt [4], as implemented in FORTRAN IV by Burnette and Roberts at AT&T Bell Laboratories. The C-language version was written by M.F. Robbins, and further developed by R.J. Morris. We have incorporated the

starting

Par.nm*rs

1 10.0000 2 0.2000 3 0.1200 4 153.0000 5 1960.0000 6 5.0000 7 1.0000 e -5.0000 9 (1.5000 Index of Fixed

xtcr. 1

: : a 7 8 9 it

:

Pqams:

R/t< 5.8972298

0

COP”W

Lmbd.

0.1000000 0.0100000 0 0010000 0.0001000 0.1000000 0.0100000 0.1000000 0.1000000 0 1000000 1.0000000

I.6946406 1.0103b93 0.9638123 0.9681b15 0.9625680

0.9562967 0.9554126 U.9544i98 0.9544505

l.GOOO

7

-0.4496 1.0000

3 4

34.42 b3 52 45.30 72.49 65.21 LT.50 53.;; 61.57 44.23 15.34

Paran

-0.1261 0.3066

0.1375 153.8076 2012.4017 3.2631 ,.OOOO -5.5041 0 6617

Std.

26.13 135 39 53.70 20 lb 19.67 12 44

-0.6012 0.2735

-0.0711 0.1600

-0.0214 0.0519

-0.5472 0.3762

1.0000

0.2638 1.0000

-0.2093 -0.6587

0.3570 -0.4273

1.0000

-0.0539 1.0000

: 7

Fznal

626.21 731 ii 187.37 98.70

0.0000 0 0000 0.0000 0.0000 0.0000 O.ODOO 1.0000

0.1506 -0.3991 0 5379 0.2417 -0.5609 0 3039 0.0000 1.0000

-0.Ob55 0.1495 0.047d 0 8642 -0.7574 -0.4022 O.OOOO 0 3202 1.0000

Error

3.2035 0.0794 0. OOP’? 0.0485 20.b213 0 2150 0.0000 0.0471 0.017')

Urzt1o9 Output 15~42 0 pea1

Fllrr 2.03.2

user

9.7

.y*

Fig. 1. Sample of the output provided by the C-version during a least-squares adjustment. The parameters represent the following (1 channel ~0.05 eV): (1) constant background, in counts; (2) Lcrentzian half-width, in channels; (3) singularity index in Doniach-Sunjic equation; (4) position of j = 312 line, as channel number; (5) height of j = 312 line, in counts; (6) Gaussian full-width due to phonon excitation, in channels; (7) ratio of Lorentzian width of j = l/2 to 3/2 line (held fixed); (8) spin-orbit splitting, in channels; (9) ratio of intensities of j = l/2 line to 3/2 line.

changes described here into both versions of the non-linear least squares algorithm. The substantive changes are confined to the function which sets up the A matrix, but all calls to the model function have been revised. A lo-parameter C version has also been run on an AT&T PC6300 using the VENIX/86 [ 51 operating system.

306

I

52

I

51

I

60 BINDING

I

49 ENERGY

Fig. 2. Fit to the Mg 2p spectrum

48 (BV)

obtained

by the run chronicled

in Fig. 1.

DISCUSSION

As a test of the efficacy of this version of the least squares algorithm we compare the cpu times required to fit the Mg 2p spin-orbit doublet taken from Citrin et al. [6]. The model function contains nine parameters, representing the background, lifetime widths of the two components, singularity index, position, amplitude, phonon width, and the spin-orbit splitting and intensity ratio. The model function contains two convolution integrals, one for the fixed 0.25 eV instrumental resolution function, and the other for phonon broadening. This is a non-trivial problem, because the spin-orbit splitting is comparable to the width of the resolution function so that the doublet is unresolved. The progress of a least squares adjustment is chronicled in Fig. 1 which shows the output obtained during a run on a VAX-11/785 [7] using a UNIX [8] operating system, in which 101 data points were fitted with a model function containing eight free parameters. As a measure of the success of the fitting procedure we use the ratio of R to the sum of the counts in the data points being fitted. Counting statistics alone lead to a ratio of unity, when the number of data points is much greater than the number of free parameters. This ratio is close to the final value after five iterations, but the program continues to iterate because we use relatively stringent convergence criteria. The rapid approach to the final parameters is characteristic of model functions without strong parameter correlations. The fit obtained is shown in Fig. 2. The average cpu time per iteration is 12 s for both

306

FORTRAN and C versions. The time depends linearly on the number of data points fitted, and roughly on the second power of the number of parameters. Many experimentalists refrain from fitting non-local model functions to their data because they believe the cost will be prohibitive without free access to a large computer. Such expectations are perhaps based on attempts to use a less efficient version of the least-squares algorithm. We have tested a lo-parameter fit on an AT&T PC6300 personal computer. While the cpu time per iteration is typically ten times what it is on a VAX, the real time required may be comparable because of the time-sharing load on the VAX. Given the proliferation of small computers in the research community, we believe that this demonstrates that this data analysis tool is accessible to essentially all experimentalists.

REFERENCES 1 2 3 4 5 6 7 8

G.K. Wertheim and S.B. DiCenzo, J. Electron Spectrosc. Relat. Phenom., 37 (1985) 57. N. Draper and H. Smith, Applied Regression Analysis, Wiley, New York, 1966. P.R. Bevington, Data Reduction and Error Analysis for the Physicsl Sciences, McGrawHill, New York, 1969. D.W. Marquardt, J. Sot. Ind. Appl. Math., 11 (1963) 431; Share Program Library SDA 3094. VENIX/86 is a trademark of VenturCom, Inc. P.H. Citrin, G.K. Wertheim and Y. Baer, Phys. Rev. B, 16 (1977) 4256. VAX is a trademark of Digital Equipment Corporation. UNIX is a trademark of AT&T Bell Laboratories.