\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,hyperref}

\newcommand{\rem}[1]{}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\ddh}[2]{\mathrm d #1/\mathrm d #2}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\ddhp}[3]{\mathrm d^{#3} #1/{\mathrm d #2}^{#3}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\opttwo}[2]{\left\{\begin{array}{l}#1\\#2\end{array}\right.}
\newcommand{\vecthree}[3]{\left[\begin{array}{c}#1\\#2\\#3\end{array}\right]}
\newcommand{\mxtwo}[1]{\left[\begin{array}{cc}#1\end{array}\right]}
\newcommand{\mxthree}[1]{\left[\begin{array}{ccc}#1\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}

\newcommand{\ii}{\mathrm i}
\newcommand{\e}{\mathrm e}

\newcommand{\B}{\mathbf B}
\newcommand{\C}{\mathbf C}
\newcommand{\F}{\mathbf F}
\newcommand{\x}{\mathbf x}
\newcommand{\p}{\mathbf p}
\newcommand{\vv}{\mathbf v}
\newcommand{\uv}{\mathbf u}
\newcommand{\bv}{\mathbf b}
\newcommand{\cv}{\mathbf c}
\newcommand{\dv}{\mathbf d}
\newcommand{\g}{\mathbf g}

\newcommand{\diag}{\mathrm{diag}}

\begin{document}

\title{Bounded Approximate Solutions of Linear Systems using SVD}
\author{Stephen Brooks}
\maketitle

\section{Definitions}
The Singular Value Decomposition (SVD) of a complex matrix is conventionally $A=U\Sigma V^*$, where $M^*$ denotes ${\bar M}^T$.  Here, $U$ and $V$ are unitary matrices with $U^{-1}=U^*$ and $\Sigma$ is diagonal with $\Sigma=\diag[\sigma_n]$.  For real matrices this is just $A=U\Sigma V^T$ and unitarity is equivalent to $U^{-1}=U^T$, i.e. orthogonality.  In fact, $V^T$ is also orthogonal since $(V^T)^{-1}=(V^{-1})^{-1}=V=(V^T)^T$, which means the simpler definition $A=U\Sigma V$ can be used for the rest of this note.

\section{Fundamental Problem}
In control systems, one often uses a linear or locally-linear model to determine the required inputs.  Suppose an input vector change $\x\in X$ produces an output reponse $A\x\in Y$ that is meant to achieve some desired change $\bv\in Y$.  The input and output spaces $X$ and $Y$ may have different dimensionalities and therefore $A$ can be a rectangular matrix.  This means that an exact solution may not be possible, particularly if $\dim Y>\dim X$.  Thus the `best' solution can be formulated as the minimisation problem of finding $\arg\min|A\x-\bv|_Y$.

However, particularly in the case of ill-conditioned matrices, the exact solution may require unacceptably large control inputs.  What is required practically is the best approximation that can be achieved while $\x$ is not too large.  This suggests casting the fundamental problem as \[\arg\min_{|\x|_X\le r}|A\x-\bv|_Y\] with $r>0$ being chosen depending on how large a solution is acceptable.  As $r\rightarrow\infty$, the value will eventually settle at the exact or optimum solution if one exists.

\section{Solution using SVD}

The SVD decomposition of $A$ gives
\[  \arg\min_{|\x|_X\le r}|A\x-\bv|_Y = \arg\min_{|\x|_X\le r}|U\Sigma V\x-\bv|_Y. \]
Here, $A$ and $\Sigma$ are possibly-rectangular matrices mapping from $X$ to $Y$, $V$ is a square orthogonal matrix mapping $X$ to itself and $U$ is another mapping $Y$ to itself.  Note that any orthogonal matrix $U$ preserves the norm as $|U\x|^2=\x^TU^TU\x=\x^TU^{-1}U\x=\x^T\x=|\x|^2$ so $|U\x|=|\x|$ as norms are non-negative.  In particular,
\[ |\x|_X=|V\x|_X \qquad\mbox{and}\qquad |U\Sigma V\x-\bv|_Y=|\Sigma V\x-U^{-1}\bv|_Y, \]
where the second equality has multiplied by the unitary matrix $U^{-1}$.  This means that
\[  \arg\min_{|\x|_X\le r}|A\x-\bv|_Y = \arg\min_{|V\x|_X\le r}|\Sigma V\x-U^{-1}\bv|_Y. \]
Defining vectors $\vv=V\x$ and $\uv=U^{-1}\bv$ this becomes
\[  \arg\min_{|\x|_X\le r}|A\x-\bv|_Y = V^{-1}\arg\min_{|\vv|_X\le r}|\Sigma\vv-\uv|_Y, \]
where the right-hand $\arg\min$ is now understood to find the value of $\vv$, so the premultiplication for $\x=V^{-1}\vv$ is required.  The problem has now been simplified into one with a diagonal matrix instead of $A$.

\subsection{Exact Minimum Solution}
If the unrestricted $\arg\min$ also satisfies $|\x|_X\le r$ then it is the solution.  The unrestricted minimum is a fixed point of the norm expression squared:
\begin{eqnarray*}
0 &=& \pdd{}{v_n}|\Sigma\vv-\uv|_Y^2 
= \pdd{}{v_n}\sum_{i=1}^{\dim Y}(\Sigma\vv-\uv)_i^2 
= \pdd{}{v_n}\sum_{i=1}^{\dim Y}(1_{i\le\dim X}\sigma_i v_i-u_i)^2 \\
&=& \pdd{}{v_n}(\sigma_n v_n-u_n)^2 
= \pdd{}{v_n}(\sigma_n^2 v_n^2-2\sigma_n v_n u_n+u_n^2) 
= 2\sigma_n^2 v_n-2\sigma_n u_n \\
&\Leftrightarrow& \sigma_n(\sigma_n v_n-u_n)=0.
\end{eqnarray*}
For each $n$, this is true if either $v_n=u_n/\sigma_n$ or $\sigma_n=0$.  In the latter case, the $\Sigma$ matrix does not range over the full dimensionality of $Y$ and any value of $v_n$ may be chosen because the minimum is non-unique.  It is usually best to choose $v_n=0$ in all such ambiguous cases, since this corresponds to the minimum with smallest $|\vv|_X=|\x|_X$.  There is also the case when $\dim Y<\dim X$, where the above equation reduces to $0=0$ for $n>\dim Y$, giving no constraint on $v_n$, which should be set to zero by the same argument.  The exact minimum can be written explicitly as
\[ \x = V^{-1} [(U^{-1}\bv)_n/^0\sigma_n],
\quad\mbox{ where }\quad
x/^0y = \opttwo{x/y \quad \mbox{ if } y\ne0}{0 \quad \mbox{ otherwise}}. \]

\subsection{Constrained Minimum}
The function $|\Sigma\vv-\uv|_Y$ does not have multiple disconnected local minima, so if the exact minimum with smallest norm found in the previous section still has $|\x|_X>r$, the constrained minimum must have $|\x|_X=r$ rather than being an interior point.  The local gradient found in the previous section
\[ \nabla_\vv|\Sigma\vv-\uv|_Y^2 = 2[\sigma_n^2 v_n-\sigma_n u_n] \]
must be a scalar multiple of the position $\vv$ because otherwise it has some component parallel to the surface of the radius $r$ hypersphere and the value of the function can be reduced.  The gradient is expected to be negative with increasing $r$, anti-parallel to $\vv$, so for some $\lambda>0$, 
\begin{eqnarray*}
\nabla_\vv|\Sigma\vv-\uv|_Y^2 &=& -2\lambda^2\vv \\
\Leftrightarrow \quad 2(\sigma_n^2 v_n-\sigma_n u_n) &=& -2\lambda^2 v_n \\
\Leftrightarrow \quad (\sigma_n^2+\lambda^2)v_n -\sigma_n u_n &=& 0 \\
\Leftrightarrow \quad v_n &=& \frac{\sigma_n u_n}{\sigma_n^2+\lambda^2}. 
\end{eqnarray*}
For the case where $n>\dim Y$, the gradient of that component is zero as before and $0=-2\lambda^2 v_n$, so $v_n=0$.   The constrained minimum can be written explicitly as
\[ \x = V^{-1} \left[\frac{\sigma_n(U^{-1}\bv)_n}{\sigma_n^2+\lambda^2}\right],
\quad\mbox{ where we set }\quad
(U^{-1}\bv)_n=0 \quad \mbox{ if } n>\dim Y.\]

The norm of $\x$ decreases monotonically with $\lambda$ because $|\x|_X=|\vv|_X$ and every element of $\vv$ decreases in magnitude with increasing $\lambda$.
As $\lambda\rightarrow0$ the constrained minimum tends towards the exact minimum.  As $\lambda\rightarrow\infty$, the constrained minimum tends towards $\mathbf 0$ but if renormalised, the limit has $v_n=\sigma_n u_n$, which is $-\frac12$ times the gradient of $|\Sigma\vv-\uv|_Y^2$ at $\vv=\mathbf 0$.  Thus the large $\lambda$ limit corresponds to a infinitesimal `steepest descent' step.

The continuity and monotonicity of $|\x|_X=r(\lambda)$ ensures a value of $\lambda$ can always be found for any value of $r$ between 0 and the norm of the exact solution point.  For example, a bisection search or root-finding algorithm can determine $\lambda$ for a given $r$, after first checking the exact solution point does not have norm less than $r$.

\subsection{Implementation Note}
Using the orthogonal property of $U$ and $V$, entries $(U^{-1}\bv)_n$ should be calculated as the much faster equivalent $(U^T\bv)_n$ and the premultiplication by $V^{-1}$ should be implemented as $V^T$.  Once the SVD is calculated, nothing slower than matrix-vector multiplication is required.

\section{Units}
Elements of the vector spaces $X$ and $Y$ can be physical quantities with units $[X]$ and $[Y]$ respectively.  By definition, $A$ has units $[Y]/[X]$.  In the SVD, the entries of $U$ and $V$ have no units as they map within the same space, leaving $\Sigma$ and its entries $\sigma_n$ with units $[Y]/[X]$.  The parameter $\lambda$ in the previous section was defined to also have units $[Y]/[X]$ but $r$ has units $[X]$.

\section{Identity with the Levenberg--Marquardt Algorithm}
The Levenberg--Marquardt algorithm involves a `damped' least squares step, which for a Jacobian matrix $J$ involves solving
\[ (J^TJ+\lambda_{LM}I)\x = J^T\bv, \]
where $\lambda_{LM}\ge0$ is called the damping factor.  If the Jacobian is decomposed via SVD as $J=U\Sigma V$, this becomes
\[ (V^T\Sigma U^TU\Sigma V+\lambda_{LM}I)\x = V^T\Sigma U^T\bv \]
and noting that $U^TU=I$ by orthogonality of U,
\[ (V^T\Sigma^2 V+\lambda_{LM}I)\x = V^T\Sigma U^T\bv. \]
Pre-multipliying both sides by $V$ and using its orthogonality $VV^T=I$ gives
\begin{eqnarray*}
 (\Sigma^2 V+\lambda_{LM}V)\x &=& \Sigma U^T\bv \\
\Rightarrow \quad (\Sigma^2+\lambda_{LM}I)V\x &=& \Sigma U^T\bv.
\end{eqnarray*}
This is starting to look vaguely familiar.  Inverting the left-hand side to give an expression for $\x$ yields
\begin{eqnarray*}
 \x &=& V^{-1}(\Sigma^2+\lambda_{LM}I)^{-1}\Sigma U^T\bv \\
&=& V^{-1}(\Sigma^2+\lambda_{LM}I)^{-1}\Sigma U^{-1}\bv.
\end{eqnarray*}
Comparing this to the constrained minimum formula with parameter $\lambda$ from a previous section:
\[ \x = V^{-1} \left[\frac{\sigma_n(U^{-1}\bv)_n}{\sigma_n^2+\lambda^2}\right] \]
and noting that $\Sigma=\diag[\sigma_n]$ reveals that these are the same formulae if $\lambda_{LM}=\lambda^2$.

\section{Constrained Maximum of a Quadratic}
As the $|\Sigma\vv-\uv|_Y^2$ minimised in the previous sections was a quadratic function of $\x$, it is natural to wonder if an arbitrary (scalar) quadratic function could be maximised using a similar method: that is, find
 \[ \arg\max_{|\x|\le r}f(\x) = \arg\max_{|\x|\le r}\left(f(\mathbf 0)+\g\cdot\x+\frac12\x^TH\x\right). \] 
$H$ is the Hessian matrix of second derivatives, so is symmetric, meaning its SVD decomposition can be written $H=U^T\Sigma U$, with $U$ orthogonal.  This permits a change of variable
\begin{eqnarray*}
f(\x) &=& f(\mathbf 0)+\g^T\x+\frac12\x^TU^T\Sigma U\x \\
 &=& f(\mathbf 0)+\g^TU^T(U\x)+\frac12(U\x)^T\Sigma(U\x) \\
\Rightarrow \qquad \arg\max_{|\x|\le r}f(\x) &=& \arg\max_{|U\x|\le r}\left(f(\mathbf 0)+(U\g)^T(U\x)+\frac12(U\x)^T\Sigma(U\x)\right).
\end{eqnarray*}
Defining $\uv=U\x$ and ignoring the constant term, this becomes 
\[ \arg\max_{|\x|\le r}f(\x) = U^T\arg\max_{|\uv|\le r}\left((U\g)^T\uv+\frac12\uv^T\Sigma\uv\right). \]
The maximised expression is a single sum as $\Sigma$ is diagonal, so its gradient vector is
\[ \nabla_\uv\left((U\g)^T\uv+\frac12\uv^T\Sigma\uv\right) = [(U\g)_n+\sigma_n u_n]. \]

\subsection{Exact Stationary Point}
If none of the $\sigma_n$ are zero, $f$ has a stationary point at $\uv=[-(U\g)_n/\sigma_n]$, which is only a maximum if all the $\sigma_n$ are negative.

\subsection{Constrained Maximum}
A constrained maximum would have, for some $\lambda>0$,
\[ [(U\g)_n+\sigma_n u_n] = [\lambda u_n] \]
and thus $\uv=[(U\g)_n/(\lambda-\sigma_n)]$.  The value of $\lambda$ must satisfy
\[ r^2 = |\x|^2 = |\uv|^2 = \sum_n \frac{(U\g)_n^2}{(\lambda-\sigma_n)^2}. \]
The expression on the right has a $+\infty$ singularity whenever $\lambda=\sigma_n$ for some $n$.  It is also not monotonic, so there could be many solutions.  However, note that $\lambda\rightarrow\infty$ still corresponds to $r\rightarrow0$, so small $r$ solutions are in the region where $\lambda>\max_n \sigma_n = \sigma_\mathrm{max}$.

What does the other end of this region, $\lambda\rightarrow\sigma_\mathrm{max}^+$ correspond to?  First note that if $\sigma_\mathrm{max}<0$ then the other end is actually $\lambda\rightarrow0$, corresponding to the exact maximum (and it really is a maximum because all the $\sigma_n$ are negative).  Otherwise, a vector element $u_n$ with $\sigma_n=\sigma_\mathrm{max}\ge0$ tends to infinity, meaning the solution is asymptotically running up the steepest parabolic ascent direction available to it, as expected of a maximum.

Finally, note that although $r^2$ is not a monotonic function of $\lambda$, it is a (locally) convex one:
\[ \dds{r^2}{\lambda} = \sum_n \frac{6(U\g)_n^2}{(\lambda-\sigma_n)^4} \ge 0. \]
Taking into account the asymptotic behaviour as $\lambda\rightarrow\infty$, this means $r^2$ in the region $\lambda>\sigma_\mathrm{max}$ is monotonically decreasing, so a value of $\lambda$ can always be found for any value of $r$ between 0 and the norm of the exact solution point (or infinity if $\sigma_\mathrm{max}\ge0$, corresponding to a saddle, ridge or minimum valley).
% This argument isn't exactly watertight, could use reflection argument (inverting sign of any u_n gives the same value, so best to stay on the same side you started - no sign changes?)

\subsection{Summary}
The locus of constrained maxima is
\[ \x(\lambda) = U^T\left[\frac{(U\g)_n}{\lambda-\sigma_n}\right] \]
for $\lambda>\max\{0,\sigma_\mathrm{max}\}$.  If $\sigma_\mathrm{max}<0$ then $\x(0)$ is the exact maximum.

\end{document}

