% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Systems of linear equations}

\section{Introduction}
A {\em system of linear equations} (or {\em linear system}) is a
collection of linear equations involving the same set of unknown
variables.  A general system of $n$ linear equations with $m$ unknowns
can be written as
	\begin{equation}
\left\{\begin{array}{ccccccccc} A_{11}x_{1} & + & A_{12}x_{2} & + &
\cdots & + & A_{1m}x_{m} & = & b_{1} \\ A_{21}x_{1} & + & A_{22}x_{2} & +
& \cdots & + & A_{2m}x_{m} & = & b_{2} \\ \vdots      &   & \vdots	&
&	 &   & \vdots	   &   &\vdots \\ A_{n1}x_{1} & + & A_{n2}x_{2}
& + & \cdots & + & A_{nm}x_{m} & = & b_{n} \end{array}\right.\;,
	\end{equation}
%or, equivalently, as
%	\begin{equation}
%\sum_{j=1}^{m}A_{ij}x_j=b_i\;\big|_{i=1\dots n} \;.
%	\end{equation}
where $x_1,x_2,\dots,x_m$ are the unknown variables,
$A_{11},A_{12},\dots,A_{nm}$ are the (constant) coefficients of the
system, and $b_1,b_2,\dots,b_n$ are the (constant) right-hand side terms.

The system can be equivalently written in the matrix form,
	\begin{equation}
\ma{A}\ve{x}=\ve{b}\;,
	\end{equation}
where $\ma{A}\doteq\{A_{ij}\}$ is the $n\times m$ matrix
of the coefficients, $\ve{x}\doteq\{x_{j}\}$
is the size-$m$ column-vector of the unknown variables, and
$\ve{b}\doteq\{b_{i}\}$ is the size-$n$ column-vector of
right-hand side terms.

A solution to a linear system is a set of values for the variables
$\ve{x}$ which satisfies all equations.

Systems of linear equations occur quite regularly in applied mathematics.
Therefore computational algorithms for finding solutions of linear
systems are an important part of numerical methods.

A system of non-linear equations can often be approximated by a linear
system -- a helpful technique (called {\em linearization}) in creating a
mathematical model of an otherwise a more complex system.

If $m=n$, the matrix $A$ is called {\em square}. A square system has a
unique solution if $A$ is invertible.

\section{Triangular systems}
An efficient algorithm to solve a square system of linear equations
numerically is to transform the original system into an equivalent {\em
triangular system}\index{triangular system},
	\begin{equation}
\ma{T}\ve{y}=\ve{c}\;,
	\end{equation}
where $\ma{T}$ is a {\em triangular matrix} -- a special kind of square
matrix where the matrix elements either below (lower triangular) or above
(upper triangular) the main diagonal are zero.

Indeed, an upper triangular system $\ma{U}\ve{y}=\ve{c}$ can be easily
solved by {\em back-substitution}\index{back-substitution},
	\begin{equation}
y_{i} = {1\over U_{ii}}\left(c_{i} - \sum_{k=i+1}^{n} U_{ik}y_{k}\right),
\; i=n,n-1,\ldots,1\; ,
	\end{equation}
where one first computes $y_n=b_n/U_{nn}$, then substitutes {\em back} into the
previous equation to solve for $y_{n-1}$, and repeats through $y_1$.

Here is a C-function implementing in-place\footnote{here
{\em in-place} means the right-hand side $\ve{c}$ is replaced by the
solution $\ve{y}$.} back-substitution\footnote{ the functions {\small\tt
vector\_get}, {\small\tt vector\_set}, and {\small\tt matrix\_get}
are assumed to implement fetching and setting the vector and matrix
elements.}:
	\lstinputlisting{c/qr/backsub.c}

For a lower triangular system $\ma{L}\ve{y}=\ve{c}$
the equivalent procedure is called {\em
forward-substitution}\index{forward-substitution},
	\begin{equation}
y_{i} = {1\over L_{ii}}\left(c_{i} - \sum_{k=1}^{i-1} U_{ik}y_{k}\right),
\; i=1,2,\ldots,n\; .
	\end{equation}

\section{Reduction to triangular form}
Popular algorithms for reducing a square system of linear equations to
a triangular form are {\em LU-decomposition} and {\em QR-decomposition}.

\subsection{QR-decomposition}
QR-decomposition\index{QR-decomposition} is a factorization of a
matrix into a product of an orthogonal matrix $\ma{Q}$, such that
$\ma{Q}^\T\ma{Q}=1$, where $^\T$ denotes transposition, and a right
triangular matrix $\ma{R}$,
	\begin{equation}
\ma{A}=\ma{Q}\ma{R} \;.
	\end{equation}
QR-decomposition can be used to convert a linear system $\ma{A}\ve{x}=\ve{b}$
into the triangular form (by multiplying with $\ma{Q}^\T$ from the left),
	\begin{equation}
\ma{R}\ve{x}=\ma{Q}^\T\ve{b} \;,
	\end{equation}
which can be solved directly by back-substitution.

QR-decomposition can also be performed on non-square matrices with few
long columns. Generally speaking a rectangular $n\times m$ matrix $\ma{A}$
can be represented as a product, $\ma{A}=\ma{Q}\ma{R}$, of an orthogonal
$n\times m$ matrix $\ma{Q}$, $\ma{Q}^\T\ma{Q}=1$, and a right-triangular
$m\times m$ matrix $\ma{R}$.

QR-decomposition of a matrix can be computed using several methods,
such as Gram-Schmidt orthogonalization, Householder
transformation~\cite{householder}\index{Householder transformation},
or Givens rotation~\cite{givens}.

\subsubsection{Gram-Schmidt orthogonalization}

{\em Gram-Schmidt orthogonalization} is an algorithm for orthogonalization
of a set of vectors in a given inner product space. It takes a linearly
independent set of vectors $\ma{A}=\{\ve{a}_1,\dots,\ve{a}_m\}$
and generates an orthogonal set $\ma{Q}=\{\ve{q}_1,\dots,\ve{q}_m\}$
which spans the same subspace as $\ma{A}$. The algorithm is given as
	\begin{code}
for $i=1$ to $m$ : \\
\hspace*{1em}$\ve{q}_i \leftarrow {\ve{a}_i}/{\|\ve{a}_i\|}$ \\
\hspace*{1em}for $j=i+1$ to $m$ :
%\hspace*{2em}
$\ve{a}_j \leftarrow \ve{a}_j-\DOT{\ve{q}_i}{\ve{a}_j}\ve{q}_i$
	\end{code}
where $\DOT{\ve{a}}{\ve{b}}$ is the inner product of two vectors,
and $\|\ve{a}\|\doteq\sqrt{\DOT{\ve{a}}{\ve{a}}}$ is the vector's norm.
This variant of the algorithm, where all remaining vectors $\ve{a}_j$
are made orthogonal to $\ve{q}_i$ as soon as the latter is calculated, is
considered to be numerically stable and is referred to as {\em stabilized}
or {\em modified}.

Stabilized Gram-Schmidt orthogonalization can be used to compute
QR-decomposition of a matrix $\ma{A}$ by orthogonalization of its
column-vectors $\ve{a}_i$ with the inner product
	\begin{equation}
\DOT{\ve{a}}{\ve{b}}=\ve{a}^\T\ve{b}\equiv
\sum_{k=1}^{n}(\ve{a})_k(\ve{b})_k\;,
	\end{equation}
where $n$ is the length of column-vectors $\ve{a}$ and $\ve{b}$, and
$(\ve{a})_k$ is the $k$th element of the column-vector,
	\begin{code}
for $i=1$ to $m$ : \\
\hspace*{1em}$R_{ii}=\sqrt{{\ve{a}_i^\T\ve{a}_i}}\;$;
	$\;\ve{q}_i=\ve{a}_i/R_{ii}$ \\
\hspace*{1em}for $j=i+1$ to $m$ : \\
\hspace*{2em}$R_{ij}=\ve{q}_i^\T\ve{a}_j\;$;
	$\;\ve{a}_j=\ve{a}_j-\ve{q}_iR_{ij} \;.$ \end{code}
After orthogonalization the matrices $\ma{Q}=\{\ve{q}_1\dots \ve{q}_m\}$ and
$\ma{R}$ are the sought orthogonal and right-triangular factors of
matrix $\ma{A}$.

The factorization is unique under requirement that the diagonal elements
of $R$ are positive. For a $n\times m$ matrix the complexity of
the algorithm is $O(m^2n)$.

Table~(\ref{lin:t:qr-gs}) shows a C implementation of the stabilized
Gram-Schmidt orthogonalization.  The corresponding function to solve
the equation
	\begin{equation}\label{lin:e:qrxb}
\ma{Q}\ma{R}\ve{x}=\ve{b}
	\end{equation}
are shown in Table~(\ref{lin:t:qrxb}).

	\begin{table}
\caption{QR-decomposition in C using stabilized Gram-Schmidt
orthogonalization.}
	\label{lin:t:qr-gs}
	\lstinputlisting{c/qr/qrdec.c}
	\end{table}

	\begin{table}
\caption{C-functions to perform QR-backsubstitution of
equation~(\ref{lin:e:qrxb})}
	\label{lin:t:qrxb}
	\lstinputlisting{c/qr/qrback.c}
	\end{table}

\subsubsection{Householder transformation}
A square matrix $\ma{H}$ of the form
	\begin{equation}\label{}
\ma{H} = 1 - \frac{2}{\ve{u}^\T \ve{u}}\ve{u}\ve{u}^\T
	\end{equation}
is called {\em Householder matrix}, where the vector $\ve{u}$ is
called a {\em Householder vector}. Householder matrices are symmetric
and orthogonal,
	\begin{equation}\label{}
\ma{H}^\T = \ma{H}\:,\;\ma{H}^\T \ma{H} = 1\:.
	\end{equation}
The transformation induced by the Householder matrix on a given vector
$a$,
	\begin{equation}\label{}
\ve{a} \rightarrow \ma{H}\ve{a} \:,
	\end{equation}
is called a {\em Householder transformation}\index{Householder
transformation} or {\em Householder reflection}\index{Householder
reflection}.  The transformation changes the sign of the affected
vector's component in the $\ve{u}$ direction, or, in other words, makes a
reflection of the vector about the hyperplane perpendicular to $\ve{u}$,
hence the name.

Householder transformation can be used to zero selected components of
a given vector $\ve{a}$. For example, one can zero all components but the
first one, such that 
	\begin{equation}\label{}
\ma{H}\ve{a} = \gamma \ve{e}_1 \:,
	\end{equation}
where $\gamma$ is a number and $\ve{e}_1$ is the unit vector
in the first direction. The factor $\gamma$ can be easily calculated,
	\begin{equation}\label{}
\|\ve{a}\|^2 \doteq \ve{a}^\T \ve{a} = \ve{a}^\T \ma{H}^\T \ma{H}\ve{a}
= (\gamma \ve{e}_1)^\T (\gamma \ve{e}_1)=\gamma^2 \:,
	\end{equation}
	\begin{equation}\label{}
 \Rightarrow\; \gamma=\pm\|a\| \:.
	\end{equation}
To find the Householder vector, we notice that
	\begin{equation}\label{}
\ve{a} = \ma{H}^\T \ma{H}\ve{a} = \ma{H}^\T \gamma \ve{e}_1
= \gamma \ve{e}_1 - \frac{2(\ve{u})_1}{\ve{u}^\T \ve{u}}\ve{u} \;,
	\end{equation}\begin{equation}\label{}
\Rightarrow\; \frac{2(\ve{u})_1}{\ve{u}^\T \ve{u}}\ve{u}
= \gamma \ve{e}^1 - \ve{a} \:,\label{eq:h:t}
	\end{equation}
where $(\ve{u})_1$ is the first component of the vector $\ve{u}$. One
usually chooses $(\ve{u})_1=1$ (for the sake of the possibility to store
the other components of the Householder vector in the zeroed elements
of the vector $\ve{a}$) and stores the factor
	\begin{equation}\label{}
\frac{2}{\ve{u}^\T \ve{u}} \equiv \tau
	\end{equation}
separately. With this convention one readily finds $\tau$ from the first
component of equation~(\ref{eq:h:t}),
	\begin{equation}\label{}
\tau = \gamma - (\ve{a})_1 \:.
	\end{equation}
where $(\ve{a})_1$ is the first element of the vector $\ve{a}$. For
the sake of numerical stability the sign of $\gamma$ has to be chosen
opposite to the sign of $(\ve{a})_1$,
	\begin{equation}\label{}
\gamma = -\mathrm{sign}\left((\ve{a})_1\right)\|\ve{a}\| \:.
	\end{equation}
Finally, the Householder reflection, which zeroes all component of a
vector $\ve{a}$ but the first, is given as
	\begin{equation}\label{eq:hh:hh}
\ma{H} = 1 - \tau \ve{u} \ve{u}^\T  \:,\;
\tau = -\mathrm{sign}((\ve{a})_1)\|\ve{a}\| - (\ve{a})_1 \:,\;
(\ve{u})_1 = 1 \:,\; (\ve{u})_{i>1} = -\frac{1}{\tau}(\ve{a})_i \:.\;
	\end{equation}

Now, a QR-decomposition of an $n\times n$ matrix $\ma{A}$ by Householder
transformations\index{Householder transformation} can be
performed in the following way:
	\begin{enumerate}

\item Build the size-$n$ Householder vector
$\ve{u}_1$ which zeroes the sub-diagonal elements of the
first column of matrix $\ma{A}$, such that
	\begin{equation}\label{li:qr:h1}
\ma{H}_1\ma{A} = \left[\begin{array}{c|ccc}
\star &\star&\dots&\star\\
\hline
0 & & & \\
\vdots & & \ma{A}_1 & \\
0 & & & \end{array}\right] \;,
	\end{equation}
where $\ma{H}_1 = 1 - \tau_1 \ve{u}_1 \ve{u}_1^\T$ and where $\star$
denotes (generally) non-zero matrix elements. In practice one does
not build the matrix $\ma{H}_1$ explicitly, but rather calculates the
matrix $\ma{H}_1\ma{A}$ in-place, consecutively applying the Householder
reflection to columns the matrix $\ma{A}$, thus avoiding computationally
expensive matrix-matrix operations. The zeroed sub-diagonal elements of
the first column of the matrix $\ma{A}$ can be used to store the elements
of the Householder vector $\ve{u}_1$ while the factor $\tau_1$ has to be
stored separately in a special array. This is the storage scheme used
by LAPACK and GSL.

\item Similarly, build the size-$(n-1)$ Householder vector $\ve{u}_2$
which zeroes the sub-diagonal elements of the first column of matrix
$\ma{A}_1$ from eq.~(\ref{li:qr:h1}). With the transformation
matrix $\ma{H}_2$ defined as
	\begin{equation}\label{}
\ma{H}_2 = \left[\begin{array}{c|ccc}
1 & 0 &\cdots& 0 \\
\hline 0 & & & \\
\vdots & & 1-\tau_2 \ve{u}_2 \ve{u}_2^\T & \\
0 & & & \end{array}\right] \;.
	\end{equation}
the two transformations together zero the sub-diagonal elements of the
two first columns of matrix $\ma{A}$,
	\begin{equation}\label{}
\ma{H}_2\ma{H}_1\ma{A} = \left[\begin{array}{cc|ccc}
\star & \star & \star & \cdots & \star \\
0 & \star & \star & \cdots & \star \\
\hline
0 & 0 &   &        &   \\
\vdots & \vdots &  & \ma{A}_3 & \\
0 & 0 &   &        &
\end{array}\right] \:,
	\end{equation}

\item Repeating the process zero the sub-diagonal elements of the
remaining columns. For column $k$ the corresponding Householder matrix
is
	\begin{equation}\label{}
\renewcommand*\arraystretch{2}
\ma{H}_k = \left[\begin{array}{c|c}
\ma{I}_{k-1} & 0 \\
\hline
0 &  1-\tau_k \ve{u}_k \ve{u}_k^\T \\
\end{array}\right] \;,
	\end{equation}
where $\ma{I}_{k-1}$ is a unity matrix of size $k-1$, $\ve{u}_k$ is the
size-(n-k+1) Householder vector that zeroes the sub-diagonal elements
of matrix $\ma{A}_{k-1}$ from the previous step. The corresponding
transformation step is
	\begin{equation}\label{}
\ma{H}_{k}\dots\ma{H}_2\ma{H}_1\ma{A} = \left[\begin{array}{c|c}
\ma{R}_k & \star \\
\hline
0  & \ma{A}_k
\end{array}\right] \:,
	\end{equation}
where $\ma{R}_k$ is a size-$k$ right-triangular matrix.

After $n-1$ steps the matrix $\ma{A}$ will be transformed into a right
triangular matrix,
	\begin{equation}\label{li:qr:ha7}
\ma{H}_{n-1}\cdots\ma{H}_2\ma{H}_1\ma{A} = \ma{R} \;.
	\end{equation}

\item Finally, introducing an orthogonal matrix
$\ma{Q}=\ma{H}_1^\T\ma{H}_2^\T\dots\ma{H}_{n-1}^\T$ and multiplying
eq.~(\ref{li:qr:ha7}) by $\ma{Q}$ from the left, we get the sought
QR-decomposition,
	\begin{equation}\label{}
\ma{A}=\ma{Q}\ma{R} \;.
	\end{equation}

In practice one does not explicitly builds the $\ma{Q}$ matrix but rather
applies the successive Householder reflections stored during the
decomposition.

	\end{enumerate}

%Indeed, setting $v=a+\beta e^1$, where $\beta$
%is a constant, gives
%	\begin{equation}\label{}
%v^\T a = a^\T a + \beta a_1 \:,
%	\end{equation}
%and
%	\begin{equation}\label{}
%v^\T v = a^\T a + 2\beta a_1 + \beta^2 \:.
%	\end{equation}
%Therefore
%	\begin{equation}\label{}
%Ha = \left( 1 - 2\frac{a^\T a+\beta a_1}{a^\T a + 2\beta a_1 + \beta^2}
%\right)a-2\beta\frac{v^\T a}{v^\T v}e^1 \;.
%	\end{equation}
%The condition that the expression in parenthesis is zero has two solutions
%	\begin{equation}\label{}
%\beta=\pm\sqrt{a^\T a} \:,
%	\end{equation}
%which gives
%	\begin{equation}\label{}
%v^\T a = a^\T a \pm\sqrt{a^\T a}a_1\:,\; v^\T v = 2a^\T a \pm
%2\sqrt{a^\T a}a_1\:,
%	\end{equation}
%and
%	\begin{equation}\label{}
%Ha = \mp \sqrt{a^\T a}e_1 \:.
%	\end{equation}
%
%To avoid round-off errors the sign of $\beta$ should be chosen equal to
%the sign of $a_1$. In practice it is useful to normalize the Householder
%vector such that $v_1 = 1$. That allows the compact storage of the
%Householder vector in place where the zeros in $a$ have been indtroduced.

%\begin{algorithmic}
%\STATE $\sigma = \sqrt{\sum_{i=1}^{n}(a_i)^2}$
%\STATE $v=\{a_1+\mathrm{sign}(a_1)\sigma,a_2,\dots,a_n\}^\T $
%\STATE $\tau = \frac{(v_1)^2}{\sigma^2+\sigma|a_1|} $
%\STATE $v/=v_1$
%\end{algorithmic}

\subsubsection{Givens rotations}
A Givens rotation is a transformation in the form
	\begin{equation}
\ma{A} \rightarrow \ma{G}(p,q,\theta)\ma{A} \:,
	\end{equation}
where $\ma{A}$ is the transformed object---matrix of vector---and
$\ma{G}(p,q,\theta)$ is the Givens rotation matrix (also known as Jacobi
rotation matrix): an orthogonal matrix in the form
	\begin{equation}
G(p, q, \theta) =
	\begin{bmatrix}
1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots &
& \vdots & & \vdots \\ 0 & \cdots & \cos\theta & \cdots & \sin\theta &
\cdots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots
& -\sin\theta & \cdots & \cos\theta & \cdots & 0 \\ \vdots & & \vdots &
& \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1
	\end{bmatrix} \begin{array}{c} \\ \\ \leftarrow \mathrm{row}\;
  p \\ \\
\leftarrow \mathrm{row}\; q \\ \\ \\
	\end{array}\:.	\end{equation}
When a Givens rotation matrix $\ma{G}(p,q,\theta)$ multiplies a vector,
only elements with indices $p$ and $q$ are affected. Restricting our
attention to these two elements, say $a$ and $b$, the Givens rotation
is given explicitely as
	\begin{equation} \begin{bmatrix}
\cos\theta & \sin\theta \\ -\sin\theta & \cos\theta
	\end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix}
