R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Solving Large-Scale Nonlinear System of Equations

sane	R Documentation

Solving Large-Scale Nonlinear System of Equations

Description

Non-Monotone spectral approach for Solving Large-Scale Nonlinear Systems of Equations

Usage

  sane(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 sane implements a non-monotone spectral residual method for finding a root of nonlinear systems. It stands for "spectral approach for nonlinear equations". It differs from the function dfsane in that it 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 and Raydan (2003).

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 and Raydan (2003); 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 equivalent steplength schemes in their EM acceleration work. Here, we have chosen method = 2 as the "default" method, as it generally performed 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 will not work. Argument control has the following components:

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 sane. 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 evaluated at `par`, 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` indicates error in function evaluation; `3` is failure due to exceeding 100 steplength reductions in line-search; `4` denotes failure due to an anomalous iteration; 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, and M Raydan (2003), Nonmonotone spectral methods for large-scale nonlinear systems, Optimization Methods and Software, 18, 583-599.

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)
    sane(par=p0, fn=trigexp)
    sane(par=p0, fn=trigexp, method=1)    
######################################
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, 10))
  sane(par=p0, fn=brent, control=list(trace=FALSE))
  sane(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/sane.Rd_%03d_medium.png", width=480, height=480)
> ### Name: sane
> ### Title: Solving Large-Scale Nonlinear System of Equations
> ### Aliases: sane
> ### 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)
>     sane(par=p0, fn=trigexp)
Iteration:  0  ||F(x0)||:  17.35215 

 iteration:  10  ||F(xn)|| =   0.0559578 

 iteration:  20  ||F(xn)|| =   8.366857e-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] 4.009437e-08

$fn.reduction
[1] 122.6983

$feval
[1] 53

$iter
[1] 26

$convergence
[1] 0

$message
[1] "Successful convergence"

>     sane(par=p0, fn=trigexp, method=1)    
Iteration:  0  ||F(x0)||:  17.35215 

 iteration:  10  ||F(xn)|| =   0.01729078 

 iteration:  20  ||F(xn)|| =   3.414079e-06 
$par
 [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
 [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[15] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[22] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[29] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[36] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[43] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[50] 0.9999999

$residual
[1] 7.965319e-08

$fn.reduction
[1] 122.6983

$feval
[1] 47

$iter
[1] 23

$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, 10))
>   sane(par=p0, fn=brent, control=list(trace=FALSE))
$par
 [1]  0.9697799  1.6932033  2.3260895  2.9062584  3.4501947  3.9669567
 [7]  4.4622043  4.9397778  5.4024377  5.8522544  6.2908326  6.7194489
[13]  7.1391405  7.5507652  7.9550431  8.3525860  8.7439200  9.1295014
[19]  9.5097296  9.8849567 10.2554952 10.6216240 10.9835936 11.3416299
[25] 11.6959379 12.0467038 12.3940980 12.7382764 13.0793826 13.4175490
[31] 13.7528980 14.0855433 14.4155904 14.7431376 15.0682770 15.3910945
[37] 15.7116707 16.0300812 16.3463974 16.6606862 16.9730112 17.2834321
[43] 17.5920056 17.8987853 18.2038224 18.5071650 18.8088593 19.1089492
[49] 19.4074763 19.7044805

$residual
[1] 9.308635e-08

$fn.reduction
[1] 325.0439

$feval
[1] 1821

$iter
[1] 752

$convergence
[1] 0

$message
[1] "Successful convergence"

>   sane(par=p0, fn=brent, control=list(M=200, trace=FALSE))
$par
 [1]  0.9697794  1.6932026  2.3260886  2.9062572  3.4501934  3.9669552
 [7]  4.4622026  4.9397760  5.4024358  5.8522524  6.2908306  6.7194468
[13]  7.1391384  7.5507630  7.9550408  8.3525837  8.7439177  9.1294991
[19]  9.5097274  9.8849545 10.2554929 10.6216217 10.9835914 11.3416278
[25] 11.6959358 12.0467018 12.3940960 12.7382744 13.0793807 13.4175472
[31] 13.7528963 14.0855416 14.4155888 14.7431361 15.0682756 15.3910931
[37] 15.7116694 16.0300800 16.3463963 16.6606852 16.9730103 17.2834313
[43] 17.5920048 17.8987847 18.2038218 18.5071645 18.8088590 19.1089489
[49] 19.4074761 19.7044804

$residual
[1] 9.639547e-08

$fn.reduction
[1] 325.0439

$feval
[1] 1649

$iter
[1] 822

$convergence
[1] 0

$message
[1] "Successful convergence"

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