# Higher-order derivatives of the QR and of the real symmetric eigenvalue
decomposition in forward and reverse mode algorithmic differentiation

Sebastian F. Walter<EMAIL_ADDRESS>Department of Mathematics,
Faculty of Mathematics and Natural Sciences II,
Humboldt-Universität zu Berlin, Rudower Chaussee 25, Adlershof, 12489 Berlin,
mail address: Unter den Linden 6, 10099 Berlin, Germany
Lutz Lehmann  Department of Mathematics, Faculty of Mathematics and Natural
Sciences II,
Humboldt-Universität zu Berlin, Rudower Chaussee 25, Adlershof, 12489 Berlin,
mail address: Unter den Linden 6, 10099 Berlin, Germany
René Lamour  Department of Mathematics, Faculty of Mathematics and Natural
Sciences II,
Humboldt-Universität zu Berlin, Rudower Chaussee 25, Adlershof, 12489 Berlin,
mail address: Unter den Linden 6, 10099 Berlin, Germany

###### Abstract

We address the task of higher-order derivative evaluation of computer programs
that contain QR decompositions and real symmetric eigenvalue decompositions.
The approach is a combination of univariate Taylor polynomial arithmetic and
matrix calculus in the (combined) forward/reverse mode of Algorithmic
Differentiation (AD). Explicit algorithms are derived and presented in an
accessible form. The approach is illustrated via examples.

Keywords: univariate Taylor polynomial arithmetic; higher-order derivatives;
QR decomposition; real symmetric eigenvalue decomposition; algorithmic
differentiation; automatic differentiation

Classification: 65D25; 12E05; 58C15; 65F15; 65F25

##  1 Introduction and related work

This paper is concerned with the efficient evaluation of higher-order
derivatives of functions $F:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{M}$ which
are implemented as computer programs that contain numerical linear algebra
functions like the QR or the real symmetric eigenvalue decomposition.

Traditionally, Algorithmic Differentiation (AD) tools like ADOL-C [GJM+99] or
CppAD [Bel10] regard the functions defined in the C header file math.h as
elementary functions. In the forward mode of AD, their approach to compute
higher-order derivatives is to generalize from real arithmetic to univariate
Taylor polynomial (UTP) arithmetic [GJM+99, GUW00, GW08]. For the reverse mode
of AD, the program evaluation is traced and stored in a computational graph or
on a sequential tape. During the so-called reverse sweep the stored
intermediate values are retrieved and used to compute derivatives (c.f.
Section 4).

As explained in Section 3, the functions in math.h suffice since all
computable functions are a concatenation of these functions. However, working
only at the expression level has also disadvantages since no global knowledge
of the function’s structure can be used. A particularly important class of
algorithms in science and engineering are numerical linear algebra (NLA)
functions. Though NLA functions are typically locally smooth, their
implementations often contain non-differentiable operations and program
branches. If no special care is taken, this may result in incorrect
computations of derivatives. Also, many NLA functions on
${\mathbb{R}}^{N\times N}$ matrices require $\mathcal{O}(N^{3})$ arithmetic
operations. Since during the reverse mode intermediate results are required,
this would yield an $\mathcal{O}(N^{3})$ memory requirement. Though it may be
possible to adapt codes to yield reduced memory requirements, as for instance
reported for the LU decomposition [Gri03], in practice it can be a cumbersome
and error-prone process. Also one would like to reuse existing, high-
performance implementations of NLA algorithms. Adding the NLA functions to the
list of elementary functions circumvents this problem. This has been realized
before [Büc02, BH96] and also UTP algorithms for some NLA functions (e.g. the
solution of linear equations) have been implemented in software [Eri03].

The contribution of this paper is to provide explicit algorithms for UTP
arithmetic applied to the QR decomposition and the real symmetric eigenvalue
decomposition. Note that our approach to the real symmetric eigenvalue
decomposition is similar to [AT98, vdAMM07] but our algorithmic result
differs. In addition, we also treat the reverse mode of AD.

The document is structured as follows: In Section 2 we give two application
examples for the algorithms presented in this document, followed by a brief
review of the underlying computational model in Section 3. We shortly describe
the basics of AD in Section 4 where we make use of the results from 3. In
Section 5 we describe the general approach of NLA functions. After that, we
apply the results from Section 4 to find extended functions for the QR and
eigenvalue decomposition in Section 6 and 7 and also provide pullback
algorithms that are necessary in in the reverse mode of AD. Finally, we
present some numerical results in Section 8.

##  2 Application examples for the proposed algorithms

The purpose of this section is to show two application examples where higher-
order derivatives of computer programs that contain the QR and the real
symmetric eigenvalue decomposition are necessary.

###  2.1 Optimum experimental design

The goal in optimum experimental design (OED) is to minimize some cost
function representing the size of the confidence region of parameters of
interest. We consider here a popular formulation where the objective function
$\Phi(C)\in{\mathbb{R}}$ depends on the covariance matrix
$C\in{\mathbb{R}}^{{N_{p}}\times{N_{p}}}$ of a constrained parameter
estimation problem, where the covariance matrix is computed by

$\displaystyle C$ $\displaystyle=$
$\displaystyle(1\\!\\!\\!\mathrm{I},0)\begin{pmatrix}J_{1}^{T}J_{1}&J_{2}^{T}\\\
J_{2}&0\\\ \end{pmatrix}^{-1}\begin{pmatrix}1\\!\\!\\!\mathrm{I}\\\
0\end{pmatrix}\;,$ (2.1)

and we assume that $J_{1}\equiv
J_{1}(q)\in{\mathbb{R}}^{{N_{m}}\times{N_{p}}}$, $J_{2}\equiv
J_{2}(q)\in{\mathbb{R}}^{{N_{\mathrm{r}}}\times{N_{p}}}$,
$1\\!\\!\\!\mathrm{I}\in{\mathbb{R}}^{N_{p}\times N_{p}}$ and
$q\in{\mathbb{R}}^{{N_{q}}}$. The notation is motivated as follows:
$p\in{\mathbb{R}}^{N_{p}}$ are model (pseudo-)parameters,
$q\in{\mathbb{R}}^{{N_{q}}}$ are control variables and $J_{1}$ and $J_{2}$ are
Jacobians of the residuals resp. of the constraint function with respect to
the parameters $p$. Typical choices for cost function $\Phi$ are the trace,
the determinant or the maximum eigenvalue of the covariance matrix $C$. Though
Eqn. (2.1) correctly describes the covariance matrix $C$, the actual
algorithmic implementation is often a code like

$\displaystyle C$ $\displaystyle=$ $\displaystyle
Q_{2}^{T}\left(Q_{2}J_{1}^{T}J_{1}Q_{2}^{T}\right)^{-1}Q_{2}\;,$ (2.2)

where $Q_{2}$ results from a QR-like decomposition of $J_{2}$, i.e.
$J_{2}^{T}=(Q_{1}^{T},Q_{2}^{T})(L,0)^{T}$. The matrices $J_{1}$ and $J_{2}$
are assumed to satisfy the constraint qualification
${\mathrm{rank}\,}(J_{2}))={N_{\mathrm{r}}}$ and the condition
${\mathrm{rank}\,}(J)={N_{p}}$, where $J=(J_{1}^{T},J_{2}^{T})^{T}$. The
matrix $Q_{2}$ spans the nullspace of $J_{2}$. For a detailed discussion we
refer to Körkel [Kör02] and to Bock and Kostina [KKSB07].

Newton-type optimization algorithms require at least the gradient
$\nabla_{q}\Phi(q)$ of the objective function $\Phi$. To obtain good
convergence near the local minimizer, it is often advantageous if exact
second-order derivatives are available. Since the number of controls ${N_{q}}$
can be large, one would like to have the possibility to compute these
derivatives in the reverse mode of AD. Robust objective functions are often
formulated in a way that require third and even higher-order derivatives, so
it is necessary to have algorithms that scale easily to arbitrary order.

Thus, this example requires the differentiation of the nullspace of a matrix,
the matrix product, matrix inversion, the QR decomposition and the objective
function evaluation, e.g. the eigenvalue decomposition.

###  2.2 Index determination of differential algebraic equations

Many dynamical problems in chemical engineering, rigid body mechanics, circuit
simulation and control theory are described by Differential Algebraic
Equations (DAEs) of the form

$0=f\left(\tfrac{{\rm d}}{{\rm d}x}d(y,x),y,x\right),\;x\in
I=[a,b]\subset{\mathbb{R}},$ (2.3)

where $y:{\mathbb{R}}\to{\mathbb{R}}^{m}$ lives in suitable function space,
$f:{\mathbb{R}}^{n}\times{\mathbb{R}}^{m}\times I\to{\mathbb{R}}^{m}$,
$d:{\mathbb{R}}^{m}\times I\to{\mathbb{R}}^{n}$ are sufficiently smooth and
typically $n$ is smaller than $m$.

Using higher–order derivatives of the functions in the DAE one can, in
general, transform the DAE system into an ODE system of order one. The
_differentiation index_ is the highest derivative order required in this
process, that is, derivatives of up to this order of the original equations
are part of any solution of the DAE. The knowledge of the index allows to
estimate the difficulty to solve the DAE.

There are many different index definitions. Here we consider the tractability
index. To compute it, the DAE is linearized along a given function
$\bar{y}(x)$ as

$\underbrace{\tfrac{\partial}{\partial
z}f(\bar{w}(x),\bar{y}(x),x)}_{=A(x)\in{\mathbb{R}}^{m\times n}}\tfrac{{\rm
d}}{{\rm d}x}\Bigl{(}\underbrace{\tfrac{\partial}{\partial
y}d(\bar{y}(x),x)}_{=D(x)\in{\mathbb{R}}^{n\times
m}}z(x)\Bigr{)}+\underbrace{\tfrac{\partial}{\partial
y}f(\bar{w}(x),\bar{y}(x),x)}_{=B(x)\in{\mathbb{R}}^{m\times
m}}z(x)=\underbrace{-f\left(\bar{w}(x),\bar{y}(x),x\right)}_{q(x)}$

with $\bar{w}(x)=\frac{{\rm d}}{{\rm d}x}d(\bar{y}(x),x)$. The coefficient
functions $A=A(x),\;D=D(x)$ and $B=B(x)$ give rise to a matrix sequence

$\displaystyle G_{0}$ $\displaystyle=AD,$ $\displaystyle B_{0}$
$\displaystyle=B,$ $\displaystyle G_{i+1}$ $\displaystyle=G_{i}+B_{i}Q_{i},$
$\displaystyle B_{i+1}$ $\displaystyle=B_{i}P_{i}-G_{i+1}D^{-}\tfrac{{\rm
d}}{{\rm d}x}(DP_{0}\cdots P_{i+1}D^{-})DP_{0}\cdots P_{i},$ (2.4)

where $Q_{i}$ describes a projector onto ${\rm ker}\,G_{i}$, $P_{i}=I-Q_{i}$
and $D^{-}$ is a generalized reflexive inverse of $D$. Now, the tractability
index is the smallest number $\mu\in{\mathbb{N}}$ where $G_{\mu}$ is
nonsingular. The projectors $Q_{i}$ can be determined mainly by use of a QR
decomposition.

A QR decomposition of the potentially singular matrix
$G\in{\mathbb{R}}^{M\times M}$ with ${\mathrm{rank}\,}G=r$ results in

$G\mathnormal{\Pi}=Q\begin{pmatrix}R_{1}&R_{2}\\\ 0&0\end{pmatrix},$

where $\mathnormal{\Pi}$ describes a column permutation matrix, $Q$ an
orthogonal matrix and $R_{1}\in{\mathbb{R}}^{r\times r}$ an upper triangular
matrix. Then a nullspace projector $Q_{G}$ onto ${\rm ker}\,G$ is given by

$Q_{G}=\mathnormal{\Pi}\begin{pmatrix}0&-R^{-1}_{1}R_{2}\\\
0&I\end{pmatrix}\mathnormal{\Pi}^{T}.$

