code: def lu_solve_ATAI(A, rho, b, lu, piv, check_finite=True):
      r
  
      N, M = A.shape
      if N >= M:
          x = linalg.lu_solve((lu, piv), b, check_finite=check_finite)
      else:
          x = (b - A.T.dot(linalg.lu_solve((lu, piv), A.dot(b), 1,
                                           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.lu_solve`.
  
      Parameters
      ----------
      A : array_like
        Matrix :math:`A`
      rho : float
        Scalar :math:`\rho`
      b : array_like
        Vector :math:`\mathbf{b}` or matrix :math:`B`
      lu : array_like
        Matrix containing U in its upper triangle, and L in its lower triangle,
        as returned by :func:`scipy.linalg.lu_factor`
      piv : array_like
        Pivot indices representing the permutation matrix P, as returned by
        :func:`scipy.linalg.lu_factor`
      check_finite : bool, optional (default False)
        Flag indicating whether the input array should be checked for Inf
        and NaN values
  
      Returns
      -------
      x : ndarray
        Solution to the linear system