# Automatic Differentiation - reverse accumulation of linear system solve

I am studying the reverse mode of automatic differentiation.

The reverse mode of automatic differentiation allows the efficient computation of a the derivative of a single dependent variable $y$ with respect to as many independent variables $x_i$ as you want. One assigns to each intermediate variable $v$ an adjoint variable $\bar{v}$ which is the derivative of a chosen dependent variable with respect to the subexpression $$\bar{v} \rightarrow \bar{v} = \frac{\partial y}{\partial v}$$

so the assignment $$y = v_1 \sin(v_2)$$ corresponds to the adjoint assignments

\begin{eqnarray} \bar{v_1} &=& \bar{y} \sin(v_2) \\ \bar{v_2} &=& \bar{y} v_1 \cos(v_2) \end{eqnarray}

where $\bar{y} = 1$ .

I am interested in the situation where you solve a linear system in the program:

$$Ax = b$$

where $y$ might be another function of $x$: $$y(x)$$

According to this tutorial of the CoDiPack software, the corresponding adjoint statements are

\begin{eqnarray} \bar{A} &=& - \lambda x^T \\ \bar{b} &=& \lambda \\ \end{eqnarray} where $\lambda$ is the the solution of the adjoint equation $$A^T \lambda = \bar{x}$$

I found the same algorithm in several other documents, for example here in section 7, Iteration and equation solving.

It is not clear to me how to arrive at these statements. I think the derivation must be similar as in the case where one wants to optimize $w(x)$ where $x$ is subject to the constraint $$Ax = b$$ See for example this document.

Starting with $b=Ax$ the forward differentiation gives $\dot b=\dot A x+A\dot x$, or, with $\dot y = \dot A x$:

$$\dot x = -A^{-1} \dot y + A^{-1} \dot b$$

Let's write this in terms of the local transformation Jacobian for the tangent linear variables $\left( \dot A, \dot y, \dot b, \dot x\right)$:

$$\left[ \begin{matrix} \dot A \\ \dot y\\ \dot b\\ \dot x\\ \end{matrix} \right] = \left[ \begin{matrix} 1 & 0 & 0 & 0 \\ x & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -A^{-1} & A^{-1} & 0 \\ \end{matrix} \right] \left[ \begin{matrix} \dot A \\ \dot y\\ \dot b\\ \dot x\\ \end{matrix} \right]$$

The adjoint operator is then the transpose matrix acting on the (local) adjoint variables $\left( \bar A, \bar y, \bar b, \bar x\right)$:

$$\left[ \begin{matrix} \bar A \\ \bar y\\ \bar b\\ \bar x\\ \end{matrix} \right] = \left[ \begin{matrix} 1 & x^T & 0 & 0 \\ 0 & 0 & 0 & -A^{-T} \\ 0 & 0 & 1 & A^{-T} \\ 0 & 0 & 0 & 0 \\ \end{matrix} \right] \left[ \begin{matrix} \bar A \\ \bar y\\ \bar b\\ \bar x\\ \end{matrix} \right],$$

or:

\begin{align} \bar y &= -A^{-T} \bar x \\ \bar A &\mathrel{+}= x^T \bar y \\ \bar b &\mathrel{+}= A^{-T} \bar x \\ \bar x &= 0. \end{align}

Reference:

Giering and Kaminski, Recipes for Adjoint Code Construction, ACM TOMS, 1998

• Thanks for pointing me to this reference, reading it was quite helpful. I noted that you can actually do the same derivation without introducing the variable $\dot{y}$. – benno May 2 '18 at 10:19
• yes. the $\dot y$ keeps the explanation a bit simpler, I think. – GoHokies May 2 '18 at 10:24

The identity connecting forward and reverse mode for $x=\text{solve}(A,b)$ is $$\bar x^⊤\dot x=\bar A(\dot A)+\bar b^⊤\dot b.$$ Where $\bar A$ is a linear functional on matrices which may be realized as the $trace(\bar A^⊤\dot A)$ which is related to the (Frobenius) scalar product for matrices $\langle X,Y\rangle=trace(X^⊤Y)=trace(XY^⊤)$, which is kind-of the euclidean scalar product as it is the sum over the element-wise product of $X$ and $Y$.

With $b=Ax$ the forward differentiation gives $\dot b=\dot A x+A\dot x$, inserting you get $$\bar x^⊤\dot x=\bar A(\dot A)+\bar b^⊤(\dot A x+A\dot x). \tag{1}$$ As that has to be valid for all choices of $\dot A, \dot x$, one gets \begin{align} 0&=trace(\bar A^⊤\dot A)+ \bar b^⊤\dot A x \\ &=trace(\bar A^⊤\dot A)+ trace(\bar b^⊤\dot A x) \\ &=trace(\bar A^⊤\dot A)+ trace(x\bar b^⊤\dot A) \\ &=trace((\bar A+\bar bx^⊤)^⊤\dot A), \\\text{and }~\bar b A\dot x&=\bar x\dot x \tag{2} \end{align} which implies $$\bar b^⊤A=\bar x^⊤\iff \bar b=\text{solve}(A^⊤,\bar x),\tag3 \\ \bar A = -\bar bx^⊤. \\$$

• Can you explain further how you arrive at the connecting identity? Is the right hand side after inserting the expression for $\dot{b}$ correct? Should it not only be $\bar{x}^T\dot{x}$? – benno May 2 '18 at 10:34
• You are correct, I started out with $b=Ax$ so that $A,x$ are the inputs and $b$ the output and did not correctly and completely change to $x=solve(A,b)$ where $A,b$ are the inputs and $x$ the output. – LutzL May 2 '18 at 10:58
• We have the identity $\bar{x}\dot{x} = \bar{y}\dot{y}$ that connects input /output of a function in reverse and forward mode. In order to be consistent, should we not choose an element-wise multiplication as the functional $\bar{A}(\cdot)$? – benno May 3 '18 at 6:57
• I did. $trace(\bar A^⊤\dot A)$ is the Euclidean scalar product for matrices. – LutzL May 3 '18 at 7:06
• Thanks for the clarification, I didn't realize that it is actually the same. I still cannot see how you go from the statement where you inserted $\dot{b} = \dot{A}x + A\dot{x}$ to the next one. Is there a dot missing in $trace(\bar{A}^T{A}) \rightarrow trace(\bar{A}^T\dot{A})$? – benno May 4 '18 at 16:24