The computation of $B_{i+1}$ via (2.4) needs the differentiation of
$DP_{0}\cdots P_{i+1}D^{-}$ with respect to $x$. Thus, higher–order
derivatives of a function that contains the QR decomposition are necessary.
For a in-depth discussion of index definition of DAEs see März [Mär02, Mär03].

##  3 Computational model

We consider functions

$\displaystyle F:{\mathbb{R}}^{N}$ $\displaystyle\rightarrow$
$\displaystyle{\mathbb{R}}^{M}$ $\displaystyle x$ $\displaystyle\mapsto$
$\displaystyle y=F(x)\;,$

that can be described by the _three-part form_

$\displaystyle v_{n-N}$ $\displaystyle=$ $\displaystyle x_{n}\quad\quad\quad
n=1,\dots,N$ $\displaystyle v_{l}$ $\displaystyle=$
$\displaystyle\phi_{l}(v_{j\prec l})\quad l=1,\dots,L$ $\displaystyle y_{M-m}$
$\displaystyle=$ $\displaystyle v_{L-m}\quad\quad m=M-1,\dots,0\;,$

where $\phi_{l}\in\\{+,-,\cdot,/,\sin,\exp,\dots\\}$ are called _elementary
functions_ , $v_{l}$ are intermediate values and $v_{i\prec l}$ denote the
tuples of input arguments of $\phi_{l}$. For instance the function
$F:{\mathbb{R}}^{2}\rightarrow{\mathbb{R}}$, $x\mapsto
y=F(x)=\sin(x_{1}+\cos(x_{2})*x_{1})$ is described by

independent | $v_{-1}$ | = | $x_{1}$ | = | $3$
---|---|---|---|---|---
independent | $v_{0}$ | = | $x_{2}$ | = | $7$
| $v_{1}$ | = | $\phi_{1}(v_{0})$ | = | $\cos(v_{0})$
| $v_{2}$ | = | $\phi_{2}(v_{1},v_{-1})$ | = | $v_{1}v_{-1}$
| $v_{3}$ | = | $\phi_{3}(v_{-1},v_{2})$ | = | $v_{-1}+v_{2}$
| $v_{4}$ | = | $\phi_{4}(v_{3})$ | = | $\sin(v_{3})$
dependent | $y$ | = | $v_{4}$ |  |

It shows a sequential representation of the computation. Alternatively, one
can describe the function evaluation as composite function

$\displaystyle F(x)$ $\displaystyle=$ $\displaystyle
P_{y}\circ\Phi_{L}\circ\Phi_{L-1}\circ\dots\circ\Phi_{1}\circ P_{x}^{T}(x)\;,$
(3.1)

where $\Phi_{l}:\mathcal{H}\rightarrow\mathcal{H}$, $s^{(l-1)}\mapsto
s^{(l)}=\Phi_{l}(s^{(l-1)})$ are called _elementary transitions_ that operate
the _state space_ $\mathcal{H}$. Each elementary transition can be written as

$\Phi_{l}=P_{l}\circ\phi_{l}\circ
Q_{l}+(1\\!\\!\\!\mathrm{I}-(1-\sigma_{l})P_{l}\circ P_{l}^{T})\;.$ (3.2)

where the functions
$\phi_{l}:\mathcal{D}_{l}\subseteq\mathcal{H}_{l}\rightarrow\mathcal{H}_{l}\in\\{+,-,*,/,\sin,\exp,\dots\\}$
are the elementary functions. The
$Q_{l}:\mathcal{H}\rightarrow\mathcal{H}_{l}$ map to the domains of the
elementary functions and the $P_{l}:\mathcal{H}_{l}\rightarrow\mathcal{H}$ map
back to the overall state space. The functions $P_{x}^{T}$ and $P_{y}$ are
used to map the independent variables $x$ to the state $s^{(0)}$ and $s^{(L)}$
to $y$. The case $\sigma_{l}=1$ corresponds to an augmented assignment
$s_{l}=s_{l}+\phi_{l}(s_{l})$ and $\sigma_{l}=0$ to the usual assignment
$s_{l}=\phi_{l}(s_{l})$. For our purposes it suffices to consider a real
vector space as state space, i.e., the mappings $P_{l}$ and $Q_{l}$ can be
written as matrices. For a more detailed discussion see Griewank [Gri03].

##  4 Algorithmic differentiation

In this section we briefly review some key results from the theory of AD that
will be necessary in Section 6 and 7. For a detailed discussion we refer to
the standard reference “Evaluating Derivatives” by Griewank and Walther
[GW08].

###  4.1 The forward mode

One can use univariate Taylor series expansions to compute higher-order
(directional) derivatives. The basic observation is that given a smooth curve
$x(t)=x_{0}+x_{1}t$ with $t\in(-\epsilon,\epsilon)$, $\epsilon>0$, and a
smooth function $F$ one obtains a smooth curve $y(t)=F(x(t))$ with the Taylor
series expansion

$y(t)=\sum_{d=0}^{D-1}y_{d}t^{d}+\mathcal{O}(t^{D})=\sum_{d=0}^{D-1}\left.\frac{1}{d!}\frac{{\rm
d}^{d}}{{\rm d}t^{d}}F(x(t))\right|_{t=0}t^{d}+\mathcal{O}(t^{D})\;.$ (4.1)

By application of the chain rule one can interpet the terms of the expansion.
The zeroth derivative is the normal function evaluation $y_{0}=F(x_{0})$ and
the first coefficient $y_{1}=\left.\frac{{\rm d}}{{\rm
d}t}F(x(t))\right|_{t=0}=\frac{{\rm d}}{{\rm d}x}F(x_{0})\cdot x_{1}$ is a
directional derivative.

In the context of AD it is advantageous to generalize the notion of Taylor
series expansions to a purely algebraic task. In other words, for arithmetic
with univariate Taylor polynomials (UTP) one extends functions
$F:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{N}$ to functions
$E_{D}(F):{\mathbb{R}}^{N}[T]/(T^{D})\rightarrow{\mathbb{R}}^{M}[T]/(T^{D})$.
We denote representing elements of the polynomial factor ring
${\mathbb{R}}^{N}[T]/(T^{D})$ as

$[x]_{D}:=[x_{1},\dots,x_{D-1}]:=\sum_{d=0}^{D-1}x_{d}T^{d}\;,$ (4.2)

where $x_{d}\in{\mathbb{R}}^{N}$ is called _Taylor coefficient_. The quantity
$T$ is an indeterminate, i.e., a formal variable. The _extended function_
$E_{D}(F)$ is defined by its action

$\displaystyle[y]_{D}$ $\displaystyle=$ $\displaystyle
E_{D}(F)([x]_{D})=\sum_{d=0}^{D-1}y_{d}T^{d}=\left.\sum_{d=0}^{D-1}\frac{1}{d!}\frac{{\rm
d}^{d}}{{\rm d}t^{d}}F(\sum_{d=0}^{D-1}x_{d}t^{d})\right|_{t=0}T^{d}\;.$ (4.3)

The notation $[x]_{D}\equiv[x]_{:D-1}$ and
$[x]_{d+1:D-1}\equiv[x]_{d+1:}\equiv[x_{d+1},\dots,x_{D-1}]$ will be useful
later on. One can show that this definition is compatible with the usual
polynomial addition and multiplication. Furthermore, any composite function
$F(x)=(H\circ G)(x)=H(G(x))$ satisfies

$\displaystyle E_{D}(F)$ $\displaystyle=$ $\displaystyle E_{D}(H)\circ
E_{D}(G)\;.$ (4.4)

I.e., the extension function $E_{D}$ is a homomorphism which preserves
function composition. An immediate consequence is that it is necessary to find
algorithms only for the very limited set of elementary functions
$\phi\in\\{+,-,*,/,\sin,\cos,\exp,\dots\\}$. Explicitly, one performs the
program transformation $E_{D}(F)=E_{D}(\Phi_{L})\circ\dots\circ
E_{D}(\Phi_{1})([x]_{D})$. We call the action of computing
$[y]_{D}=E_{D}(F)([x]_{D})$, i.e., the resolution of the symbolic dependence
to obtain the numerical value $[y]_{D}$, the _pushforward_ of the function
$E_{D}(F)$.

Many functions are implicitly defined by equations of the type
$0=F(x,y)\in{\mathbb{R}}^{M}\;$, where $x\in{\mathbb{R}}^{N}$ are the inputs
and $y\in{\mathbb{R}}^{M}$ the outputs. The idea is to demand that the
_defining equations of order $D$_

$0{\stackrel{{\scriptstyle D}}{{=}}}E_{D}(F)([x]_{D},[y]_{D})$ (4.5)

should be satisfied. By ${\stackrel{{\scriptstyle D}}{{=}}}$ it is meant that
$[x]{\stackrel{{\scriptstyle D}}{{=}}}[y]$ if $x_{d}=y_{d}$ for
$d=0,\dots,D-1$. This is also often written either as
$[x]=[y]+\mathcal{O}(T^{D})$ or $[x]=[y]\mod T^{D}$. The defining equations
lead directly to an algorithmic approach to compute $[y]_{D}$, the so-called
_Newton-Hensel lifting_. In the literature it is often also just called
Hensel-lifting or Newton’s method [GW08]. Assuming $[y]_{D}$ is already known
and satisfies $0{\stackrel{{\scriptstyle D}}{{=}}}E_{D}(F)([x],[y]_{D})$, one
can lift the computation to a higher degree. Explicitly, one tries to solve
$0{\stackrel{{\scriptstyle D+E}}{{=}}}E_{D+E}(F)([x],[y]_{D+E})$. Splitting
$[y]_{D+E}=[y]_{D}+[\Delta y]_{E}T^{D}$ and performing a first order Taylor
expansion of $F$ about $[y]_{D}$ yields after a short calculation

$\displaystyle[\Delta y]_{E}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle-[F_{y}]_{E}^{-1}[\Delta F]_{E}\;,$ (4.6)

where $E_{D+E}(F)([x],[y]_{D}){\stackrel{{\scriptstyle D+E}}{{=}}}[\Delta
F]_{E}T^{D}$ and $[F_{y}]_{E}:=E_{E}(\frac{{\rm d}F}{{\rm d}y})([x],[y]_{E})$.
Setting $E=D$ means that at each step the number of correct coefficients is
doubled. In this case we call it _Newton’s method_. In the case $E=1$ only the
next coefficient is computed. We call the special case $E=1$ _sequential
Hensel lifting_ which is also the formula that is often given as part of the
implicit function theorem. The difference is that Newton-Hensel lifting is a
purely algebraic task. For a discussion on how to obtain asymptotically fast
algorithms and for the nomenclature see e.g. Bernstein [Ber01, Ber08].

###  4.2 The reverse mode

The basic idea of the reverse mode of AD is to pullback linear forms $\alpha$
to obtain an explict mapping $\bar{y}\mapsto\bar{x}$. I.e., given
$F:{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}^{M},\;y=F(x)$ one has

$\displaystyle\alpha(\bar{y},y)$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{M}\bar{y}_{m}{\rm
d}y_{m}=\sum_{m=1}^{M}\bar{y}_{m}\sum_{n=1}^{N}\frac{\partial F_{m}}{\partial
x_{n}}{\rm d}x_{n}=\sum_{n=1}^{N}\bar{x}_{n}{\rm d}x_{n}=\alpha(\bar{x},x)\;,$
(4.7)

where $\bar{x}_{n}=\sum_{m=1}^{M}\bar{y}_{m}\frac{\partial F_{m}}{\partial
x_{n}}$. For notational reasons one uses $\sum_{n=1}^{N}\bar{x}_{n}{\rm
d}x_{n}\equiv\bar{x}^{T}{\rm d}x$. We call the action of going back one level
of the symbolic dependence the _pullback_ of the linear form
$\alpha(\bar{y},y)$. For a more detailed discussion on calculations with
differentials see Magnus and Neudecker [MN99].

It is also possible to compute higher-order derivatives by combining UTP
arithmetic and the reverse mode of AD. For that, the UTP algorithms are
interpreted as functions mapping $D$ coefficients $x_{d}$, $0\leq d<D$ to $D$
coefficients $y_{d}$, $0\leq d<D$, i.e., a mapping from ${\mathbb{R}}^{N\times
D}\rightarrow{\mathbb{R}}^{M\times D}$ with a special structure. One can
formally define a linear form by

