R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Solving Large-Scale Nonlinear System of Equations

dfsane

R Documentation

Solving Large-Scale Nonlinear System of Equations

Description

Derivative-Free Spectral Approach for solving nonlinear systems of equations

Usage

  dfsane(par, fn, method=2, control=list(),
         quiet=FALSE, alertConvergence=TRUE, ...)

Arguments

`fn`	a function that takes a real vector as argument and returns a real vector of same length (see details).
`par`	A real vector argument to `fn`, indicating the initial guess for the root of the nonlinear system.
`method`	An integer (1, 2, or 3) specifying which Barzilai-Borwein steplength to use. The default is 2. See Details.
`control`	A list of control parameters. See Details.
`quiet`	A logical variable (TRUE/FALSE). If `TRUE` warnings and some additional information printing are suppressed. Default is `quiet = FALSE` Note that the `control` variable `trace` and `quiet` affect different printing, so if `trace` is not set to `FALSE` there will be considerable printed output.
`alertConvergence`	A logical variable. With the default `TRUE` a warning is issued if convergence is not obtained. When set to `FALSE` the warning is suppressed.
`...`	Additional arguments passed to `fn`.

Details

The function dfsane is another algorithm for implementing non-monotone spectral residual method for finding a root of nonlinear systems, by working without gradient information. It stands for "derivative-free spectral approach for nonlinear equations". It differs from the function sane in that sane requires an approximation of a directional derivative at every iteration of the merit function F(x)^t F(x).

R adaptation, with significant modifications, by Ravi Varadhan, Johns Hopkins University (March 25, 2008), from the original FORTRAN code of La Cruz, Martinez, and Raydan (2006).

A major modification in our R adaptation of the original FORTRAN code is the availability of 3 different options for Barzilai-Borwein (BB) steplengths: method = 1 is the BB steplength used in LaCruz, Martinez and Raydan (2006); method = 2 is equivalent to the other steplength proposed in Barzilai and Borwein's (1988) original paper. Finally, method = 3, is a new steplength, which is equivalent to that first proposed in Varadhan and Roland (2008) for accelerating the EM algorithm. In fact, Varadhan and Roland (2008) considered 3 similar steplength schemes in their EM acceleration work. Here, we have chosen method = 2 as the "default" method, since it generally performe better than the other schemes in our numerical experiments.

Argument control is a list specifing any changes to default values of algorithm control parameters. Note that the names of these must be specified completely. Partial matching does not work.

M: A positive integer, typically between 5-20, that controls the monotonicity of the algorithm. M=1 would enforce strict monotonicity in the reduction of L2-norm of fn, whereas larger values allow for more non-monotonicity. Global convergence under non-monotonicity is ensured by enforcing the Grippo-Lampariello-Lucidi condition (Grippo et al. 1986) in a non-monotone line-search algorithm. Values of M between 5 to 20 are generally good, although some problems may require a much larger M. The default is M = 10.
maxit: The maximum number of iterations. The default is maxit = 1500.
tol: The absolute convergence tolerance on the residual L2-norm of fn. Convergence is declared when sqrt(sum(F(x)^2) / npar) < tol. Default is tol = 1.e-07.
trace: A logical variable (TRUE/FALSE). If TRUE, information on the progress of solving the system is produced. Default is trace = TRUE.
triter: An integer that controls the frequency of tracing when trace=TRUE. Default is triter=10, which means that the L2-norm of fn is printed at every 10-th iteration.
noimp: An integer. Algorithm is terminated when no progress has been made in reducing the merit function for noimp consecutive iterations. Default is noimp=100.
NM: A logical variable that dictates whether the Nelder-Mead algorithm in optim will be called upon to improve user-specified starting value. Default is NM=FALSE.
BFGS: A logical variable that dictates whether the low-memory L-BFGS-B algorithm in optim will be called after certain types of unsuccessful termination of dfsane. Default is BFGS=FALSE.

Value

A list with the following components:

`par`	The best set of parameters that solves the nonlinear system.
`residual`	L2-norm of the function at convergence, divided by `sqrt(npar)`, where "npar" is the number of parameters.
`fn.reduction`	Reduction in the L2-norm of the function from the initial L2-norm.
`feval`	Number of times `fn` was evaluated.
`iter`	Number of iterations taken by the algorithm.
`convergence`	An integer code indicating type of convergence. `0` indicates successful convergence, in which case the `resid` is smaller than `tol`. Error codes are `1` indicates that the iteration limit `maxit` has been reached. `2` is failure due to stagnation; `3` indicates error in function evaluation; `4` is failure due to exceeding 100 steplength reductions in line-search; and `5` indicates lack of improvement in objective function over `noimp` consecutive iterations.
`message`	A text message explaining which termination criterion was used.

