Merge branch 'master' of tea.maggioni.xyz:maggicl/ICS

notes/lecture6.tex

@@ -0,0 +1,147 @@
+% vim: set ts=2 sw=2 et tw=80:
+
+\documentclass[12pt,a4paper]{article}
+
+\usepackage[utf8]{inputenc} \usepackage[margin=2cm]{geometry}
+\usepackage{amstext}
+\usepackage{amsmath}
+\usepackage{array}
+\usepackage[utf8]{inputenc}
+\usepackage[margin=2cm]{geometry}
+\usepackage{amstext}
+\usepackage{array}
+
+\newcommand{\lra}{\Leftrightarrow}
+\newcolumntype{L}{>{$}l<{$}}
+\DeclareMathOperator*{\argmax}{arg\,max}
+\DeclareMathOperator*{\argmin}{arg\,min}
+
+\title{Lecture notes 6 -- Introduction to Computational Science}
+\author{Micheal Multerer \\ Copied by: Claudio Maggioni}
+
+\begin{document}
+\maketitle
+
+\section*{2.5: Partial pivoting}
+Obviously, Algorithm 2.8 fils if one of the pivots becomes zero. In this case, we need to choose a different pivot element.
+
+\paragraph{Simple approach:} Column pivoting choose $|a^{(i)}_{k,i}|=\max_{i \leq l \leq n}|a^{(i)}_{l,i}|$
+In order to move the pivot element and the  corresponding row we switch row k and row i by a \textit{primitive matrix}:
+
+\[
+\bar{P}_i :=  \begin{bmatrix} 
+    1_1 \\
+    & \ddots \\
+    && 1_{i-1}  \\
+    &&& 0 & \hdots & 1_k \\    
+    &&&\vdots &1_{i+1} \\
+    &&& 1_i && 0 \\
+    &&&&&& \ddots \\
+    &&&&&&& 1_n  \\
+    \end{bmatrix}
+\]
+
+The following rules apply:
+\begin{enumerate}
+\item Multiplication by $P_i$ from the left $\Rightarrow$ interchange rows i and k
+\item Multiplication by $P_i$ from the right $\Rightarrow$ interchange columns i and k
+\item $(\bar{P}_i)^2 = \bar{P}_i \cdot \bar{P}_i = \bar{I}$
+\end{enumerate}
+
+Then performing an LU decompisition with pivoting can be
+written in matrix notation as:
+$$A_{i+1} = L_i \cdot \bar{P}_i \cdot A_i$$
+
+\paragraph{Note:} For $j < i$, it holds $\bar{P}_i \bar{L}_j = \widetilde{L}_j \bar{P}_i$ where $\widetilde{L}$ is the same matrix as $\bar{L}_j$ except that $[\widetilde{l}_j]_i$ and $[\widetilde{L}_j]_k$ are interchanged:
+
+\[
+\bar{L}_i :=  \begin{bmatrix} 
+    1\\
+    &\ddots\\
+    &&1\\
+    &&&\hat{l}^{(j)}_2\\
+    &&&\hat{l}^{(j)}_i&&1\\
+    &&&\hat{l}^{(j)}_k&&&\ddots\\
+    &&&\hat{l}^{(j)}_n&&&&1\\
+    \end{bmatrix}
+    \text{\hspace{1cm}}
+\widetilde{L}_i :=  \begin{bmatrix} 
+    1\\
+    &\ddots\\
+    &&1\\
+    &&&\hat{l}^{(j)}_2\\
+    &&&\hat{l}^{(j)}_k&&1\\
+    &&&\hat{l}^{(j)}_i&&&\ddots\\
+    &&&\hat{l}^{(j)}_n&&&&1\\
+    \end{bmatrix}
+\]
+
+Resolving (*) then yields:
+
+$$\bar{A}_{i+1} = \bar{L}_i \bar{P}_i \bar{L}_{i-1} \bar{P}_{i-1} \ldots \bar{L}_1 \bar{P}_1 A$$
+or
+$$\bar{L}_{n-1} \bar{P}_{n-1} \bar{L}_{n-2} \bar{P}_{n-2} \ldots \bar{L}_1 \bar{P}_1 \bar{A} = \bar{U}$$
+
+Now we can exploit $P_2 \bar{L}_1 \bar{P}_1 = \widetilde{L}_1 \bar{P}_2 \bar{P}_1$ and so on. This yields:
+
+$$\bar{P}\bar{A} = \bar{L}\bar{U}$$
+with
+$$\bar{P} = \bar{P}_{n-1} \bar{P}_{n-2} \ldots \bar{P}_1$$ 
+and
+$$\widetilde{L} = \widetilde{L}^{-1}_1 \ldots \widetilde{L}^{-1}_{n-1}$$
+and $$ \widetilde{L}_{n-1} = \bar{L}_{n-1} $$$$
+\widetilde{L}_{n-2} = \bar{P}_{n-1} \bar{L}_{n-2} \bar{P}_{n-1} $$$$
+\vdots $$$$
+\widetilde{L}_1 = \bar{P}_{n-1} \bar{P}_{n-2} \ldots \bar{P}_{2} \bar{L}_1 \bar{P}_2 \ldots \bar{P}_{n-2} \bar{P}_{n-1} $$
+
+\paragraph{Note:} If $\bar{A} \in R^{n \times n}$ is non-singular, the pivoted LU decomposition $\bar{P}\bar{A} = \bar{L}\bar{U}$ always exists.
+\\\\
+We can easily add column priority in Algorithm 2.8:
+
+\paragraph{Algorithm 2.9}(Outer product LU decomposition with column pivoting) \\
+input:  matrix $\bar{A} = [a_{i,j}]^n_{i,j = 1} \in R^{n\times n}$ \\
+output: pivoted LU decomposition $ \bar{L}\bar{U} = \bar{P}\bar{A}$
+
+\begin{enumerate}
+\item Set $\bar{A}_1 = \bar{A}, \bar{p} = [1,2,\ldots,n]$
+\item For $i = 1,2,\ldots,n$ \begin{itemize}
+\item compute: k = $\argmax_{1 \leq j \leq n} | a^{(i)}_{p_j,i} |$ \% find pivot
+\item swap: $ p_i \gets\to p_k$
+\item $\bar{l}_i := \bar{a}^{(i)}_{:,i} / a^{(i)}_{p_i,p_i}$
+\item $\bar{u}_i := a_{p_i,:}$
+\item compute: $\bar{A}_{i+1} = \bar{A}_i - \bar{l}_i \cdot \bar{u}_i$
+\end{itemize}
+\item set $\bar{P} := [\bar{l}_{p_1},\bar{l}_{p_2},\ldots, \bar{l}_{p_n}]^T$ \% $\bar{e}_i$ is i-th unit vector
+\item set $bar{L} = bar{P}[\bar{l}_1,\bar{l}_2,\cdots,\bar{l}_n]$
+\end{enumerate}
+
+\paragraph{Example 2.10} \textit{(omitted)}
+
+\section*{2.6: Cholesky decomposition}
+If $\bar{A}$ is symmetric and positive definite, i.e. all eigenvalues of $\bar{A}$ are bigger than zero or equivalently $bar{x}^T\bar{A}\bar{x} > 0$ for all $\bar{x} \neq 0$, we can compute a symmetric decomposition of $\bar{A}$.
+
+\paragraph{Note:} if $\bar{A}$ is symmetric and positive definite, then the \textit{Schur complement} $\bar{S}:=\bar{A}_{2:n,2:n} - (\bar{a}_{2:n},1/a_{1,1})\bar{a}^T_{2:n,1}$, is symmetric and positive definite as well. In particular, it holds $s_{i,i} > 0$ and $a_{i,i} > 0$ !
+ 
+\paragraph{Definition 2.11} A decomposition $\bar{A} = \bar{L}\bar{L}^T$ with a lower triangular matrix $\bar{L}$ with positive diagonal elements is called \textit{Cholesky decomposition of $\bar{A}$}.
+ 
+\paragraph{Note:} A Cholesky decomposition exists, if $\bar{A}$ is symmetric and positive definite.
+
+\paragraph{Algorithm 2.12} (outer product of Cholesky decomposition) \\
+input: matrix $\bar{A}$ symmetric and positive definite \\
+output: Cholesky decomposition $\bar{A} = \bar{L}\bar{L}^T =
+[\bar{l}_1,\bar{l}_2,\ldots,\bar{l}_n] [\bar{l}_1,\bar{l}_2,\ldots,\bar{l}_n]^T$
+
+\begin{enumerate}
+\item set: $\bar{A}_1 := \bar{A}$
+\item for $i = 1,2,\ldots,n$ \begin{itemize}
+\item set: $\bar{l}_i := a^{(i)}_{:,i} / \sqrt{a^{(i)}_{i,i}}$
+\item set: $\bar{A}_{i+1} := \bar{A}_{i} - \bar{l}_i\bar{l}^T_i$ 
+\end{itemize}
+\item set: $\bar{L} = [\bar{l}_1,\bar{l}_2,\ldots,\bar{l}_n]$
+\end{enumerate}
+
+The computational cost is $\frac{1}{6}n^{3}+O(n^{2})$ and thus only half the cost of LU decomposition.
+
+\paragraph{Example 2.13} \textit{(omitted)}
+
+\end{document}