$E_{D}(\alpha)([\bar{y}]_{D},[y]_{D}):=[y]_{D}^{T}{\rm d}[y]_{D}\;.$ (4.8)

Here, ${\rm d}[y]_{D}=\sum_{d=0}^{D-1}{\rm d}y_{d}T^{d}$ is a formal
polynomial where each coefficient is a differential and $[\bar{y}]_{D}^{T}{\rm
d}[y]_{D}=\sum_{m=1}^{M}[\bar{y}_{m}]_{D}{\rm d}[y_{m}]_{D}$ computes the
polynomial multiplication of formal polynomials. One can show that

$\displaystyle E_{D}(\alpha)([\bar{y}]_{D},[y]_{D})$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[\bar{y}]_{D}^{T}E_{D}(\frac{\partial F}{\partial
x})([x]_{D}){\rm d}[x]_{D}{\stackrel{{\scriptstyle
D}}{{=}}}[\bar{x}_{n}]_{D}^{T}{\rm
d}[x_{n}]_{D}=E_{D}(\alpha)([\bar{x}]_{D},[x]_{D})$ (4.9)

holds [Chr91]. One can interpret this result as follows: If
$[\bar{y}]_{D}=w\in{\mathbb{R}}^{M}$ then
$[\bar{x}]_{D}=E_{D}(w^{T}\frac{\partial F}{\partial x})([x]_{D})$. Setting
$w=e_{i}$ a Cartesian basis vector would yield the Taylor expansion of the
$i$’th row of the Jacobian. The interpretation of the Taylor coefficients as
derivatives yields higher-order derivatives. If $M=1$ and $w=1$ one obtains
the Taylor expansion of the gradient $[\bar{x}]_{D}=E_{D}(\nabla F)([x]_{D})$.
E.g., propagating the UTP $[x]_{2}=x_{0}+x_{1}T$ would yield
$[\bar{x}]_{2}=\bar{x}_{0}+\bar{x}_{1}T$ where $\bar{x}_{0}=\nabla_{x}F(x)$
and $\bar{x}_{1}=\nabla_{x}^{2}F(x)\cdot x_{1}$, i.e., a Hessian-vector
product.

##  5 Defining equations of numerical linear algebra functions

As briefly mentioned in the introduction, Numerical Linear Algebra (NLA)
functions can be viewed as algorithms representing a concatenation of
functions like $+,-,*,/,\sin,\cos,\dots$ and thus it is possible to apply the
AD techniques described in the previous section directly to the algorithm.
However, there is also the possibility to regard the problem from a more
abstract point of view. Many NLA functions are implicitly defined by a system
of equations.

For instance the QR decomposition is defined by the defining equations

$\displaystyle 0$ $\displaystyle=$ $\displaystyle QR-A$ (5.1) $\displaystyle
0$ $\displaystyle=$ $\displaystyle Q^{T}Q-1\\!\\!\\!\mathrm{I}$ (5.2)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle P_{L}\circ R\;,$ (5.3)

where $A,R\in{\mathbb{R}}^{M\times N}$ with $M\geq N$ and
$Q\in{\mathbb{R}}^{M\times M}$. The functional dependence of the defining
equations is denoted

$\displaystyle Q,R={\rm qr\,}(A)\;.$ (5.4)

Only the first $N$ rows $R_{:N,:}\in{\mathbb{R}}^{N\times N}$ of $R$ are
nonzero. For convenience reasons we use the slicing notation
$i:j=(i,i+1,\dots,j)$.

The defining equations of the symmetric eigenvalue decomposition are given by

$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}AQ-\Lambda$ (5.5)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}Q-1\\!\\!\\!\mathrm{I}$
(5.6) $\displaystyle 0$ $\displaystyle=$
$\displaystyle(P_{L}+P_{R})\circ\Lambda\;,$ (5.7)

where $A\in{\mathbb{R}}^{M\times M}$ is symmetric. The functional dependence
is denoted $\Lambda,Q={\rm eigh\,}(A)$. We call the matrices
$(P_{L})_{ij}=\delta_{j<i}$ and $(P_{R})_{ij}=\delta_{i<j}$ skeletal
projectors since their elementwise product with a matrix returns strictly
lower resp. strictly upper triangular matrices.

##  6 The QR decomposition

Before we derive algorithms based on the defining equations, we briefly
investigate what can go wrong if a typical implementation of the QR
decomposition using Householder reflections is evaluated in UTP arithmetic.
Consider Algorithm 5.1.1 from the book “Matrix Computations” by Golub and Van
Loan [GVL96] which we adapted to our notation in Algorithm 1. From the AD
point of view, the problematic part in the code is the check $\sigma=0$. Since
a paradigm of AD tools is that the control flow must remain unchanged, the
check $\sigma=0$ only considers the zeroth coefficient $x_{0}$ of a UTP.
Hence, if $[x]_{2}=e_{1}+x_{1}T$ is given as input and $x_{1}\neq 0$, the
algorithm will simply evaluate $\beta=0$ and return. As final result, one
obtains a matrix $[R]_{2}$ where $R_{1}$ is not upper triangular. The LAPACK
implementation (LAPACK-3.2.2) of DGEQRFP.f calls the subroutine DLARFGP.f
which contains a similar check. Hence, automatic augmentation based on AD
principles can go wrong in such cases.

As a side remark, note that additionally the function realized by this
algorithm has a pole at $\sigma=0$, producing numerical overflow for
$\sigma\approx 0$.

input : $x\in{\mathbb{R}}^{N}$

output : $v\in{\mathbb{R}}^{N}$ with $v_{1}=1$

output : $\beta\in{\mathbb{R}}$

$\sigma=x_{2:}^{T}x_{2:}$

$v=\begin{pmatrix}1\\\ x_{2:}\end{pmatrix}$

if _$\sigma=0$_ then $\beta=0$ else $\mu=\sqrt{x_{1}^{2}+\sigma}$

if _$x_{1}\leq 0$_ then $v_{1}=x_{1}-\mu$ else $v_{1}=-\sigma/(x_{1}+\mu)$ end
$\beta=2v_{1}^{2}/(\sigma+v_{1}^{2})$

$v=v/v_{1}$ end

Algorithm 1 Householder Vector. The reflector is $P_{v}=I-\beta vv^{\top}$,
with $v_{1}=1$.

###  6.1 Pushforward in Taylor arithmetic

We now come to the higher–level approach that is based on the defining
equations given in Section 5. To compute $[Q]_{D},[R]_{D}=E_{D}({\rm
qr\,})([A]_{D})$ one can apply Newton-Hensel lifting to solve

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}[R]_{D}-[A]_{D}$ (6.1) $\displaystyle 0$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I}$ (6.2) $\displaystyle 0$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$ $\displaystyle
P_{L}\circ[R]_{D}\;.$ (6.3)

The direct application of Eqn. (4.6) should be avoided since $F_{y}$ is sparse
and has a lot of structure. Rather, one assumes that one has already computed
$[Q]_{D}$ and $[R]_{D}$ and computes the next $1\leq E\leq D$ coefficients by
performing a first order Taylor expansion $[Q]_{D+E}=[Q]_{D}+[\Delta
Q]_{E}T^{D}$ and $[R]_{D+E}=[R]_{D}+[\Delta R]_{E}T^{D}$ and tries to solve
for the yet unknown $[\Delta R]_{E}$ and $[\Delta Q]_{E}$. As result one
obtains Proposition 1. For convenience, we use the convention that $R_{d;i,j}$
is the $i$’th row and $j$’th column of the $d$’th coefficient of $[R]_{D}$.

######  Proposition 1.

Let $[A]_{D+E}\in{\mathbb{R}}^{M\times N}[T]/(T^{D+E})$ with $M\geq N$ and
$1\leq E\leq D$, $[R]_{D}\in{\mathbb{R}}^{M\times N}[T]/(T^{D})$ where
$[R_{:N,:}]_{D}$ is upper triangular with nonsingular $R_{0;:N,:}$ and
$[Q]_{D}\in{\mathbb{R}}^{M\times M}[T]/(T^{D})$ orthogonal be given and
satisfy the defining equations of order $D$. Then $[\Delta
R_{:N,:}]_{E}\equiv[R_{:N,:}]_{D:D+E-1}$ and $[\Delta
Q]_{E}\equiv[Q]_{D:D+E-1}$ are given by

$\displaystyle\;[\Delta F]_{E}T^{D}$ $\displaystyle{\stackrel{{\scriptstyle
D+E}}{{=}}}$ $\displaystyle-[Q]_{D}[R]_{D}+[A]_{D+E}$ (6.4)
$\displaystyle\;[\Delta G]_{E}T^{D}$ $\displaystyle{\stackrel{{\scriptstyle
D+E}}{{=}}}$ $\displaystyle-[Q]^{T}_{D}[Q]_{D}+1\\!\\!\\!\mathrm{I}$ (6.5)
$\displaystyle\;[S]_{E}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle\frac{1}{2}[\Delta G]_{E}$ (6.6)
$\displaystyle\;P_{L}\circ([X_{:,:N}]_{E})$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle
P_{L}\circ\left([Q]^{T}_{E}[\Delta
F]_{E}[R_{:N,:}]_{E}^{-1}\right)-P_{L}\circ[S_{:,:N}]_{E}$ (6.7)
$\displaystyle\;[\Delta R]_{E}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle[Q]_{E}^{T}[\Delta F]_{E}-([S]_{E}+[X]_{E})[R]_{E}$
(6.8) $\displaystyle\;[\Delta Q]_{E}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle[Q]_{E}\left([S]_{E}+[X]_{E}\right)\;,$ (6.9)

where $P_{L}\in{\mathbb{R}}^{M\times N}$ with $(P_{L})_{ij}=\delta_{j<i}$.

###### Proof.

In the Appendix A.1.1. ∎

###  6.2 Pullback

######  Proposition 2.

Let $A,R,\bar{R}\in{\mathbb{R}}^{M\times N}$ resp.
$Q,\bar{Q}\in{\mathbb{R}}^{M\times M}$ be given and it holds $M\geq N$,
${\mathrm{rank}\,}(A)=N$, $Q,R={\rm qr\,}(A)$. Then
$\bar{A}\in{\mathbb{R}}^{M\times N}$ can be computed by

$\displaystyle\bar{A}$ $\displaystyle=$ $\displaystyle
Q\left(\bar{R}+\left(P_{L}\circ\left(R\bar{R}^{T}-\bar{R}R^{T}+Q^{T}\bar{Q}-\bar{Q}^{T}Q\right)\right)R^{+T}\right)\;.$
(6.10)

Here, $R^{+}$ denotes the Moore-Penrose pseudoinverse of $R$. That means it
satisfies $RR^{+}R=R$ and since $R$ has full column rank also
$R^{+}R=1\\!\\!\\!\mathrm{I}$.

###### Proof.

In Appendix A.1.2. ∎

###  6.3 Explicit algorithms

One can use Proposition 1 to derive an explicit algorithm as shown in
Algorithm 2, where at each step $E=1$ is used.

input : $[A]_{D}=[A_{0},\dots,A_{D-1}]$, where $A_{d}\in{\mathbb{R}}^{M\times
N}$, $d=0,\dots,D-1$ and ${\mathrm{rank}\,}(A_{0})=N$, $M\geq N$.

output : $[Q]_{D}=[Q_{0},\dots,Q_{D-1}]$ matrix with orthonormal column
vectors, where $Q_{d}\in{\mathbb{R}}^{M\times N}$ , $d=0,\dots,D-1$

output : $[R]_{D}=[R_{0},\dots,R_{D-1}]$ upper triangular, where
$R_{d}\in{\mathbb{R}}^{N\times N}$ , $d=0,\dots,D-1$

$Q_{0},R_{0}={\rm qr\,}(A_{0})$

for _$d=1$ to $D-1$_ do $\Delta F=A_{d}-\sum_{k=1}^{d-1}Q_{d-k}R_{k}$