= \begin{bmatrix} a\cos\theta+b\sin\theta \\ -a\sin\theta+b\cos\theta
	\end{bmatrix} \:.  \end{equation}
Apparently, the rotation can zero the lowest element, if the angle
$\theta$ is chosen as
	\begin{equation}
\tan\theta = \frac{b}{a} \;\;\Rightarrow\;\;\theta=\mathrm{atan2}(b,a) \:.
	\end{equation}

A sequence of Givens rotations can zero all elements of a matrix below
the main diagonal.  Which obviously amounts to QR-decomposition of the
matrix with $\ma{Q}$ being the product of all the applied Givens matrices.

Since each Givens rotation only affects two rows of the matrix it is
possible to apply a set of rotations in parallel.  Givens rotations are also
more efficient on sparse matrices.

\subsection{LU-decomposition}
LU-decomposition\index{LU-decomposition} is a factorization of a square
matrix $\ma{A}$ into a product of a lower triangular matrix $\ma{L}$
and an upper triangular matrix $\ma{U}$,
	\begin{equation}\label{le:eq:alu}
\ma{A}=\ma{L}\ma{U}\;.
	\end{equation}

The linear system $\ma{A}\ve{x}=\ve{b}$ after LU-decomposition of
the matrix $\ma{A}$ becomes $\ma{L}\ma{U}\ve{x}=\ve{b}$ and can be solved
by first solving $\ma{L}\ve{y}=\ve{b}$ for $\ve{y}$ and then
$\ma{U}\ve{x}=\ve{y}$ for $\ve{x}$ with two runs of forward and
backward substitutions.

