# Large NLS example
# (https://www.gnu.org/software/gsl/doc/html/nls.html#large-nonlinear-least-squares-example)
## number of parameters
p <- 250
## model function
f <- function(theta) {
c(sqrt(1e-5) * (theta - 1), sum(theta^2) - 0.25)
}
## jacobian function
jac <- function(theta) {
rbind(diag(sqrt(1e-5), nrow = length(theta)), 2 * t(theta))
}
## dense Levenberg-Marquardt
# \donttest{
gsl_nls_large(
fn = f, ## model
y = rep(0, p + 1), ## (dummy) responses
start = 1:p, ## start values
algorithm = "lm", ## algorithm
jac = jac, ## jacobian
control = list(maxiter = 250)
)
# }
## dense Steihaug-Toint conjugate gradient
# \donttest{
gsl_nls_large(
fn = f, ## model
y = rep(0, p + 1), ## (dummy) responses
start = 1:p, ## start values
jac = jac, ## jacobian
algorithm = "cgst" ## algorithm
)
# }
## sparse Jacobian function
jacsp <- function(theta) {
rbind(Matrix::Diagonal(x = sqrt(1e-5), n = length(theta)), 2 * t(theta))
}
## sparse Levenberg-Marquardt
gsl_nls_large(
fn = f, ## model
y = rep(0, p + 1), ## (dummy) responses
start = 1:p, ## start values
algorithm = "lm", ## algorithm
jac = jacsp, ## sparse jacobian
control = list(maxiter = 250)
)
## sparse Steihaug-Toint conjugate gradient
gsl_nls_large(
fn = f, ## model
y = rep(0, p + 1), ## (dummy) responses
start = 1:p, ## start values
jac = jacsp, ## sparse jacobian
algorithm = "cgst" ## algorithm
)
Run the code above in your browser using DataLab