nleqslv  R Documentation 
The function solves a system of nonlinear equations with either a Broyden or a full Newton method. It provides line search and trust region global strategies for difficult systems.
nleqslv(x, fn, jac=NULL, ...,
method = c("Broyden", "Newton"),
global = c("dbldog", "pwldog",
"cline", "qline", "gline", "hook", "none"),
xscalm = c("fixed","auto"),
jacobian=FALSE,
control = list()
)
x 
A numeric vector with an initial guess of the root of the function. 
fn 
A function of 
jac 
A function to return the Jacobian for the 
... 
Further arguments to be passed to 
method 
The method to use for finding a solution. See ‘Details’. 
global 
The global strategy to apply. See ‘Details’. 
xscalm 
The type of x scaling to use. See ‘Details’. 
jacobian 
A logical indicating if the estimated (approximate) jacobian in the solution should be returned. See ‘Details’. 
control 
A named list of control parameters. See ‘Details’. 
The algorithms implemented in nleqslv
are based on Dennis and Schnabel (1996).
Method Broyden
starts with a computed Jacobian of the function and
then updates this Jacobian after each successful iteration using the socalled Broyden update.
This method often shows super linear convergence towards a solution.
When nleqslv
determines that it cannot
continue with the current Broyden matrix it will compute a new Jacobian.
Method Newton
calculates a Jacobian of the function fn
at each iteration.
Close to a solution this method will usually show quadratic convergence.
Both methods apply a socalled (backtracking) global strategy to find a better (more acceptable) iterate. The function criterion used by the algorithm is half of the sum of squares of the function values and “acceptable” means sufficient decrease of the current function criterion value compared to that of the previous iteration. A comprehensive discussion of these issues can be found in Dennis and Schnabel (1996). Both methods apply an unpivoted QRdecomposition to the Jacobian as implemented in Lapack. The Broyden method applies a rank1 update to the Jacobian at the end of each iteration and is based on a simplified and modernized version of the algorithm described in Reichel and Gragg (1990).
When applying a full Newton or Broyden step does not yield a sufficiently smaller
function criterion value nleqslv
will attempt to decrease the steplength using
one of several socalled global strategies.
The global
argument indicates which global strategy to use or to use no global strategy
cline
a cubic line search
qline
a quadratic line search
gline
a geometric line search
dbldog
a trust region method using the double dogleg method as described in Dennis and Schnabel (1996)
pwldog
a trust region method using the Powell dogleg method as developed by Powell (1970).
hook
a trust region method described by Dennis and Schnabel (1996) as The locally constrained optimal (“hook”) step. It is equivalent to a LevenbergMarquardt algorithm as described in Moré (1978) and Nocedal and Wright (2006).
none
Only a pure local Newton or Broyden iteration is used. The maximum stepsize (see below) is taken into account. The default maximum number of iterations (see below) is set to 20.
The double dogleg method is the default global strategy employed by this package.
Which global strategy to use in a particular situation is a matter of trial and error. When one of the trust region methods fails, one of the line search strategies should be tried. Sometimes a trust region will work and sometimes a line search method; neither has a clear advantage but in many cases the double dogleg method works quite well.
When the function to be solved returns nonfinite function values for a parameter vector x
and the algorithm is not evaluating a numerical Jacobian, then any nonfinite
values will be replaced by a large number forcing the algorithm to backtrack,
i.e. decrease the line search factor or decrease the trust region radius.
The elements of vector x
may be scaled during the search for a zero of fn
.
The xscalm
argument provides two possibilities for scaling
fixed
the scaling factors are set to the values supplied in
the control
argument and remain unchanged during the iterations.
The scaling factor of any element of x
should be set
to the inverse of the typical value of that element of x
,
ensuring that all elements of x
are approximately equal in size.
auto
the scaling factors are calculated from the euclidean norms of the columns of the Jacobian matrix. When a new Jacobian is computed, the scaling values will be set to the euclidean norm of the corresponding column if that is larger than the current scaling value. Thus the scaling values will not decrease during the iteration. This is the method described in Moré (1978). Usually manual scaling is preferable.
When evaluating a numerical Jacobian, an error message will be issued on detecting nonfinite function values. An error message will also be issued when a user supplied jacobian contains nonfinite entries.
When the jacobian
argument is set to TRUE
the final Jacobian or Broyden matrix
will be returned in the return list.
The default value is FALSE
; i.e. to not return the final matrix.
There is no guarantee that the final Broyden matrix resembles the actual Jacobian.
The package can cope with a singular or illconditioned Jacobian if needed
by setting the allowSingular
component of the control
argument.
The method used is described in Dennis and Schnabel (1996);
it is equivalent to a LevenbergMarquardt type adjustment with a small damping factor.
There is no guarantee that this method will be successful.
Warning: nleqslv
may report spurious convergence in this case.
By default nleqslv
returns an error
if a Jacobian becomes singular or very illconditioned.
A Jacobian is considered to be very illconditioned when the estimated inverse
condition is less than or equal to a specified tolerance with a default value
equal to 10^{12}
; this can be changed and made smaller
with the cndtol
item of the control
argument.
There is no guarantee that any change will be effective.
The control
argument is a named list that can supply any of the
following components:
xtol
The relative steplength tolerance.
When the relative steplength of all scaled x values is smaller than this value
convergence is declared. The default value is 10^{8}
.
ftol
The function value tolerance.
Convergence is declared when the largest absolute function value is smaller than ftol
.
The default value is 10^{8}
.
btol
The backtracking tolerance.
When the relative steplength in a backtracking step to find an acceptable point is smaller
than the backtracking tolerance, the backtracking is terminated.
In the Broyden
method a new Jacobian will be calculated if the Jacobian is outdated.
The default value is 10^{3}
.
cndtol
The tolerance of the test for ill conditioning of
the Jacobian or Broyden approximation. If less than the machine precision it will
be silently set to the machine precision.
When the estimated inverse condition of the (approximated) Jacobian matrix is less than or equal to
the value of cndtol
the matrix is deemed to be illconditioned,
in which case an error will be reported if the allowSingular
component is set to FALSE
.
The default value is 10^{12}
.
sigma
Reduction factor for the geometric line search. The default value is 0.5.
scalex
a vector of scaling values for the parameters. The inverse of a scale value is an indication of the size of a parameter. The default value is 1.0 for all scale values.
maxit
The maximum number of major iterations. The default value is 150 if a global strategy has been specified. If no global strategy has been specified the default is 20.
trace
Nonnegative integer. A value of 1 will give a detailed report of the
progress of the iteration. For a description see Iteration report
.
chkjac
A logical value indicating whether to check a user supplied Jacobian, if
supplied. The default value is FALSE
. The first 10 errors are printed.
The code for this check is derived from the code in Bouaricha and Schnabel (1997).
delta
Initial (scaled) trust region radius.
A value of 1.0
or "cauchy"
is replaced by the length of the Cauchy step in the initial point.
A value of 2.0
or "newton"
is replaced by the length of the Newton step in the initial point.
Any numeric value less than or equal to 0 and not equal to 2.0
, will be replaced by 1.0
;
the algorithm will then start with the length of the Cauchy step in the initial point.
If it is numeric and positive it will be set to the smaller of the value supplied or the maximum stepsize.
If it is not numeric and not one of the permitted character strings then an error message will be issued.
The default is 2.0
.
stepmax
Maximum scaled stepsize.
If this is negative then the maximum stepsize is set to the largest positive representable number.
The default is 1.0
, so there is no default maximum stepsize.
dsub
Number of non zero subdiagonals of a banded Jacobian.
The default is to assume that the Jacobian is not banded.
Must be specified if dsuper
has been specified and must be larger than zero when dsuper
is zero.
dsuper
Number of non zero super diagonals of a banded Jacobian.
The default is to assume that the Jacobian is not banded.
Must be specified if dsub
has been specified and must be larger than zero when dsub
is zero.
allowSingular
A logical value indicating if a small correction
to the Jacobian when it is singular or too illconditioned is allowed.
If the correction is less than 100*.Machine$double.eps
the correction
cannot be applied and an unusable Jacobian will be reported.
The method used is similar to a LevenbergMarquardt correction and
is explained in Dennis and Schnabel (1996) on page 151.
It may be necessary to choose a higher value for cndtol
to enforce the correction.
The default value is FALSE
.
A list containing components
x 
final values for x 
fvec 
function values 
termcd 
termination code as integer. The values returned are

