1 #ifndef STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP
2 #define STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP
3
4 #include <stan/math/prim/fun/add_diag.hpp>
5 #include <stan/math/rev/core.hpp>
6 #include <stan/math/prim/err.hpp>
7 #include <stan/math/prim/fun/Eigen.hpp>
8 #include <stan/math/rev/fun/inverse.hpp>
9
10 namespace stan {
11 namespace math {
12
13 namespace internal {
14 /*
15 * Reverse mode specialization of calculating the generalized inverse of a
16 * matrix.
17 * <ul><li> Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses
18 * and Nonlinear Least Squares Problems Whose Variables Separate. <i>SIAM
19 * Journal on Numerical Analysis</i>, Vol. 10, No. 2 (Apr., 1973), pp.
20 * 413-432</li></ul>
21 */
22 template <typename T1, typename T2>
generalized_inverse_lambda(T1 & G_arena,T2 & inv_G)23 inline auto generalized_inverse_lambda(T1& G_arena, T2& inv_G) {
24 return [G_arena, inv_G]() mutable {
25 G_arena.adj()
26 += -(inv_G.val_op().transpose() * inv_G.adj_op()
27 * inv_G.val_op().transpose())
28 + (-G_arena.val_op() * inv_G.val_op()
29 + Eigen::MatrixXd::Identity(G_arena.rows(), inv_G.cols()))
30 * inv_G.adj_op().transpose() * inv_G.val_op()
31 * inv_G.val_op().transpose()
32 + inv_G.val_op().transpose() * inv_G.val_op()
33 * inv_G.adj_op().transpose()
34 * (-inv_G.val_op() * G_arena.val_op()
35 + Eigen::MatrixXd::Identity(inv_G.rows(), G_arena.cols()));
36 };
37 }
38 } // namespace internal
39
40 /*
41 * Reverse mode specialization of calculating the generalized inverse of a
42 * matrix.
43 *
44 * @param G specified matrix
45 * @return Generalized inverse of the matrix (an empty matrix if the specified
46 * matrix has size zero).
47 *
48 * @note For the derivatives of this function to exist the matrix must be
49 * of constant rank.
50 * Reverse mode differentiation algorithm reference:
51 *
52 * <ul><li> Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses
53 * and Nonlinear Least Squares Problems Whose Variables Separate. <i>SIAM
54 * Journal on Numerical Analysis</i>, Vol. 10, No. 2 (Apr., 1973), pp.
55 * 413-432</li></ul>
56 *
57 * Equation 4.12 in the paper
58 *
59 * See also
60 * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse
61 *
62 */
63 template <typename VarMat, require_rev_matrix_t<VarMat>* = nullptr>
generalized_inverse(const VarMat & G)64 inline auto generalized_inverse(const VarMat& G) {
65 using ret_type = promote_var_matrix_t<VarMat, VarMat>;
66
67 if (G.size() == 0)
68 return ret_type(G);
69
70 if (G.rows() == G.cols()) {
71 arena_t<VarMat> G_arena(G);
72 Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd>
73 complete_ortho_decomp_G
74 = G_arena.val().completeOrthogonalDecomposition();
75 if (!(complete_ortho_decomp_G.rank() < G.rows())) {
76 return ret_type(inverse(G_arena));
77 } else {
78 arena_t<ret_type> inv_G(complete_ortho_decomp_G.pseudoInverse());
79 reverse_pass_callback(
80 internal::generalized_inverse_lambda(G_arena, inv_G));
81 return ret_type(inv_G);
82 }
83 } else if (G.rows() < G.cols()) {
84 arena_t<VarMat> G_arena(G);
85 arena_t<ret_type> inv_G((G_arena.val_op() * G_arena.val_op().transpose())
86 .ldlt()
87 .solve(G_arena.val_op())
88 .transpose());
89 reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G));
90 return ret_type(inv_G);
91 } else {
92 arena_t<VarMat> G_arena(G);
93 arena_t<ret_type> inv_G((G_arena.val_op().transpose() * G_arena.val_op())
94 .ldlt()
95 .solve(G_arena.val_op().transpose()));
96 reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G));
97 return ret_type(inv_G);
98 }
99 }
100
101 } // namespace math
102 } // namespace stan
103 #endif
104