Extended Newton-Raphson's Method

Bela Palancz

Department of Photogrammetry and Geoinformatics
Budapest University of Technology and Economy,H-1521,Hungary
e-mail:palancz@epito.bme.hu

Abstract

The standard Newton - Raphson' s method has been extended for under and overdetermined systems or  in more general, for situations when the Jacobian matrix is not regular, has deficient rank, i.e. it is singular or ill-conditioned, namely its condition number is high. Instead of using the inverse of the Jacobian, the pseudoinverse of the Jacobian computed by singular value decomposition is employed. The efficiency of the method is illustrated via many examples.

Keywords : Newton-Raphson method, over and underdetermined systems, singular Jacobian, Jacobian with deficient rank,  isolated multiple roots, low convergence, pseudoinverse, singular value decomposition.

Introduction

For solving  a set of nonlinear equations

f (x) = 0

where f:((ℜ )^n)->(ℜ^m) , the  Newton-Raphson procedure  is a well-known method.

If  n = m  and the Jacobi matrix  has a full rank everywhere, with other words the system of equations is regular, and if  in addition, the initial value of the iteration  is close enough to the solution, the method   ensures quadratic convergence.
One of these conditions fails, for example the system is over or under determined, or the the Jacobi matrix is singular, one can use extended Newton-Raphson method see. Haselgrove (1961), Ben-Israel (1966), Fletcher (1970), Ben-Israel and Greville (1974), Quoc-Nam Tran (1994) etc. 
In addition isolated multiple roots may cause low convergency. To restore quadratic convergence, deflation method can be used, see  Zhao (2007). 
Let us see some examples illustrating these problems.

Example 1

Let us consider  a  simple monomial system, Sommese and Wampler (2005)

(Subscript[f, 1])=(x^2);

(Subscript[f, 2])=x y;

(Subscript[f, 3])= (y^2);

This system is a "monomial ideal" and trivial for computer algebra.

Ideal=GroebnerBasis[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{x,y}]

That is why the solution is

Off[Solve::"svars"]

