A computer program for non-linear curve fitting

A computer program for non-linear curve fitting J. P. OTTOY and G. C. VANSTEENKISTE Department of Applied Mathematics, University of Ghent, Coupure Links 533, 9000 Ghent, Belqium

Recently several techniques for non-linear curve fitting have been developed. The implementation of a non-linear curve fitting procedure is treated for mathematical models in which the linear and the nonlinear parameters are separable. The technique of Golub and Pereyra is used so that a minimization algorithm only for the non-linear parameters is needed. The minimization algorithm of Marquardt has been completed with an eigenvalue analysis. In order to reduce the computation steps the inverses of matrices of the form A + 21 are calculated with the eigenvalues and eigenvectors of the matrix A. Of particular interest is the obtained convergence speed and the ease with which the method can be applied.

y=(p(~,Ex)= L ajq~(~;x)+,;Oo(Ex) j=l with ~ R " , ~ R k and n + k ~
The functions ~oj(~ x) are not linear in ~ One has then to minimize the following sum of squares of deviations: ,',(a-;b-3= .=

y~-~oo(~ x,)-

~°, aj%(~;

x~.)

(2)

If: q~ = the vector ~.Rm with/-component y~-~~o0(b~xi) and = the matrix ~R .... with (i-j)-element ~oj(b;x~), the nonlinear functional (2) can be written as:

,.1(~,~ = I1@-~ II2

INTRODUCTION During the last decade the techniques of non-linear curve fitting have emerged as an important subject for study and research. The increasingly widespread application of this subject has been stimulated by the availability of digital computers and the necessity of using them in complicated systems. The intention of curve fitting can be: (a) to verify the correspondence between a mathematical model and some experimental data (x~, Yi, i = 1, m); (b) to determine certain unknown parameters by means of a valid supposed mathematical model. As long as all the parameters a~(i= 1, n) are linear in the model, e.g. a model of the form:

(i)

(3)

First, it is proven that it is possible to transform such a separable problem to a minimization problem involving the non-linear parameters only.

THEOREM Suppose rl(~,b~) has a simple isolated minimum for b~eS c R k and that the matrix ~(b~ is of rank n with continue derivatives in this region for ~ If ~ is constant, the minimum of rt(d,~ is attained for ~'being the solutions of the set of normal equations:

(,:I)T(1))~= ,:I)T~

(4)

y = a o + a Iq~i(x) +... + a,,cp,,(x) the determination of the unknown parameters is not difficult. They are solutions of the linear normal equations. However, if the model depends also on non-linear parameters, the fitting is more difficult and needs iterative methods. If there are several non-linear parameters, the computer time can therefore increase very rapidly. In this paper a computer program is described to fit in an efficient way a non-linear mathematical model. Only one independent variable x is considered, but the method can easily be extended for several independent variables. It is supposed that the linear and the non-linear parameters are separable. A least squares problem is called separable if the fitting function can be written as a linear combination of functions ~0~(b~x) involing further parameters in a non-linear manner. Suppose the data (x~, Yi, i = 1, m) has to be fitted the model:

0141-1195/81/020055-0752.00

© 1981 CML Publications

for the corresponding linear regression. This equation can be solved for d as:

z,(~')=(OTO)-'OTq'

(5)

because of the rank of(OrO) is the same as the rank of O and thus equal n. After putting this result in equation (3) we get the non-linear functiofial: r2(~ = II0P- ~(~WO)- 1~@112

(6)

which is only function of the non-linear parameters ~.. The following theorem is now proven: If o~S is the minimum of r2(b-') and if ~ is given by equation (5)in which O _ Orb') is substituted, then (~b') is the minimum of rl(~b~).

Adv. Eng. Software, 1981, Vol. 3, No. 2

55

Proof One can always find an orthogonal m by m matrix Q, which reduce • to his 'trapezoidal' form:

,7, with U a triangular n by n matrix of rank n, and QjeR "'m and Q2eR "-"'r" (see Actona). As a result of the orthogonality of Q, the following relation between Q~ and Q2 is obtained:

QTQ = QTQ, + QV2Q2= I

(8)

and from equation (7) is found: Fn7 QT[oJ and q~r= [UT0]Q

~=

(9)

Using these relations one proves easily after some calculations that:

Otherwise one has always: min q(~,b-}~< min r ~( a-',b3 ~,b~s /;ss a'=aTb3 or with equation (16): min ,'I(K~< min r2(~) ,T.~s ~s

(19)

The inequalities (18) and (19) prove the theorem. The minimization problem of the non-linear sum of squares in the n + k parameters of the model (1) is thus reduced to a minimization problem in the k purely nonlinear parameters only. Models of the form (1) are used frequently and this reduction of the number of iteration parameters can sometimes improve enormously in computer-time and numerical convergence. Consider now the minimization algorithm for the remaining nonlinear parameters. MINIMIZATION ALGORITHM

O,Oo, oT0 9

"0'

To find iteratively the minimum of r2(b-) a so-called second order method is used. The non-linear functional rz( ~ in equation (14) is written as:

and I

- -

(~((I)T(I}) - I (I)T = QTQ2

r2(b)=~T~ with ~=Q2i~.R "-n

For the identity matrices I, we have denoted the rank n as I,, only where it was necessary. Out of the last equations it can even be proven that the matrices QTQI and Qv2Q2 are idempotent with power 2:

(QTO,)2 QTQ, (Q~Q2) 2 = Q2TQ2 =

To approximate the minimum value of the objective function r2(b--)from points ffnear to the minimum b~+ A--b~ the Taylor expansion series with second order terms is used:

,.2(b*q- A-b)= ,.2(b-)q- ~*. A-~q_ ~ T H A b *

(12)

Let us now put: g'= ~ - ~ d

(13)

~ - for2 "~eR 1.k 9-~-0~) (14)

r, (~,~ = -U~-- -~TQTQ~ = ~TQTQ, ~ + ~TQ~Q2~

(16)

The three terms in the above equation are certainly positive, so that: (17)

Choosing for h*= a~b-)= (q)T~)-I (1)Ti~/' the second term in the second member of equation (17) is zero, which results in: min rl(~,~>~ min r2(b--) (18) •.~s ~s

56

Adv. Eng. Software, 1981, Vol. 3, No. 2

(23)

(15)

following relation between these two functionals:

min r,(a-,b-)>1 min r2(b--)+ min 2TQTQ,-~ •.~s ~s ~s

" "~ERk'k

H = \c3b,Obj (i,j = 1,k)]

A combination of equations (14) and (15) results in the I'l (a~b--~= 1"2(b-~ "[- UQTQ, ~

(22)

and

( O2r

and the primary functional (3) is transformed as:

(21)

In this expansion 0"and H are respectively the Jacobian gradient vector and the Hessian matrix. They are defined by:

then the modified functional (6) is given by: r2tb~ = ~TQr2Q2~ = "~TQ2TQ2~"

(20)

(11)

Our objective is to determine the vector Ab of the moqement, required to approximate the minimum from To determine Ab-"approximately, consider c-)'and H as fixed and differentiate partially the increment:

0'.Ab +

T.H.Ab

(24)

with respect to Ab. Setting this result to zero gives: 9--'(b~+ H(b~. AT = 0

(25)

and after solving for Ab yields:

A--b= -(H(b~)- 'g-lb~

(26)

as the approximation for therequired movement to the minimum b~, from a point b near to the minimum (the current point). Equation (26) is fundamental to all second order solutions for a minimization problem. When it is used directly to generate successive movements toward a minimum from a given initial value bo, the method is known as the Gauss-Newton algorithm. Direct use of equation (26) is limited, however, because the Hessian matrix H must be computed and inverted at each steE of any iterative procedure. If the partial derivation of r2(b) is analytically too difficult to perform one can have recourse to numerical derivations, and then several methods to approximate H-x are available. Only one is mentioned. A first partial differentiation equation (20) gives:

o,=Or2=2z~02.'i = 1,k) 02,'2

0~ 02'

c7~2=

= 1,k)

o y 02.

Pk(A + 2klk)-1C~

0 ] Qm~

0

;;

(32)

o

Calculation of Oz/ab i

(29)

(30)

Finally, we observe the numerical convergence of the method. It can be proven that the matrix H, evaluated at the minimum is positive definite. However, H in equation (26) is not necessarily positive definite, since it is evaluated at a point other than the minimum, so that the process may not converge. This situation is most likely to occur at some distance from the minimum. Moreover, the Hessian H is approximated by the matrix A in equation (30). This is the reason why in some situations it is important firstly to limit the step size (taken as a fraction pk < 1) SO that a solution is not predicted outside the range of a valid first order approximation to H and secondly to add a positive scalar 2k to the diagonal elements of A so that A + 2kl is certainly positive definite. Taking this into account, the iteration scheme becomes finally: "ffk+ a =~k--

0

[10

(28)

where A is a symmetric matrix with/-j-element:

A~.~- Ob~Obj (i,j = 1,k)

=

0 ] [10_ 2 Q,?_ ]

(The zeros in the second member mean zero-matrices of different kind.) Applying the same sequence of transformations on the vector i~, we get the vector 2.as the last ran elements of this result.

The Newton-Raphson least squares procedure assumes that the second term in equation (28) can be neglected. Therefore we take: A--b= - ~A - 'O.

ixX]o

QO=[Io_ a

(27)

and a second partial differentiation yields:

Hi.- -= 9 "~ - - + 2~ r ~ - ( i , j "~ Ob~db~ -Ob iobj coicoj

formation Q"! transforms • to a matrix with on the first column all zeros except the first element, a second transformation reduces the second column to all zeros, except the two first numbers etc. This can be visualized as:

For the purpose of calculating gi and Ai,.i following equations (27) and (30) are needed the derivatives 02./0bv They can be calculated analytically if the derivatives &pj/c3bi are determined. It is of course also possible to determine them numerically. The last procedure is less laborious for the user of the program, because no supplementary program is needed to calculate &pj/Sb v Yet, it introduces a supplementary inaccuracy in the calculations. I]' the functions rpj are sufficiently smooth, it appears that numerical derivation yields no serious risk to lose convergency, but one has to be careful if the functions ~oj are liable to error noise (e.g. results of numerical integrations )

Scaling of the matrix A in equation (31) In minimization problems of the sum of squares of deviations the covariance matrix A in the normal equations is usually scaled to obtain a correlation matrix. In this way the normal equations are better conditioned. If A*= T A T

(33)

with T a k by k diagonal matrix with diagonal elements A?.y 2, then A* is the scaled correlation matrix. Using equation (33) the increments Ab from equation (29) are obtained: .____..+

a-if= - ~T(A*)- ' T0'

(34)

(31)

The practical determination of the scalars Pk and 2 k will be treated in the next paragraph. (We use the index k to mean the kth iteration step. There will be no confusion with the number of non-linear parameters also denoted as k.) COMPUTATIONAL PROCEDURE

In order to obtain also for A* a positive definite matrix, we add to A* the diagonal matrix 2kl, and finally equation (31) is reduced to:

~k+l =~k--pkT(A* + Afl)-ITO.

(35)

Determination of Pk and 2 k Following equation (29) once can choose Pk 1/2 and if necessary a more appropriate Pk can be determined with the method of successive division by two, or another line search method. Some authors 2'7 have treated this problem in detail. Here we will try to reduce the number of =

Calculation of 2.= Q2~ The determination a of the orthogonal matrix Q, which reduces • to its trapezoidal form, can best be done using a sequence of Householder transformations. A first trans-

Adv. Eng. Software, 1981, Vol. 3, No. 2

57

evaluations of r2(~ 'rather than the number of iteration steps. In order that A* +2kl should be positive definite, it is necessary and sufficient to choose 2k>maximum (0, --e), with e the smallest eigenvalue of A*. It is, however, no_~tnecessary to take 2k tOO great, because if 2k--' O0 then IIAbll--,0 and the convergence can be slowed down. In his algorithm, Marquardt 3 has not used an eigenvalue analysis. He staets with a small value for 2 k and takes Pk =0.5. If the kth iteration step successful then 2k+, is set equal 2k/V (V> 1), if not, the angle between Ab and ~'is calculated. If this angle is lower than 7r/4, then Pk +1 = Pk/2 2k + 1 = 2k

(36)

is taken. In the other case: Pk+ 1 =Pk

(37)

"/'k+ I ='~'k ~'

In our program the eigenvalues are determined with Pk =1 and 2k=maximum (0, --e) with e the smallest eigenvalue. If one iteration step is not successful, e.g. the kth step, Pk and 2k are transformed following equation (36). But if Pk becomes smaller than 0.03, the direction is changed, with Pk÷~= 1 and 2 is transformed as follows: 2k + 1 = 22k(if 2k 4: 0) 2k+ I = e (if 2k =0)

(38)

The inverse of (A* + 2kI ) Let the eigenvalues and the eigenvectors of A* be denoted by ej and vj (j = 1.... k) respectively, then we have:

Stopcriterium If r2(b~k)is the value of r 2 after k iterations, and if k + 1 is a successful iteration than b~k÷l is taken as a local minimum if:

,

2(b'k÷ 1) ,.:(b-.k)
(42)

with e a predefined small number, given by the user of the program. In order to take action in a case of nonconvergence, the user of the program must also give a value for the maximum number of iterations. LISTING OF THE COMPUTER PROGRAMS The above algorithm can be implemented in a computer program as follows. First, a main subroutine is needed which determines the minimum of r2(~ and the corresponding parameters fi' and b" of the regression model. For this subroutine which is named 'BBO3M2' a flowchart can be drawn as shown in Fig. 1. The variables have been labeled as in the text as much as possible. This main subroutine requires two other subroutines. A first one which is named 'BBO3EI' determines the eigenvalues and the eigenvectors of a real symmetric matrix, and a second one which is named 'BBO3TH' reduces a general m by n matrix to his trapezoidal form using a series oforthogonal transformations. In the flowchart the subroutine 'BBO3TH' is needed in the boxes denoted by '*'. This subroutine 'BBO3TH' as well as the subroutine'BBO3EI' are, of course, standard routines; they have been copied

k

(A*)-' = ~ e f 1viva.

(39)

Calculate kk

Initialization

j=l

The eigenvectors of A* +2kl are the same as those of A*, but the eigenvalues of this matrix are given by ej+ 2 k, so that:

Ak • max (O.-e)

k=0

e smallest eigenvalue of A*

Calculationof the initial sum of squares

,,"usinOeouat,on,,O, '^ I I

I

*Modify the non-linear parametersconlorm (35) and calculate the corresponding sum of sq. . . . .

I I [

!

k

(A* + ~.kI)- 1 = ~ (ej + )~k)- 1VjuT

(40)

j=l

~1~ Printinitial / )arameters and l initial sum of / squares

If it is necessary to test several values for the parameter )~k (see above), this procedure seems to be very interesting, because we don't need any matrix inversion.

/

I.

~ I

I"

~

. /2

r k " ~'kl

/

TYES

~YES Calculation of the linear parameters After the minimization of r2(b) the linear parameters ~'(b) given by equation (4) can be calculated in an easy way 7. Using equation (9), equation (4) is transformed in:

E T01Q TE I ;E T0 Q

I

*Numerical

Oetermination of the

YES =l linearparameters equation41)

A* (equation 33}

(41)

Adv. Enq. Software, 1981, Vo/. 3, No. 2

Figure 1

Printfinal I \ parametersand / \ calcu,ated data/ points / \

NO

Because U is a triangular matrix, this set of equations is easily solved for ~..

58

1

matrix A tO obtain

or

U~=Q,~

Print newparemlers /

calculation of

end A (equations 27 and 30)

from the Scientific Subroutine Package (SSP-library) of IBM. All programs are written in Fortran IV.

1 2 C 3 L A L 5 L 6 C 7 C b L 9 L 10 C 11 C 12 L 13 L 1A C 1SO lc, ~. 17 C 1~ C 1~; L 2(. L 21 L •2 C 23 L 2A C Z~ C 2L L 27 L

3O 31 3~ 33 3A 35 3D 37 3~ 3S i,~,

C C C C L L L L L C C C

SUEkBUIlNE

Bbb~h~{A,ALINeF,Ye~eNeR,PIleEPS)

SUBkTUIINF

ObU3H~

PUDIOSE 10 PEPFOR~ A GENERAL LEAST SBUARES E S I I N A T I O N PARAqEEERS USAGE L4LL

OF N O N - L I N E A R

OESCRIPFIO~ ~F ThE pARARETERS A : DOUSLE I ' k E C I S I O M VECTOR OF D I q E N S I O N ~e ONEEH CONTA|NS THE I N I T I A L NONLINEAR pARAMETER E S T I R A T I D V S . ON RETURN T H I S VECTOR C O h T A I H b THE F I N A L N G N L I N E A k PARAHETER E S T I N A I Z O N S o A L I N : DOUbLL P R E C I S I O N OUTPUT VECTbR OF D I q E N S I O N H , C O N T A I N I N G THE L I N E A R PARAHCTERS X : DOUBLE P / E L I S I O N VECTOR Of D I R E N ~ I U N H . OHICH CONTAINS THE G A l A UI THL INDEPENDENT V A R I A B L E T : DOUBLE EULCZbIOM VECTOR OF OZqENS|ON q.MHZCH CONTAINS THE DATA UE 1 H I DEPENDENT VARIABLE M: NUMBER Ot O~T~ P O | N T S ( e A R I ~ U ~ 5 0 ) h : NUNBER OE L I N E A R P A R A N E I E R S I R A X I q u M 1 : ) K : HUNGER bE N O H - L I N E A R PARAMETERS N I T : R A X I N U h NUROER OT I T E R A T I O N S ALLORED LPS : DOUBLE F R L L I S | O M I N P U T PAR~RETER. 1RE I T E R A T I O N STOPS IF THE R E L k T I V L CHANGE~ENT OF ThE SUM OF SQUARES I S LOVER THAN L P E . SUBROUTINES k L U b | R L u bSO3UF : h~S~ UL FURNISHEO NY THE U~Ek b303El CALLULATES TNE EIGENVALUES OF A REAL SYMMETRIC MATRIX BOO3TH : T R A N b l U R M S A q A T Q I R TO ~ZS TRAPEZOZDAL FORN°USING hOUSEHOLDER TRANSFORqAEIONS PLHARKS 1 ) FOR PROULEHS bHERE q | S GREATEk THAN SO OR N IS GREATER IHAN I ~ . O N L T THL R | q E M S I O H S T k l E R E N T S HAVE TO BE ADAPTED 2) THE PRObKAh Lb V R I T T E N I N ROUBLE P R E C I S I O N

IMPLICIT REAL*b(A-H.O-Z) QIWC~SIOH A ( 1 ) . k L I N ( 1 ) . R ( 1 ) . Y ( 1 ) QI~LNSIU~ PSI(~)oFbIA($O)°AJ(SO#IO) RI~LVSIBh F((5Go11).FINUL(~O).A~(60) GG(lO) oE(10).Ah(1U).k(lOO).DA(10) EOUIV~LEOCE(FII].ll)oFINUL(1)).(R(1)oAJ(1.1)) :

D1 ~2 g3 L4 DD ~ ~? D~ ?C ;'1 ?2 73 ;'A 7S 7D ??

¢ C C ~

tC ~3 EN ['. Ee LL N b~ L 9C L

~3 S~

1L~ L 1~,1 7 1~.3 1L;4 TU~ 1LO 1b~

I.

111

11

113 114 115 11(. 117 11~ 119 1;.L IZ1

L

123

Of

THE R E R | V A T Z V E S ( E Q U A T I O N S

(3.8)

FSIA(I):Y(I)-FIH~L(I) CALL B ~ O T T H I M . N ° I I . F S I A ) DUb |UkloM |lm|-N AJ(IE.J)=IPSIA(I)-FSI(I))/DEL G(JIs~(J)~AJ (IT,J)tFSI(I) LbhllNU£ OF T h L HATkLX A(SEE

EVUATION C A . 2 ) )

LL=~ BO 1 0 6 I = 1 , K LI=LI~] 1F(AR(LI).RE.0.)

ANb ( 3 . 1 1 )

OF THE RATRIR

THE

1E3

15

1~? 1~ 1oS 17~ 171 17L 173 17& 17~ 17C 177 17b

10 C C C L C L

40 27 ,~

~1

2LG C 2~1L ZL~ L 2 9 3 23

3E ~3 ~9 C L 39 9~

3S c~ 95 3G 107

PROCESS HAS 5 T O P P E D I , I ,

bETERHINAVT

Of

AN

INVERTING

MATRIX

IS

ZEROI*/B

INITIAL

VALUE OF L A H B O A ( V A R I A G L E 8 L )

BLm-ENIN

CALCULATION ~ f THL LNVERSE OF THE SCALEQ MATRIA Aq USING THE EIGEVVALUES nhO THE EIGENVECTORS(SEE EQUATION ( E . 9 ) ) O0 12 [ = I . R SO~g0. DO 15 J = I . R SUHaO. 00 14 L : I . K LI:¢L-1)*K+I LJm(L-1)*K+J LLILeIL~I)I2 SUNEEUM#R(LIIaR(LJJ/(Aq(LL)+RL) SUH:SUM*G(J) SOmwSTm~SUn GALl)e-SON IF(E(I).EQ.S.) bO TO 12 bA(I)=-SO~I~(I) COkTIVUE H O L I F I C A T I O N 01 THE NONLINEAR PARANET~RS &NO CALCULATION Of CORRESPONDENt ~bF. LF S O , A R E S ( S E E EQUATION ( 4 o L ) )

THE

00 1~ l u S , g AN(I)sAflI+DA(I)*STLP CALL U b ~ T U F ( A N ° X . t I . R o H ) DO 15 I B l . q PSI(I)'¥(I)-ELNULII) CALL 3 B ~ S T ~ ( R ° f l , F I , F S I ) SS~b=3. b0 15 I ~ N 1 , P SSkbuSSGU~PSI(L)*t~IEI) IL$1|H6 I f ThE I T L k A T I O N IS SUCCESSFUL OR NOT ANO Z f THE XAA|NUM NUNbEQ OF I T L R ~ T I O H S OR TOE N I N I q U ~ I S NOT TEE REACHED. I f kECESSARY k h ~ ( V A R I A B L E S I E P ) ANN L A q S O A ( V ~ R I A N L E E L ) ARI C~AR~£D ZF(SEQB.LT*SSu) bUTU 17 |FISIEP.r.T.~.U3) GLTO ~6 SPIRe1. BL:bL~bE IF(uL.EN.O.) bL~¢hxH BOTh 3E. STL|:SEEP~b.~ GGT~ t ~ VkITE(3.2~)S~b. (AN(1).ISloC) F O P R t T ( I . ~ S L C C L ~ S l b L I T S R & T I O N o° 11°1~X. ~ SUq~bt S~AkES~°°ElS.~.IOX. eI°I~.'NOH-~IHLAK I~RA~ETERS:'~SEIS.~,I°32Yo~ElS.S) DO 22 1 ~ 1 , ¢ A(1)SAN(I)

IF(SSOH.LT.I.Q-7) ~ TO ~0 IF(D4GS(SSU-bSbU)I~$R.GT.EPS) GO TO ~O VR|TE(3.?7) FORMAT( t I T E N A I I ~ H I, ORMSLLY EMO£D I) GOTL ~3 NIT:NIT~I IF(HIT.bE.PI1) b b l O 21 SSOuS~uS STEFul. 60TO 2~ URZTE(~.cE) F O R ~ L T ( t THE H A E I N L h NUMBER OF I T E R A T I O N S DETEEMIhATIOh

OE TPL

FINAL

LINFAR

IG REAIHED

RAkAMET[RS(SEE

o)

EMUATION E L . T 0 ) )

IF(N.EO~Q) GbTD 3~ DO 31 l ~ l , H QL~IZII.I) PSI(I)=~SI(I)IOL DO 51 J ~ I . N I|(I.J)~FI(I°J)/DL 00 3~ I I a l . N I:N-||~I

&LIN(I)ePSI(I) IF(I.EQ.I) GOT~ 33 BO 52 J ~ I , I M 1 PSL(J)~PSI(J)-ALIN(1)iFZ(J.I) CONTINUE URITE(3°IV) |ALLH(L).Iul.N) FORHAT( ~ L I H L A k P A R A ~ E T E R S S • o O E 1 5 . S * I . 1 9 1 . 6 E 1 3 . A ) OOT|UT OF THL b A | A - t 0 I k T S CALL U U O 3 U F ( H ° X ° E | . h ~ N ) NRlTEI3.99) FORHAT(III,~)o I X - V A L U ~ S l o l A X o | T - V A L U E S ° , A ] • 1 IESE|~ATED T-NALULSI,12X,tDEVIATIOqS°,I) DO ~3 I ~ l . N SOM~FINUL(I) lF(h*EO.O) GbTO ~5 DO 35 J ~ l , q SORuSGM~ALIN(J)*FLCI,J) OEVmT(L)-SON vRITE(3.30) A(IIoT(X)oSOq.OEV TORhAT(E13.3.3122.3) RE|bEN ENR

1

SUbl.0UTINE RL~3L|(~*R*N,qV)

3 C 4 t.

.....................................................

., ......

. ....

E E (

SUbROUIINL LbL~bL

CALL P B @ S F I f A N , k ° A I U ) bLICkHINARI

1E

226 22? 2Zb 229 23Q E31 Z3, 233 23N 23S

L A L L J L A T I O N Of T h L ~IGENVALUES AND E|GEMVECTTRS OF THE SCPLED N A T R I X ~R

THf

l~v

21e 217 218 219 22C 2~1L 222 223 22E

IS=C DO 7 I = I o K DO E J = l , ! IS©lS+l SUM:~. DO ~ I N = I , H N EUK=SUM÷AJIIH.|)~AJ(INoJ) ~(ISI-SUP EII)zDSGRI(Sbh) lS=0 DO 11 I = l o K D0 1~ J = l , l IS=IS÷I IF(L(I)*E(J)*LQ.U.) GO TO 1 0 AM(IS):AM(IS)IIE(1)*E(J)) CbHll~uE IF(E(I).(Q.~.) ~0 T~ 11 G(I)'~IJ)I~(*) COH~INUE

1E$11E

1~ C C L C

~15

L L C C

13

~OS 2U~ 2C? 2 0 8 31 ZC9 Elb 211 21L 213

UO 4 J : I . K O:~(J) 6(J)aO. DEL:S*O*O~S lt(OEL.EU.O.~) LELuG.OOS A(J)IBtD~L CALL B B O T U F ( A . X o F I . h o N ) A(J)IO

SLALIYb

C C C C 3~

1~3 15A TBS lS~ 157 1SB

1GO 19L 1S1 19~ 193 19~ 19S IS~ 197 19b

L A L C U L A I l O N OI ThE L N Z ] I A L SUN OF SQUARES IALL ~HG3UF[~.~.EIoH.N) BO 1 I x l . H FSI(|I=Y(1)-FLhUL(|) CALL ~ o ~ 3 1 H ( h . h . l / . P S l ) ssu=n. NlzH÷l STEFI1, DO 2 I = ~ I . M SS~:SS~+PSIEI)~-2 bRIIE(3oT) S~.(k(II°lal#K) F D k H A T ( I I . I L H L I ~ A L S U , OF S Q U A R E S I I . E l S . S . */.I~ko'NON*LII;LkR EARAqETERS~.SE15.5.1~32SeSElS.5)

o

TRY 0 1 f l L R START/NG V A L U E S ° ) GO 1 ) 1 0 7 COkTINUE •