$S=-\frac{1}{2}\sum_{k=1}^{d-1}Q^{T}_{d-k}Q_{k}$

$X_{:,:N}=P_{L}\circ(Q_{0}^{T}\Delta FR_{0;:N,:N}^{-1}-S_{:,:N})$

$X_{:,N+1:}=0$

$X=X-X^{T}$

$R_{d}=Q_{0}^{T}\Delta F-(S+X)R_{0}$

$Q_{d}=Q_{0}(S+X)$ end

Algorithm 2 Sequential Hensel lifting for the QR decomposition.

The pullback can be computed in Taylor arithmetic. In the global derivative
accumulation it is necessary to update the value of $[\bar{A}]_{D}$. This
happens if $[A]_{D}$ is input of more than one function. The algorithm for the
pullback takes this into consideration.

input : $[A]_{D}=[A_{0},\dots,A_{D-1}]$, where $A_{d}\in{\mathbb{R}}^{M\times
N}$, $d=0,\dots,D-1$, $M\geq N$.

input : $[Q]_{D}=[Q_{0},\dots,Q_{D-1}]$, where $Q_{d}\in{\mathbb{R}}^{M\times
M}$ , $d=0,\dots,D-1$

input : $[R]_{D}=[R_{0},\dots,R_{D-1}]$, where $R_{d}\in{\mathbb{R}}^{M\times
N}$ , $d=0,\dots,D-1$

input/output : $[\bar{A}]_{D}=[\bar{A}_{0},\dots,\bar{A}_{D-1}]$, where
$\bar{A}_{d}\in{\mathbb{R}}^{M\times N}$, $d=0,\dots,D-1$, $M\geq N$.

input : $[\bar{Q}]_{D}=[\bar{Q}_{0},\dots,\bar{Q}_{D-1}]$, where
$\bar{Q}_{d}\in{\mathbb{R}}^{M\times M}$ , $d=0,\dots,D-1$

input : $[\bar{R}]_{D}=[\bar{R}_{0},\dots,\bar{R}_{D-1}]$, where
$\bar{R}_{d}\in{\mathbb{R}}^{M\times N}$ , $d=0,\dots,D-1$

$\displaystyle[\bar{A}]_{D}$ $\displaystyle=$
$\displaystyle[\bar{A}]_{D}+[Q]_{D}\cdot$
$\displaystyle\cdot\left([\bar{R}]_{D}+\left(P_{L}\circ\left([R]_{D}[\bar{R}]_{D}^{T}-[\bar{R}]_{D}[R]_{D}^{T}+[Q]_{D}^{T}[\bar{Q}]_{D}-[\bar{Q}]_{D}^{T}[Q]_{D}\right)\right)[R]_{D}^{+T}\right)$

Algorithm 3 Pullback of the QR decomposition in Taylor arithmetic. The inputs
$[A]_{D},[Q]_{D},[R]_{D}$ must satisfy the defining equations.

##  7 The real symmetric eigenvalue decomposition

The problem of finding eigenvalues and eigenvectors arises in a wide variety
of practical applications. As for the QR decomposition, we want to have
algorithms that compute the real symmetric eigenvalue decomposition in UTP
arithmetic as well as pullback algorithms. The symmetric eigenvalue
decomposition is also important since the Singular Value Decomposition (SVD)
of real matrices is closely related to it. More explicitly, one can compute
the SVD of a matrix $A\in{\mathbb{R}}^{M\times N}$ of rank $r$., i.e.,
$A=U\Sigma V^{T}$, where $\Sigma=\mathop{\rm diag}(\Sigma_{1},0)$,
$U=(U_{1},U_{2})$, $U_{1}\in{\mathbb{R}}^{M\times r}$, $V=(V_{1},V_{2})$,
$V_{1}\in{\mathbb{R}}^{N\times r}$ as

$\displaystyle C=\begin{pmatrix}0&A\\\
A^{T}&0\end{pmatrix}=P^{T}\begin{pmatrix}\Sigma_{1}&0&0\\\ 0&-\Sigma_{1}&0\\\
0&0&0\\\ \end{pmatrix}P\;,$

where

$\displaystyle
P=\frac{1}{\sqrt{2}}\begin{pmatrix}U_{1}&U_{1}&\sqrt{2}U_{2}&0\\\
V_{1}&-V_{1}&0&\sqrt{2}V_{2}\end{pmatrix}^{T}$

is orthogonal [Bjö96].

###  7.1 Pushforward in Taylor arithmetic

Given the symmetric polynomial matrix $[A]_{D}\in{\mathbb{R}}^{N\times
N}[T]/(T^{D})$. The eigenvalue decomposition is the solution
$[\Lambda]_{D},[Q]_{D}\in{\mathbb{R}}^{N\times N}[T]/(T^{D})$ of the implicit
system

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[A]_{D}[Q]_{D}-[\Lambda]_{D}$ (7.1) $\displaystyle 0$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I}$ (7.2) $\displaystyle 0$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle(P_{L}+P_{R})\circ[\Lambda]_{D}\;,$ (7.3)

which is called the _defining equations of order $D$_. We also assume that the
eigenvalues are sorted as
$[\Lambda_{11}]_{D}\leq[\Lambda_{22}]_{D}\leq\dots\leq[\Lambda_{NN}]_{D}$. The
functional dependence is denoted

$\displaystyle[\Lambda]_{D},[Q]_{D}$ $\displaystyle={\rm eigh\,}([A]_{D})\;.$
(7.4)

Let $\Lambda,Q={\rm eigh\,}(A)$ be the usual symmetric eigenvalue
decomposition. We denote the diagonal of $[\Lambda]_{D}$ as
$[\lambda]_{D}=\mathop{\rm diag}([\Lambda]_{D})$. If eigenvalues are repeated,
i.e., multiple, the eigenvectors generalize to eigenspaces and the columns of
$Q$, that are associated to such a multiple eigenvalue, are not unique.
Rather, any orthonormal basis could be the result. This has consequences for
the Hensel-Newton lifting approach, because given $[Q]_{D}$ and $[R]_{D}$ that
satisfy the defining equations of order $D$ it is generally not possible to
find a $[\Delta Q]_{E}$ and $[\Delta R]_{E}$ such that
$[Q]_{D+E}=[Q]_{D}+[\Delta Q]_{E}T^{D}$ and $[R]_{D+E}=[R]_{D}+[\Delta
R]_{E}T^{D}$ satisfy the defining equations of order $D+E$. The higher-order
coefficients $[\Delta A]_{E}$ enforce additional conditions on the chosen
basis of the eigenspaces. A wrong choice of $[Q]_{D}$ means that
$0{\stackrel{{\scriptstyle D+E}}{{=}}}(P_{L}+P_{R})\circ[\Lambda]_{D+E}$
cannot be satisfied. However, $0{\stackrel{{\scriptstyle
D}}{{=}}}P_{b}^{D}\circ[\Lambda]_{D+E}$ can be satisfied. The matrix
$P_{b}^{D}$ is a skeletal projector with zero blocks on the main diagonal
whose size corresponds to the multiplicity of an eigenvalue $[\lambda]_{D}$
and all other entries are ones. The _multiplicity_ $m^{d}([\lambda_{j}]_{D})$
of an eigenvalue $[\lambda_{j}]_{D}$ of _level_ $d$ is defined to be the
number of $i\in{\mathbb{N}}$ s.t. $[\lambda_{j}]_{D}\stackrel{{\scriptstyle
d}}{{=}}[\lambda_{i}]_{D}$. I.e.,

$\displaystyle\mathop{\rm diag}([\Lambda]_{d})$ $\displaystyle=$
$\displaystyle(\underbrace{[\lambda_{1}]_{d},\dots,[\lambda_{1}]_{d}}_{m^{d}([\lambda_{1}]_{D})\mbox{
times}},\dots,\underbrace{[\lambda_{N_{\mathrm{b}}^{d}}]_{d},\dots,[\lambda_{N_{\mathrm{b}}^{d}}]_{d}}_{m^{d}([\lambda_{N_{\mathrm{b}}^{d}}]_{D})\mbox{
times}}),$

where $N_{\mathrm{b}}^{d}$ is the number of different eigenvalues at level
$d$. We define $b^{d}\in{\mathbb{N}}^{N_{\mathrm{b}}^{d}+1}$ to be a vector
satisfying
$m^{d}([\lambda_{{n_{\mathrm{b}}}}]_{D})=b_{{n_{\mathrm{b}}}+1}^{d}-b_{{n_{\mathrm{b}}}}^{d}$.
The symbol $b$ is used because it relates to blocks in the matrix. The
elements of $P_{b}^{d}$ satisfy
$(P^{d}_{b})_{ij}=1-\sum_{{n_{\mathrm{b}}}=1}^{N_{\mathrm{b}}^{d}+1}\delta_{b_{n_{\mathrm{b}}}^{d}\leq
i<b_{{n_{\mathrm{b}}}+1}^{d}}\delta_{b_{n_{\mathrm{b}}}^{d}\leq
j<b_{{n_{\mathrm{b}}}+1}^{d}}$. This notation is a little cumbersome but turns
out to be helpful. One defines $b^{0}=[0,N+1]$. The vector $b^{1}$ represents
the multiplicities of the usual symmetric eigenvalue decomposition. E.g., for
$N=3$ and $b^{d}=[1,3,4]$ one has

$\displaystyle P_{b}^{d}=\begin{pmatrix}0&0&1\\\ 0&0&1\\\
1&1&0\end{pmatrix}\;.$

We reformulate the overall problem as a sequence of subproblems. We call the
implicit system

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q^{d}]_{D}^{T}[A]_{D}[Q^{d}]_{D}-[\Lambda^{d}]_{D}$ (7.5)
$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q^{d}]_{D}^{T}[Q^{d}]_{D}-1\\!\\!\\!\mathrm{I}$ (7.6)
$\displaystyle 0$ $\displaystyle\stackrel{{\scriptstyle d}}{{=}}$
$\displaystyle(P_{L}+P_{R})\circ[\Lambda]_{d}$ (7.7) $\displaystyle 0$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$ $\displaystyle
P_{b}^{d}\circ[\Lambda^{d}]_{D}\;,$ (7.8)

the _relaxed problem of level $d$ and order $D$_. I.e., it is assumed that up
to order $d$ the original problem is solved but only block diagonalized for
the higher order coefficients.

To give an illustrative example consider this relaxed problem of order $3$ and
level $2$. At this point of the algorithm, one has potentially obtained a
matrix polynomial $[\Lambda]_{3}=\sum_{d=0}\Lambda_{d}T^{d}$ with coefficients
of the form

$\displaystyle\Lambda_{0}=\begin{pmatrix}1\\\ &1\\\ &&1\\\ &&&2\\\ &&&&2\\\
&&&&&3\end{pmatrix},\quad\Lambda_{1}=\begin{pmatrix}2\\\ &2\\\ &&3\\\ &&&2\\\
&&&&2\\\ &&&&&2\end{pmatrix},\quad\Lambda_{2}=\begin{pmatrix}1&3\\\ 3&5\\\
&&7\\\ &&&1&2\\\ &&&2&3\\\ &&&&&7\end{pmatrix}\;.$

I.e., $\Lambda_{0}$ and $\Lambda_{1}$ are already diagonal. Since there are
two eigenvalues with multiplicity $m^{2}([\lambda]_{3})=2$ it follows that
$\Lambda_{2}$ is only block diagonal. Note that the eigenvalues are not
globally sorted by value in the higher coefficients but only in the subblocks
defined by the lower order coefficients. In this example, the repeated
eigenvalues in the first block split at the lift from $d=0$ to $d=1$. The
blocks are defined by $b^{1}=[1,4,6,7]$ and $b^{2}=[1,3,4,6,7]$. The blocks in
$\Lambda_{2}$ are defined by $b^{2}$.

The function that solves the relaxed problem of order $D$ and level $d$ is
denoted

$\displaystyle[\Lambda^{d}]_{D},[Q^{d}]_{D}$ $\displaystyle=$
$\displaystyle{\rm eigh\,}_{d}([A]_{D})\;.$ (7.9)