If $\ma{A}$ is an $n\times n$ matrix, the condition (\ref{le:eq:alu}) is a
set of $n^2$ equations,
	\begin{equation}\label{le:eq:lua}
\sum_{k=1}^{n}L_{ik}U_{kj}=A_{ij} \;\big|_{i,j=1\dots n} \;,
	\end{equation}
for $n^2+n$ unknown elements of the triangular matrices $\ma{L}$ and
$\ma{U}$. The decomposition is thus not unique.

Usually the decomposition is made unique by providing extra $n$ conditions
e.g. by the requirement that the elements of the main diagonal of
the matrix $\ma{L}$ are equal one,
	\begin{equation}\label{le:eq:lii}
L_{ii}=1\:,\;i=1\dots n \:.
	\end{equation}

The system (\ref{le:eq:lua}) with the extra conditions (\ref{le:eq:lii})
can then be easily solved row after row using the {\em Doolittle's
algorithm},
	\begin{code}
for $i=1 \dots n$ :\\
\hspace*{1em} $L_{ii}=1$ \\
\hspace*{1em} for $j=i\dots n$ : $U_{ij}=A_{ij}-\sum_{k<i}L_{ik}U_{kj}$\\
\hspace*{1em} for $j=i+1\dots n$ :
$L_{ji}=\frac{1}{U_{ii}}\left(A_{ji}-\sum_{k<j}L_{jk}U_{ki}\right)$
	\end{code}

