multiphysicsL6.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL6}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%\usepackage{kbordermatrix}
%
%\beginArtNoToc
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 6: Singular value decomposition, and conditioning number.  Taught by Prof.\ Piero Triverio}
%\chapter{Matrix norm, singular decomposition, and conditioning number}
\label{chap:multiphysicsL6}
%
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
\section{Singular value decomposition.}
\index{singular value decomposition}
%\section{Matrix norm}
\index{matrix norm}
%
Recall that the matrix norm of \( M \), for the system \( \By = M \Bx \) was defined as
%
\begin{equation}\label{eqn:multiphysicsL6:21}
\Norm{M} = \max_{\Norm{\Bx} = 1} \Norm{ M \Bx }.
\end{equation}
%
The \( L_2 \) norm will typically be used, resulting in a matrix norm of
%
\begin{equation}\label{eqn:multiphysicsL6:41}
\Norm{M}_2 = \max_{\Norm{\Bx}_2 = 1} \Norm{ M \Bx }_2.
\end{equation}
%
It can be shown that
%
\begin{equation}\label{eqn:multiphysicsL6:61}
\Norm{M}_2 = \max_i \sigma_i(M),
\end{equation}
%
where \( \sigma_i(M) \) are the (SVD) singular values.
%
\index{SVD|see {singular value decomposition}}
\index{singular value decomposition}
\makedefinition{Singular value decomposition (SVD).}{dfn:multiphysicsL6:1}{
        Given \( M \in \text{\R{m \times n}} \), a factoring of \( M \) can be found with the form
%
\begin{equation}\label{eqn:multiphysicsL6:81}
M = U \Sigma V^\T,
\end{equation}
%
where \( U \) and \( V\) are orthogonal matrices such that \( U^\T U = 1 \), and \( V^\T V = 1 \), and
%
\begin{equation}\label{eqn:multiphysicsL6:101}
\Sigma =
\begin{bmatrix}
  \sigma_1 &          &        &          &         &         &\\
           & \sigma_2 &        &          &         &         &\\
           &          & \ddots &          &         &         &\\
           &          &        & \sigma_r &         &         &\\
           &          &          &        &       0 &         &         \\
           &          &          &        &         & \ddots  &         \\
           &          &          &        &         &         & 0       \\
\end{bmatrix}
\end{equation}
%
The values \( \sigma_i, i \le \min(n,m) \) are called the singular values of \( M \).  The singular values are subject to the ordering
%
\begin{equation}\label{eqn:multiphysicsL6:261}
\sigma_{1} \ge \sigma_{2} \ge \cdots \ge 0
\end{equation}
%
If \(r\) is the rank of \( M \), then the \( \sigma_r \) above is the minimum non-zero singular value (but the zeros are also called singular values).
%
%The Matlab command to compute this is \underline{svd}.
% O(4 n^3) div O(26 n^3) ???
}
%
Observe that the condition \( U^\T U = 1 \) is a statement of orthonormality.  In terms of column vectors \( \Bu_i \), such a product written out explicitly is
%
\begin{equation}\label{eqn:multiphysicsL6:301}
\begin{bmatrix}
\Bu_1^\T \\ \Bu_2^\T \\ \vdots \\ \Bu_m^\T
\end{bmatrix}
\begin{bmatrix}
\Bu_1 & \Bu_2 & \cdots & \Bu_m
\end{bmatrix}
=
\begin{bmatrix}
1 &   & & \\
  & 1 & & \\
  &   & \ddots & \\
  &   &        & 1
