docs: def _pseudoinverse(self, A, tol=1.0e-10):
          """
          Compute the Moore-Penrose pseudoinverse.
  
          REQUIRED ARGUMENTS
            A (np KxK matrix) - the square matrix whose pseudoinverse is to be computed
  
          RETURN VALUES
            Ainv (np KxK matrix) - the pseudoinverse
  
          OPTIONAL VALUES
            tol - the tolerance (relative to largest magnitude singlular value) below which singular values are to not be include in forming pseudoinverse (default: 1.0e-10)
  
          NOTES
            This implementation is provided because the 'pinv' function of np is broken in the version we were using.
  
          TODO
            Can we get rid of this and use np.linalg.pinv instead?
  
          """
  
          # DEBUG
          # TODO: Should we use pinv, or _pseudoinverse?
          # return np.linalg.pinv(A)
  
          # Get size
          [M, N] = A.shape
          if N != M:
              raise DataError("pseudoinverse can only be computed for square matrices: dimensions were %d x %d" % (
                  M, N))
  
          # Make sure A contains no nan.
          if(np.any(np.isnan(A))):
              print("attempted to compute pseudoinverse of A =")
              print(A)
              raise ParameterError("A contains nan.")
  
          # DEBUG
          diagonal_loading = False
          if diagonal_loading:
              # Modify matrix by diagonal loading.
              eigs = linalg.eigvalsh(A)
              most_negative_eigenvalue = eigs.min()
              if (most_negative_eigenvalue < 0.0):
                  print("most negative eigenvalue = %e" % most_negative_eigenvalue)
                  # Choose loading value.
                  gamma = -most_negative_eigenvalue * 1.05
                  # Modify Theta by diagonal loading
                  A += gamma * np.eye(A.shape[0])
  
          # Compute SVD of A.
          [U, S, Vt] = linalg.svd(A)
  
          # Compute pseudoinverse by taking square root of nonzero singular
          # values.
          Ainv = np.matrix(np.zeros([M, M], dtype=np.float64))
          for k in range(M):
              if (abs(S[k]) > tol * abs(S[0])):
                  Ainv += (1.0/S[k]) * np.outer(U[:, k], Vt[k, :]).T
  
          return Ainv