ERINsAH(IS) BLwO. ITIEqZN.LE.O.O)

I~S l b ~ S2 1B? 2g

MUUNER OE 1 1 L R A T I O V S

NUHLRICAL LALLUL~TLbR

ld6 C C C

1LU ~

klT:~ S~ L '-7 L 3g S~

1 2

1E2 17 1~3 ~

~C hi1

150 lsl

NRITE(3°IOO) EURRAT( ° TN~ L I E R A I Z O N

1U0

1A3 1A{

HR~TkII~.~LZNIX.YoRIN.R.MJEoEPG)

4~ C A3 L DS L A¢ A? AP 4~ ~C

I~A 125 12o 127 128 123 130 131 13c 133 13A 133 136 137 138 139 14L 1A1

AH.[S

MOT ZERO

t~ 7 E $ 1L 11 1; 13 14

L L c C C C L L L

E

PURFUSE CU~PUTE ELGLhVACUES AHD E I ; E N V E C T O R S OF A RE~L SYMMETRIC MAIRIY

E E

USAGE LALL U B O T C L ( A . R * h . RV)

E

£ESCRIPTION

Of

| !

PARAMETERS E

GO TO 106

lc, L 17 L

A -

O R I ~ I A k L P.ATRJX ( S Y U M E T R I C ) e DESTROYED IN CORPUTAT(ON° RESULTAH1 EIGEHVALUES ARE REVELOPEO I N DIAGONAL OF

Adv. Enq. So/hvare, 1981, Vol. 3, No. 2 59

| E

IG ( 1V C L~ t ~Z L ~3 C 2; C ~ C 27 C ~E C ~$ C 3~ C 31L ~ C 33 34 L 3~ C 3~ 39 C &C

R q ~V-

~ A T k l k A XN DESCENDING ORDER. RESULTANT ~ A T k i A OF EIGENVECTORS (STORED C O L U R N H I S E , I~ SAhE SERUENCE AS E I S E N V k L U E S ) ORGLE OI MATRICES A AND R IN]El EbbL U c o q P U T E EIGENVALUES ANR EIGENVEETORS 1 LU~PUTE EIGENVALUES ONLY (R NEED NO1 BE UEnENSIOREO OUT RUST S T I L L APPEAR | N C A L L I N G SkqUENCE)

kEUA~KS O R I G I N A L F;ATRI) A MUST RE REAL EYqNGTRIE (STORAGE RODEOS) M&TR|X A LAhNOT BE I N THE SAME LOCATION AS ~ A T R I S R SUUROUTINLS Ahb NONE

FUNCTION SUBPRRGRAMS

I"

RERUIRER

t.ETHOR b l A U R N N b I L A T I U M PETHOq O R I G I N A T E U bY J k C O B I AkR ADAPTED HY VON ~LUKAhA FOR LARGE COMPUTERS AS FOUNR I N tNATNEMATICALE qETHORG F b k G / G I T A L C O N P U T E R S ' , ~OITED 8Y R . RALSTON ANO [ H.S. blLF, J O h n b I L E T tNO SONS, ~EH YORE, 1 9 6 2 , CHAPTER ? [ . . . . . . . . . . . . . . . . . . . . . . . . . . . .

***.*

.............

*...

...........

. ....

F ! F

IMPLICIT RFAL*b(A-h.O-I) OIMLqE|Oh A(1).k(1)

4~ kS A~ 4? ~L

~ C E C

I !

.......

. ..............

~EHFRAT[

. .................

IGLt*TtTT

.*° ..........

161 lS~ lb3 154 1S~ lEE 157 158 15~ 1DU 1~1 16~ 1O~ 164 lOS lEE 1~7 lag 16S 170 171 ~TZ 173 17G 17S 17b 177 17b

C £

COMPARE THRLShULO t I T H 166

IE(T~R*ANRqR)

1~$

IO=-~ O0 l q S

C C C

F I N A L qORq

1D~.1¢5.45

S3RT E I G E h S A b U L E RMO EIGENVECTORS lml,N

IV=IR÷N

ITb

175

LLnI~(I*I-I)IZ JOmA*(I'2) RO ~ES J s I . N JQ~JQ~R NNmJtIJ*J*J)12 IF(A(LL)-A(~H)I 17~.18S.18~ X:A(LI) A(LL)mA(~M) A(~M)mX /F(MV-1) 17S,I~5017~ bO 1SO N = 1 . N

ILPulO+[ IHR~JR+K XeflIILR) H(ILRI=R(|MR) IE~ R(IMRIsX 1bE CONTINUE RETURN END

.***.****.i

WATRIk [

RAh~Eml.0R-I~ II(h¢-l) 10.~.1~ IG:-~ DU ~ J J l . N

~3

SUBNOUTINE B b U S T H ( R , N , F I , P S I )

1 2 C 3 £ A C

..* .........

.........

. ....

...*..*.*.*.*..*

.................

* ....

,

SUBROUTINE U b Q S T h 5~

GO ~0 I : l . q

D S 7 C

IJml}*l If(I'J) 2~,15,c, $(IJ):1.~ ~L CUNTIRI,~

cI 01

¢

~2 C 63 ¢ o4 66 b8 ~9

70 71 72 73 T4 ¢ 7S C 7~ L 77 7B 79 8e bl 8Z

COMPUTE I N I T I A L

AND F I N A L

NORMS (ANORR AND ANORRE)

2S AHORR=O.Ü DO 35 I m l . N DO 35 J u | , N IE(l-~) 30.3~.30 39 IASI*(J*J-J)I~ ANORReAMORM~A{|A)tAIIA) 35 CONTIHUE |F(RNORN) 16S.16~.AG ~0 A N O k q u I . A 1 E * U S G N I ( A N O R R ) ANG~XUANORMeEAkbE/FLOAT(N) I N I T I A L I Z E |hOICAFORS

AND COMPUTE THRESHOLR. TNR

8O

G7 ~L ~9 91

93 94 95 ~ 97 98 9V lbl C luZ C 1L3 C 104 1bE 1;7 1Cb 1C¥ 11~ Ill 11E 113 114 I~S

C L c t E

TESIS

FOR C G h F b ~ T l O N

TEST Tog R •

LAST COLUMN

13U ] F i b - k ) 135.140.13S 1~ quN+1 GO TO 6 3 . t C C

TEST FOR L • 14U I F ( L - I N - I ) ) 1 4 5 LmLe~ GO TO 55 1 5 ~ IF L I N D - l ) 1 5 5 IRRw~ GG TO 5~

D E S C R I P T I O N Of THE pARAMETERS F, : RUMGEG U] kObS OF RATRIR E l h : HUNGER bF COLUMNS Of RATRIR F l FI : OOUBLL P R b C I S I O N SO BY 11 INPUT M A T R I X * ON RETURN E l CONTAINS THE TNAPEEOIOAL FORH P S I : DOUBLE F R L C I S I O N VECTOR OF DOMENSiOq SOnON ERICH THE SAqE |RANSEURMATLQNS ARE PERFORMED

C

L ¢

L L £ C

SUEkOUTIHES kLGUIRLG NONE

REMARKS 1) FOR PRUNL(~S NNERE q I S GREqTER THAN 5 0 OR N I S GREATER THAN 1 1 , ONLY ThE D I q E N S I O N STATEMENTS HAVE TO BE ADAPTER ~) T~E PRObkAP I S N R I T T E N I N ROUELE P R E C I S I O N