References

J Barzilai, and JM Borwein (1988), Two-point step size gradient methods, IMA J Numerical Analysis, 8, 141-148.

L Grippo, F Lampariello, and S Lucidi (1986), A nonmonotone line search technique for Newton's method, SIAM J on Numerical Analysis, 23, 707-716.

W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual mathod without gradient information for solving large-scale nonlinear systems of equations, Mathematics of Computation, 75, 1429-1448.

R Varadhan and C Roland (2008), Simple and globally-convergent methods for accelerating the convergence of any EM algorithm, Scandinavian J Statistics.

R Varadhan and PD Gilbert (2009), BB: An R Package for Solving a Large System of Nonlinear Equations and for Optimizing a High-Dimensional Nonlinear Objective Function, J. Statistical Software, 32:4, http://www.jstatsoft.org/v32/i04/

Examples

  trigexp <- function(x) {
# Test function No. 12 in the Appendix of LaCruz and Raydan (2003)
    n <- length(x)
    F <- rep(NA, n)
    F[1] <- 3*x[1]^2 + 2*x[2] - 5 + sin(x[1] - x[2]) * sin(x[1] + x[2])
    tn1 <- 2:(n-1)
    F[tn1] <- -x[tn1-1] * exp(x[tn1-1] - x[tn1]) + x[tn1] * ( 4 + 3*x[tn1]^2) +
        2 * x[tn1 + 1] + sin(x[tn1] - x[tn1 + 1]) * sin(x[tn1] + x[tn1 + 1]) - 8 
    F[n] <- -x[n-1] * exp(x[n-1] - x[n]) + 4*x[n] - 3
    F
    }

    p0 <- rnorm(50)
    dfsane(par=p0, fn=trigexp)  # default is method=2
    dfsane(par=p0, fn=trigexp, method=1)
    dfsane(par=p0, fn=trigexp, method=3)
    dfsane(par=p0, fn=trigexp, control=list(triter=5, M=5))
######################################
 brent <- function(x) {
  n <- length(x)
  tnm1 <- 2:(n-1)
  F <- rep(NA, n)
  F[1] <- 3 * x[1] * (x[2] - 2*x[1]) + (x[2]^2)/4 
  F[tnm1] <-  3 * x[tnm1] * (x[tnm1+1] - 2 * x[tnm1] + x[tnm1-1]) + 
              ((x[tnm1+1] - x[tnm1-1])^2) / 4   
  F[n] <- 3 * x[n] * (20 - 2 * x[n] + x[n-1]) +  ((20 - x[n-1])^2) / 4
  F
  }
  
  p0 <- sort(runif(50, 0, 20))
  dfsane(par=p0, fn=brent, control=list(trace=FALSE))
  dfsane(par=p0, fn=brent, control=list(M=200, trace=FALSE))

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(BB)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/BB/dfsane.Rd_%03d_medium.png", width=480, height=480)
> ### Name: dfsane
> ### Title: Solving Large-Scale Nonlinear System of Equations
> ### Aliases: dfsane
> ### Keywords: multivariate
> 
> ### ** Examples
> 
>   trigexp <- function(x) {
+ # Test function No. 12 in the Appendix of LaCruz and Raydan (2003)
+     n <- length(x)
+     F <- rep(NA, n)
+     F[1] <- 3*x[1]^2 + 2*x[2] - 5 + sin(x[1] - x[2]) * sin(x[1] + x[2])
+     tn1 <- 2:(n-1)
+     F[tn1] <- -x[tn1-1] * exp(x[tn1-1] - x[tn1]) + x[tn1] * ( 4 + 3*x[tn1]^2) +
+         2 * x[tn1 + 1] + sin(x[tn1] - x[tn1 + 1]) * sin(x[tn1] + x[tn1 + 1]) - 8 
+     F[n] <- -x[n-1] * exp(x[n-1] - x[n]) + 4*x[n] - 3
+     F
+     }
> 
>     p0 <- rnorm(50)
>     dfsane(par=p0, fn=trigexp)  # default is method=2
Iteration:  0  ||F(x0)||:  27.0933 
iteration:  10  ||F(xn)|| =   0.08661234 
iteration:  20  ||F(xn)|| =   1.536949e-05 
$par
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 1 1 1 1 1 1 1 1 1 1 1 1

$residual
[1] 2.325577e-08

$fn.reduction
[1] 191.5786

$feval
[1] 28

$iter
[1] 27

$convergence
[1] 0

$message
[1] "Successful convergence"