message 
a string describing the termination code 
scalex 
a vector containing the scaling factors, which will be the final values when automatic scaling was selected 
njcnt 
number of Jacobian evaluations 
nfcnt 
number of function evaluations, excluding those required for calculating a Jacobian and excluding the initial function evaluation (at iteration 0) 
iter 
number of outer iterations used by the algorithm. This excludes the initial iteration.
The number of backtracks can be calculated as the difference between the 
jac 
the final Jacobian or the Broyden approximation if 
You cannot use this function recursively.
Thus function fn
should not in its turn call nleqslv
.
Bouaricha, A. and Schnabel, R.B. (1997), Algorithm 768: TENSOLVE: A Software Package for Solving Systems of Nonlinear Equations and Nonlinear Leastsquares Problems Using Tensor Methods, Transactions on Mathematical Software, 23, 2, pp. 174–195.
Dennis, J.E. Jr and Schnabel, R.B. (1996), Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Siam.
Moré, J.J. (1978), The LevenbergMarquardt Algorithm, Implementation and Theory, In Numerical Analysis, G.A. Watson (Ed.), Lecture Notes in Mathematics 630, SpringerVerlag, pp. 105–116.
Golub, G.H and C.F. Van Loan (1996), Matrix Computations (3rd edition), The John Hopkins University Press.
Higham, N.J. (2002), Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, pp. 10–11.
Nocedal, J. and Wright, S.J. (2006), Numerical Optimization, Springer.
Powell, M.J.D. (1970), A hybrid method for nonlinear algebraic equations, In Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz (Ed.), Gordon & Breach.
Powell, M.J.D. (1970), A Fortran subroutine for solving systems nonlinear equations, In Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz (Ed.), Gordon & Breach.
Reichel, L. and W.B. Gragg (1990), Algorithm 686: FORTRAN subroutines for updating the QR decomposition, ACM Trans. Math. Softw., 16, 4, pp. 369–377.
If this function cannot solve the supplied function then it is a good idea to try the function testnslv in this package. For detecting multiple solutions see searchZeros.
# Dennis Schnabel example 6.5.1 page 149
dslnex < function(x) {
y < numeric(2)
y[1] < x[1]^2 + x[2]^2  2
y[2] < exp(x[1]1) + x[2]^3  2
y
}
jacdsln < function(x) {
n < length(x)
Df < matrix(numeric(n*n),n,n)
Df[1,1] < 2*x[1]
Df[1,2] < 2*x[2]
Df[2,1] < exp(x[1]1)
Df[2,2] < 3*x[2]^2
Df
}
BADjacdsln < function(x) {
n < length(x)
Df < matrix(numeric(n*n),n,n)
Df[1,1] < 4*x[1]
Df[1,2] < 2*x[2]
Df[2,1] < exp(x[1]1)
Df[2,2] < 5*x[2]^2
Df
}
xstart < c(2,0.5)
fstart < dslnex(xstart)
xstart
fstart
# a solution is c(1,1)
nleqslv(xstart, dslnex, control=list(btol=.01))
# Cauchy start
nleqslv(xstart, dslnex, control=list(trace=1,btol=.01,delta="cauchy"))
# Newton start
nleqslv(xstart, dslnex, control=list(trace=1,btol=.01,delta="newton"))
# final Broyden approximation of Jacobian (quite good)
z < nleqslv(xstart, dslnex, jacobian=TRUE,control=list(btol=.01))
z$x
z$jac
jacdsln(z$x)
# different initial start; not a very good final approximation
xstart < c(0.5,2)
z < nleqslv(xstart, dslnex, jacobian=TRUE,control=list(btol=.01))
z$x
z$jac
jacdsln(z$x)
## Not run:
# no global strategy but limit stepsize
# but look carefully: a different solution is found
nleqslv(xstart, dslnex, method="Newton", global="none", control=list(trace=1,stepmax=5))
# but if the stepsize is limited even more the c(1,1) solution is found
nleqslv(xstart, dslnex, method="Newton", global="none", control=list(trace=1,stepmax=2))
# Broyden also finds the c(1,1) solution when the stepsize is limited
nleqslv(xstart, dslnex, jacdsln, method="Broyden", global="none", control=list(trace=1,stepmax=2))
## End(Not run)
# example with a singular jacobian in the initial guess
f < function(x) {
y < numeric(3)
y[1] < x[1] + x[2]  x[1]*x[2]  2
y[2] < x[1] + x[3]  x[1]*x[3]  3
y[3] < x[2] + x[3]  4
return(y)
}
Jac < function(x) {
J < matrix(0,nrow=3,ncol=3)
J[,1] < c(1x[2],1x[3],0)
J[,2] < c(1x[1],0,1)
J[,3] < c(0,1x[1],1)
J
}
# exact solution
xsol < c(.5, 5/3 , 7/3)
xsol
xstart < c(1,2,3)
J < Jac(xstart)
J
rcond(J)
z < nleqslv(xstart,f,Jac, method="Newton",control=list(trace=1,allowSingular=TRUE))
all.equal(z$x,xsol)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.