[R] Non-Linear Optimization - Query

[Ravi Varadhan]

Hi Paul and Berend,

Here is an example for Broyden's tridiagonal system:

broydt <- function(x, h=2) {
n <- length(x)
f <- rep(NA, n)
f[1] <- ((3 - h*x[1]) * x[1]) - 2*x[2] + 1
tnm1 <- 2:(n-1)
f[tnm1] <- ((3 - h*x[tnm1]) * x[tnm1]) - x[tnm1-1] - 2*x[tnm1+1] + 1
f[n] <- ((3 - h*x[n]) * x[n]) - x[n-1] + 1
f
}
 
	require(BB)
	require(nleqslv)

	set.seed(123)
	
	p0 <- -runif(1000)  # 1000 variables

# Timings on R 2.8.0, windows XP, 1 GB Ram and 3.40 GHz Pentium processor
 	
	system.time(ans.bb <- dfsane(par=p0, fn=broydt,
control=list(trace=FALSE)))[1]

	system.time(ans.nl <- nleqslv(x=p0, fn=broydt))[1] 	

	max(abs(ans.nl$x - ans.bb$par))  # comparing the two solutions	 

> system.time(ans.bb <- dfsane(par=p0, fn=broydt,
control=list(trace=FALSE)))[1]
user.self 
     0.02 
> 
> system.time(ans.nl <- nleqslv(x=p0, fn=broydt))[1] 
user.self 
     8.17 
>
> max(abs(ans.nl$x - ans.bb$par))  # comparing the two solutions 
[1] 7.179915e-08
>
>

For large systems, O(10^3), such as arising in discretized solutions of
nonlinear, partial differential equations, methods that do not require
numerical jacobians, or can handle specialized structure will be certainly
faster.  For smaller problems, nleqslv() would be faster than dfsane(), but
computational effort is not a serious issue for smaller problems.  Note that
dfsane() does not use the numerical Jacobian at all - it is derivative free.
Another important point is that dfsane() is coded entirely in R, and is
still quite fast. This is because of the efficiency of the underlying
spectral gradient (Barzilai-Borwein steplength) algorithms. 

In any case, it is nice to have two very different approaches for solving
nonlinear systems, one implementing the more recent, derivative-free,
spectral gradient method, and the other based on the more classical,
Broyden's and Newton's methods.  When I first submitted the BB package to
CRAN in April 2008, there was not a single option available to the R users
for "directly" solving nonlinear systems.  Now, there are two solid
alternatives.

Best,
Ravi.  
  


----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Paul Smith
Sent: Wednesday, March 18, 2009 10:19 AM
To: r-help at r-project.org
Subject: Re: [R] Non-Linear Optimization - Query

On Tue, Mar 17, 2009 at 7:55 PM, Berend Hasselman <bhh at xs4all.nl> wrote:
> You can also try my package "nleqslv" for solving systems of non 
> linear equations (using Broyden or Newton with a selection of global
strategies).
>
> library(nleqslv)
>
> xinit <- rep(1,3)               # or rep(0,3) for a singular start
> nleqslv(xinit,f2,control=list(trace=1))  # f2 defined as above
>
> (You don't need the trace=1; but especially when doing first 
> explorations of a system it can come in handy).

Is your package also able to deal with large systems, Berend?

Paul

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.