The idea is to implement an algorithm that successively increases $d$ by one.
For convenience we define $[Q^{0}]_{D}:=1\\!\\!\\!\mathrm{I}$ and
$[\Lambda^{0}]_{D}:=[A]_{D}$.

######  Theorem 3.

Let $[A]_{D}$ be given. Then the solution of

$\displaystyle[Q^{d+1}]_{D},[\Lambda^{d+1}]_{D}{\stackrel{{\scriptstyle
D}}{{=}}}{\rm eigh\,}_{d+1}([A]_{D})$ (7.10)

can be computed from the solution
$[Q^{d}]_{D},[\Lambda^{d}]_{D}{\stackrel{{\scriptstyle D}}{{=}}}{\rm
eigh\,}_{d}([A]_{D})$ by computing

$\displaystyle[\hat{\Lambda}_{s,s}]_{D-d},[\hat{Q}_{s,s}]_{D-d}$
$\displaystyle\stackrel{{\scriptstyle D-d}}{{=}}$ $\displaystyle{\rm
eigh\,}_{1}([\Lambda^{d}_{s,s}]_{d:})\;,$ (7.11)

where $s={b^{d}_{{n_{\mathrm{b}}}}:b^{d}_{{n_{\mathrm{b}}}+1}-1}$ are slice
indices and ${n_{\mathrm{b}}}=1,\dots,N_{\mathrm{b}}^{d}$. All other elements
of $[\hat{Q}]_{D-d}$ and $[\hat{\Lambda}]_{D-d}$ are zero. I.e.,
$[\hat{Q}]_{D-d}$ and $[\hat{\Lambda}]_{D-d}$ are block diagonal. It holds
that

$\displaystyle\;[\Lambda^{d+1}]_{D}$ $\displaystyle{\stackrel{{\scriptstyle
D}}{{=}}}$ $\displaystyle[\Lambda^{d}]_{d}+[\hat{\Lambda}]_{D-d}T^{d}$ (7.12)
$\displaystyle\;[Q^{d+1}]_{D}$ $\displaystyle{\stackrel{{\scriptstyle
D}}{{=}}}$ $\displaystyle[Q^{d}]_{D}[Q]_{D}\;,$ (7.13)

where $[Q]_{D}=[\hat{Q}]_{D-d}+[\Delta Q]_{d}T^{D-d}$ for some $[\Delta
Q]_{D-d}$ that satisfies

$\displaystyle 0{\stackrel{{\scriptstyle
D}}{{=}}}[Q]^{T}_{D}[Q]_{D}-1\\!\\!\\!\mathrm{I}\;.$ (7.14)

###### Proof.

We need to show that $[\Lambda^{d+1}]_{D}$, $[Q^{d+1}]_{D}$ is a solution to
the relaxed equations of level $d+1$ and order $D$. From the definition of
${\rm eigh\,}_{1}$ it follows that $0=(P_{L}+P_{R})\circ[\Lambda^{d+1}]_{d+1}$
and $0=P_{b}^{d+1}\circ[\Lambda^{d+1}]_{D}$ is satisfied. We also know that
$0{\stackrel{{\scriptstyle
D}}{{=}}}[Q^{d+1}]_{D}^{T}[Q^{d+1}]_{D}-1\\!\\!\\!\mathrm{I}{\stackrel{{\scriptstyle
D}}{{=}}}[Q]_{D}^{T}[Q^{d}]_{D}^{T}[Q^{d}]_{D}[Q]_{D}-1\\!\\!\\!\mathrm{I}$ is
fulfilled because $0{\stackrel{{\scriptstyle
D}}{{=}}}[Q^{d}]_{D}^{T}[Q^{d}]_{D}-1\\!\\!\\!\mathrm{I}$ and
$0{\stackrel{{\scriptstyle
D}}{{=}}}[Q]_{D}^{T}[Q]_{D}^{T}-1\\!\\!\\!\mathrm{I}$. Hence, it only remains
to show that the third defining equation is satisfied which is shown by the
following straight-forward calculation:

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[Q^{d}]_{D}^{T}[A]_{D}[Q^{d}]_{D}[Q]_{D}-[\Lambda^{d+1}]_{D}$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[\Lambda^{d}]_{D}[Q]_{D}-[\Lambda^{d+1}]_{D}$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}([\Lambda^{d}]_{d}+[\Lambda^{d}]_{d:}T^{d})[Q]_{D}-[\Lambda^{d}]_{d}-[\hat{\Lambda}]_{D-d}T^{d}$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[Q]_{D}^{T}[\Lambda^{d}]_{d}[Q]_{D}+[Q]_{D}^{T}[\Lambda^{d}]_{d:}[Q]_{D}T^{d}-[\Lambda^{d}]_{d}-[\hat{\Lambda}]_{D-d}T^{d}$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[\Lambda^{d}]_{d}[Q]_{D}^{T}[Q]_{D}+[\hat{Q}]_{D-d}^{T}[\Lambda^{d}]_{d:}[\hat{Q}]_{D-d}T^{d}-[\Lambda^{d}]_{d}-[\hat{\Lambda}]_{D-d}T^{d}$
$\displaystyle{\stackrel{{\scriptstyle D}}{{=}}}$
$\displaystyle[\hat{Q}]_{D-d}^{T}[\Lambda^{d}]_{d:}[\hat{Q}]_{D-d}T^{d}-[\hat{\Lambda}]_{D-d}T^{d}$
$\displaystyle\stackrel{{\scriptstyle D-d}}{{=}}$
$\displaystyle[\hat{Q}]_{D-d}^{T}[\Lambda^{d}]_{d:}[\hat{Q}]_{D-d}-[\hat{\Lambda}]_{D-d}\;.$

In the fifth line it has been used that the diagonalization has only to be
performed for block diagonal matrices. If the eigenvalues are already distinct
there is nothing to diagonalize and the step can be skipped. It also means
that one may interchange $[\Lambda^{d}]_{d}$ with $[Q]_{D}$. ∎

The following proposition gives us the means to diagonalize a matrix in the
zeroth degree and block diagonalize w.r.t. the blocks defined by the repeated
eigenvalues. I.e., it gives the justification that the solution of Eqn. (7.11)
can be found. In the case of distinct eigenvalues the application of this
algorithm already solves the original problem.

######  Proposition 4.

Let $[A]_{D+E}=[A]_{D}+[\Delta A]_{E}T^{D}$ be given and $[\Lambda^{d}]_{D}$,
$[Q^{d}]_{D}$ be a solution of the relaxed problem of level $d=1$ and order
$D$. Then it exist $[\Delta\Lambda^{d}]_{E}$ and $[\Delta Q^{d}]_{E}$ such
that $[\Lambda^{d}]_{D+E}=[\Lambda^{d}]_{D}+[\Delta\Lambda^{d}]_{E}T^{D}$ and
$[\Delta Q^{d}]_{D+E}=[\Delta Q^{d}]_{D}+[\Delta Q^{d}]_{E}T^{D}$ are a
solution of the relaxed problem of level $d=1$ and order $D+E$. A closed form
solution is

$\displaystyle\;[\Delta\Lambda^{d}]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle\bar{P}_{b}^{d}\circ[K]_{E}$ (7.15) $\displaystyle\;[\Delta
Q^{d}]_{E}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle[Q^{d}]_{E}\left([\Delta
G]_{E}+P_{b}^{d}\circ\left([K]_{E}/[E]_{E}\right)\right)$ (7.16)

where $[\Delta F]_{E}T^{D}{\stackrel{{\scriptstyle
D+E}}{{=}}}[Q^{d}]^{T}_{D}[A]_{D}[Q^{d}]_{D}-[\Lambda^{d}]_{D}$ and $[\Delta
G]_{E}T^{D}{\stackrel{{\scriptstyle
D+E}}{{=}}}-\frac{1}{2}\left([Q^{d}]^{T}_{D}[Q^{d}]_{D}-1\\!\\!\\!\mathrm{I}\right)$,
$[K]_{E}{\stackrel{{\scriptstyle E}}{{=}}}[\Delta F]_{E}+([\Lambda]_{E}[\Delta
G]_{E}+[\Delta G]_{E}[\Lambda]_{E})+[Q^{d}]^{T}_{E}[\Delta A]_{E}[Q^{d}]_{E}$
and $[E_{ij}]_{E}{\stackrel{{\scriptstyle
E}}{{=}}}[\Lambda^{d}_{jj}]_{E}-[\Lambda^{d}_{ii}]_{E}$. The expression
$[K]_{E}/[E]_{E}$ denotes an element-wise division. $P_{b}^{d}$ is a matrix
with only ones on the diagonal blocks defined by the multiplicity of
eigenvalues in $\Lambda_{0}$. $\bar{P}_{b}^{d}$ is defined s.t.
$\bar{P}_{b}^{d}+P_{b}^{d}$ is a matrix full of ones. One can see here that if
the eigenvalues are distinct, then $\bar{P}_{b}^{d}$ is the identity matrix
$1\\!\\!\\!\mathrm{I}$.

###### Proof.

We set $Q^{d}\equiv Q$ etc. for notational simplicity. Applying Newton-Hensel
lifting to the defining equations yields

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle([Q]_{D}+[\Delta Q]_{E}T^{D})^{T}([Q]_{D}+[\Delta
Q]_{E}T^{D})-1\\!\\!\\!\mathrm{I}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle-2[\Delta G]_{E}+[\Delta
Q]_{E}^{T}[Q]_{E}+[Q]_{E}^{T}[\Delta Q]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle-2[\Delta
G]_{E}+2[S]_{E}\;,$ $\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle
D+E}}{{=}}}$ $\displaystyle([Q]_{D}+[\Delta Q]_{E}T^{D})^{T}([A]_{D}+[\Delta
A]_{E}T^{D})([Q]_{D}+[\Delta
Q]_{E}T^{D})-([\Lambda]_{D}+[\Delta\Lambda]_{E}T^{D})$ (7.17)
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle[\Delta
F]_{E}+[Q]_{E}^{T}[\Delta A]_{E}[Q]_{E}+[\Delta
Q]_{E}^{T}[Q]_{E}[\Lambda]_{E}+[\Lambda]_{E}[Q]_{E}^{T}[\Delta
Q]_{E}-[\Delta\Lambda]_{E}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle[K]_{E}+[X]_{E}[\Lambda]_{E}-[\Lambda]_{E}[X]_{E}-[\Delta\Lambda]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle[K]_{E}+[E]_{E}\circ[X]_{E}-[\Delta\Lambda]_{E}\;.$

Thus $[\Delta\Lambda]_{E}{\stackrel{{\scriptstyle
E}}{{=}}}\bar{P}_{b}^{d}\circ[K]_{E}$ and $[X]_{E}^{T}{\stackrel{{\scriptstyle
E}}{{=}}}P_{b}^{d}\circ([K]_{E}/[E]_{E})$. Above, $[\Delta
Q]_{E}^{T}[Q]_{E}{\stackrel{{\scriptstyle E}}{{=}}}[S]_{E}+[X]_{E}$, $[S]_{E}$
symmetric and $[X]_{E}$ antisymmetric (Lemma 15) has been used. ∎

It remains to show that Eqn. (7.14) can be satisfied.

######  Lemma 5.

Let $[Q]_{D}$ be given and it satisfies the defining equation
$0{\stackrel{{\scriptstyle D}}{{=}}}[Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I}$.
Then the solution can be lifted to $D+E$ with $E\leq D$. I.e., it is possible
to find $[Q]_{D+E}:=[Q]_{D}+[\Delta Q]_{E}T^{D}$ s.t.
$0{\stackrel{{\scriptstyle
D+E}}{{=}}}[Q]_{D+E}^{T}[Q]_{D+E}-1\\!\\!\\!\mathrm{I}$. A closed form
solution for $[\Delta Q]_{E}$ is given by

$\displaystyle\;[\Delta Q]_{E}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle[Q]_{E}[S]_{E}\;,$ (7.18)

where $[S]_{E}T^{D}\stackrel{{\scriptstyle
D+E}}{{=}}-\frac{1}{2}\left([Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I}\right)$.

