code: def cho_solve_ATAI(A, rho, b, c, lwr, check_finite=True):
      r
  
      N, M = A.shape
      if N >= M:
          x = linalg.cho_solve((c, lwr), b, check_finite=check_finite)
      else:
          x = (b - A.T.dot(linalg.cho_solve((c, lwr), A.dot(b),
                                            check_finite=check_finite))) / rho
      return x
docs: r"""
      Solve the linear system :math:`(A^T A + \rho I)\mathbf{x} = \mathbf{b}`
      or :math:`(A^T A + \rho I)X = B` using :func:`scipy.linalg.cho_solve`.
  
      Parameters
      ----------
      A : array_like
        Matrix :math:`A`
      rho : float
        Scalar :math:`\rho`
      b : array_like
        Vector :math:`\mathbf{b}` or matrix :math:`B`
      c : array_like
        Matrix containing lower or upper triangular Cholesky factor,
        as returned by :func:`scipy.linalg.cho_factor`
      lwr : bool
        Flag indicating whether the factor is lower or upper triangular
  
      Returns
      -------
      x : ndarray
        Solution to the linear system