l I

2~ 3~ 31 3~ 33 34 35 36 ST ~E 39 4C 41 42 A] 44 A~ 4~ 47 4E 4~ ~ ~I

!

53

¢ ;

! I I ! I

£ C ( t C C £

C L IMPLICIT REAL*~(A-N.O-E] R1MENSION F l ( 5 ~ , 1 1 ) , V ( 5 0 ) , P S l ( S O ) IF(N.EQ.O) RLEURN bO 1 I = l * N SUNmO.

bO 2 J l l , ~ SUPmSUMAFI(J,I)tTI(J*II SgGSORTISUMI IF(FI(I,I)*LT.O*I Sw-S TkIIS*(SAfI(l,LJ) IF(TKI.EQ*C*) bL Tb I DO 10 J = I , M V(J)~FItJ,IJ V(l)lw(I)~S T~G. GO 5 J ' l . v TITAW(J)*PSI(J) TmTITK2 RO ~ J m J . ~ ISI(J)=PSI(J)-kIJ)*T DR ~ J = l , ~

1~

~

IdG I L R n I L O + I IqFmlqO+| XmR(ILR)eCOSX-RIIRk)eSINR R(IMR)=RIILRI*SIhL~A(IPR)*COSR R(ILR)=X 1 2 5 CONTINUE XmI.0~A(LM)eSINLS TmA(LL)*CUSXI÷A(Nk)eSINA~-~ XmA(LL)*SIHX~*A(~N)tCOSRE~R A(LMIm(P(LLI-A(NR))eSINCS*A(LNIa(COSXE-SIHRE) A(LL)U¥ A(MR) mR

126 127 IZE

(

USAGE CALL a O O 3 T H t R . N , F I . P S l )

C

~5 ~C ~7

1U-(l*l-l)12

1L1 122 1 23 1L4

60

M CCLUMN5

IF¢I'L) EOo115.EL b~,115,Gb 5~ IflL-q) 65 l q : l + M b GO TO 9~ SO IM=M*Id lqS,TLG.1L~ S~ l f f l - L ) lug IL~I*LR GR TO 1 1 0 ILmLtIQ 11L Isl(IL)*COSA'A(L~)*~INI A(IR)=A(|L)*SIhk*A(IR)*COEX A(IL)SX IF(KV-1) 12~,12~.12C

117 11E 119

130 131 13~ 133 134 13b 13e 137 13~ 13~ 14C 161 162 143 144 1¢S |AD 147 l¢b 1;9 l~b

ROTATE L k h b ]L~mR*(L-F) IMUwN*(M-1) O0 1~5 l m l . N

1U 19 ~U Z1 ~2 ~3 ~4

PURPOSE TO PERFOEh A GLhERAL N BY N M ~ T R | I f l TO H I S TRAPEEOIDAL FORM, USING HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS* THE SAME TRANSFORNATIONS ARE PERFORN[R ON THE VECTOR P S I *

hETh0D HOUSE~OLRER URThOGOHAL TRANSFORMITIONS

COMPUTE S I N AND SOS 60 NR=(qtM-M)I2 LOm(LeL*L)I2 LMOL~MQ 130.65.65 62 I F f b R B S ( A ( L M ) ) - T h R ) 6~ IMC©I LL~LtLO NN=N*RR Xm~.S*(A(LL)-A(NFI) 58 Y = - A ( L M ) / R S O E T ( ~ I L R ) ' A ( L H ) + E ~ P ) 1FIX) 70,75,75 70 TU-T 7S S I N X m T I D S b R I ( L . U * ( 1 . 0 ~ f R S O R T ( 1 . 0 - Y e Y ) ) ] ) SINAITSIHX*S~NX 78 C O S X a R S U R T ( I . O - S I h E ~ ) COSk~mCOS)*CbG~ SINCE ~ S I N k * C G S ~

C C C C

2~ L 27 C ~G L

INRmO ThNmSNORN 45 T N R m T H R I F L O A T ( ~ ) SO L = I SS N m L * I

84 L 8~

9 1~ 11 12 13 14 15 Ib 17

SEL~ND FROM LAST EOLUNA

165,15L,1A$ 16u.lS~.lbO

C

Adv. Eng. Software, 1981, Vol. 3, No. 2

~9 &~ &l 6Z

bOb iwl.. T=T~d(K)*FI(K,J) TmTITR2 DO 7 R U l , n F|(L,J)mFI(K,JI-k(N)*T CONTINUE CONTINUE RETURN EN~

63 64 7 65 66 1 ~7 U bE

/

AN EXAMPLE

For each problem the user of the above subroutines has to furnish himself two Fortran programs. The first one is the main program, in which the subroutine 'BBO3M2' is called after that the data-points (x;, Yi) are read in or generated. The second one is a subroutine BBO3UF (A,X,FI,M,N), in which the terms ~oj(~ x) of the fitting function q) are generated for given non-linear parametervalues ~.. In this subroutine the first argument A contains these non-linear parameters, the second argument X contains the x-values of the data-points, while the last two ones M and N respectively denote the number of datapoints and the number of non-linear parameters. In the elements FI(I,J) of the matrix FI (50,11) the terms (#j(~ xi)

with corresponding indices have to be stored. The special values ~oo(b; x~) have always to be programmed in the elements FI (I,11). If there is no function ~Oo(~; x) the elements F I (I,11) must be put zero. The above described program has been tested for several regression models. We mention only one. As example 50 data-points (x~, y~) have been generated from the equation: y = a~ + a2th[b~(log x -

b2)-]

(43)

with parameter-values: a~=200,

a2=150,

b~=3,

b2=l

and with: xi=0.2, 0.4, 0.6 . . . . . 10 This regression model contains two non-linear parameters CUt,b2) and two linear parameters (a~,a~) with corresponding functions:

x-b2)']

~P2(~ x ) = t h [ b i ( l o g

(44)

There is no function ~Oo(~ x). Starting from the initial estimations (7., 2.) for the non-linear parameters (b~,bz) the proposed program has been used to fit the model (43) to the generated data-points. The main program and the subroutine B B O 3 U F which are needed for this problem are listed below. 3 4 5 8 10 1~ 11 12 13 14 1 2 3 4 5 ¢ 7 1 V

PROGRAN uRO]vk INPLICLT RERL*b(A-H~O*Z) DIkENSION X ( S U ) , T ( 5 0 ) , A L I N I Z ) , A I Z ) ~5~ Kl2 N~2 O0 10 I ' l , R X(I):I*O.~ T(I)=ZOO.U*I~.b*DTkNH(S.ODO*(RLOG(R(I))-I.0DO)) A(1)'7. A(2)~2. CALL BBFlS~2(A.ALIh.X#Y*N*N*K*30*O*O~05) STOP ENP

RESULTS OF THE PROGRAM 01

0.~07360

O0

SUCCESSFUL ITERRlION SUN OF SOUARFS= 0.~11440-0~ NON-LINEAR PARANFTFRS= ~.29998D ~1

0.100)00

01

SUCCFSSFUL IIERRT|ON SU~ OF SQUARES~ 0.?~10-07 NON-LINEAR ~RUR~FTiPS s C*3~300R Ol ITERATION NORNRLL¥ ENDFR LINEAR PARAqETERS: O.?~uCO0 O~ ~.15~0C0 C~

I-VALUFS d.2Cb 0.400 0.600 0.8~0 1.bOO 1.200 1.400 1.600 1.8U0 2.riO0 2.200 2.400 2.~00 2.80U ].O~O 3.2UU 3.40~ 3.6U0 3.qUu 6.000 4.ROU 4.4uU 4.60~ 4.~U • .~GU ~.20~ ~.40') S.6~C s.RO0 6.009 ~.2u0 0.400 ~.600 6.~00 7.000 7.2b0 7.~00 7.609 7.800 8.000 0.200 R.4UO ~.6~0 R.UCU o.O~O 9.20~ 9.&00 QoQOG 9.~;0 10.000

0.1UOOGD P1

Y-VALUES

ESTI'ATkO Y-VlLUES

~.9~0 ~0.003 sC.~35 ~0.195 ~.742 52.2U& 5~.4~7 h1.978 73.~?~ 91.07o 11¶.P15 14~.NTp l~C.00~ 213.294 44~.1~5 ~6P.C67 2P7.878 J~.~Q3 ~14,555 323.131 ~?Q.EAtP v~4.199 ~37.747 ~C.42~ ~2.449 ~44.C01 34~.107 ~4F .12o ~4~.ES4 ~&7.&28 ]47.P84 x4~.~4Q ~AF.54~ x6~.781 ~.o75 ~4o.13& ~&9.265 ~49.~73 x~Q.46& ~49.539 ]49.~02 $4Q.~56 X4Q.701 ~ ~Ao.74L ~4~.772 349.2C1 ~4Q.R25 ]~Q.E45 349.e63 ~9.N79