###### Proof.

In Appendix A.1.3. ∎

###  7.2 Pullback

The eigenvalue decomposition is non-differentiable at points where eigenvalues
are repeated and hence the defining equations do not define a _well behaved
implicit mapping_ as described by Christianson [Chr98]. However, the
eigenvalue decomposition is typically used within a global context where the
non-uniqueness and non-differentiability can be worked around. Here, we give
only the pullback algorithm that is correct for unique eigenvalues.

######  Proposition 6 (Pullback of the Symmetric Eigenvalue Decomposition with
Distinct Eigenvalues:).

Given $A,Q,\Lambda,\bar{Q},\bar{\Lambda}$, where all eigenvalues are distinct,
one can compute $\bar{A}$ by

$\displaystyle\;H_{ij}$ $\displaystyle=$
$\displaystyle(\lambda_{j}-\lambda_{i})^{-1}\quad\mbox{if}\quad i\neq j,\quad
0\quad\mbox{ else}$ (7.19) $\displaystyle\bar{A}$ $\displaystyle=$
$\displaystyle Q\left(\bar{\Lambda}+H\circ(Q^{T}\bar{Q})\right)Q^{T}$ (7.20)

###### Proof.

In Appendix A.1.4. ∎

###  7.3 Explicit algorithms

input : $[Q]_{d}=[Q_{0},\dots,Q_{d-1}]$ with $0\stackrel{{\scriptstyle
d}}{{=}}[Q]_{d}^{T}[Q]_{d}-1\\!\\!\\!\mathrm{I}$

input : $D\in{\mathbb{N}}$

output : $[Q]_{D}=[Q_{0},\dots,Q_{D-1}]$, where $0\stackrel{{\scriptstyle
D}}{{=}}[Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I}$

for _$k=d$ to $D-1$_ do
$Q_{k}=-\frac{1}{2}Q_{0}\sum_{i=1}^{k-1}Q_{i}^{T}Q_{k-i}$ end

Algorithm 4 This algorithm computes $[Q]_{D}=\mbox{qlift}([Q]_{d},D)$ as
described in Proposition 5 using sequential Hensel-lifting ($E=1$).

input : $[A]_{D}=[A_{0},\dots,A_{D-1}]$, where $A_{d}\in{\mathbb{R}}^{N\times
N}$ symmetric positive definite, $d=0,\dots,D-1$

output : $[\Lambda]_{D}=[\Lambda_{0},\dots,\Lambda_{D-1}]$, where
$\Lambda_{0}\in{\mathbb{R}}^{N\times N}$ diagonal and
$\Lambda_{d}\in{\mathbb{R}}^{N\times N}$ block diagonal $d=1,\dots,D-1$.

output : $[Q]_{D}=[Q_{0},\dots,Q_{D-1}]$ orthogonal, where
$Q_{d}\in{\mathbb{R}}^{N\times N}$

output : $b\in{\mathbb{N}}^{N_{b}+1}$, array of integers defining the blocks.
The integer $N_{b}$ is the number of blocks. Each block has the size of the
multiplicity of an eigenvalue $\lambda_{n_{b}}$ of $\Lambda_{0}$ s.t. for
$s=b_{n_{b}}:b_{n_{b}+1}-1$ one has
$(Q_{0;:,s})^{T}A_{0}Q_{0;:,s}=\lambda_{n_{b}}1\\!\\!\\!\mathrm{I}$.

$\Lambda_{0},Q_{0}={\rm eigh\,}(A_{0})$

compute $b\in{\mathbb{R}}^{N_{b}+1}$

$E_{ij}=\Lambda_{0;jj}-\Lambda_{0;ii}$

$H=P_{b}\circ(1/E)$

for _$d=1$ to $D-1$_ do $\Delta F=\sum_{|i|=d}Q_{i_{1}}^{T}A_{i_{2}}Q_{i_{3}}$

$S=-\frac{1}{2}\sum_{k=1}^{d-1}Q^{T}_{d-k}Q_{k}$

$K=\Delta F+Q_{0}^{T}A_{d}Q_{0}+S\Lambda_{0}+\Lambda_{0}S$

$Q_{d}=Q_{0}(S+H\circ K)$

$\Lambda_{d}=\bar{P}_{b}\circ K$ end

Algorithm 5 This algorithm computes $[\Lambda]_{D},[Q]_{D},b={\rm
eigh\,}_{1}([A]_{D})$ as specified by 4 using sequential Hensel-lifting
($E=1$). I.e., the zeroth coefficient is diagonalized and the higher order
coefficients are block diagonalized. The symbol $i\in{\mathbb{N}}_{0}^{3}$
denotes a multi-index, i.e., the summation $\sum_{|i|=d}$ goes over all
possible $i$ such that $|i|\equiv\sum_{k=1}^{3}i_{k}=d$.

input : $[A]_{D}=[A_{0},\dots,A_{D-1}]$ symmetric with
$A_{d}\in{\mathbb{R}}^{N\times N}$

output : $[\Lambda]_{D}=[\Lambda_{0},\dots,\Lambda_{D-1}]$, where
$\Lambda_{d}\in{\mathbb{R}}^{N\times N}$ diagonal for $d=0,\dots,D-1$.

output : $[Q]_{D}=[Q_{0},\dots,Q_{D-1}]$ orthogonal, where
$Q_{d}\in{\mathbb{R}}^{N\times N}$

$[\Lambda^{0}]_{D}=[A]_{D}$, $[Q^{0}]_{D}=1\\!\\!\\!\mathrm{I}$ and
$b^{0}=[1,N+1]$

for _$d=0$ to $D-1$_ do for _$n_{b}=1$ to $N_{b}^{d}$_ do
$s=b^{d}_{n_{b}}:b^{d}_{n_{b}+1}-1$ (slice index)

$[\hat{\Lambda}_{s,s}]_{D-d},[\hat{Q}_{s,s}]_{D-d},b^{d+1}={\rm
eigh\,}_{1}([\Lambda^{d}_{s,s}]_{d:})$

$[Q_{s,s}]_{D}=\mbox{qlift}([\hat{Q}_{s,s}]_{D-d},D)$ end
$[\Lambda^{d+1}]_{D}=[\Lambda^{d}]_{d}+[\hat{\Lambda}]_{D-d}T^{d}$

$[Q^{d+1}]_{D}=[Q^{d}]_{D}[Q]_{D}$ end

Algorithm 6 This algorithm computes $[\Lambda]_{D},[Q]_{D}={\rm
eigh\,}([A]_{D})$ as described in Theorem 3. The algorithm uses internally
Algorithm 5 and 4.

##  8 Numerical tests and examples

###  8.1 Taylor polynomial arithmetic on real symmetric eigenvalue problem

As an example to test the validity of the pushforward in UTP arithmetic we
consider the following system [AT98]:

$\displaystyle Q(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{3}}\begin{pmatrix}\cos(x(t))&1&\sin(x(t))&-1\\\
-\sin(x(t))&-1&\cos(x(t))&-1\\\ 1&-\sin(x(t))&1&\cos(x(t))\\\
-1&\cos(x(t))&1&\sin(x(t))\\\ \end{pmatrix}$ $\displaystyle\Lambda(t)$
$\displaystyle=$ $\displaystyle\mathop{\rm
diag}(x^{2}-x+\frac{1}{2},4x^{2}-3x,\delta(-\frac{1}{2}x^{3}+2x^{2}-\frac{3}{2}x+1)+(x^{3}+x^{2}-1),3x-1)\;,$

where $x\equiv x(t):=1+t$. The constant $\delta$ is some predefined constant.
In Taylor arithmetic one obtains

$\displaystyle\Lambda_{0}$ $\displaystyle=$ $\displaystyle\mathop{\rm
diag}(1/2,1,1+\delta,2)$ $\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\mathop{\rm diag}(1,5,5+\delta,3)$ $\displaystyle\Lambda_{2}$
$\displaystyle=$ $\displaystyle\mathop{\rm diag}(2,8,8+\delta,0)$
$\displaystyle\Lambda_{3}$ $\displaystyle=$ $\displaystyle\mathop{\rm
diag}(0,0,6-3\delta,0)$ $\displaystyle\Lambda_{d}$ $\displaystyle=$
$\displaystyle\mathop{\rm diag}(0,0,0,0),\quad\forall d\geq 4\;.$

One can see that in the case $\delta=0$ one obtains one repeated eigenvalue
that splits at $d=3$. We apply Algorithm 6 to reconstruct $[\Lambda]_{D}$. The
reconstructed values are denoted $[\tilde{\Lambda}]_{D}$. The numerical
results are shown in Fig. (8.1).

Figure 8.1:  One left side one can see that error between the true
$\lambda_{1}$ and reconstructed $\tilde{\lambda}_{1}$ is close to machine
precision. On the right side one can see that the absolute error
$|\tilde{\lambda}_{2}-\lambda_{2}|$ has a jump at $\delta\approx 10^{-7}$.
This is due to the fact that that the algorithm treats eigenvalues
$|\lambda_{i}-\lambda_{j}|<10^{-7}$ as repeated eigenvalues. One can see that
when $\delta$ approaches $10^{-16}$ the error gets smaller. The eigenvalue
$\lambda_{4}$ shows the same qualitative behavior as $\lambda_{1}$ and
$\lambda_{3}$ the same as $\lambda_{2}$.

###  8.2 Covariance matrix computation

To test the validity of the covariance matrix computation of Eqn. (2.1) we
first check that the first directional derivatives of the covariance matrix
$C$ w.r.t. $J_{1}$ and $J_{2}$ coincide with the results from the complex-step
derivative approximation, abbreviated here for convenience as CSDA. The CSDA
computes directional derivatives of a real-valued function $y=f(x)$ as
$\dot{y}\approx\frac{\Im(f(x+i\epsilon\dot{x}))}{\epsilon}=\frac{f(x+i\epsilon\dot{x})-f(x-i\epsilon\dot{x})}{2i\epsilon}$,
i.e., $\Im$ extracts the imaginary part and $i\equiv\sqrt{-1}$. The number
$\epsilon\in{\mathbb{R}}$ can be made very small. For a detailed discussion
that also shows the relation to AD see Martins et al. [MSA01, MSA03]. Having
verified the first order derivatives by UTP arithmetic we can check if the UTP
arithmetic on Eqn. (2.2) yields the same result. Unfortunately, it is not
possible to use the CSDA in a straight-forward fashion since for complex
matrices the QR decomposition does not yield an orthogonal but a unitary $Q$.
For reproducibility we define $J_{1}$ and $J_{2}$ rather arbitrarily as