Map[Solve[#==0,{x,y}]&,Ideal]

We can see that the origin is an isolated singular root with multiplicity of  3. General polynomial solver using Jenkins-Traub algorithm can solve it

NSolve[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{x,y}]

However,  standard Newton - Raphson method fails,

FindRoot[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{{x,0.1},{y,0.1}}]

while the system is overdetermined.

Although we can transform the overdetermined system into a determined one in sense of least squares. The objective function is the sum of the square of the residium of the equations

obj=((Subscript[f, 1])^2)+((Subscript[f, 2])^2)+((Subscript[f, 3])^2)

Considering the necessary condition for the minimum

eqx=D[obj,x]

eqy=D[obj,y]

Then using the Newton - Raphson method we get

FindRoot[{eqx,eqy},{{x,0.1},{y,0.1}}]

Or we can use global minimization with genetic algorithm, too

NMinimize[obj, {x,y}]

The convergence is slow and the accuracy of the result is low because of  the multiple root,

{solxy,listxy}=Reap[FindRoot[{eqx,eqy},{{x,1},{y,1}},EvaluationMonitor:>Sow[{x,y}]]];

The solution is

solxy

The steps of the iteration are

listxy

the error of the solution

Norm[Last[Flatten[listxy,1]]]

The number of steps

Length[Flatten[listxy,1]]

P1=ListPlot[Flatten[listxy,1],PlotRange->All]

Now the multiplicity of the root (0, 0) has increased from 3 up to 9,

NSolve[{eqx,eqy},{x,y}]

Length[%]

Example 2

A usual test problem for parameter estimation procedures is the Bard problem, Bard (1974)

The measured data are

data={{0.14,1,15,1},
{0.18,2,14,2},
{0.22,3,13,3},
{0.25,4,12,4},
{0.29,5,11,5},
{0.32,6,10,6},
{0.35,7,9,7},
{0.39,8,8,8},
{0.37,9,7,7},
{0.5,10,6,6},
{0.73,11,5,5},
{0.96,12,4,4},
{1.34,13,3,3},
{2.1,14,2,2},
{4.39,15,1,1}};

The system of the nonlinear equations

eqs=Map[(p1+((#[[2]])/(p2 #[[3]]+p3 #[[4]]))-#[[1]])&,data]

Length[eqs]

We have 15 equations and 3 unknown parameters. The system is overdetermined. In that case general polynomial solver (NSolve[#]), Newton-Raphson method (FindRoot[#]) as well as elimination technique of computer algebra (GroebnerBasis[#]) fail.

Again, one can consider the minimum of the sum of square of  errors of the equations. The objective function

obj=Apply[Plus,Map[#^2&,eqs]]

Then one can find its minimum via global minimization

NMinimize[obj,{p1,p2,p3}]//Timing

Or transform the problem into a determined system, see above in Example 1. The necessary condition of the minimum

eq1=D[obj,p1]

eq2=D[obj,p2]

eq3=D[obj,p3]

Then we can use Newton method

FindRoot[{eq1,eq2,eq3},{{p1,1},{p2,1},{p3,1}}]//Timing

Although this can also fail, when the  starting values are far from the solution and the Jacobian is singular

FindRoot[{eq1,eq2,eq3},{{p1,-1},{p2,1},{p3,1}}]//Timing

General polynomial solver with Jenkins -Traub algorithm provides solution, but besides the correct solution in sense of least square, we have other real solutions, too, even more positive real solutions,

NSolve[{eq1,eq2,eq3},{p1,p2,p3}]//Timing

This is because the determined system may have more solutions than the original overdetermined system has.

Example 3

Now let us consider a determined system of polynomial equations, which is a somewhat modified version of the example of Ojika (1987),

f1[x_,y_]:=x^2+ y-3

f2[x_,y_]:=x+1/8 y^2-1

Show[{ContourPlot[{f1[x,y]==0,f2[x,y]==0},{x,-10,10},{y,-10,10},FrameLabel->{"x","y"}],Graphics[Text["(Subscript[f, 1])(x,y)= 0",{0,-9}]],Graphics[Text["(Subscript[f, 2])(x,y)=0",{-8,6.5}]]}]

First, let us  try to solve the problem using Newton Rapson method, starting with  (Subscript[x, 0])= - 1 and (Subscript[y, 0])=  -1,

FindRoot[{f1[x,y],f2[x,y]},{x,-1},{y,-1}]

This result is not correct

{f1[x, y], f2[x, y]} /.%

or with (Subscript[x, 0])=  1 and (Subscript[y, 0])=  -1,

FindRoot[{f1[x,y],f2[x,y]},{x,1},{y,-1}]

again, this is not a solution

{f1[x, y],f2[x, y]}/.%

These mean that none of them is a correct solution, Newton - Raphson method fails.

Example 4

Let us consider the following underdetermined system, see Quoc-Nam Tran (1998),

g1=(((x-u))^2)+(((y-v))^2)-1//Expand

g2=2 v (x-u)+3 (u^2)(y-v)//Expand

g3=(3 w (u^2)-1) (2 w v-1)//Expand

It goes without saying, that we have infinite roots. The general polynomial solver gives some solutions,

NSolve[{g1,g2,g3},{x,y,u,v,w}]

To make the solution unique, one may consider the solution with minimal norm. It means the problem can be transformed into a problem of  minimization with constrains

NMinimize[{(Sqrt[(x^2)+(y^2)+(u^2)+(v^2)+(w^2)]),g1==0,g2==0,g3==0},{x,y,u,v,w}]

In order to overcome these problems one can extend the Newton - Raphson method using pseudoinverse of the Jacobian matrix, which can be computed by  singular value decomposition.

Singular Value Decomposition

Every A matrix m × n, m ≥ n can be decomposed as

A = ((U Σ V)^T)

where (((.))^(T )) denotes the transposed matrix and  U is m × n matrix, V is n × n matrix  satifying

(U^T) U =(Subscript[(V^T) V =V (V^T) = I, n])

and Σ =<(Subscript[σ, 1]),...,(Subscript[σ, n])> a diagonal matrix.

These (Subscript[σ, i])'s,   (Subscript[σ, 1])≥(Subscript[σ, 2])≥,...,(Subscript[σ, n])≥0are the square root of the non negative eigenvalues  of  (A^T)A

and are called as the singular values of  matrix A.

As it is well - known from linear algebra, see i.e.Press et al (1995) singular value decomposition is a technique to compute pseudoinverse for singular or ill-conditioned  matrix of linear systems. In addition this method provides least square solution for overdetermined system and minimal norm solution in case of undetermined system.

Pseudoinverse

The  pseudoinverse of a matrix A  m × n is a matrix (SuperPlus[A])  n × m satisfying

A (SuperPlus[A]) A=A , (SuperPlus[A]) A (SuperPlus[A]) = (SuperPlus[A]),   (SuperStar[((SuperPlus[A]) A)]) =(SuperPlus[A]) A, (SuperStar[(A (SuperPlus[A]))])=A (SuperPlus[A])

where (SuperStar[(.)]) denotes the conjugate transpose of the matrix.

Always exists  a unique  (SuperPlus[A])  whic can be computed using SVD:

a) If m≥n and  A=U Σ (V^T)  then

(SuperPlus[A]) = V (Σ^(-1))(U^T)

where (Σ^(-1)) = <1/(Subscript[σ, 1]),...,1/(Subscript[σ, n])>

b) If m <  n then compute the (SuperPlus[((A^T))]), pseudoinverse of (A^T) and then

(SuperPlus[A]) = ((((SuperPlus[((A^T))])))^T)

Newton - Raphson Method with Pseudoinverse

The idea of using pseudoinverse in order to generalize of Newton method is not new but has been suggested in Haselgrove (1961), Ben-Israel (1966), Fletcher (1970), Ben-Israel and Greville (1974), Quoc-Nam Tran (1994) etc. 
It means that in the iteration formula, the pseudoinverse of the Jacobian matrix will be employed,

(Subscript[x, i+1])=(Subscript[x, i])-(SuperPlus[J])((Subscript[x, i])) f((Subscript[x, i]))

In Mathematica pseudoinverse can be computed in symbolic as well as in numeric form. For example considering Example 1, the Jacobi matrix is

jac=Outer[D,{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{x,y}];MatrixForm[jac]

its pseudoinverse

pinvjac=Simplify[PseudoInverse[J],Element[{x,y},Reals]];MatrixForm[pinvjac]

and at the point  (x, y) = (1., 1.)

jacN=jac/.{x->1.,y->1.};MatrixForm[jacN]

pinvjacN=PseudoInverse[jacN];MatrixForm[pinvjacN]

or

(pinvjac/.{x->1.,y->1.})//MatrixForm

The new values  ((Subscript[x, i+1]),(Subscript[y, i+1])) in the next iteration step are

{1,1}-pinvjacN.({(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])}/.{x->1.,y->1.})

Implementation in Mathematica

The implementation is a modified adoptation of  the function written by Ruskeepää (2004).

NewtonExtended[f_List,x_List,x0_List,eps_:10^-12,n_:100]:=With[{jac=N[Outer[D,f,x]]},FixedPointList[(#+PseudoInverse[jac/.Thread[x->N[#]]].(-f/.Thread[x->N[#]]))&,N[x0],n,SameTest->(Sqrt[Abs[(#1-#2).(#1-#2)]]<eps&)]]

where

f-    the list of the equations of the system, {(Subscript[f, 1]),...,(Subscript[f, m])}

x-     the variables of the system, x ={(Subscript[x, 1]),...(Subscript[x, n])}

x0-   the numerical start value of the iteration,  x0 ={(Subscript[x0, 1]),...(Subscript[x0, n])}

eps- the limit of the error , employing Frobenius norm, the default value is 10^-12

n -    the maximum number of the iterations, the default value is 100

In the following section we shall solve the Examples 1 - 4. with this method.

Solutions of the Examples with Extended Newton's Method

Example 1

sol=NewtonExtended[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{x,y},{1,1},5 10^-8]

The solution is

Last[sol]

Its error

Norm[%]

The number of steps of the iteration

Length[sol]

Show[{P1,ListPlot[sol,PlotStyle->{PointSize[0.015],RGBColor[1,0,0]},PlotRange->All]}]

In case of the same magnitude of error, this method requires less steps, which means it has faster convergence than the standard Newton-Raphson method.

Example 2

Now, we can solve directly the original overdetermined system

eqs

sol=NewtonExtended[eqs,{p1,p2,p3},{-1,1,1}]//Timing

Again, Extended Newton - Raphson Method is robust, Jacobian singularity is not a problem as well as considerable faster than NSolve(#).

Example 3

sol=NewtonExtended[{f1[x,y],f2[x,y]},{x,y},{-1,-1}]

sol=NewtonExtended[{f1[x,y],f2[x,y]},{x,y},{1,-1}]//Last

This is really a solution

{f1[x,y],f2[x,y]}/.{x->sol[[1]],y->sol[[2]]}

The method is very robust, its convergence does not depend on the start values,

t1=NewtonExtended[{f1[x,y],f2[x,y]},{x,y},{-1,-1}];

Last[t1]

t2=NewtonExtended[{f1[x,y],f2[x,y]},{x,y},{1,-1}];

Last[t2]

t3=NewtonExtended[{f1[x,y],f2[x,y]},{x,y},{-1,1}];

Last[t3]

GraphicsArray[{ListPlot[Map[#[[1]]&,t1],Joined->True,PlotRange->All,Mesh->All],ListPlot[Map[#[[2]]&,t1],Joined->True,PlotRange->All,Mesh->All]}]

GraphicsArray[{ListPlot[Map[#[[1]]&,t2],Joined->True,PlotRange->All,Mesh->All],ListPlot[Map[#[[2]]&,t2],Joined->True,PlotRange->All,Mesh->All]}]

GraphicsArray[{ListPlot[Map[#[[1]]&,t3],Joined->True,PlotRange->All,Mesh->All],ListPlot[Map[#[[2]]&,t3],Joined->True,PlotRange->All,Mesh->All]}]

Example 4

Let us solve the underdetermined system with different random intial values

ig=Table[Random[Integer,{-1,1}],{5}]

sol=NewtonExtended[{g1,g2,g3},{x,y,u,v,w},ig]//Last

This really a solution

{g1,g2,g3}/.{x->sol[[1]],y->sol[[2]],u->sol[[3]],v->sol[[4]],w->sol[[5]]}

The Froebenius norm of this solution

Norm[sol]

which is quite close the the minimal norm. Let us generate an other random initial set of values

ig=Table[Random[Integer,{-1,1}],{5}]

sol=NewtonExtended[{g1,g2,g3},{x,y,u,v,w},ig]//Last

{g1,g2,g3}/.{x->sol[[1]],y->sol[[2]],u->sol[[3]],v->sol[[4]],w->sol[[5]]}

Norm[sol]

Conclusions

Extended Newton - Raphson method using pseudoinverse of the Jacobian can overcome problems caused by deficient rank in case of  over and underdetermined systems, singular or ill-conditioned Jacobian as well as can restore quadratic  convergence at multiple roots. The implemented method in Mathematica is fast and robust concerning the starting values of the iteration procedure.

References

Bard Y. (1974) Nonlinear Parameter Estimation, Academic Press,New York.

Ben-Israel, A  (1966)  A Newton-Raphson Method for the solution of systems of equations. J.Math.Anal.Applic. 15. pp.243-252

Ben Israel, A. and Greville T.N.E (1974) Generalized Inverse Matrices, Wiley-Interscience, New York.

Fletcher, R.(1970) Generalized inverses for nonlinear eqautions and optimization. In Rabinowitz, P. ed. Numerical ethods for Nonlinear Algebraic Equations, pp.75-85, Gordon and Breach, London

Haselgrove, C.B. (1961)  The solution of nonlinear equations and of diferntial equations with two-point boundary conditions. Computing J. 4., pp.255-259

Ojika, T.  (1987) Modified deflation algorithm for the solution of singular problems. I. A system of nonlinear algebraic equations. J.Math.Anal.Appl. 123., pp. 199-221

Press W.H, Teukolsky S.A., Vetterling W.T. and Flannery B.P. (1995) Numerical Recipes in C, Second Edition, Cambridge, University Press

Ruskeepää H. (2004). Mathematical Navigator, Second Edition, Elsevier Academic Press, Amsterdam


Quoc-Nam Tran (1994)  Extended Newton's method  for  finding the roots of an arbitrary system of nonlinear equations. In Hamza, M.H. ed. Proc. 12th IASTED Int. Conf.on Applied Informatics, IASTED, Anaheim, CA. 


Quoc-Nam Tran (1998), A Symbolic-Numerical Method for Finding a Real Solution of an Arbitrary System of Nonlinear Algebraic Equations, J.Symbolic Computation 26. pp. 739-760


Sommese, A.J. and Wampler, C.W. (2005) The Numerical Solution of Systems of Polynomials arising in Engineering and Sciences, Word Science Press.


Zhao A. (2007) Newton's Method with Deflation for Isolated Singularities of Polynomial Systems, Ph.D. Thesis. University of Illinois at Chicago, USA