0

I am trying to solve a penalty function utilizing a newton method algorithm. It uses the hessian and a transposed gradient to solve a linear system that will give me approximate values.

The problem with my code starts at the tenth iteration (but it may vary depending on the initial guess, which I call 'chute'), because the hessian turns out to be a singular matrix. What should be done in this case? Is there a problem with my algorithm?

I'm not actually pasting the function to be optimized because it is very extensive, but if otherwise solicited I'll do it.

I have tried different functions, but the only correct answer I got from the algorithm was for the function F(x, y) = x^2 + y^2. (x, y) = (0,0).

for(i in 1:n){
  hessian <- hessian(f = f, x = chute, centered = TRUE)
  gradient.transposed <- -1 * t(gradient(f = f, x = chute))

  # solving the linear system
  sk <- qr.solve(hessian, gradient.transposed)
  xk.1 <- xk + sk
  k[[i]] <- xk.1

  # calculating a test to verify the solution
  first.test <- (norm(xk - xk.1) < epsilon * (1 + norm(xk)))
  g[[i]] <- first.test
  if(first.test == TRUE){
    # calculating a test to verify the solution
    second.test <- (norm(gradient.transposed) <= delta * (1 + abs(f(xk))))
    h[[i]] <- second.test
    if(second.test == TRUE){
      root.approx <- tail(k)[[1]]
      res <- list('aprox' = root.approx, 'iter' = k)
      return(res)
    }
  }

  xk <- xk.1

  chute <- c(xk.1)
}

The error message I try to get past certain iterations is:

"Error in qr.solve(hessian, gradient.transposed) : singular matrix 'a' in solve"

Thank you so much for the attention and future answers.

Edit: I'm utilizing the rootSolve package to find numeric hessians and gradients. Also fixed some typos.

Bruno Nascimento
  • 15
  • 5
  • singularity means no inverse of the matrix due to 0 determinant. Maybe try jittering (adding a bit of noise) to your points...? – Sotos May 17 '19 at 13:12
  • You might want to skip the error and continue iterating. [This question](https://stackoverflow.com/questions/14748557/skipping-error-in-for-loop#14749552) might help. In your case, just check whether your hessian matrix (or hessiana?) is inversable, if not, use a `skip` statement. You could also use the `tryCatch()` function. – Trusky May 17 '19 at 14:45

0 Answers0