5G.gUO 55.0~3 50.035 SO.19] 5U.T&2 5Z.ZUE 55.497 61.978 ~3.32~ 01.076 115.~15 146.410 18C.098 ~13.29A 2~3.125 258.067 287.NTR J~].OQ] 31&.$~5 323.101 329.460 334.199 ~32.747 340.42~ 342.449 ]44.001 345.197 346.126 ]~6.E54 ]&7.428 347.884 ]4b,~49 ~48.543 $48.7~1 3N~.975 3&Q.134 ~&9.~65 340.375 349.464 349.$39 349.611~ T49.656 ~49.?01 ~49.740 349.772 ]49.801 ~&9.O2~ 349.845 349.863 349.O79

OEUIATIONS 0.000 0.000 d . ~O 0 O.UO0 0.090 O.O00 0.000 -0.C00 -0.000 -0.000 -O.OOC -0.000 -O.O00 "0.000 0.000 0.OUC 0.000 0.000 0.~00 O.O00 0.000 0.000 0.000 0.000 C.qO0 O.O00 11.000 0.000 -0.0~0 -0.000 -~.GO0 -o.no0 -O.OUO -0.000 -O.OO0 -0.000 -0.000 -0.F00 -0.000 -0.000 -0.000 -~*000 -0.000 -0.000 -~.000 -0.000 -0.000 -OoGO0 -0.000 -~.000

