Maybe bugs in Jacobian function

[ctfysh]

Below is the code:

library(pnd)
library(numDeriv)

f <- function(x) c(sin(x[1]) + cos(x[2]), exp(x[1]*x[2]))
x0 <- c(1, 1)

cat("=== Jacobian Computation Comparison ===\n\n")

cat("1. numDeriv::jacobian (Richardson method):\n")
J_numDeriv <- numDeriv::jacobian(f, x0, method = "Richardson")
print(J_numDeriv)
cat("\n")

cat("2. pnd::Jacobian:\n")
# pnd::Jacobian prints output but may not return a matrix properly
# Capture output and try to extract numeric values
J_pnd <- tryCatch({
  capture.output(result <- pnd::Jacobian(f, x0), type = "output")
  if (is.numeric(result)) {
    as.matrix(result)
  } else {
    NA
  }
}, error = function(e) {
  cat("  Error:", e$message, "\n")
  NA
})

if (is.matrix(J_pnd) && is.numeric(J_pnd) && !anyNA(J_pnd)) {
  print(J_pnd)
} else {
  cat("  pnd::Jacobian could not compute valid Jacobian\n")
  J_pnd <- NULL
}
cat("\n")

vectorized_Jacobian <- function(f, x0, ...) {
  if (!is.numeric(x0)) {
    stop("x0 must be a numeric vector")
  }
  if (!is.function(f)) {
    stop("f must be a function")
  }

  f_x0 <- f(x0)
  if (!is.numeric(f_x0)) {
    stop("f(x0) must return a numeric vector")
  }
  out_dim <- length(f_x0)

  in_dim <- length(x0)

  J <- vapply(1:out_dim, function(i) {
    f_i <- function(x) f(x)[i]
    pnd::Grad(f_i, x0, ...)
  }, numeric(in_dim))

  return(t(J))
}

cat("3. vectorized_Jacobian (optimized with vapply):\n")
J_vectorized <- vectorized_Jacobian(f, x0)
print(J_vectorized)
cat("\n")

cat("=== Verification ===\n")
cat("Max absolute difference between numDeriv and vectorized:",
    max(abs(J_numDeriv - J_vectorized)), "\n")

if (!is.null(J_pnd) && is.matrix(J_pnd) && !anyNA(J_pnd)) {
  cat("Max absolute difference between numDeriv and pnd:",
      max(abs(J_numDeriv - J_pnd)), "\n")
  cat("Max absolute difference between pnd and vectorized:",
      max(abs(J_pnd - J_vectorized)), "\n")
} else {
  cat("pnd::Jacobian result not available for comparison\n")
}

Results summary:
Output:

1. numDeriv::jacobian (Richardson method):
   [,1] [,2]
   [1,] 0.5403023 -0.841471
   [2,] 2.7182818 2.718282
2. pnd::Jacobian:
   pnd::Jacobian could not compute valid Jacobian
3. vectorized_Jacobian (optimized with vapply):
   [,1] [,2]
   [1,] 0.5403023 -0.841471
   [2,] 2.7182818 2.718282

[Fifis]

@ctfysh I have pushed a commit that should fix this undesirable behaviour. Thank you for the report. Pull the latest dev version from GitHub (in the code below) and see if this works (should work!):

devtools::install_github("Fifis/pnd")
library(pnd)
library(numDeriv)

f <- function(x) c(sin(x[1]) + cos(x[2]), exp(x[1]*x[2]))
x0 <- c(1, 1)
jac.true <- matrix(c(cos(1), exp(1), -sin(1), exp(1)), ncol = 2)

J_numDeriv <- numDeriv::jacobian(f, x0, method = "Richardson")
print(J_numDeriv)

J_pnd     <- matrix(as.numeric(pnd::Jacobian(f, x0)), ncol = 2)
J_pnd_acc <- matrix(as.numeric(pnd::Jacobian(f, x0, acc.order = 8)), ncol = 2)
print(J_pnd)
print(J_pnd_acc)

jac.true - J_numDeriv
jac.true - J_pnd
jac.true - J_pnd_acc  # Much lower error than other methods

If this works fine on your machine, I shall push a new version to CRAN.