\documentclass[12pt]{article}

\usepackage{fullpage,amsmath,amsfonts,enumitem}

\def\R{\mathbb{R}}
\def\d{\mathrm{d}}

\DeclareMathOperator{\Tr}{Tr}

\begin{document}

\title{Root-finding in $d$ dimensions -- Newton method and Broyden method}
\author{Advanced Programming}
\date{\today}
\maketitle

\section{Continuous differentiability in $d$ dimensions}

Recall that if $f \colon \R^d \to \R^d$ is at least once continuously
differentiable, then for every $x \in \R^d$ there exists a
\emph{Jacobian} matrix $Df(x)$ depending continuously on $x$ such that
\begin{equation}
  f(x + \Delta x)
  = f(x) + Df(x) \, \Delta x + o(\lVert \Delta x \rVert) \,.
  \label{e.diff}
\end{equation}
Here, $o(\lVert \Delta x \rVert)$ is shorthand for some function
$\phi_x \colon \R^d \to \R^d$, defined for every $x\in \R^d$, such
that
\begin{equation}
  \lim_{\Delta x \to 0} \frac{\phi_x(\Delta x)}{\lVert \Delta x
  \rVert} = 0 \,.
\end{equation}
The Jacobian is the matrix of partial derivatives,
\begin{equation}
  Df(x) =
  \begin{pmatrix}
    \partial_1 f_1 (x) & \cdots & \partial_d f_1 (x) \\
    \vdots & & \vdots \\
    \partial_1 f_d (x) & \cdots & \partial_d f_d (x) 
  \end{pmatrix} \,.
\end{equation}