From the computer results one can verify the very rapid convergency from the chosen initial estimates for bl and b 2 towards the exactly optimal values. N o t e also that th.e linear parameters are updated only at the very last iteration step, so that an initial choice is not needed for t h e procedure. Of particular interest is the obtained convergency speed and the ease with which the method can be applied. This has certainly two reasons: first the linear and the non-linear parameters have been separated and second the number of computation steps are reduced due to the fact that the inverses of matrices of the form A + 21 have been calculated with the eigenvalues and eigenvectors.

SUBQOUT|NE BbC3UF(A.R.FI~N.N) 1NNLICJT HEAL*8(A--H.O*E) RINENSION A f 1 ) o X I 1 ) , F I ( 5 0 , 1 1 ) PO 1 I 0 1 . u II(l.11)=U. F1(1.1)=1. FICIoZ)BOTANN(k(1)*(RLOE(XI[))-A(2))) RETURN END

I N I T I A L SUN OF SQUARESs G.51216D 06 NON-LINEAR PARAaETSRSu 0.70000D

SUCCESSFUL I I F R R T I ~ SUq OF SOUARFS= 0.12970D 02 NON-LINEAR pARAqFTrPS= O*~9&blD 01

O.ZOOUOR 01

REFERENCES

SUCCESSFUL ITERATION SUN OF S~URRES= 0 . 5 5 1 2 0 0 O6 NON-LINEAR PARAqETERSm -O.Z1439D CI