\end{bmatrix}.
\end{equation}
%
This is both normality \( \Bu_i^\T \Bu_i = 1 \), and orthonormality \( \Bu_i^\T \Bu_j = 1, i \ne j \).
%
\makeexample{A \( 2 \times 2 \) case.}{example:multiphysicsL6:1}{
(for column vectors \( \Bu_i, \Bv_j \)).
%
\begin{equation}\label{eqn:multiphysicsL6:281}
M =
\begin{bmatrix}
\Bu_1 & \Bu_2
\end{bmatrix}
\begin{bmatrix}
        \sigma_1 &  \\
        & \sigma_2
\end{bmatrix}
\begin{bmatrix}
\Bv_1^\T \\
\Bv_2^\T
\end{bmatrix}
\end{equation}
%
Consider \( \By = M \Bx \), and take an \( \Bx \) with \( \Norm{\Bx}_2 = 1 \)
%
Note: I've chosen not to sketch what was drawn on the board in class.  See instead the animated gif of the same in \citep{wiki:svd}, or the Wolfram online SVD demo in \citep{wolframdemo:svd}.
}
%
A very nice video treatment of SVD by Prof Gilbert Strang can be found in \citep{ocw:svd}.
%
\section{Conditioning number.}
\index{conditioning number}
%
Given a perturbation of \( M \Bx = \Bb \) to
%
\begin{equation}\label{eqn:multiphysicsL6:121}
\lr{ M + \delta M }
        \lr{ \Bx + \delta \Bx } = \Bb,
\end{equation}
%
or
\begin{equation}\label{eqn:multiphysicsL6:141}
\cancel{ M \Bx - \Bb } + \delta M \Bx + M \delta \Bx + \delta M \delta \Bx = 0.
\end{equation}
%
This gives
%
\begin{equation}\label{eqn:multiphysicsL6:161}
M \delta \Bx = - \delta M \Bx - \delta M \delta \Bx,
\end{equation}
%
or
\begin{equation}\label{eqn:multiphysicsL6:181}
\delta \Bx = - M^{-1} \delta M \lr{ \Bx + \delta \Bx }.
\end{equation}
%
Taking norms
%
\begin{equation}\label{eqn:multiphysicsL6:201}
\Norm{ \delta \Bx}_2 = \Norm{
        M^{-1} \delta M \lr{ \Bx + \delta \Bx }
}_2
\le
\Norm{ M^{-1} }_2 \Norm{ \delta M }_2 \Norm{ \Bx + \delta \Bx }_2,
\end{equation}
%
or
\begin{equation}\label{eqn:multiphysicsL6:221}
\mathLabelBox
[
   labelstyle={yshift=1.2em},
   linestyle={}
]
{
\frac{\Norm{ \delta \Bx}_2 }{
\Norm{ \Bx + \delta \Bx}_2 }
}
{relative error of solution}
\le
\mathLabelBox
[
   labelstyle={below of=m\themathLableNode, below of=m\themathLableNode}
]
{
\Norm{M}_2
\Norm{M^{-1}}_2
}
{conditioning number of \(M\)}
\mathLabelBox
[
   labelstyle={xshift=2cm},
   linestyle={out=270,in=90, latex-}
]
{
\frac{ \Norm{ \delta M}_2 }
{\Norm{M}_2}
}
{
relative perturbation of \( M \)
}.
\end{equation}
%
The \textAndIndex{conditioning number} can be shown to be
%
\begin{equation}\label{eqn:multiphysicsL6:241}
        \cond(M) =
        \frac
        {\sigma_{\mathrm{max}}}
        {\sigma_{\mathrm{min}}}
        \ge 1
\end{equation}
%
FIXME: justify.
%
\paragraph{sensitivity to conditioning number}
\index{conditioning number}
%
Double precision relative rounding errors can be of the order \( 10^{-16} \sim 2^{-54} \), which allows the relative error of the solution to be gauged
%
\captionedTable{Relative errors.}{tab:multiphysicsL6:1}{
\begin{tabular}{|llll|}
	\hline
	relative error of solution & \( \le \) & \( \cond(M) \) & \( \frac{ \Norm{ \delta M } }{\Norm{M} } \) \\
	\hline
	\( 10^{-15} \)             & \( \le \) & 10             & \( \sim 10^{-16} \) \\
	\( 10^{-2} \)              & \( \le \) & \( 10^{14} \)  & \( 10^{-16} \) \\
	\hline
\end{tabular}
}
%
%\EndArticle
%\EndNoBibArticle