In a slightly different {\em Crout's algorithm} it is the matrix $\ma{U}$
that has unit diagonal elements,
	\begin{code}
for $i=1 \dots n$ :\\
\hspace*{1em} $U_{ii}=1$ \\
\hspace*{1em} for $j=i\dots n$ : $L_{ji}=A_{ji}-\sum_{k<i}L_{jk}U_{ki}$\\
\hspace*{1em} for $j=i+1\dots n$ :
$U_{ij}=\frac{1}{L_{ii}}\left(A_{ji}-\sum_{k<j}L_{jk}U_{ki}\right)$
	\end{code}

Without a proper ordering (permutations) in the matrix, the factorization
may fail. For example, it is easy to verify that $A_{11}=L_{11}U_{11}$. If
$A_{11}=0$, then at least one of $L_{11}$ and $U_{11}$ has to be zero,
which implies either $\ma{L}$ or $\ma{U}$ is singular, which is impossible
if $\ma{A}$ is non-singular. This is however only a procedural problem. It
can be removed by simply reordering the rows of $\ma{A}$ so that the
first element of the permuted matrix is nonzero (or, even better, the
largest in absolute value among all elements of the column below the
diagonal).  The same problem in subsequent factorization steps can be
removed in a similar way.  Such algorithm is refered to as {\em partial
pivoting}\index{partial pivoting}.  It requires an extra integer array
to keep track of row permutations.

\section{Determinant of a matrix}
LU- and QR-decompositions allow $O(n^3)$ calculation of the
determinant of a square matrix. Indeed, for the LU-decomposition,
	\begin{equation}
\det{\ma{A}}=\det \ma{L}\ma{U} = \det{\ma{L}}\det{\ma{U}}
=\det \ma{U} = \prod_{i=1}^n U_{ii}\;.
	\end{equation}
For the QR-decomposition
	\begin{equation}
\det \ma{A} =\det \ma{Q}\ma{R} =\det \ma{Q} \det \ma{R} \;.
	\end{equation}
Since $\ma{Q}$ is an orthogonal matrix $(\det \ma{Q})^2=1$ and therefore
	\begin{equation}
|\det \ma{A}|=|\det \ma{R}|=\left|\prod_{i=1}^n R_{ii}\right|\;.
	\end{equation}

\section{Matrix inverse}
The inverse $\ma{A}^{-1}$ of a square $n\times n$ matrix $\ma{A}$ can
be calculated by solving $n$ linear equations
	\begin{equation}
\ma{A}\ve{x}_i=\ve{e}_i \:\Big|_{i=1 \dots n}\;, 
	\end{equation}
where $\ve{e}_i$ is a column where all elements are equal zero except
for the element number $i$, which is equal one; that is, columns $\ve{e}_i$
form a unity matrix.  The matrix made of columns $\ve{x}_i$ is apparently
the inverse of $\ma{A}$.

A C implementation of this algorithm using QR-decomposition is shown
in Table~(\ref{lin:t:i}).

	\begin{table}[htb]
\caption{Matrix inverse with QR-decomposition in C}
\label{lin:t:i}
\lstinputlisting{c/qr/qrinverse.c}
   \end{table}
