code: def Gf(T, ff):
      """
      Subroutine for the gradient of f using numerical derivatives.
      """
      k = T.shape[0]
      ep = 1e-4
      G = np.zeros((k, k))
      for r in range(k):
          for s in range(k):
              dT = np.zeros((k, k))
              dT[r, s] = ep
              G[r, s] = (ff(T + dT) - ff(T - dT)) / (2 * ep)
      return G
  
file_path: statsmodels/statsmodels/multivariate/factor_rotation/_gpa_rotation.py
imports: import numpy
label_desc: Use image node_id 2 for calling a global function with example usage: Gf(T, ff) and returns: G