$\displaystyle J_{1}(x)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\sin(x_{1})x_{2}&\cos(x_{2})\\\
\exp(x_{1})&x_{1}x_{2}\\\ x_{1}\log(x_{2})&\log(1+\exp(\cos(x_{1})))\\\
x_{2}+x_{1}&x_{1}(x_{2}+\cos(x_{1})\\\
\end{pmatrix}\;,\quad\;J_{2}(x)^{T}=\begin{pmatrix}x_{1}\log(x_{2}+3\sin(x_{1}x_{2}))\\\
x_{2}\exp(\sin(x_{1})+\cos(x_{1}x_{2}))\end{pmatrix}\;.$

The numerical results are shown in Figure 8.2. Note that $x_{1}$ and $x_{2}$
are here elements of the vector $x$ and not coefficients.

Figure 8.2:  This plot shows the maximum absolute differences of the
directional derivatives at $x=t(3,1)^{T}$, where $t\in[0.1,1]$ in direction
$\dot{x}=(5,7)^{T}$. The circles show the difference between the CSDA solution
and the first order UTP solution using Eqn. (2.1). The diamonds show the
difference between the UTP solution of Eqn. (2.1) and Eqn. (2.2). One can see
that the difference is close to machine precision of IEEE 754 float64, which
is approximately $10^{-16}$.

##  9 Summary and outlook

We have shown how computer codes containing real symmetric eigenvalue
decompositions and QR decompositions can be evaluated in univariate Taylor
polynomial arithmetic. Furthermore, the reverse mode of AD has been treated.
Explicit algorithms have been presented that can be used in combination with
existing AD software, e.g. general purpose AD tools like ADOL-C [GJM+99] or
CppAD [Bel10] but also differentiated DAE solvers like SolvIND [AK]. Numerical
tests have been used to check the algorithms.

Other algorithms that contain the the QR decomposition and the real symmetric
eigenvalue decompostion can be differentiated using the shown algorithms. In
particular, we think of the differentiation of the SVD or the computation of
pseudoinverses. We believe that these algorithms, in modified form, may also
be valuable for the eigenvalue optimization problem where eigenvalues are
repeated in the solution point.

## Acknowledgements

The authors wish to thank Bruce Christianson for his comments that helped to
gain a deeper understanding of the matter and also greatly helped to improve
the readability of this work.

This research was partially supported by the Bundesministerium für Bildung und
Forschung (BMBF) within the project NOVOEXP (Numerische Optimierungsverfahren
für die Parameterschätzung und den Entwurf optimaler Experimente unter
Berücksichtigung von Unsicherheiten für die Modellvalidierung
verfahrenstechnischer Prozesse der Chemie und Biotechnologie) (03GRPAL3),
Humboldt Universität zu Berlin.

##  Appendix A Additional proofs

###  A.1 Proofs of QR decomposition

####  A.1.1 Proof of Proposition 1

###### Proof.

We look at the first defining equation and try to separate the known from the
unknown quantities:

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle[Q]_{D+E}[R]_{D+E}-[A]_{D+E}$ (A.1)
$\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle([Q]_{D}+[\Delta Q]_{E}T^{D})([R]_{D}+[\Delta
R]_{E}T^{D})-[A]_{D+E}$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle[Q]_{D}[R]_{D}-[A]_{D+E}+([\Delta Q]_{E}[R]_{D}+[Q]_{D}[\Delta
R]_{E})T^{D}$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle-[\Delta F]_{E}T^{D}+([\Delta Q]_{E}[R]_{D}+[Q]_{D}[\Delta
R]_{E})T^{D}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-[\Delta F]_{E}+[\Delta Q]_{E}[R]_{E}+[Q]_{E}[\Delta R]_{E}\;.$

Similarly for the second defining equation

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle[Q]^{T}_{D+E}[Q]_{D+E}-1\\!\\!\\!\mathrm{I}$
$\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle[Q]^{T}_{D}[Q]_{D}-1\\!\\!\\!\mathrm{I}+([Q]^{T}_{D}[\Delta
Q]_{E}+[\Delta Q]^{T}_{E}[Q]_{D})T^{D}$ $\displaystyle\Rightarrow 0$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle-[\Delta
G]_{E}+[Q]^{T}_{E}[\Delta Q]_{E}+[\Delta Q]^{T}_{E}[Q]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle-[\Delta
G]_{E}+[S]_{E}+[X]_{E}+[S]_{E}-[X]_{E}$ $\displaystyle\Rightarrow\quad S$
$\displaystyle=$ $\displaystyle\frac{1}{2}[\Delta G]_{E}\;,$

where $[S]_{E}+[X]_{E}=[Q]_{E}^{T}[\Delta Q]_{E}$ and it has been used that
every matrix can be written as the sum of a symmetric and an antisymmetric
matrix. Now multiply (A.1) by $[Q]^{T}_{E}$ from the left to obtain

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-[Q]^{T}_{E}[\Delta F]_{E}+[Q]^{T}_{E}[\Delta
Q]_{E}[R]_{E}+[\Delta R]_{E}$ (A.2) $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle-[Q]^{T}_{E}[\Delta
F]_{E}+([S]_{E}+[X]_{E})[R]_{E}+[\Delta R]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-[Q]^{T}_{E}[\Delta F]_{E}+[S]_{E}[R]_{E}+[X]_{E}[R]_{E}+[\Delta
R]_{E}\;.$

Multiplication of $[R_{:N,:}]_{E}^{-1}$ from the right yields

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-[Q]^{T}_{E}[\Delta
F]_{E}[R_{:N,:}]_{E}^{-1}+[S]_{E}[R]_{E}[R_{:N,:}]_{E}^{-1}+[X]_{E}[R]_{E}[R_{:N,:}]_{E}^{-1}+[\Delta
R]_{E}[R_{:N,:}]_{E}^{-1}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-[Q]^{T}_{E}[\Delta
F]_{E}[R_{:N,:}]_{E}^{-1}+[S_{:,:N}]_{E}+[X_{:,:N}]_{E}+[\Delta
R]_{E}[R_{:N,:}]_{E}$ $\displaystyle\Rightarrow P_{L}\circ([X_{:,:N}]_{E})$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$ $\displaystyle
P_{L}\circ\left([Q]_{E}^{T}[\Delta
F]_{E}[R_{:N,:}]_{E}^{-1}-[S_{:,:N}]_{E}\right)\;.$

The coefficients of $X_{:,N+1:}$ are not specified and can for instance be set
to zero. Since $X$ is antisymmetric it is already defined by the above
equation. Since $[S]_{E}+[X]_{E}{\stackrel{{\scriptstyle
E}}{{=}}}[Q]_{E}^{T}[\Delta Q]_{E}$ one can obtain $[\Delta Q]_{E}$ as

$\displaystyle[\Delta Q]_{E}$ $\displaystyle=$
$\displaystyle[Q]_{E}([S]_{E}+[X]_{E})$

because for quadratic $Q$ one has the identity $QQ^{T}=1\\!\\!\\!\mathrm{I}$.
∎

####  A.1.2 Proof of Proposition 2

###### Proof.

We differentiate the implicit system

$\displaystyle 0$ $\displaystyle=$ $\displaystyle A-QR$ $\displaystyle 0$
$\displaystyle=$ $\displaystyle Q^{T}Q-1\\!\\!\\!\mathrm{I}$ $\displaystyle 0$
$\displaystyle=$ $\displaystyle P_{L}\circ R$

and obtain

$\displaystyle 0$ $\displaystyle=$ $\displaystyle{\rm d}A-{\rm d}QR-Q{\rm
d}R\quad(*)$ $\displaystyle 0$ $\displaystyle=$ $\displaystyle{\rm
d}Q^{T}Q+Q^{T}{\rm d}Q\quad(**)\;.$

We define the antisymmetric “matrix” $X:=Q^{T}{\rm d}Q$. Multiplication of
Eqn. (*) from the left with $Q^{T}$ yields

$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}{\rm d}A-XR-{\rm d}R$
$\displaystyle\mbox{hence}\quad{\rm d}R$ $\displaystyle=$ $\displaystyle
Q^{T}{\rm d}A-XR\;.$

The multipication of this last equation from the right with the Moore-Penrose
pseudoinverse $R^{+}=(R_{:N,:}{}^{-1},0)$ yields the equivalent equation

$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}{\rm
d}AR^{+}-XRR^{+}-{\rm d}RR^{+}$ $\displaystyle\mbox{ and thus }\;P_{L}\circ X$
$\displaystyle=$ $\displaystyle P_{L}\circ(Q^{T}{\rm d}AR^{+})\;,$

where we have chosen arbitrarily that $X_{:,N+1:}=0$. Since $X$ is
antisymmetric we have

$\displaystyle X$ $\displaystyle=$ $\displaystyle(P_{L}\circ X)-(P_{L}\circ
X)^{T}\;.$

We can use these results to compute the pullback:

$\displaystyle\mathop{\rm tr\,}(\bar{R}^{T}{\rm d}R)+\mathop{\rm
tr\,}(\bar{Q}^{T}{\rm d}Q)$ $\displaystyle=\mathop{\rm tr\,}(Q\bar{R}{\rm
d}A^{T})-\mathop{\rm tr\,}(R\bar{R}^{T}X)+\mathop{\rm
tr\,}(\bar{Q}^{T}QQ^{T}{\rm d}Q)$ $\displaystyle=\mathop{\rm
tr\,}(Q\bar{R}{\rm d}A^{T})+\mathop{\rm
tr\,}(\underbrace{(\bar{Q}^{T}Q-R\bar{R}^{T})}_{=:K}X)$
$\displaystyle=\mathop{\rm tr\,}(Q\bar{R}{\rm d}A^{T})+\mathop{\rm
tr\,}((K-K^{T})(P_{L}\circ X))$ $\displaystyle=\mathop{\rm tr\,}(Q\bar{R}{\rm
d}A^{T})+\mathop{\rm tr\,}(R^{+T}{\rm d}A^{T}Q(P_{L}\circ(K^{T}-K)))$
$\displaystyle=\mathop{\rm
tr\,}(Q[\bar{R}+\\{P_{L}\circ(Q^{T}\bar{Q}-\bar{Q}^{T}Q+R\bar{R}^{T}-\bar{R}R^{T})\\}R^{+T}]{\rm
d}A^{T})$ $\displaystyle=\mathop{\rm tr\,}(\bar{A}{\rm d}A^{T})\;.$

In the above derivation we have used Lemmas 7, 8 and 9.

∎

####  A.1.3 Proof of Lemma 5

###### Proof.

$\displaystyle 0$ $\displaystyle{\stackrel{{\scriptstyle D+E}}{{=}}}$
$\displaystyle([Q]_{D}+[\Delta Q]_{E}T^{D})^{T}([Q]_{D}+[\Delta
Q]_{E}T^{D})-1\\!\\!\\!\mathrm{I}$ $\displaystyle{\stackrel{{\scriptstyle
D+E}}{{=}}}$
$\displaystyle([Q]_{D}^{T}[Q]_{D}-1\\!\\!\\!\mathrm{I})+([Q]_{D}^{T}[\Delta
Q]_{E}+[\Delta Q]_{E}^{T}[Q]_{D})T^{D}$ $\displaystyle{\stackrel{{\scriptstyle
E}}{{=}}}$ $\displaystyle[\Delta G]_{E}+[Q]_{E}^{T}[\Delta Q]_{E}+[\Delta
Q]_{E}^{T}[Q]_{E}$ $\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle[\Delta G]_{E}+2[S]_{E}$ $\displaystyle\;[\Delta Q]_{E}$
$\displaystyle{\stackrel{{\scriptstyle E}}{{=}}}$
$\displaystyle-\frac{1}{2}[Q]_{E}[\Delta G]_{E}\;,$

where $[\Delta Q]_{E}^{T}[Q]_{E}=[S]_{E}+[X]_{E}$, $[S]_{E}$ symmetric and
$[X]_{E}$ antisymmetric and $[\Delta G_{E}]_{E}T^{D}\stackrel{{\scriptstyle
D+E}}{{=}}\left(Q^{T}Q-1\\!\\!\\!\mathrm{I}\right)$ . Since no condition
defines constraints on $[X]_{E}$ it has been set to zero. ∎

####  A.1.4 Proof of Proposition 6

###### Proof.

We want to compute $\mathop{\rm tr\,}(\bar{A}^{T}{\rm d}A)=\mathop{\rm
tr\,}(\bar{\Lambda}^{T}{\rm d}\Lambda)+\mathop{\rm tr\,}(\bar{Q}^{T}{\rm
d}Q)$. We differentiate the implicit system

$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}AQ-\Lambda$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle Q^{T}Q-1\\!\\!\\!\mathrm{I}$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle(P_{L}+P_{R})\circ\Lambda$

and obtain

$\displaystyle{\rm d}\Lambda$ $\displaystyle=$ $\displaystyle Q^{T}{\rm
d}AQ+{\rm d}Q^{T}AQ+Q^{T}A{\rm d}Q$ $\displaystyle=$ $\displaystyle Q^{T}{\rm
d}AQ+{\rm d}Q^{T}Q\Lambda+\Lambda Q^{T}{\rm d}Q$ $\displaystyle 0$
$\displaystyle=$ $\displaystyle{\rm d}Q^{T}Q+Q^{T}{\rm d}Q\;.$

A straight forward calculation shows:

