  This paper is concerned with the efficient evaluation of higher-order
derivatives of functions $f$ that are composed of matrix operations. I.e., we
want to compute the $D$-th derivative tensor $\nabla^D f(X) \in \mathbb
R^{N^D}$, where $f:\mathbb R^{N} \to \mathbb R$ is given as an algorithm that
consists of many matrix operations. We propose a method that is a combination
of two well-known techniques from Algorithmic Differentiation (AD): univariate
Taylor propagation on scalars (UTPS) and first-order forward and reverse on
matrices. The combination leads to a technique that we would like to call
univariate Taylor propagation on matrices (UTPM). The method inherits many
desirable properties: It is easy to implement, it is very efficient and it
returns not only $\nabla^D f$ but yields in the process also the derivatives
$\nabla^d f$ for $d \leq D$. As performance test we compute the gradient
$\nabla f(X)$ % and the Hessian $\nabla_A^2 f(A)$ by a combination of forward
and reverse mode of $f(X) = \trace (X^{-1})$ in the reverse mode of AD for $X
\in \mathbb R^{n \times n}$. We observe a speedup of about 100 compared to
UTPS. Due to the nature of the method, the memory footprint is also small and
therefore can be used to differentiate functions that are not accessible by
standard methods due to limited physical memory.
