R version 4.4.0 alpha (2024-03-26 r86209 ucrt)
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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.

> ### BFGSR-related tests
> 
> ## minimize quadratic form t(D) %*% W %*% D with p.d. weight matrix
> ## (ie unbounded problems).
> set.seed(0)
> options(digits=4)
> quadForm <- function(D) {
+    C <- c(1, 2, 3)  # 3-dimensional case
+    return( - t(D - C) %*% W %*% (D - C) )
+ }
> ## a) test quadratic function t(D) %*% D
> library(maxLik)
Loading required package: miscTools

Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.

If you have questions, suggestions, or comments regarding the 'maxLik' package, please use its GitHub page and the issues there:
https://github.com/arne-henningsen/maxLik
> W <- diag(3)
> D <- rep(1/3, 3)
> res <- maxBFGSR(quadForm, start=D)
> coef(res)
[1] 1 2 3
> all.equal(coef(res), c(1, 2, 3), tolerance=1e-4)
[1] TRUE
> all.equal(gradient(res), rep(0, 3), tolerance=1e-3)
[1] TRUE
> all.equal(nIter(res) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(res) < 4, TRUE)
[1] TRUE
> 
> ## Next, optimize hat function in non-concave region.  Does not work well.
> hat <- function(param) {
+    ## Hat function.  Hessian negative definite if sqrt(x^2 + y^2) < 0.5
+    x <- param[1]
+    y <- param[2]
+    exp(-(x-2)^2 - (y-2)^2)
+ }
> 
> hatNC <- maxBFGSR(hat, start=c(1,1), tol=0, reltol=0)
> all.equal(coef(hatNC), rep(2,2), tolerance=1e-4)
[1] TRUE
> all.equal(gradient(hatNC), rep(0,2), tolerance=1e-3)
[1] TRUE
> all.equal(nIter(hatNC) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(hatNC) < 4, TRUE)
[1] TRUE
> 
> ## Test BFGSR with fixed parameters and equality constraints
> ## Optimize 3D hat with one parameter fixed (== 2D hat).
> ## Add an equality constraint on that
> hat3 <- function(param) {
+    ## Hat function.  Hessian negative definite if sqrt((x-2)^2 + (y-2)^2) < 0.5
+    x <- param[1]
+    y <- param[2]
+    z <- param[3]
+    exp(-(x-2)^2-(y-2)^2-(z-2)^2)
+ }
> sv <- c(x=1,y=1,z=1)
> ## constraints: x + y + z = 8, z = 1
> ## solution: x = 3.5, y = 3.5, z = 1
> A <- matrix(c(1,1,1), 1, 3)
> B <- -8
> constraints <- list(eqA=A, eqB=B)
> hat3CF <- maxBFGSR(hat3, start=sv, constraints=constraints, fixed=3)
> all.equal(coef(hat3CF), c(x=3.5, y=3.5, z=1), tolerance=1e-4)
[1] TRUE
> all.equal(nIter(hat3CF) < 100, TRUE)
[1] TRUE
> all.equal(returnCode(hat3CF) < 4, TRUE)
[1] TRUE
> all.equal(sum(coef(hat3CF)), 8, tolerance=1e-4)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   1.34    0.15    1.48