code: def stable_cholesky(M, add_to_diag_till_psd=True):
      """Returns L, a 'stable' cholesky decomposition of M. L is lower triangular and
      satisfies L*L' = M.
      Sometimes nominally psd matrices are not psd due to numerical issues. By adding a
      small value to the diagonal we can make it psd. This is what this function does.
      Use this iff you know that K should be psd. We do not check for errors.
      """
      _printed_warning = False
      if M.size == 0:
          return M  # if you pass an empty array then just return it.
      try:
          # First try taking the Cholesky decomposition.
          L = np.linalg.cholesky(M)
      except np.linalg.linalg.LinAlgError as e:
          # If it doesn't work, then try adding diagonal noise.
          if not add_to_diag_till_psd:
              raise e
          diag_noise_power = -11
          max_M = np.diag(M).max()
          diag_noise = np.diag(M).max() * 1e-11
          chol_decomp_succ = False
          while not chol_decomp_succ:
              try:
                  diag_noise = (10**diag_noise_power) * max_M
                  L = np.linalg.cholesky(
                      M + diag_noise * np.eye(M.shape[0])
                  )
                  chol_decomp_succ = True
              except np.linalg.linalg.LinAlgError:
                  if diag_noise_power > -9 and not _printed_warning:
                      from warnings import warn
  
                      warn(
                          (
                              "Could not compute Cholesky decomposition despite adding %0.4f to the "
                              "diagonal. This is likely because the M is not positive semi-definite."
                          )
                          % ((10**diag_noise_power) * max_M)
                      )
                      _printed_warning = True
                  diag_noise_power += 1
              if diag_noise_power >= 5:
                  raise ValueError(
                      (
                          "Could not compute Cholesky decomposition despite adding"
                          + " %0.4f to the diagonal. This is likely because the M is not "
                          + "positive semi-definite or has infinities/nans."
                      )
                      % (diag_noise)
                  )
      return L
  
imports: from __future__ import print_function
  from __future__ import division
  import numpy
  from scipy.sparse import dok_matrix
  from scipy.linalg import solve_triangular
label_desc: Use image node_id 13 for calling a global function with example usage: stable_cholesky(M, add_to_diag_till_psd) and returns: L, M