code: def batch_cholesky_numerical_grad(x, vgrads, upper, *args, **kwargs):
      # Our implementation of cholesky decomposition returns symmetric grad instead of lower triangular grad
      from scipy.optimize import approx_fprime
  
      def func(vector):
          matrix = np.reshape(vector, x.shape)
          L = batch_cholesky(matrix, upper) * vgrads
          return np.sum(L)
  
      x_c = np.concatenate([i.flatten() for i in x if i is not None])
      epsilon = 1e-4
      n_grad = approx_fprime(x_c, func, epsilon)
      n_grad = np.reshape(n_grad, x.shape)
      n_grad = (n_grad + np.transpose(n_grad, axes=(0, 2, 1))) * 0.5
      return n_grad.flatten()
  
imports: import pytest
  import numpy
  import nnabla.functions
  from nbla_test_utils import list_context, function_tester
label_desc: Use image node_id 2 for calling a global function with example usage: batch_cholesky_numerical_grad(x, vgrads, upper) and returns: n_grad, np