0.194690

fll

SUCCESSFUL ITERATION SUN OF SRUARFS~ ~.115150 06 NON-LINEAR PARANETERS~ 0.206730-01

0.2057&~

01

SUCCESSFUL ITERATION SUu OF S~UkRSS: L . ~ 7 2 ~ 3 0 05 NON-LINEAR FARR~ET~RS~ 0.310190

O0

0.2bS4~0 ~1

SUCCESSFUL ITERATION SUN OF SQURRFS~ C.845910 05 NON-LINEAR PARAqETERS= 0.428230

1 ALton,F. S. Numerical Methods that Work, Harper,New York, 1970 2 Bard,Y. Comparison of gradient methods for the solution of nonlinear parameterestimation problems,SIAM J. Num. Anal. 1970, 7, 157

OO

0.1062~0

01

SUCCESSFUL ITERATION SON OF SQUARFS= C*695580 05 NON.~.INEAR PARA~E]ERSa 0.771460 00

0.630260

00

SUCCESSFUL ITERATION SO~ OF SVUAN|S= ~.$16360 Ob NON-LINEAR PRRA~ETENS= 0.160300

01

0.126090

~1

SUCCESSFUL ITERATION SUN OF SQUARES= C.116280 05 NON-LINEAR P*RAUETFPS ~ 0.190840

~1

0.~8N06~

Q~

SIJCCESSFUL ITF~ATZOK SU~ OF SQUAkFS= ~.127410 ~ NON-L|N~AR PARAUFTFRS: 0.Z41570 01

0.1~23~0

01

3 Golub, G. H. and Pereyra,V. The differentiationof pseudo-inverses and non-linear least squares problems whose variables separate, SIAM J. Num. Anal. 1973, 10, 413 4 Marquardt,D. W. An algorithmfor least squares estimation of nonlinear parameters,J. Soc. Ind. Appl. Math. 1963, 11, 431 5 Osborne, M. R. Some special non-linear least squares problems, SIAM J. Num. Anal. 1975, 12, 571 6 Pereyra, V. Iterative methods for solving non-linear least squares problems, SIAM J. Num. Anal. 1967, 4, 27 7 Rao, R. C. and Mitra, S. R. Generalized Inverse of Matrices and its Applications, John Wiley, New York, 1971 8 Shanno, D. F. Parameterselection for modified Newton methods for function minimization, SIAM J. Num. Anal. 1970, 7, 366

Adv. Eng. Software, 1981, Vol. 3, No. 2 61