  In scientific computing and machine learning applications, matrices and more
general multidimensional arrays (tensors) can often be approximated with the
help of low-rank decompositions. Since matrices and tensors of fixed rank form
smooth Riemannian manifolds, one of the popular tools for finding the low-rank
approximations is to use the Riemannian optimization. Nevertheless, efficient
implementation of Riemannian gradients and Hessians, required in Riemannian
optimization algorithms, can be a nontrivial task in practice. Moreover, in
some cases, analytic formulas are not even available. In this paper, we build
upon automatic differentiation and propose a method that, given an
implementation of the function to be minimized, efficiently computes Riemannian
gradients and matrix-by-vector products between approximate Riemannian Hessian
and a given vector.