$\displaystyle\mathop{\rm tr\,}(\bar{\Lambda}^{T}{\rm d}\Lambda)$
$\displaystyle=$ $\displaystyle\mathop{\rm tr\,}(Q\bar{\Lambda}Q^{T}{\rm
d}A)+\mathop{\rm tr\,}(\Lambda\bar{\Lambda}{\rm d}Q^{T}Q)+\mathop{\rm
tr\,}(\bar{\Lambda}\Lambda Q^{T}{\rm d}Q)$ $\displaystyle=$
$\displaystyle\mathop{\rm tr\,}(Q\bar{\Lambda}Q^{T}{\rm d}A)\;,$
$\displaystyle\mathop{\rm tr\,}(\bar{Q}^{T}{\rm d}Q)$ $\displaystyle=$
$\displaystyle\mathop{\rm tr\,}(\bar{Q}^{T}Q(H\circ(Q^{T}{\rm d}AQ)))$
$\displaystyle=$ $\displaystyle\mathop{\rm
tr\,}(Q(H^{T}\circ(\bar{Q}^{T}Q))Q^{T}{\rm d}A)\;,$ $\displaystyle\mathop{\rm
tr\,}(\bar{A}^{T}{\rm d}A)$ $\displaystyle=$ $\displaystyle\mathop{\rm
tr\,}\left((Q(\bar{\Lambda}+H^{T}\circ(\bar{Q}^{T}Q))Q^{T}){\rm d}A\right)$

where we have used

$\displaystyle{\rm d}\Lambda$ $\displaystyle=$ $\displaystyle Q^{T}{\rm
d}AQ-(Q^{T}{\rm d}Q)^{T}\Lambda+\Lambda Q^{T}{\rm d}Q$ $\displaystyle=$
$\displaystyle Q^{T}{\rm d}AQ-K\circ(Q^{T}{\rm d}Q)$ $\displaystyle\implies
Q^{T}{\rm d}Q$ $\displaystyle=$ $\displaystyle H\circ(Q^{T}{\rm d}AQ-{\rm
d}\Lambda)$ $\displaystyle=$ $\displaystyle H\circ(Q^{T}{\rm d}AQ)$

where we have defined $K_{ij}:=\Lambda_{jj}-\Lambda_{ii}$ and
$H_{ij}=(K_{ij})^{-1}$ for $i\neq j$ and $H_{ij}=0$ otherwise and used the
property $X\Lambda-\Lambda X=K\circ X$ for all $X\in{\mathbb{R}}^{N\times N}$
and diagonal $\Lambda\in{\mathbb{R}}^{N\times N}$. ∎

###  A.2 Basic results used in the proofs

######  Lemma 7.

Let $X\in{\mathbb{R}}^{N\times N}$ be an antisymmetric matrix, i.e.,
$X^{T}=-X$ and $P_{L}$ defined as above. We then can write

$\displaystyle X$ $\displaystyle=$ $\displaystyle P_{L}\circ X-(P_{L}\circ
X)^{T}\;.$ (A.4)

###### Proof.

$X=P_{L}\circ X+P_{R}\circ X=P_{L}\circ X+(P_{L}\circ X^{T})^{T}=P_{L}\circ
X-(P_{L}\circ X)^{T}$ ∎

######  Lemma 8.

Let $A\in{\mathbb{R}}^{N\times N}$ and $P_{L}$ resp. $P_{R}$ defined as above.
Then

$\displaystyle(P_{L}\circ A)^{T}$ $\displaystyle=$ $\displaystyle P_{R}\circ
A^{T}\;.$ (A.5)

###### Proof.

$B_{ij}:=(P_{L}\circ A)_{ij}=A_{ij}(i>j)$ and
$B_{ij}^{T}=B_{ji}=A_{ji}(j>i)=A_{ij}^{T}P_{R}=P_{R}\circ A$ ∎

######  Lemma 9.

Let $A,B,C\in{\mathbb{R}}^{M\times N}$. We then have

$\displaystyle\mathop{\rm tr\,}\left(A^{T}(B\circ C)\right)$ $\displaystyle=$
$\displaystyle\mathop{\rm tr\,}\left(C^{T}(B\circ A)\right)$ (A.6)

###### Proof.

$\mathop{\rm tr\,}(A^{T}(B\circ
C))=\sum_{i=1}^{N}\sum_{j=1}^{M}A_{ij}B_{ij}C_{ij}=\mathop{\rm
tr\,}(C^{T}(B\circ A))$ ∎

######  Lemma 10.

Let $A,B$ be lower triangular matrices. Then the following expression holds:

$\displaystyle P_{D}\circ(AB)$ $\displaystyle=$ $\displaystyle(P_{D}\circ
A)(P_{D}\circ B)\;.$ (A.7)

###### Proof.

$AB$ is also lower trinangular and thus
$P_{D}\circ(AB)=diag(a_{ii}b_{ii})=diag(a_{ii})diag(b_{ii})=(P_{D}\circ
A)(P_{D}\circ B)$ ∎

######  Lemma 11.

The formula

$\displaystyle P_{D}\circ(A^{T})$ $\displaystyle=$ $\displaystyle P_{D}\circ
A$ (A.8)

holds for all quadratic matrices $A$.

###### Proof.

$(P_{D}\circ(A^{T}))_{ij}=\delta_{ij}A_{ji}=\delta_{ij}A_{ij}=(P_{d}\circ
A)_{ij}$ ∎

######  Lemma 12.

Let $A$ be a nonsingular quadratic lower triangular matrix. Then the formula

$\displaystyle P_{D}\circ(A^{-1})$ $\displaystyle=$ $\displaystyle(P_{D}\circ
A)^{-1}$ (A.9)

holds.

###### Proof.

Using Lemma 10 we obtain $(P_{D}\circ(A^{-1}))(P_{D}\circ A)=P_{D}\circ
1\\!\\!\\!\mathrm{I}=1\\!\\!\\!\mathrm{I}$. Since the quadratic matrices form
a group, the inverse is unique. Therefore, equality between
$(P_{D}\circ(A^{-1}))=(P_{D}\circ(A))^{-1}$ must hold. ∎

######  Lemma 13.

Let $A\in{\mathbb{R}}^{N\times N}$ be strictly lower triangular and
$B\in{\mathbb{R}}^{N\times N}$ lower triangular. Then their product $C=AB$ is
strictly lower triangular.

###### Proof.

$C$ is lower triangular and the diagonal entries are $C_{ii}=A_{ii}B_{ii}=0$
since $B$ has a zero diagonal. ∎

######  Corollary 14.

Let $A\in{\mathbb{R}}^{N\times N}$ be strictly lower triangular and
$D\in{\mathbb{R}}^{N\times N}$ diagonal. Then their product $C=AD$ is strictly
lower triagonal.

######  Lemma 15.

Every quadratic matrix $A$ can be written as the sum of a symmetric matrix
$S=\frac{1}{2}(A+A^{T})$ and an antisymmetric matrix $X=\frac{1}{2}(A-A^{T})$,
i.e.

$\displaystyle A$ $\displaystyle=$ $\displaystyle S+X$ (A.10)

###### Proof.

$A=\frac{1}{2}(A+A^{T}+A-A^{T})=\frac{1}{2}(A+A^{T})+\frac{1}{2}(A-A^{T})=S+X$.
∎

## References

  * [AK] Jan Albersmeyer and C. Kirches.  The SolvIND webpage.  http://www.solvind.org.
  * [AT98] Alan L. Andrew and Roger C. E. Tan.  Computation of derivatives of repeated eigenvalues and the corresponding eigenvectors of symmetric matrix pencils.  SIAM Journal on Matrix Analysis and Applications, 20(1):78–100, 1998.
  * [Bel10] B. Bell.  Cppad: a package for c++ algorithmic differentiation.  http://www.coin-or.org/CppAD, 2010.
  * [Ber01] Daniel J. Bernstein.  Multidigit multiplication for mathematicians, 2001.
  * [Ber08] Daniel J. Bernstein.  Fast multiplication and its applications.  Algorithmic Number Theory, Lattices, Number Fields, Curves and Cryptography, 2008.
  * [BH96] Christian H. Bischof and Mohammad R. Haghighat.  Hierarchical approaches to automatic differentiation.  In Martin Berz, Christian Bischof, George Corliss, and Andreas Griewank, editors, Computational Differentiation: Techniques, Applications, and Tools, pages 83–94. SIAM, Philadelphia, PA, 1996.
  * [Bjö96] Å. Björck.  Numerical Methods for Least Squares Problems.  SIAM, Philadelphia, 1996.
  * [Büc02] H. M. Bücker.  Hierarchical Algorithms for Automatic Differentiation.  Habilitationsschrift, Aachen, nov 2002.
  * [Chr91] Bruce Christianson.  Reverse accumulation and accurate rounding error estimates for Taylor series coefficients.  Optimization Methods and Software, 1(1):81–94, 1991.  Also appeared as Tech. Report No. NOC TR239, The Numerical Optimisation Centre, University of Hertfordshire, U.K., July 1991.
  * [Chr98] Bruce Christianson.  Reverse accumulation and implicit functions.  Optimization Methods and Software, 9(4):307–322, 1998.
  * [Eri03] Eric Todd Phipps.  Taylor Series Integration of Differential-Algebraic Equations: Automatic Differentation as a Tool For Simulationg Rigid Body Mechanical Systems.  PhD thesis, Cornell University, February 2003.
  * [GJM+99] Andreas Griewank, David Juedes, H. Mitev, Jean Utke, Olaf Vogel, and Andrea Walther.  ADOL-C: A package for the automatic differentiation of algorithms written in C/C++.  Technical report, Institute of Scientific Computing, Technical University Dresden, 1999.  Updated version of the paper published in ACM Trans. Math. Software 22, 1996, 131–167.
  * [Gri03] Andreas Griewank.  A mathematical view of automatic differentiation.  In Acta Numerica, volume 12, pages 321–398. Cambridge University Press, 2003.
  * [GUW00] Andreas Griewank, Jean Utke, and Andrea Walther.  Evaluating higher derivative tensors by forward propagation of univariate Taylor series.  Mathematics of Computation, 69:1117–1130, 2000.
  * [GVL96] Gene H. Golub and Charles F. Van Loan.  Matrix computations (3rd ed.).  Johns Hopkins University Press, Baltimore, MD, USA, 1996.
  * [GW08] Andreas Griewank and Andrea Walther.  Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation.  Number 105 in Other Titles in Applied Mathematics. SIAM, Philadelphia, PA, 2nd edition, 2008.
  * [KKSB07] S. Körkel, E. Kostina, J. P. Schlöder, and H. G. Bock.  Robustness Aspects in Parameter Estimation, Optimal Design of Experiments and Optimal Control.  Springer Berlin Heidelberg, 2007.
  * [Kör02] Stefan Körkel.  Numerische Methoden für Optimale Versuchsplanungsprobleme bei nichtlinearen DAE-Modellen.  PhD thesis, Universitat Heidelberg, 2002.
  * [Mär02] R. März.  The index of linear differential algebraic equations with properly stated leading terms.  Result. Math., 42:308–338, 2002.
  * [Mär03] R. März.  Differential Algebraic Systems with Properly Stated Leading Term and MNA Equations.  In Modeling, Simulation and Optimization of Integrated Circuits, International Series of Numerical Mathematics, volume 146, pages 135–151, Basel, 2003. Birkhauser Verlag.
  * [MN99] Jan R. Magnus and Heinz Neudecker.  Matrix differential calculus with applications in statistics and econometrics.  John Wiley & Sons, 2nd edition, 1999.
  * [MSA01] Joaquim R. R. A. Martins, Peter Sturdza, and Juan J. Alonso.  The connection between the complex-step derivative approximation and algorithmic differentiation, 2001.
  * [MSA03] Joaquim R. R. A. Martins, Peter Sturdza, and Juan J. Alonso.  The complex-step derivative approximation.  ACM Trans. Math. Softw., 29(3):245–262, 2003.
  * [vdAMM07] N.P. van der Aa, H.G. Ter Morsche, and R.R.M. Mattheij.  Computation of eigenvalue and eigenvector derivatives for a general complex-valued eigensystem.  Electronic Journal of Linear Algebra ISSN 1081-3810, 16:300–314, 2007.