Iterative methods for solving nonlinear least squares problem

Pereyra, V.

English

NASA Technical Reports Server

MATHEMATICS RESEARCH CENTER, UNITED STATES ARMY 
THE UNIVERSITY OF WISCONSIN . 
Contract No. : DA-11-022-ORD-2059 
ITERATIVE METHODS FOR SOLVING 
NONLINEAR LEAST SQUARES PROBLEMS 
Victor Pereyra 
MRC Technical  Summary Report #622  
February 1966 
Madison, Wisconsin 
https://ntrs.nasa.gov/search.jsp?R=19670015871 2020-03-16T18:21:50+00:00Z
ABS TRAC T 
In many technical applications it is desired to f i t  a nonlinear model t o  
a s e t  of observations. 
t o  determine a bes t  s e t  of parameters in the l ea s t  squares  sense .  
Several iterative techniques have been devised in  order 
In t h i s  paper we  d i scuss  conditions for convergence, and give error 
es t imates  for a c l a s s  of methods, which includes as particular c a s e s  some 
well  known techniques. 
modified Newton's i terations for a suitable functional equation, and then a 
general  theorem, f i rs t  indicated by Bartle, is proved and applied to this  
particular case.  The hypotheses are s e t  in such a way that their checking by 
an  automatic computer i s  made possible. 
I t  is shown t h a t  those methods can  be considered as  
Some numerical examples are given. The main aim is to show that  t he  
automatic error estimation procedure works, rather than attempting t o  optimize 
the computational scheme.  
ITERATIVE METHODS FOR SOLVING NONLINEAR LEAST SQUARES PROBLEMS 
Victor Pereyra 
1. Introduction. 
In many technical applications it is desired t o  f i t  a nonlinear model t o  a 
s e t  of observations. When the bes t  fit i s  sought in the l e a s t  squares  s e n s e  the 
problem can  be s ta ted  a s  follows: 
Given the nonlinear transformation F( x) = y between the  f in i te  dimensional . ~ - -  -- 
n 
Euclidean spaces  E 
b E E m ,  find a vector x'" E E 
and Em ( n  m )  , and the vector of "observations" -
-" n which minimizes the  L -norm of F( xJ - b. -2 - _ -  
We consider in  t h i s  paper a general class of i terative methods for finding 
2 
stat ionary points of ( 1  F( x) -bI[ 
of the  form: 
If we call f(  x) = F( x) - - b ,  these methods are  of - -  - -  - 2 '  - -
( v  = O , l ,  . . . ) ,  
-1 T 
X = x  - [ T  1 fpV) f (x)  Y t 1  -7/ Y 
n 
where t h e  T a re  linear, nonsingular transformations of E in  i tself ,  and 
[2v] 
-V 
T is the t ranspose of the Jacobian matrix of the transformation - f . In 
components 
By choosing T appropriately we can  obtain some well  known methods 
-V 
used  i n  the solution of nonlinear l ea s t  squares  problems. For instance,  i f  
Sponsored in part by the Mathematics Research Center, United S ta tes  Army under 
Contract  No. : DA-11-022-ORD-2059, and in  part by NASA Research Grant 
NGR-50-0 22-0 28. 
then the resulting i teration corresponds t o  the Gaus s-Newton T T = f g ( x  ) L'((x,)  
-I, - - v  
method ( see Moore and Zeigler [ 71 or Hartley [ 41 and references there in) .  
T 
If -v T = - f g (  -0 x ) - f 9 (  -0 x ) for all v ,  then we could call th i s  the simplified 
Gauss-Newton method. 
Jakovlev [ 51, which replace f 9 (  x ) f 9 (  x ) by approximate expressions could 
In general, schemes l ike  those  of Powell [ 91, and 
T 
-v - -I, - 
be considered into the  class of methods described by ( 1.1).  
Sufficient conditions for the convergence of the Gaus s-Newton method 
were obtained i n  Zadunaisky and Pereyra [ 111 by applying a standard fixed point 
t h e  ore m, 
2 
2 Observing that  to find the s ta t ionary points of 11 f (  x)I/ is equivalent t o  
finding the zeros  of i t s  gradient, t h e  problem is reduced to the  solut ion of the 
n X n system of nonlinear equations 
T 
= f ' (  x) l f (  x) = - 0 . - -  - -  
If we  now regard the i teration (1. 1) as a method for solving ( 1. 2) then a 
general  theorem by Bartle [ 11 c a n  be applied in  order t o  obtain suff ic ient  conditions 
for convergence, and error bounds for the s u c c e s s i v e  approximations. We include 
a proof of th i s  theorem s ince  in  Bartless paper it is only vaguely indicated. 
Considering the solution of overdetermined systems of nonlinear equations 
n m  
of t h e  form - -  f (  x) - 7  0 - f : E -+ E 
gence of generalized versions of Newton's method and the  simplified ( or modified) 
, Ben-Israel [ 21, [ 31 h a s  s tud ied  the  conver- 
Newton method. The general izat ions c o n s i s t  in  replacing the inverses  of the 
Jacobian matrix of - -  f( x) , which appear in the standard nonsingular n X  n case, 
by their  pseudoinverses. This in  particular permits considerations of the  case 
in  which rank f * (  x) < n . - -  
- L- #622 
In this  paper we will  only be concerned with t h e  c a s e  of full rank ( = n )  
Since i n  tha t  case 
( 1 . 3 )  [ f s (  x) T f '(  x) ]- ' f ' (  x) pseudoinverse of f ' (  x) E [ f(5) ] t , - -  - - - - -  - -  
we see tha t  the Gauss-Newton method coincides  with Ben-Israels s generalized 
Newton-Raphson's method. 
The conditions obtained from Bartle's theorem for the general  i teration (1.1) 
are spec ia l ized  t o  t h e  Gauss-Newton method 
(1. 3) for the pseudoinverse of f ' (  x) a l lows 
cedure which can  be implemented on a digita 
- -  
in Theorem 4. 1, and the representation 
u s  t o  give an error estimation pro- 
computer. Ben-Israel's theorem and 
our resu l t  do  not seem t o  be comparable. On one s ide  we require conditions on 
the second derivatives of t h e  function f (  x) 
we require only the boundedness of 11 [ - fs( xo) fL(x,)J I( while he  uses the 
condition 11 [ - -  f '  ( x) ] 
while he  dDes not. On the other hand - -
T 1 
- [_f' (z)] '11 - < N(I - x-y(I which is Dbviously more difficult 
t o  verify. 
An implementation of the error estimation procedure is briefly explained 
in  Section 5. Finally we present in  Section 6 two numerical examples. 
#622 -3  - 
2. Approximate Newton type iteration. 
A s  we  said in the Introduction, the problem of finding the stationary values  
of IIf(5) - b(/Z is equivalent to that  of finding the solutions of the n X  n system 2 
of nonlinear equations: 
T q(2) = f ' ( x )  f ( x )  = o  -- -  - -  
T 
where f (  x) = - -  F( x) - - b ,  and - -  f ' (  x) is i t s  Jacobian matrix. If w e  put - f ' (  - -  x) f ' ( - x)' - -  
N( - -  x) and assume tha t  - -  f (  x) is twice Frgchet differentiable in a certain region 
S2 c En , then the Fre'chet derivative ( Jacobian) of - -  q( x) can  be written as 
This is e a s i l y  obtained by applying the operational calculus  with Fre'chet deriva- 
t ives  ( cf. Dieudonnd [ 121). 
I In the benefit of those readers not familiar with this  calculus  we wil l  obtain 
th i s  formula by operating on the components of the vector functions involved. We 
I will  u s e  tensor  notation with the summation convention. On the f i rs t  place 
and 
Thus, going back to  the matrix notation w e  obtain ( 2. 2 ) .  
The standard Newton-Kantorovich method for solving ( 2. 1) is: 
( 2 . 3 )  
-1- #622 
and the approximate Newton iteration ( s e e  Bartle [ 11) i s  obtained i f  - -  c p ’ (  x v )  is 
replaced by linear, nonsingular operators T which a re  c lose  to cp’( x ) i n  
some sense .  
-V - -0 
But this is essent ia l ly  what is s ta ted  in  equation (1.1) .  Thus we 
can apply Bartle’s theorem t o  that iteration obtaining sufficient conditions for its 
convergence. We can write those conditions in  the following form: 
n 
Theorem 2.1. &t T be a sequence of linear nonsingular - _ _  operators from E 
into itself, such that, for x E En and p > 0 the sphere S( xo ~ p)  c R ,  
-V 
- -0 ~- - 
-k< 
( e )  k = A ( ~ + E )  < 1 ,  r -  l - k - p  
Under these  conditions t h e  iteration (1.1) is wel l  defined and converges 
_. -. 
.I. ‘6%J- 
t o  a solution - x”. of g( 5) = - 0 .  
solution contained in this  sphere. 
1lC-x 11 < kvr  . 
Proof: First of a l l  from ( 1. 1 ) ,  ( a)  ,and ( d) i t  follows tha t  
Furthermore, 115 - xo ( 1  - < r and is the only 
The  rapidity of the convergence is given by 
.(r 
-v - 
Furthermore, by ( 1. 1) , ( b) , and  ( c) 
#622 - 5-  
- 6- 
With formulas ( 2. 4) and ( 2. 5) we  have s tar ted an  induction argument. 
Assume that, for u = 1 , .  . . , n 
( i u )  I I X  - x II 5 P , 
-1, -0 
( ii ) II x - x I1 < l/I'p(xv-l) /I 
(i i i  
U -v -u-1 - 
ll(p(x > l l  < ( P I - E ) I I x  -u  - x  -u-1 II . U -u  - 
From ( 1.1) and ( a )  we obtain 
II x - x \I (&+ -n ) I1 , -n t l  -n ( ii ) n t l  
and by ( iii ) , ( ii ) ,  ( iiiv-l), . . . 
U V 
Furthermore 
which is ( i  ) . n t l  
Having proved th is  w e  c a n  u s e  (c) on 
#622 
t h u s  obtaining (iii  ) which completes the induction argument. n t l  
With this  we can es tab l i sh  that the sequence {x 1 generated by (1. 1) 
-1, 
is a Cauchy sequence and simultaneously we can obtain the error estimation. 
In fact, for any p > 0 
.b 
Consequently, there ex i s t s  x“: E S( x , r )  such that x --f x *,. , and -0 -V - 
fur  therm ore 
can  write 
.L ‘8 .  J- 4, -,,. , 1  
which is impossible unless  x = x . - - 
# 6 2 2  -7-  
3 .  Another s e t  of sufficient conditions. 
Sometimes it is computationally advantageous t o  replace hypothesis ( c) 
in  Theorem 2 . 1  by other conditions which, though giving a s l ight ly  less general 
result ,  a re  more e a s y  to handle. 
Lemma 3.1. Suppose that y ,  K > 0 ex i s t  s u c h  that 
2 
( c2) I l g g l l  ( K 2 ,  zc S ( X 0 ,  P )  , 
are  sat isf ied,  then ( c )  follows with p = K p t 
Proof: We will f i rs t  calculate  a bound for !I - 'p' 
( 2. 2)  and the hypotheses 
2 
If we call  K p -t y = p then by Bartless L e m m a  1 : 2 
ll&) -&I - ' p s ( x g ) ( l ( - - y )  I 1  ( P l l l L - y I l  
which is ( c ) ,  
-8- #622 
4. The Gauss-Newton method. 
By taking T = N( x ) i n  ( 1.1) we obtain the wel l  known Gauss-Newton -v --v 
method fo r  solving nonlinear l ea s t  squares problems. We will  express  now the 
conditions ( a )  and ( b )  of Theorem 2.1  in terms of some of the quant i t ies  used in  
Section 3 .  If ( cl) and ( c2)  hold then 
implies condition ( b) of Theorem 2. 1 with E = p . This  follows e a s i l y  from 
We want to  show now that the ex is tence  and uniform boundedness of 
-1 E() 
-1 N_(x,) . 
(ZE g(xo ,  p))  is a consequence of the existence and boundedness of 
To do s o  we wil l  s ta te  without proof a wel l  known resul t  of the theory 
of matrices: 
Lemma 4.1. Let B and C be n X n  matrices. Assume that  - - - - -  
( i) E is nonsingular and I \  - B-lIl  - < CY ; 
( iii) CY 6 < 1 , 
then - -  C is nonsingular and 
From th i s  we can  prove the following: 
#6 22 -9 - 
ex i s t s  and I( N( - -  x)-l(( < 1. 
Proof: Take in  Lemma 4. 1 
1 
2 B = N ( x ) ,  C = N ( x ) ,  a = - A  and 6 = K  2 p .  - - -  - -0 - 
1 1 
2 2  4 Since a6 =-K p we have by ( 3 .  1) that a 6  < - and thus N( - x )  is nonsingular. 
Furthermore 
-10- 
Collecting these  results together we can  s t a t e  the following theorem. 
Theorem 4. 3 .  With the same notation as above and for x E R ,  l e t  u s  assume that 
-0 
and that  N ( x  ) is nonsingular. .- Define - - - 0  -1- 
A E 2 \I N(xO)- ' / l  - -  , 
and assume tha t  k = 2 A (  K p t y) < 1 ,  Define -____ 2 -- 
(4.4) 
( 4 . 5  
( 4 . 6 )  
( 4 .  7 )  
Assume further that r < p ,  then the sequence (x 1 defined by -v ------ -
(4.  8 )  
-1 
- v t l  X = x  -v -N_((xv) p(Ijv) 
converges to  the unique solution x" of - - -  cp( x) = - 0 in the sphere S( xo, r )  . 
Moreover, the rate of convergence is estimated by _--~ ____-- 
1 1 ~ ' ' -  x 11 < kvr  . ( 4 . 9 )  -v - - 
#622 
__ 
I 
5. Automatic error estimation. 
A Fortran 63 program has  been written for the CDC 1604 computer a t  the 
University of Wisconsin Computing Center, which implements the Gauss-Newton 
method with the error estimation procedure given in Theorem 4. 3. 
The hypotheses a re  checked automatically a s  the i teration proceeds, and 
a s  soon a s  they are sat isf ied it can  be ensured that  the process  converges and 
the bound ( 4 .  9 )  used to es t imate  the norm of the error. 
A s  s een  in  s t ep  V of the procedure described below: a r e l axa t im  technique 
is used  in  order to  prevent divergence in  the earlier s tages  of the iteration 
( cf. Hartley [ 41 ) . 
An interesting feature i s  
Reiter [ l o ] )  to  compute y and 
a s  a n  n-dimensional hypercube 
the use of interval arithmetic ( cf. Moore [ 61 ~ 
K 2  in Theorem 4. 3. Given the region L? 
C ( 5 ,  p)  with edge 2p and center  - x ~ we compute 
N 7- 
-I N 
( 4. 3 )  and ( 4. 4) in interval arithmetic with the argument x = ( x 
and from there we can  obtain y and K immediately. 
- p ~ x i  -t p )  i - 
2 
A further sophistication would be t o  use a pragram to generate the code 
for the partial derivatives which a re  needed in the discussion.  
this  s ince  our t e s t  c a s e s  were very simple, but such a program is avai lable  
( Reiter [ 101 ) and it has  been successfu l ly  applied to  the solution of some 
complicated s y s  tems of nonlinear equations by Newton's method. 
We have not done 
We will  briefly describe n3w the computational scheme.  We u s e  in this  
description informal Algol. 
Naur [ 81 . 
For details  on the Algorithmic Language Algol, cf. 
The notation is t h e  same as  in Section 4. The narm used  is the Id@ norm 
#622 -11- 
2 
except  for the residual  11 fl12 where, of course,  w e  u s e  the L norm. 2 
Let x and p be given. -0 
begin procedure error; 
i f  N( x ) is singular then go to  error 1 ;  ----0 - 
- -1 T 
11: A X  :=N_(so) T ( x  ) f ( x  ) ;  p : = p  ; -0 -0 --0 
IV: i f  2y X < 1 then go to  Rest of error procedure -- e l s e  --- -0 0  - 
V: comment this is the relaxation routine; 
for i I = 0 s t ep  1 until  p d o  begin - -- 
comment p is a given integer; 
i -i x = x  - 2  A x  ; -1 -0 -0  
i x : = x  ; g o  to  I end e n d ;  -0 -1 -- --
Rest of error procedure: K 2  a -  max ( ( I  N_Vx)ll) :
X € S ( X  , p )  -0 - 
k : = 2 X ( K 2 p S y  0 ) :  
-12- f i622  
VI : i f  - 
i f  - 
II 'p( so) I1 
k < 1 then begin R O -  - ; 
0 .- 1-  k --
R > p then go t o  V ;  --- 0 
comment E is the desired accuracy; 
for i : = O  s t e p  1 until q do - --
x :=x -N_(xi) -1fs(xi) T r ( x . )  end - i t1  -i -1 -
go to v ;  -- 
error 1 :  e n d .  
#622 
else 
-13 - 
6. Test  c a s e s .  - 
The following t e s t  c a s e s  have been run on t h e  CDC 1604: 
x t  
1) Fi( xl, x2, x3) = g( t i ,  2) = x2 e l i + X  3 ’  
1 0 < m < 2 0 ,  I t . \  5 1 0 ,  yi taken from tab les  with different accurac ies  
1 - -  
( f o r  given - x::’ ) . 
t 
i 
Sample results are  given in Table 1 for y = 0. 1 e - 5 truncated a t  the i 
f i f th  significant figure; m = 10, - 2  < t .  < 2. 5 .  
1 -  - 
The conditions were usually fulfilled when the i terates  were quite c lose  
to the exac t  solution. 
g ( t i , x  1’ x 2’ x )  3 = x  2 s i n ( x t )  l i  + x  3 ’  
Taking y = s in  t t o  3D for t .  = ( 0 .  105, 0. 25, 0 . 4 ,  0. 55, 0. 7,  0 . 9 ,  
We have obtained the resu l t s  
i i 1 
1. 1, 1. 25, 1, 35, 1. 45, 1. 55, 1. 57, 1. 6 ) ,  m = 13. 
shown in Table 2. 
In t h e  second example we see that  before the 3rd i teration the condition 
k < 1 is not fulfilled. Then 2Xy becomes l e s s  than one and simultaneously 
k < 1 is a l s o  fulfilled. Thus r can  be calculated and i t  comes t o  be < p .  The 
error estimations a re  then read in the corresponding column. In this  example we 
have, instead of using the second part of VI, chosen to  proceed as i f  every 
i teration were the first ,  s ince  in  this way the error estimation is much better 
than the one obtained directly from the  theorem. 
have two error estimations: 7. 9 X 10 
as  x and 3. 5 X  10 i f  we recompute everything anew. For the  5th i terate 
Thus, in the 4th i teration we 
-5 
obtained by considering the third i terate  
-10 
-0 ’ 
-14- #622 
we give in  parentheses the error estimation obtained i f  the 4th i terate  is taken 
as  x . 
-0 
The disparity showed by these  estimations stem from the fac t  that in this  
2 i tself  is quite s m a l l .  2 11 f ( x") ( 1  c a s e  the convergence is fas te r  than l inear s ince 
This does not have to  be the  c a s e  in more real  problems in which the model cannot 
be expected to reproduce the observations very accurately.  
- -
The conclusion we can  draw from these  and other experiments we have 
carried out is that the most difficult part i n  t h i s  error procedure i s  the choice of 
a n  appropriate p . 
Also i t  is clear  that in i t s  present form the conditions a re  fulfilled only 
when we are  quite c lose  to  the exac t  solution. 
t he  definition of r that  this  quantity can be obtained with very l i t t le  effort i f  k 
is ignored and, from a practical  point of view, small values  of r could be enough 
assurance  of convergence, and they could be used as error estimators without any 
further check. 
s t ep  I of t he  error procedure. 
On the other hand, we see from 
In this c a s e  the only extra computation would be that of t h e  
Acknowledgement. 
comments and suggestions.  
The author is grateful to  the referees for their  most valuable 
#622 -15- 
I Q  
N 
X 
X" 
I 
, 
i 
(L, 
iJ .ri
. 
* L n m  
co 
0 
* 
Ln 
cr) 
co 
0 
N 
Ln 
In 
m 
0 
H 
H 
N 
M co 
N) r- 
M 
In 
0 
rs, 
0 
3 
" 
x 
N 
Q 
e- 
X 
N 
X 
X" 
I_ 
i 
2 
-i 
I 
cc 
O N  
m N  
m o  
0 0  
0 0  
0 0 0 0  
0 0 0 0  
0 0 0 0  
Ln 
Q D  
+ \ D  
d o  
0 0  
__ 
H M m w m  
0 0 0 0 0  
0 0 0 0 0  
" 0 0 0 0 0  
0 0 0 0 0 0  
I I I I I 
r\ 
0 
4 
Ln 
M 
0 
0 
0 
3 
REFERENCES 
Bartle R. G. ,  "Newton's method in  Banach spaces" ,  Proc. A. M. S. - 6 ( 1955) ~ 
pp. 827-831. 
Ben-Israel A. ,  "A modified Newton-Raphson me thod  for the solution of 
equations", I s rae l  J. Math. - 3 ( 19 65) ,  pp. 94-9 8. 
"On Newton's method for the solution of systems of equations", 
to appear i n  J. Math. Anal. Appl. 
Hartley H. O. ,  "The modified Gauss-Newton method for the fi t t ing of 
nonlinear regression functions by l e a s t  squares'!, Technometrics 3 - 
( 1961), pp. 269-280. 
Jakovlev M. N. ,  "On the solution of nonlinear equations by i terations",  
Dokl. Akad. Nauk SSSR 156 ( 19 64) ,  pp. 522-524. -
Moore R. E . ,  "The automatic ana lys i s  and  control of error in digi ta l  
computation based on the u s e  of interval  numbers", in  Error in  Digital 
Computation, Vol. I ( ed. by L. B. Rall) ( 19 6 5 ) ,  Wiley, New York. 
Moore R. H. and Zeigler R. K . ,  "The solution of the general least squares  
problem with spec ia l  reference to high speed  computers", Los Alamos 
Sc. Lab. LL-2367 ( 19 59) ,  3 8 pp. 
Naur P. ,  "Revised report on the algorithmic language Algol 60", Comm. ACM 
- 6 ( 19 6 3 ) ,  pp. 1-17. 
Powell M. J. D., ' I  A method for minimizing a sum of squares of nonlinear 
functions without calculating derivatives", Comp. J. - 7 ( 19 65) ,  
pp. 303-307. 
Reiter A . ,  "Compiler of differentiable expressions ( CODEX) ' I ,  MRC 
Computer Program #l; I'Interval arithmetic package ( INTERVAL) 1 1 ,  
MRC Computer Program #2. 
( 1965) . 
University of Wisconsin,  Madison 
Zaduanisky P. and  Pereyra V.,  "On the  convergence and precision of a 
process  of successive differential  corrections", to appear in  Proc. 
IFIPS 65. 
Dieudonn6 J. , "Foundations of Modern Analysis", Academic Press, 
New York ( 1960). 
#622 -17- 


Iterative methods for solving nonlinear least squares problems

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server