gen_inv.hpp

#ifndef STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP
#define STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP

auto f = [](const auto& a) {
  return stan::math::generalized_inverse(a);
};
Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10);
test_ad(f, A);
#include <stan/math/rev/core.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/rev/fun/cholesky_decompose.hpp>
#include <stan/math/rev/fun/transpose.hpp>
#include <stan/math/rev/fun/tcrossprod.hpp>
#include <stan/math/rev/fun/crossprod.hpp>
#include <stan/math/rev/fun/inverse_spd.hpp>
#include <stan/math/rev/fun/mdivide_left_spd.hpp>
#include <stan/math/rev/fun/mdivide_right_spd.hpp>
#include <stan/math/rev/fun/inverse.hpp>

namespace stan {
namespace math {

/**
 * Returns the Moore-Penrose generalized inverse of the specified matrix.
 *
 * The method is based on the Cholesky computation of the transform as specified
 in
 *
 * <ul><li> Courrieu, Pierre. 2008.  Fast Computation of Moore-Penrose Inverse
 Matrices.
 * <i>arXiv</i> <b>0804.4809</b> </li></ul>
 *
 * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase)
 *
 * @param A specified matrix
 * @return Aeneralized inverse of the matrix (an empty matrix if the specified
 * matrix has size zero).
 * @note Because the method inverts a SPD matrix internally that interal matrix
 may result in small numerical issues that result in a non-SPD error. There are
 two
 * <code>generalized_inverse</code> functions, one that uses one input matrix
 (this one)
 * and another that works with an input matrix and a small jitter to the
 diagonal of the internal SPD
 * matrix.
 */
template <typename EigMat, require_eigen_vt<is_var, EigMat>* = nullptr>
inline auto generalized_inverse(const EigMat& A) {
  using value_t = value_type_t<EigMat>;
  if (A.size() == 0) {
    return {};
  }

  const auto n = A.rows();
  const auto m = A.cols();

  if (A.rows() == A.cols()) {
    return inverse(A);
  }
  if (n < m) {
    arena_t<plain_type_t<EigMat>> A_arena(A);
    arena_t<EigMat> inv_A(mdivide_left_spd(tcrossprod(A_arena.val()),
      A_arena.val()));
    reverse_pass_callback([A_arena, inv_A]() mutable {
      A_arena.adj() += -inv_A *  * inv_A + 
    });
    return ret;
  } else {
    return transpose(mdivide_right_spd(A, crossprod(A)));
  }
}

/**
 * Returns the Moore-Penrose generalized inverse of the specified matrix.
 *
 * The method is based on the Cholesky computation of the transform as specified
 in
 *
 * <ul><li> Courrieu, Pierre. 2008.  Fast Computation of Moore-Penrose Inverse
 Matrices.
 * <i>arXiv</i> <b>0804.4809</b> </li></ul>
 *
 * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase)
 *
 * @param A specified matrix
 * @param a real constant
 * @return Aeneralized inverse of the matrix (an empty matrix if the specified
 * matrix has size zero).
 * @note Because the method inverts a SPD matrix internally that interal matrix
 may result in small numerical issues that result in a non-SPD error. There are
 two
 * <code>generalized_inverse</code> functions, one that uses one input matrix
 * and another (this one) that works with an input matrix and a small jitter to
 the diagonal of the internal SPD
 * matrix.
 */
template <typename EigMat, typename Scal, require_eigen_t<EigMat>* = nullptr,
          require_stan_scalar_t<Scal>* = nullptr>
inline Eigen::Matrix<value_type_t<EigMat>, EigMat::RowsAtCompileTime,
                     EigMat::ColsAtCompileTime>
generalized_inverse(const EigMat& A, const Scal& a) {
  using value_t = value_type_t<EigMat>;
  if (A.size() == 0) {
    return {};
  }

  const auto n = A.rows();
  const auto m = A.cols();

  if (A.rows() == A.cols()) {
    return inverse(A);
  }

  if (n < m) {
    Eigen::Matrix<value_t, Eigen::Dynamic, Eigen::Dynamic> A = tcrossprod(A);
    A.diagonal().array() += a;
    return transpose(mdivide_left_spd(A, A));
  } else {
    Eigen::Matrix<value_t, Eigen::Dynamic, Eigen::Dynamic> A = crossprod(A);
    A.diagonal().array() += a;
    return transpose(mdivide_right_spd(A, A));
  }
}

}  // namespace math
}  // namespace stan

#endif