>     dfsane(par=p0, fn=trigexp, method=1)
Iteration:  0  ||F(x0)||:  27.0933 
iteration:  10  ||F(xn)|| =   0.05559639 
iteration:  20  ||F(xn)|| =   2.043314e-05 
$par
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 1 1 1 1 1 1 1 1 1 1 1 1

$residual
[1] 6.70401e-08

$fn.reduction
[1] 191.5786

$feval
[1] 28

$iter
[1] 27

$convergence
[1] 0

$message
[1] "Successful convergence"

>     dfsane(par=p0, fn=trigexp, method=3)
Iteration:  0  ||F(x0)||:  27.0933 
iteration:  10  ||F(xn)|| =   0.07684265 
iteration:  20  ||F(xn)|| =   4.18081e-05 
$par
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 1 1 1 1 1 1 1 1 1 1 1 1

$residual
[1] 5.903011e-08

$fn.reduction
[1] 191.5786

$feval
[1] 28

$iter
[1] 27

$convergence
[1] 0

$message
[1] "Successful convergence"

>     dfsane(par=p0, fn=trigexp, control=list(triter=5, M=5))
Iteration:  0  ||F(x0)||:  27.0933 
iteration:  5  ||F(xn)|| =   12.76111 
iteration:  10  ||F(xn)|| =   0.08661234 
iteration:  15  ||F(xn)|| =   0.001736451 
iteration:  20  ||F(xn)|| =   1.536949e-05 
iteration:  25  ||F(xn)|| =   4.263627e-06 
$par
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 1 1 1 1 1 1 1 1 1 1 1 1

$residual
[1] 2.325577e-08

$fn.reduction
[1] 191.5786

$feval
[1] 28

$iter
[1] 27

$convergence
[1] 0

$message
[1] "Successful convergence"

> ######################################
>  brent <- function(x) {
+   n <- length(x)
+   tnm1 <- 2:(n-1)
+   F <- rep(NA, n)
+   F[1] <- 3 * x[1] * (x[2] - 2*x[1]) + (x[2]^2)/4 
+   F[tnm1] <-  3 * x[tnm1] * (x[tnm1+1] - 2 * x[tnm1] + x[tnm1-1]) + 
+               ((x[tnm1+1] - x[tnm1-1])^2) / 4   
+   F[n] <- 3 * x[n] * (20 - 2 * x[n] + x[n-1]) +  ((20 - x[n-1])^2) / 4
+   F
+   }
>   
>   p0 <- sort(runif(50, 0, 20))
>   dfsane(par=p0, fn=brent, control=list(trace=FALSE))
$par
 [1]  0.9697796  1.6932029  2.3260890  2.9062578  3.4501940  3.9669560
 [7]  4.4622034  4.9397769  5.4024368  5.8522534  6.2908316  6.7194479
[13]  7.1391395  7.5507642  7.9550420  8.3525849  8.7439189  9.1295003
[19]  9.5097286  9.8849557 10.2554941 10.6216229 10.9835925 11.3416289
[25] 11.6959369 12.0467028 12.3940970 12.7382754 13.0793816 13.4175481
[31] 13.7528971 14.0855424 14.4155895 14.7431369 15.0682763 15.3910938
[37] 15.7116700 16.0300806 16.3463968 16.6606857 16.9730107 17.2834316
[43] 17.5920052 17.8987850 18.2038221 18.5071648 18.8088591 19.1089490
[49] 19.4074762 19.7044805

$residual
[1] 8.16586e-08

$fn.reduction
[1] 134.9049

$feval
[1] 768

$iter
[1] 602

$convergence
[1] 0

$message
[1] "Successful convergence"

>   dfsane(par=p0, fn=brent, control=list(M=200, trace=FALSE))
$par
 [1]  0.9697796  1.6932029  2.3260889  2.9062576  3.4501939  3.9669558
 [7]  4.4622032  4.9397766  5.4024364  5.8522531  6.2908312  6.7194475
[13]  7.1391391  7.5507637  7.9550415  8.3525844  8.7439184  9.1294998
[19]  9.5097281  9.8849552 10.2554936 10.6216224 10.9835920 11.3416284
[25] 11.6959364 12.0467024 12.3940966 12.7382750 13.0793813 13.4175477
[31] 13.7528968 14.0855421 14.4155892 14.7431366 15.0682760 15.3910935
[37] 15.7116698 16.0300804 16.3463966 16.6606855 16.9730105 17.2834315
[43] 17.5920050 17.8987849 18.2038220 18.5071647 18.8088591 19.1089490
[49] 19.4074761 19.7044805

$residual
[1] 6.58967e-08

$fn.reduction
[1] 134.9049

$feval
[1] 635

$iter
[1] 631

$convergence
[1] 0

$message
[1] "Successful convergence"

> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>