\section{Newton's method in $d$ dimensions}

We proceed as in the case of Newton's method for functions of a single
variable.  Suppose that $x_n \in \R^d$ is the current approximation of
a root of $f$.  Writing $\Delta x_{n+1} = x_{n+1} - x_n$, we seek to
determine $x_{n+1}$ as the root of the linear approximation to $f$ at
$x=x_n$.  Dropping the higher-order terms in \eqref{e.diff}, this
linear approximation reads
\begin{equation}
  f(x_{n+1}) = f(x_n + \Delta x_{n+1})
  \approx f(x_n) + Df(x_n) \, \Delta x_{n+1} \,.
\end{equation}
Setting the right hand side to zero, we obtain the $d$-dimensional
Newton update $\Delta x_{n+1}$ by solving the linear system
\begin{equation}
  Df(x_n) \, \Delta x_{n+1} = -f(x_n) \,.
  \label{e.newton-system}
\end{equation}
This can be written in the form
\begin{equation}
  x_{n+1} = x_n - Df(x_n)^{-1} f(x_n)
\end{equation}
to emphasize that it is a proper generalization of the familiar
formula for Newton's method in one variable.  However, when $d$ is
large, it is better to solve \eqref{e.newton-system} directly as there
are iterative linear solvers that are much faster than the $O(d^3)$
complexity required to compute the inverse Jacobian.

Newton's method can be shown to converge quadratically in a
sufficiently small neighborhood of a root.  Even though, as in one
dimension, there is no guarantee that it converges at all for
arbitrary starting points $x_0$, in practice it is fairly robust and
fast.

The drawback is that a Jacobian is required, which may not be readily
available: for example, $f$ may not be available in the form of a
simple mathematical expression, but is defined implicitly as the
output of a complex piece of computer code.  A second drawback is that
the Jacobian often requires $O(d^2)$ memory which may be too much for
very large problems.  An example of such a dense Jacobian is the
discrete integral equation in Section~\ref{s.integral} below; on the
other hand, the Jacobian in a second example, a discrete boundary
value problem described in Section~\ref{s.bvp}, is sparse and can be
stored more efficiently.


\section{Secant-like methods in $d$ dimensions}

To generalize the derivation of the secant method from one dimension,
let us write $J_n$ to denote the approximation to the Jacobian at
iteration $n$.  Writing $f_n = f(x_n)$, setting
$\Delta f_n = f_n - f_{n-1}$, and dropping, once again, higher order
terms in \eqref{e.diff}, we obtain the so-called \emph{secant
condition}
\begin{equation}
  \Delta f_n = J_n \, \Delta x_n \,.
  \label{e.secant}
\end{equation}
This again is a linear system of equations, but the unknowns are the
$d^2$ elements of $J_n$.  On the other hand, we only have $d$
equations, so the system is highly underdetermined when $d>1$.

Broyden's idea (in a modern interpretation) is to start with a
guess for the initial Jacobian $J_0$.  Then, at each step, first
update $x_n$, then use the resulting $\Delta x_n$ and $\Delta f_n$
solve an optimization problem for $J_n$, where the secant condition
\eqref{e.secant} appears as a constraint.  In particular, the ``good
Broyden method'' results from
\begin{equation}
\begin{gathered}
  \text{minimize } \lVert J_n - J_{n-1} \rVert_F^2 \\
  \text{subject to } J_n \, \Delta x_n = \Delta f_n 
  \end{gathered}
\end{equation}
which, as shown in Appendix~\ref{a.update}, implies the update
\begin{gather}
  J_n = J_{n-1}
        + \frac{\Delta f_n - J_{n-1} \Delta x_n}%
               {\lVert \Delta x_n \rVert^2} \, \Delta x_n^T \,.
  \label{e.good-update}
\end{gather}
To avoid the need for matrix inversion (or for solving a linear
system) at every step, the Sherman--Morrison formula, see equation
\eqref{e.sm} in Appendix~\ref{a.sm}, can be used to convert this
update rule into an update for the inverse matrix,
\begin{equation}
  J_n^{-1} = J_{n-1}^{-1}
             + \frac{\Delta x_n - J_{n-1}^{-1} \Delta f_n}%
                    {\Delta x_n^T \, J_{n-1}^{-1} \Delta f_n} \,
               \Delta x_n^T \, J_{n-1}^{-1} \,.
  \label{e.good-inverse}
\end{equation}

Alternatively, the ``bad Broyden method'' results from 
\begin{equation}
\begin{gathered}
  \text{minimize } \lVert J_n^{-1} - J_{n-1}^{-1} \rVert_F^2 \\
  \text{subject to } J_n \, \Delta x_n = \Delta f_n 
  \end{gathered}
\end{equation}
which implies the update
\begin{gather}
  J_n^{-1} = J_{n-1}^{-1}
             + \frac{\Delta x_n - J_{n-1}^{-1} \Delta f_n}%
                    {\lVert \Delta f_n \rVert^2} \, \Delta f_n^T \,.
  \label{e.bad}
\end{gather}
While the ``good'' method appears to converge better in many simple
test cases, this is not universally true; there are also problems
where the ``bad'' method is better.  Moreover, the ``bad'' method
tends to be computationally cheaper.  Thus, this historical
distinction is not always justified \cite{griewank}.

This gives the following algorithm:
\begin{enumerate}
\item Choose an initial guess $x_0$ for the root and $J_0^{-1}$ for
the inverse Jacobian.

\item At step $n= 1, 2, \dots$, first update the location of the root
via
\begin{equation}
  x_{n} = x_{n-1} - J_{n-1}^{-1} \, f(x_{n-1}) \,.
\end{equation}

\item If $f(x_n)$ is smaller than some tolerance (e.g., by testing the
maximal absolute value of its components), terminate with success.

\item If not, update the inverse Jacobian via \eqref{e.good-inverse}
(``good Broyden'') or \eqref{e.bad} (``bad Broyden'') and go to Step~2.

\end{enumerate}
\emph{Note:} Without loss of generality, we may set $J_0=I$, the
identity matrix.  For if another (nonsingular) initial Jacobian $J_0$
proves to be better, we can replace the original problem $f(x)=0$ by
the equivalent problem $g(x)=0$ where
\begin{equation}
  g(x) = J_0^{-1} f(x) \,.
\end{equation}
Then the sequence of Broyden updates for $f$ with initial matrix $J_0$
is identical to the sequence of Broyden updates of $g$ with initial
matrix $I$.  Computationally, this inverse can often be computed
efficiently by solving a linear system (e.g.\ for the discrete
boundary value problem in Section~\ref{s.bvp} below).

\section{Limited memory Broyden methods}

In the Broyden method, each update adds a rank-$1$ matrix to the
approximate Jacobian or to the inverse Jacobian.  The update of the
inverse formulation of the Broyden method can be written
\begin{equation}
   J_n^{-1}
   = J_{n-1}^{-1} + u_{n} \, v_{n}^T
  \label{e.good-inverse2}
\end{equation}
or, by recursive insertion with $J_0=I$, 
\begin{equation}
   J_n^{-1}
   = I + \sum_{k=1}^{n} u_k \, v_{k}^T
  \label{e.good-inverse3}
\end{equation}
where, e.g.\ for the ``good'' Broyden method, see
\eqref{e.good-inverse},
\begin{equation}
  u_k = \frac{\Delta x_{k} - J_{k-1}^{-1} \Delta f_{k}}%
             {\Delta x_{k}^T \, J_{k-1}^{-1} \Delta f_{k}}
  \quad \text{and} \quad
  v_k^T = \Delta x_k^T \, J_{k-1}^{-1} \,.
  \label{e.update-vectors}
\end{equation}
Thus, if $d \gg n$, it is more efficient to store the $u_k$ and $v_k$.
In particular, when $J_0=I$, matrix-vector products can be written
\begin{equation}
  J_n^{-1} x = x + \sum_{k=1}^{n} u_k \, (v_{k}^T x)
  \label{e.matrix-vector}
\end{equation}
where the term in parentheses is an inner product, thus can be
evaluated in $O(d)$ elementary arithmetic operations, so that the
overall complexity of evaluating \eqref{e.matrix-vector} is $O(dn)$
which, when $d \gg n$, is much better than the $O(d^2)$ operations
required for a full matrix-vector product.

Since both memory requirements and computation time increase with each
Broyden iteration, it is often important to keep the rank of the
update matrix bounded, say, with maximal rank $p+1$.  This can be done
in two ways.

\subsection{Simple pruning}

When $n>p$, keep only the most recent $p$ vectors in the sum
\eqref{e.good-inverse3}.  To be precise, for the limited memory
Broyden (LMB) method, perform of Steps~1--3 as before, followed by
\begin{itemize}
\item[4'.]  Retain only $u_{n-p}, \dots, u_{n-1}$ and
$v_{n-p}, \dots, v_{n-1}$.  Then compute a new pair of update vectors
$u_n$ and $v_n$ (using \eqref{e.update-vectors} for the ``good''
Broyden method or analogous expressions for the ``bad'' Broyden
method), where the required matrix-vector products use
\eqref{e.matrix-vector}, and matrix-vector products use an analogous
expression, summing over update vectors with index $k=n-p, \dots, n-1$
only.  Then go to Step~2, retaining update vectors
$u_{n-p}, \dots, u_{n}$ and $v_{n-p}, \dots, v_{n}$.
\end{itemize}

\subsection{SVD-based pruning}

It is often better not to simply discard the oldest update vectors,
but rather compute a best rank-$p$ approximation of the update.  In
principle, a best rank-$p$ approximation with respect to the Frobenius
norm is known to be given by the singular value decomposition of the
update matrix, retaining only the $p$ largest singular values.  To
implement this efficiently when $d \gg p$, we proceed as follows (see
\cite{rotten}).

First, notice that \eqref{e.good-inverse3} can be written as
\begin{equation}
   J_n^{-1}
   = I + U V^T 
\end{equation}
where $U$ is $d \times q$ and $V^T$ is $q \times d$ (with $q = n$ if
no pruning was applied in previous steps) with
\begin{equation}
  U =
  \begin{pmatrix}
    \mid && \mid \\
    u_1 &  \dots &  u_q \\
    \mid && \mid 
  \end{pmatrix}
  \quad \text{and} \quad
  V^T =
  \begin{pmatrix}
    \text{---} & v_1^T & \text{---} \\
    & \vdots & \\
    \text{---} & v_q^T & \text{---}
  \end{pmatrix} \,.
\end{equation}
Then compute the $RQ$-decomposition of $V^T$,
\begin{equation}
  RQ = V^T \,,
\end{equation}
where $Q$ is an orthogonal $q \times d$-matrix and $R$ is $q \times
q$.  Next, compute the singular value decomposition
\begin{equation}
  U \Sigma W^T = UR \,.
\end{equation}
Assuming that the singular values in $\Sigma$ are ordered in
decreasing magnitude (as guaranteed in standard library
implementations), a new rank-$p$ pair of update matrices is given by
the first $p$ columns of $U\Sigma$ and the first $p$ rows of $W^TQ$,
respectively.

Using these reduced update matrices in the expressions for the Broyden
update, a new pair of update vectors is generated just as before, and
appended to the pruned update matrices.



\section{Test problems}

The following problems have been used historically as test problems
for Broyden-type methods (see, e.g., \cite[Appendix~A]{rotten}).

\subsection{Discrete boundary value problem}
\label{s.bvp}

On the interval $I = (0,1)$, the boundary value problem
\begin{gather}
  u''(t) = a \, (u(t) + t + 1)^3 \,, \\
  u(0) = u(1) = 0
\end{gather}
can be discretized on $d$ subintervals of equal length $h=1/d$.  Here,
$a$ is a parameter with default value $a=\tfrac12$.  Setting
$t_i = ih$ and writing $u_i$ as the discrete approximation to
$u(t_i)$, we can write the standard second order difference
approximation as
\begin{gather}
  \frac{u_{i-1} - 2 \, u_i + u_{i+1}}{h^2}
  - a \, (u_i + t_i + 1)^3 = 0
  \label{e.discrete-bvp}
\end{gather}
for $i=1, \dots, d-1$, with $u_0=u_d=0$.  Thus, the true dimension of
this system of nonlinear equations is $d-1$.  Note that if the vector
of unknowns is denoted $u = (u_1, \dots, u_{d-1})^T$, the nonlinear
system \eqref{e.discrete-bvp} takes the form
\begin{equation}
  B u - N(u) = 0 \,,
  \label{e.discrete-bvp2}
\end{equation}
where
\begin{equation}
  B = \frac1{h^2}
  \begin{pmatrix}
  -2 & 1  &  & 0\\
  1 & \ddots & \ddots  \\
    & \ddots & \ddots & 1 \\
  0 & & 1 & -2
  \end{pmatrix}
  \quad \text{and} \quad
  N(u) = a \, 
  \begin{pmatrix}
    (u_1 + t_1 + 1)^3 \\ \vdots \\ (u_{d-1} + t_{d-1} + 1)^3 
  \end{pmatrix} \,.
\end{equation}
Rearranging terms, \eqref{e.discrete-bvp2} can be written in fixed
point form,
\begin{equation}
  u = B^{-1} N(u) \equiv \Phi(u) \,.
  \label{e.discrete-bvp-fp}
\end{equation}
Computationally, the inverse of $B$ should not be computed explicitly,
as that would destroy the sparse, tri-diagonal structure of $B$, but
rather be implemented by solving a system of linear equations.


\subsection{Discrete integral equation}
\label{s.integral}

The nonlinear integral equation
\begin{gather}
  u(t) + a \int_0^1 H(s,t) \, (u(s) + s + 1)^3 \, \d s = 0
\end{gather}
kernel
\begin{gather}
  H(s,t) =
  \begin{cases}
    s \, (1-t) & \text{for } s \leq t \,, \\
    t \, (1-s) & \text{for } s > t \,,
  \end{cases}
\end{gather}
is approximated, on the same partition as before, by the discrete
Riemann sum
\begin{equation}
  u_i + a \, h \,
  \biggl(
    \sum_{j=1}^i (1-t_i) \, t_j \, (u_j + t_j + 1)^3
    + \sum_{j=i+1}^d t_i \, (1-t_j) \, (u_j + t_j + 1)^3
  \biggr) = 0
  \label{e.discrete-int}
\end{equation}
for $i = 1, \dots, d$.  As before, the parameter $a$ takes the default
value $a=\tfrac12$.  Writing $u = (u_1, \dots, u_d)$ to denote the
vector of unknowns, \eqref{e.discrete-int} can be written as
\begin{equation}
  u - H N(u) = 0
\end{equation}
where $H=(h_{ij})$ is the $d \times d$-matrix with matrix elements
\begin{equation}
  h_{ij} =
  \begin{cases}
    t_j \, (1-t_i) & \text{for } i \geq j \\
    (1-t_j) \, t_i & \text{for } i < j
  \end{cases}
\end{equation}
and
\begin{equation}
  N(u) = - a \, h
  \begin{pmatrix}
    (u_1 + t_1 + 1)^3 \\ \vdots \\ (u_{d} + t_{d} + 1)^3 
  \end{pmatrix} \,.
\end{equation}
Thus, this test problem also has a fixed point form similar to
\eqref{e.discrete-bvp-fp}, namely
\begin{equation}
  u = HN(u) \equiv \Phi(u) \,.
\end{equation}

\appendix
\section{Technical details}

\subsection{The trace of a matrix}

Matrix optimization problems involving the Frobenius norm can often be
written in a coordinate-free way using the trace.  This leads to
shorter and more concise computations, and will be the strategy
pursued here.

Given a square matrix $A = (a_{ij})_{ij} \in \R^{d \times d}$, its
\emph{trace} is defined as the sum of its diagonal elements, i.e.
\begin{equation}
  \Tr A = \sum_{i=1}^d a_{ii} \,.
\end{equation}
The following properties, which are easy to check and will be stated
without proof, are used:
\begin{enumerate}[label={\upshape(\roman*)}]
\item The trace is a linear map,
\item $\Tr A = \Tr A^T$ (invariance under transposition),
\item\label{iii} $\Tr (AB) = \Tr (BA)$ for all $A \in \R^{d\times k}$
and $B \in \R^{k \times d}$.
\end{enumerate}
Property \ref{iii} implies immediately that the trace is invariant
under all cyclic permutations of matrix products and that the inner
product of two vectors $u, v \in \R^d$ can be written
\begin{equation}
  u^Tv = \Tr(u^Tv) = \Tr(v u^T) \,.
  \label{e.inner}
\end{equation}
\begin{enumerate}[label={\upshape(\roman*)},resume]
\item\label{iv} $\langle A, B \rangle_F=\Tr (AB^T)$ defines an inner
product on the vector space of $d \times k$ matrices.
\end{enumerate}
This inner product is called the \emph{Frobenius inner product}, and
the associated norm the \emph{Frobenius norm},
\begin{equation}
  \lVert A \rVert_{F}^2
  = \Tr (A A^T)
  = \sum_{i=1}^d \sum_{j=1}^k a_{ij}^2 \,.
  \label{e.fnorm}
\end{equation}

\subsection{Lagrange multipliers}

Recall the Lagrange multiplier theorem from Analysis II: Let
$f \in C^1(\R^d, \R)$ and $g \in C^1 (\R^d, \R^k)$.  Suppose $f(x)$
has a local extremum at $x=x^* \in \R^d$ subject to the constraint
$g(x)=0$, and the Jacobian matrix $Dg(x^*)$ has full rank.  Then there
exists a unique $\lambda \in \R^k$ such that the Lagrangian
\begin{equation}
  L(x;\lambda) = f(x) + \lambda^T g(x)
\end{equation}
has a critical point at $x=x^*$, i.e.\ $\nabla L(x^*) = 0$.

In practice, the theorem is always applied in reverse order: we
compute solutions of $\nabla L(x^*) = 0$ as \emph{candidate}
solutions to the constrained extremal problem, then use other means
(e.g.\ convexity, compactness) to argue for sufficiency.

In our setting, the objective function $f$ is the Frobenius norm of a
matrix expression.  In other words, the domain of $f$ is not $\R^d$,
but some space of matrices.  While it is certainly possible to
identify a space of matrices with some $\R^d$ of matching dimension
and use the formulation of the Lagrange multiplier theorem as stated
above, this leads to terrible notation and is best avoided in favor of
a slightly more abstract approach, which is explained now.

Let $M$ be a normed vector space.  Let $x \colon \R \to M$ be a smooth
curve in $M$.  Then $x'(t)$ is a tangent vector to the curve for any
$t \in \R$.  Since $M$ is a vector space, its tangent space is
isomorphic to $M$ itself.  In particular,
\begin{equation}
  \delta x \equiv x'(0)
  \label{e.tangent}
\end{equation}
is a tangent vector and as such, is an element from $M$.  Vice versa,
any tangent vector is of the form \eqref{e.tangent}.  More generally,
for a function $f \in C^1(M,\R)$, we define the \emph{variational
derivative}
\begin{equation}
  \delta f = \frac{\d f(x(t))}{\d t} \bigg|_{t=0}
           = Df(x(0)) \, \delta x \,.
\end{equation}
Then $f$ has a critical point at $x=x^*$ if and only if
$\delta f(x^*)=0$ for any $\delta x \in M$.

In this setting, the Lagrange multiplier theorem reads: Let
$f \in C^1(M, \R)$ and $g \in C^1 (M, \R^k)$.  Suppose $f(x)$ has a
local extremum at $x=x^* \in M$ subject to the constraint $g(x)=0$,
and the Jacobian operator $Dg(x^*)$ has full rank.  Then there exists
a unique $\lambda \in \R^k$ such that the Lagrangian
\begin{equation}
  L(x;\lambda) = f(x) + \lambda^T g(x)
  \label{e.lagrangian}
\end{equation}
has a critical point at $x=x^*$, i.e.\ $\delta L(x^*) = 0$ for any
$\delta x \in M$.

\emph{Important note:} In typical applications in the calculus of
variations, this setting is further generalized by taking $M$ to be
``admissible sets'' that are not vector spaces, but, for example,
affine spaces or manifolds.  In this case, tangent vectors $\delta x$
are not elements from $M$.  In this case, one has to think more
carefully about the concept of tangent space.  Everything above still
stands, and the notation is useful, except that the membership of
$\delta x$ needs to be adapted to the correct notion of tangent space.

\subsection{Derivation of the ``good'' Broyden method}
\label{a.update}

Recall that the ``good'' Broyden method assumes that $J_{n-1}$ is
given and $J_{n}$ is determined by solving the constrained optimization
problem 
\begin{equation}
\begin{gathered}
  \text{minimize } \lVert J_n - J_{n-1} \rVert_F^2
  = \Tr [(J_n - J_{n-1})(J_n - J_{n-1})^T] \\
  \text{subject to } J_n \, \Delta x_n = \Delta f_n  \,.
  \end{gathered}
\end{equation}
(Cf.\ \eqref{e.fnorm} for expressing the Frobenius norm in terms of
the trace.)  Rewriting the inner product in \eqref{e.lagrangian} in
terms of the trace, cf.\ \eqref{e.inner}, we can write out the
Lagrangian for this problem as
\begin{equation}
  L(J_n; \lambda)
  = \tfrac12 \, \Tr [(J_n - J_{n-1})(J_n - J_{n-1})^T]
  + \Tr [(J_n \, \Delta x_n - f_n) \lambda^T] \,.
\end{equation}
(The factor $\tfrac12$ is a convenience factor that does not change
the problem because linear scaling of the objective function does not
change the location of the extrema.)

We compute
\begin{align}
  \delta L
  & = \tfrac12 \, \Tr [\delta J_n \, (J_n - J_{n-1})^T]
      + \tfrac12 \, \Tr [(J_n - J_{n-1}) \, \delta J_n^T]
      + \Tr [\delta J_n \, \Delta x_n \, \lambda^T]
      \notag \\
  & = \Tr [(J_n - J_{n-1}) \, \delta J_n^T]
      + \Tr [\lambda \, \Delta x_n^T \, \delta J_n^T]
      \notag \\
  & = \Tr [(J_n - J_{n-1} + \lambda \, \Delta x_n^T) \,
           \delta J_n^T]
\end{align}
where, in the second equality, we have used the invariance of the
trace under transposition.  Since we can interpret the final
expression as an inner product between two matrices, cf.\
property~\ref{iv} above, we see that $\delta L=0$ for all $\delta J_n
\in \R^{d \times d}$ if and only if
\begin{equation}
  J_n - J_{n-1} + \lambda \, \Delta x_n^T = 0 \,.
  \label{e.update0}
\end{equation}
To eliminate the Lagrange multiplier $\lambda$ from this expression,
we multiply \eqref{e.update0} with $\Delta x_n$ from the right and use
the constraint $J_n \, \Delta x_n = \Delta f_n$:
\begin{equation}
  \Delta f_n - J_{n-1} \, \Delta x_n 
    + \lambda \, \Delta x_n^T \, \Delta x_n = 0 \,,
  \label{e.update1}
\end{equation}
so that
\begin{equation}
  \lambda = \frac{J_{n-1} \, \Delta x_n - f_n}%
                 {\Delta x_n^T \, \Delta x_n} \,.
  \label{e.lambda}
\end{equation}
Inserting \eqref{e.lambda} into \eqref{e.update0} and solving for
$J_n$, we obtain the ``good'' Broyden update \eqref{e.good-update}.
The ``bad'' Broyden update \eqref{e.bad} can be derived similarly.


\subsection{Derivation of the Sherman--Morrison formula}
\label{a.sm}

Suppose we have an invertible matrix $A \in \R^{d \times d}$ and
perform a rank-$1$ update, setting
\begin{equation}
  B = A + uv^T
\end{equation}
for some given vectors $u, v \in \R^d$.  We ask the question whether
the inverse of $A$ is related to the inverse of $B$ also by a rank-$1$
update.  In other words, do there exist $x, y \in \R^d$ such that the
ansatz
\begin{equation}
  B^{-1} = A^{-1} + x y^T
  \label{e.ansatz}
\end{equation}
holds true?

To answer this question, let us verify the properties of the inverse
matrix.  In order for \eqref{e.ansatz} to hold, we must have
\begin{equation}
  I = B B^{-1}
    = (A + uv^T)(A^{-1} + x y^T)
    = I + A x y^T + u v^T A^{-1} + uv^T x y^T 
\end{equation}
or, canceling $I$ and multiplying by $A^{-1}$ from the left,
\begin{equation}
  x y^T + A^{-1} u v^T A^{-1} + A^{-1} uv^T x y^T = 0 \,.
  \label{e.firstcondition}
\end{equation}
We can thus determine the \emph{direction} of $x$ by multiplying
\eqref{e.firstcondition} with some column vector from the right.  Then
most factors in the resulting expression will be scalar, except for $x$
in the first term and $A^{-1} u$ in the second and third term.  Thus,
$x$ must be in the direction of $A^{-1} u$.

Similarly, we must have
\begin{equation}
  I = B^{-1} B
    = (A^{-1} + x y^T)(A + uv^T)
    = I + A^{-1} u v^T + x y^T A + x y^T uv^T 
\end{equation}
or, canceling $I$ and multiplying by $A^{-1}$ from the right,
\begin{equation}
  A^{-1} u v^T A^{-1} + x y^T + x y^T uv^T A^{-1} = 0 \,.
  \label{e.secondcondition}
\end{equation}
We can thus determine the \emph{direction} of $y$ by multiplying
\eqref{e.secondcondition} with some row vector from the left.  Then
most factors in the resulting expression will be scalar, except for
$y^T$ in the second term and $v^T A^{-1}$ in the first and third term.
Thus, $y^T$ must be in the direction of $v^T A^{-1}$.

Thus, we have established that \eqref{e.ansatz} can only hold true if
\begin{equation}
  x y^T = \alpha \, A^{-1} u v^T A^{-1}
  \label{e.xyt}
\end{equation}
for some value of $\alpha$.  To determine $\alpha$, we multiply
\eqref{e.firstcondition} or, equivalently, \eqref{e.secondcondition}
by $v^T$ from the left and by $u$ from the right, to obtain
\begin{equation}
  v^T x y^T u + (v^T A^{-1} u)^2 + v^T A^{-1} uv^T x y^T u = 0 \,.
  \label{e.alpha1}
\end{equation}
Likewise, we multiply \eqref{e.xyt} by $v^T$ from the left and by $u$
from the right, to obtain
\begin{equation}
   v^T x y^T u = \alpha \, (v^T A^{-1} u)^2 \,.
  \label{e.alpha2}
\end{equation}
Inserting \eqref{e.alpha2} into \eqref{e.alpha1}, we find that
\begin{equation}
  \alpha = - \frac{1}{1 + v^T A^{-1} u} \,.
\end{equation}

This derivation is only formal, as it assumes invertibility of $B$,
but suggests the following theorem: Suppose $A \in \R^{d \times d}$ is
invertible and $u, v \in \R^d$.  Then $A + uv^T$ is invertible if and
only if $1 + v^T A^{-1} u \neq 0$.  In this case,
\begin{equation}
  (A + uv^T)^{-1}
  = A^{-1} - \frac{A^{-1} u v^T A^{-1}}{1 + v^T A^{-1} u} \,.
  \label{e.sm}
\end{equation}

A full proof is an exercise in Linear Algebra and left to the reader.
Expression \eqref{e.sm} is known as the \emph{Sherman--Morrison
formula}.




\begin{thebibliography}{99}
\bibitem{griewank} A. Griewank, 
\textit{Broyden updating, the Good and the Bad!}
Doc. Math. (2012), 301--315. 

\bibitem{rotten} B. van de Rotten,
\textit{A limited memory Broyden method to solve
high-dimensional systems of nonlinear equations}, PhD thesis,
Universiteit Leiden, 2003.

\end{thebibliography}

\end{document}

