chapterAppend.tex

\appendix
\rohead{\headmark}
\lehead{\headmark}

\definecolor{prussianblue}{rgb}{0.0, 0.19, 0.33}

\renewcommand{\chapterformat}{\raggedleft \colorbox{prussianblue}{%
\centering\textit{\textcolor{white}{{\Large Appendix} {\Huge \thechapter}}}}}

\chapter{Supplementary Information for the Main Text}

\section{Chapter \ref*{chapter:invdet}}
\label{section:invdetappend}

\paragraph{Properties \ref*{proper:zerodet}}
Consider an $n \times n$ square matrix $A$ that has two identical and adjacent rows with indices $i_1$ and $i_2$, where $i_2 = i_1 + 1$ (hence one of the indices is odd and another is even), then cofactor expansion along the odd row (let's say $i_1$) will give
\begin{subequations}
\begin{align}
\abs{A} &= \sum_{k=1}^{n} A_{i_1k}C_{i_1k} \nonumber \\
&= \sum_{k=1}^{n} (-1)^{i_1+k}A_{i_1k}\det(M_{i_1k})
\end{align}
by Properties \ref{proper:cofactorex} and Definition \ref{defn:cofactor}. Similarly, by considering the even row, we have
\begin{align}
\abs{A} &= \sum_{k=1}^{n} A_{i_2k}C_{i_2k} \nonumber \\
&= \sum_{k=1}^{n} (-1)^{i_2+k}A_{i_2k}\det(M_{i_2k})
\end{align}
\end{subequations}
But since the $i_1$-th and $i_2$-th row are identical, $A_{i_1k} = A_{i_2k}$. Furthermore, as these two identical rows are also adjacent, the minors $M_{i_1k} = M_{i_2k}$ are also equal. The only difference between the two expressions for $|A|$ above is the $(-1)^{i+j}$ factor. And because $i_1$ is odd and $i_2$ is even, they are differed by a negative sign only. Explicitly, we have
\begin{align}
\abs{A} &= \sum_{k=1}^{n} (-1)^{i_1+k}A_{i_1k}\det(M_{i_1k}) \nonumber \\
&= \sum_{k=1}^{n} (-1)^{(i_2-1)+k}A_{i_2k}\det(M_{i_2k}) \nonumber \\
&= (-1) \sum_{k=1}^{n} (-1)^{i_2+k}A_{i_2k}\det(M_{i_2k}) \nonumber \\
&= -\abs{A}
\end{align}
Therefore $\abs{A} = -\abs{A}$ and $\abs{A} = 0$ must equal to zero. Now we can generalize the results to non-adjacent, proportional rows by doing the first and third kind (multiplication and swapping) of elementary row operations when appropriate with the aid of Properties \ref{proper:elementaryopdet}, and the second case in Properties \ref{proper:zerodet} is completed. Subsequently, for the addition/subtraction type of elementary row operations, let's say $R_{p} + cR_{q} \to R_{p}$ is applied on some matrix $A$ (this $A$ is not the same one as in the first part) to produce $A'$, then
\begin{align*}
A' = 
\begin{bmatrix}
\vdots & \vdots & & \vdots\\
A_{p1} + cA_{q1} & A_{p2} + cA_{q2} & \cdots & A_{pn} + cA_{qn} \\
\vdots & \vdots & & \vdots\\
A_{q1} & A_{q2} & \cdots & A_{qn} \\
\vdots & \vdots & & \vdots
\end{bmatrix}
\end{align*}
where we have only written out the rows $R_p$ and $R_q$. By applying cofactor expansion along $R_p$ following Properties \ref{proper:cofactorex}, we have
\begin{align}
\abs{A'} &= \sum_{k=1}^{n} [(A_{pk} + cA_{qk}) C_{pk}] \nonumber \\
&= \sum_{k=1}^{n} A_{pk}C_{pk} + c\sum_{k=1}^{n} A_{qk}C_{pk}
\end{align}
We identify the first term with $\abs{A}$ that is computed using cofactor expansion on the row $R_p$ of $A$. The second term can be thought of as the determinant of a matrix $\tilde{A}$ that is formed by replacing the $p$-th row $R_p$ by the $q$-th row $R_q$ in $A$ and subsequently expanding along that row. So $\tilde{A}$ practically has two identical rows $R_p = R_q$, and by the previous result, the value of $\abs{\tilde{A}}$ is zero. Therefore
\begin{align}
\abs{A'} = \abs{A} + c\abs{\tilde{A}} = \abs{A} + c(0) = \abs{A}    
\end{align} implying that the addition/subtraction type of elementary row operations does not affect the value of the determinant.

\section{Chapter \ref*{chap:vec_space}}
\label{section:vecspaceappend}

\paragraph{Properties \ref*{proper:subspace_n_span}}  

$\mathcal{W}$ as a subspace of $\mathbb{R}^n$ contains at most $n$ linearly independent vectors. Assume $\mathcal{W}$ is non-empty and take any non-zero vector in $\mathcal{W}$, denote it by $\vec{v}^{(1)}$. The span of $\mathcal{\beta} = \{\vec{v}^{(1)}\}$ is itself a subspace $\mathcal{B} = \text{span}(\mathcal{\beta})$ of the subspace $\mathcal{W}$ by noting that it is a subset of $\mathcal{W}$ and using the first part of Properties \ref{proper:subspace_n_span}. If this subspace $\mathcal{B}$ of the subspace $\mathcal{W}$ is exactly $\mathcal{W}$, in the sense that no other $\vec{v} \in \mathcal{W}$ is not included by $\mathcal{B} = \text{span}(\mathcal{\beta})$, then we are done because this subspace is a span by construction. Otherwise, take another non-zero vector $\vec{v}^{(2)} \in \mathcal{W}$ which is linearly independent of $\vec{v}^{(1)}$ and add it to $\mathcal{\beta}$. Now $\mathcal{B} = \text{span}(\mathcal{\beta}) = \text{span}(\{\vec{v}^{(1)}, \vec{v}^{(2)}\})$ is enlarged by one dimension but still a subspace of $\mathcal{W}$ and we can check if it coincides with $\mathcal{W}$, otherwise, repeat the procedure with $\vec{v}^{(3)}, \vec{v}^{(4)}, \ldots$ until $\text{span}(\mathcal{B}) = \mathcal{W}$ then we can stop. Remember we can at most add up to $\vec{v}^{(n)}$ since $\mathcal{W}$ contains at most $n$ linearly independent vectors, so in the middle of somewhere we would done with $\vec{v}^{(p)}$, where $p \leq n$, and hence we know that $\mathcal{W}$ is some span (and a $p$-dimensional subspace).

\paragraph{Dependence relations}

We will show this for the case of row addition/subtraction only because the other two types of elementary row operations are easy to check. Without loss of generality, take a dependence relation in the form of $\vec{v}^{(r+1)} = c_1\vec{v}^{(1)} + c_2\vec{v}^{(2)} + \cdots + c_r\vec{v}^{(r)}$ where $A=[\vec{v}^{(1)}|\vec{v}^{(2)}|\cdots|\vec{v}^{(r)}|\vec{v}^{(r+1)}]$, i.e. \begin{align*}
\begin{bmatrix}
\vdots \\
\vec{v}_p^{(r+1)} \\
\vdots \\
\vec{v}_q^{(r+1)} \\
\vdots
\end{bmatrix}
= c_1
\begin{bmatrix}
\vdots \\
\vec{v}_p^{(1)} \\
\vdots \\
\vec{v}_q^{(1)} \\
\vdots
\end{bmatrix}
+ c_2
\begin{bmatrix}
\vdots \\
\vec{v}_p^{(2)} \\
\vdots \\
\vec{v}_q^{(2)} \\
\vdots
\end{bmatrix} + \cdots + 
c_r
\begin{bmatrix}
\vdots \\
\vec{v}_p^{(r)} \\
\vdots \\
\vec{v}_q^{(r)} \\
\vdots
\end{bmatrix}
\end{align*}
The elementary row operation of adding $c_q$ times row $R_q$ to row $R_p$ then produces a new matrix $A'$ with column vectors
\begin{align*}
\vec{v}^{(j)'}
=
\begin{bmatrix}
\vdots \\
\vec{v}_p^{(j)} + c_q\vec{v}_q^{(j)} \\
\vdots \\
\vec{v}_q^{(j)} \\
\vdots
\end{bmatrix}
\end{align*}
for all $j$. Therefore,
\begin{align*}
\vec{v}^{(r+1)'} &= \begin{bmatrix}
\vdots \\
\vec{v}_p^{(r+1)} + c_q\vec{v}_q^{(r+1)} \\
\vdots \\
\vec{v}_q^{(r+1)} \\
\vdots
\end{bmatrix} \\
&=
\begin{bmatrix}
\vdots \\
(c_1\vec{v}_p^{(1)} + c_2\vec{v}_p^{(2)} + \cdots + c_r\vec{v}_p^{(r)}) + c_q(c_1\vec{v}_q^{(1)} + c_2\vec{v}_q^{(2)} + \cdots + c_r\vec{v}_q^{(r)}) \\
\vdots \\
(c_1\vec{v}_q^{(1)} + c_2\vec{v}_q^{(2)} + \cdots + c_r\vec{v}_q^{(r)}) \\
\vdots
\end{bmatrix}  
\end{align*}
by the original dependence relation, which is then equal to
\begin{align*}
&\quad \begin{bmatrix}
\vdots \\
c_1(\vec{v}_p^{(1)} + c_q\vec{v}_q^{(1)}) \\
\vdots \\
c_1\vec{v}_q^{(1)} \\
\vdots
\end{bmatrix}
+
\begin{bmatrix}
\vdots \\
c_2(\vec{v}_p^{(2)} + c_q\vec{v}_q^{(2)}) \\
\vdots \\
c_2\vec{v}_q^{(2)} \\
\vdots
\end{bmatrix} 
+ \cdots +
\begin{bmatrix}
\vdots \\
c_r(\vec{v}_p^{(r)} + c_q\vec{v}_q^{(r)}) \\
\vdots \\
c_r\vec{v}_q^{(r)} \\
\vdots
\end{bmatrix} \\
&= c_1\vec{v}^{(1)'} + c_2\vec{v}^{(2)'} + \cdots + c_r\vec{v}^{(r)'}
\end{align*}
This shows that the same dependence relation holds between the new column vectors $\vec{v}^{(j)'}$, $j = 1,2,\ldots,r+1$.

\paragraph{Steinitz Replacement Theorem}

\begin{thm}[Steinitz Replacement Theorem]\index{Steinitz Replacement Theorem}
\label{thm:Steinitz}
%G = \gamma
%S = \eta
%H = \gamma'
Given a vector space $\mathcal{V}$ generated by a set $\mathcal{\gamma}$ consisting of $n$ vectors (not necessarily linearly independent), and another set $\mathcal{\eta}$ containing $m$ linearly independent vectors from $\mathcal{V}$. Then $m \leq n$ and there exists a subset $\mathcal{\gamma'}$ of $\mathcal{\gamma}$ with exactly $n-m$ vectors such that $\mathcal{\eta} \cup \mathcal{\gamma'}$ spans $\mathcal{V}$.
\end{thm}
\begin{proof}
We proceed with mathematical induction on $m$. The base case $m = 0$ is trivial as $\mathcal{\eta} = \varnothing$ and $\mathcal{\gamma'} = \mathcal{\gamma}$. Now assume that the statement is true for some integer $m = k \geq 0$ and we have to prove it for $m = k+1$. Let $\mathcal{\eta} = \{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}, \vec{v}^{(k+1)}\}$ be a subset of $\mathcal{V}$ with $m = k+1$ linearly independent vectors. It is apparent that after removing $\vec{v}^{(k+1)}$ from $\mathcal{\eta}$, $\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}\}$ is still linearly independent. Then we can use the induction hypothesis to obtain $k \leq n$ and a subset $\{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}$ of $\mathcal{\gamma}$ so that $\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}\} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}$ generates $\mathcal{V}$. So $\vec{v}^{(k+1)} \in \mathcal{V}$ can be written as the linear combination of
\begin{align*}
\vec{v}^{(k+1)} = a_1\vec{v}^{(1)} + a_2\vec{v}^{(2)} + \cdots + a_k\vec{v}^{(k)} + b_1\vec{u}^{(1)} + b_2\vec{u}^{(2)} + \cdots + b_{n-k}\vec{u}^{(n-k)} 
\end{align*}
It must be true that $n - k \geq 1$, so that some $b_j$ are present for otherwise $\vec{v}^{(k+1)}$ will be reduced to a linear combination of $\vec{v}^{(j)}$ and contradicts the assumption that $\mathcal{\eta}$ is linearly independent. Hence $m = k + 1 \leq n$. By part (b) of Theorem \ref{thm:plusminus}, $\text{span}(\mathcal{\eta} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) = \text{span}(\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}, \vec{v}^{(k+1)}\} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) = \text{span}(\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}\} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) \\ = \mathcal{V}$ (the last equality follows from the induction hypothesis). For the same reason, in particular, some $b_j$ has to be non-zero, let's say $b_1$. Therefore we can write
\begin{align*}
\vec{u}_1 = -\frac{a_1}{b_1}\vec{v}^{(1)} -\frac{a_2}{b_1}\vec{v}^{(2)} - \cdots - \frac{a_k}{b_1}\vec{v}^{(k)} + \frac{1}{b_1}\vec{v}^{(k+1)} - \frac{b_2}{b_1}\vec{u}^{(2)} - \cdots - \frac{b_{n-k}}{b_1}\vec{u}^{(n-k)}
\end{align*}
as a linear combination of $\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}, \vec{v}^{(k+1)}$ and $\vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}$. Let $\mathcal{\gamma}' = \{\vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}$. Again by part (b) of Theorem \ref{thm:plusminus}
\begin{align*}
&\quad \text{span}(\mathcal{\eta} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) \\
&= \text{span}(\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}, \vec{v}^{(k+1)}\} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) \\
&= \text{span}(\{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(k)}, \vec{v}^{(k+1)}\} \cup \{\vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) \\
&= \text{span}(\mathcal{\eta} \cup \mathcal{\gamma'})
\end{align*}
But we just have $\text{span}(\mathcal{\eta} \cup \{\vec{u}^{(1)}, \vec{u}^{(2)}, \ldots, \vec{u}^{(n-k)}\}) = \mathcal{V}$ from above, which readily shows that $\text{span}(\mathcal{\eta} \cup \mathcal{\gamma'}) = \mathcal{V}$. As $\mathcal{\gamma'}$ is a subset of $\mathcal{\gamma}$ with $n - k-1 = n - (k+1)$ vectors, the theorem is true for $m = k+1$, and the induction is completed.
\end{proof}
The key point of the theorem is that some vectors in a generating set of a vector space can be replaced by the same number of linearly independent vectors from that vector space, hence comes the name of the replacement theorem.

\section{Chapter \ref*{chap:DFT}}
\label{section:DFTappend}

\paragraph{Bessel's Inequality}
By part (d) of Theorem \ref{thm:spectralinner}, we can expand a function $f$ as
\begin{align}
f = \lim_{p \to \infty} \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} 
\end{align}
where $\varphi^{(j)}$ are the orthonormal basis vectors/functions and the partial sum $S_p[f]$ will be just $\sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}$. Particularly, the completeness of the orthonormal basis in Properties \ref{proper:hilbertorthosys} means that (see the footnote below)\footnotemark
\begin{align}
\norm{f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}^2 = 0 \label{eqn:fcomplete}
\end{align}
On the other hand,\footnotetext{Otherwise, if $\smash{\norm{\smash{f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}}^2 > 0}$, then consider $\smash{\tilde{\varphi} = f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}$ which is now a non-zero vector. It will be orthogonal to any of the original basis vectors $\smash{\varphi^{(j')}}$ as 
\begin{align*}
\langle \tilde{\varphi}, \varphi^{(j')} \rangle &= \langle f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} , \varphi^{(j')} \rangle \\
&= \langle f, \varphi^{(j')} \rangle - \langle \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} , \varphi^{(j')} \rangle \\
&= \langle f, \varphi^{(j')} \rangle - (\cdots + (0) +\langle f, \varphi^{(j')} \rangle (1) + (0) + \cdots) \\
& \quad \text{(Orthonormality of the basis: $\langle \varphi^{(j)}, \varphi^{(j')} \rangle = 1$ } \\ 
& \quad \text{only when $j = j'$ and $0$ if $j \neq j'$)} \\
&=\langle f, \varphi^{(j')} \rangle - \langle f, \varphi^{(j')} \rangle = 0
\end{align*}
which violates the premise of completeness as it will become another vector that is linearly independent of all the other basis vectors and can be added to the basis.}
\begin{align} 
\norm{f - S_p[f]}^2 &= \norm{f - \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}^2 \nonumber \\
&= \left\langle f - \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}, f - \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \right\rangle \nonumber \\
&= \langle f, f \rangle - \left\langle f, \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \right\rangle - \left\langle \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}, f \right\rangle \nonumber \\
&\quad + \left\langle \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}, \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \right\rangle \nonumber \\
&= \norm{f}^2 - \left\langle f, \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \right\rangle \nonumber \\
&\quad - \overline{\left\langle f, \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \right\rangle} \quad \text{((1) of Definition \ref{defn:innerprod})} \nonumber \\
&\quad + \sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 \quad \text{(Orthonormality of the basis)} \nonumber \\
&= \norm{f}^2 - 2\Re{\langle f, \sum_{j=1}^{p} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} \rangle} + \sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 \nonumber \\
&= \norm{f}^2 - 2\Re{\sum_{j=1}^{p} (\overline{\langle f, \varphi^{(j)} \rangle} \langle f, \varphi^{(j)} \rangle)} \nonumber \\
& \quad \text{((4) of Properties \ref{proper:innerprod2})}\nonumber \\
&\quad + \sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 \nonumber \\
&= \norm{f}^2 - 2\sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 + \sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 \nonumber \\
&= \norm{f}^2 - \sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2
\end{align}
and since $\norm{f - S_p[f]}^2 \geq 0$ is always non-negative, we have the \textit{Bessel's Inequality} as
\begin{align}
\sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2 \leq \norm{f}^2  \label{eqn:Bessel}
\end{align}
From Exercise \ref{ex:triangular2}, we can apply the Triangle Inequality to get
\begin{align}
&\quad \norm{\left(f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}\right) + \left(\sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} - S_p[f]\right)}^2 \nonumber \\
&= \norm{f - S_p[f]}^2 \leq \norm{f - \sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}^2 + \norm{\sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} - S_p[f]}^2
\end{align}
By Equation (\ref{eqn:fcomplete}), the first term on R.H.S. of the equality is zero. And the second term
\begin{align}
&\quad \norm{\sum_{j=1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)} - S_p[f]}^2 \nonumber \\
&= \norm{\sum_{j=p+1}^{\infty} \langle f, \varphi^{(j)} \rangle \varphi^{(j)}}^2 = \sum_{j=p+1}^{\infty} \abs{\langle f, \varphi^{(j)} \rangle}^2
\end{align}
is a remainder term and when $p \to \infty$, must tend to zero, because $\sum_{j=1}^{p} \abs{\langle f, \varphi^{(j)} \rangle}^2$ is a convergent sequence, seen by applying the \textit{Monotone Convergence Theorem} from elementary Analysis on Bessel's Inequality (\ref{eqn:Bessel}). Therefore, the L.H.S.
\begin{align}
\lim_{p \to \infty}\norm{f - S_p[f]}^2 = 0
\end{align}
also tends to zero as we take the same limit, which implies the $L^2$ convergence for a complete orthonormal basis.

\section{Chapter \ref*{chapter:Markov}}
\label{section:Markovappend}

\paragraph{Properties \ref*{proper:positivestoceig}}

\textit{Reference Materials: \cite{markov}}

Let $\vec{q}$ be any eigenvector of $A^T$ that has an eigenvalue of modulus $\abs{\lambda} = 1$. Then $A^T\vec{q} = \lambda \vec{q}$ by the definition of an eigenvalue-eigenvector problem, i.e.\
\begin{align}
\sum_{k=1}^{n} a_{ki}\vec{q}_k = \lambda \vec{q}_i
\end{align}
where $a_{ji}$ is the $(j,i)$ [$(i,j)$] entry of the matrix $A$ [$A^T$]. Then
\begin{align}
\abs{\lambda \vec{q}_i} = \abs{\lambda}\abs{\vec{q}_i} &= \abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k} \nonumber \\
\abs{\vec{q}_i} &\leq \sum_{k=1}^{n} \abs{a_{ki}}\abs{\vec{q}_k} = \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k} \quad (\abs{\lambda} = 1)\\
&\text{(Triangle Inequality for complex numbers and $a_{ij}$ are real)} \nonumber
\end{align}
Now assume that for some fixed $i$, $\abs{\vec{q}_i} \geq \abs{\vec{q}_k}$ for $1 \leq k \leq n$. Then the inequality above becomes
\begin{align}
\abs{\vec{q}_i} \leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k} &\leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_i} = \abs{\vec{q}_i}  & \text{($\sum_{k=1}^{n} a_{ki}$ sums to $1$)}
\end{align}
By squeezing using $\abs{\vec{q}_i}$ at the both ends, the part $\sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k} \leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_i}$ forces that $\abs{\vec{q}_k} = \abs{\vec{q}_i}$ for all $k$ and it also holds for any $i$. This is where the positiveness of $A$ is needed, otherwise $a_{ki}$ can be $0$ and the squeeze becomes vacuous ($0\abs{\vec{q}_k} = 0\abs{\vec{q}_i} = 0$). Moreover, as $\abs{\vec{q}_i} = \abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k}$, incorporating it into the inequality
\begin{align}
\abs{\vec{q}_i} = \abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k} \leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k} &\leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_i} = \abs{\vec{q}_i}
\end{align}
and apply squeezing again, the part of Triangle Inequality $\abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k} \leq \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k}$ becomes an equality $\abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k} = \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k}$, which means that all the components $\vec{q}_k$ have to lie along the same direction in the complex plane.\footnote{For any two complex number $z_1$ and $z_2$, if $\abs{z_1 + z_2} = \abs{z_1} + \abs{z_2}$, then
\begin{align*}
\abs{z_1 + z_2}^2 &= (z_1 + z_2) \overline{(z_1 + z_2)} \\
&= z_1 \overline{z_1} + z_1 \overline{z_2} + z_2 \overline{z_1} + z_2 \overline{z_2} \\
&= \abs{z_1}^2 + z_1 \overline{z_2} + z_2 \overline{z_1} + \abs{z_2}^2
\end{align*}
but also
$(\abs{z_1} + \abs{z_2})^2 = \abs{z_1}^2 + 2\abs{z_1}\abs{z_2} + \abs{z_2}^2$. Hence we have
\begin{align*}
z_1 \overline{z_2} + z_2 \overline{z_1} = 2\abs{z_1}\abs{z_2}
\end{align*}
Assume that $z_1$ and $z_2$ points in some directions so that $z_1 = \abs{z_1}e^{i\theta_1}$ and $z_2 = \abs{z_2}e^{i\theta_2}$, then
\begin{align*}
z_1 \overline{z_2} + z_2 \overline{z_1} &= \abs{z_1}e^{i\theta_1} \abs{z_2}e^{-i\theta_2} + \abs{z_2}e^{i\theta_2} \abs{z_1}e^{-i\theta_1} \\
&= \abs{z_1}\abs{z_2} e^{i(\theta_1-\theta_2)} + \abs{z_1}\abs{z_2} e^{-i(\theta_1-\theta_2)} \\
&= 2 \abs{z_1}\abs{z_2} (\frac{e^{i(\theta_1-\theta_2)} + e^{-i(\theta_1-\theta_2)}}{2}) \\
&= 2 \abs{z_1}\abs{z_2} \cos(\theta_1-\theta_2) & \text{(Properties \ref{proper:sincoscomplex})}
\end{align*}
but this has to equal to $2\abs{z_1}\abs{z_2}$, thus $\cos(\theta_1-\theta_2) = 1$ and the arguments $\theta_1$ and $\theta_2$ will be the same. Now apply this logic repetitively to $\abs{\sum_{k=1}^{n} a_{ki}\vec{q}_k} = \sum_{k=1}^{n} a_{ki}\abs{\vec{q}_k}$.} Therefore, $\vec{q}$ must be be in the form of $c(1,1,1,\ldots,1)^T$ where $c$ will be any complex constant and this shows that over $\mathbb{C}$, the only linearly independent eigenvector for $A^T$ corresponding to $\abs{\lambda} = 1$ will be $(1,1,1,\ldots,1)^T$, matching the derivation in Properties \ref{proper:markoveigen1} where the eigenvalue is exactly $\lambda = 1$ and the geometric multiplicity is hence $1$ as well. Again, by Properties \ref{proper:eigentransinv}, this result will then also hold for the original matrix $A$ (but the form of the eigenvector will be different).

\section{Chapter \ref*{chapter:Tensor}}
\label{section:tensorappend}

\paragraph{Uniqueness of rank-$2$ and $3$ Isotropic Tensors} The uniqueness of the Kronecker delta as the only isotropic rank-$2$ tensor can be shown as follows. By a \SI{90}{\degree} positive rotation about the $3$-axis, we have
\begin{subequations}
\begin{align}
\textbf{e}'_1 &= \textbf{e}_2 \\
\textbf{e}'_2 &= -\textbf{e}_1 \\
\textbf{e}'_3 &= \textbf{e}_3
\end{align}    
\end{subequations}
and according to (\ref{eqn:aij}), the change of coordinates matrix is
\begin{align}
a_{ij} = 
\begin{bmatrix}
\textbf{e}_1 \cdot \textbf{e}'_1 & \textbf{e}_1 \cdot \textbf{e}'_2 & \textbf{e}_1 \cdot \textbf{e}'_3 \\
\textbf{e}_2 \cdot \textbf{e}'_1 & \textbf{e}_2 \cdot \textbf{e}'_2 & \textbf{e}_2 \cdot \textbf{e}'_3 \\
\textbf{e}_3 \cdot \textbf{e}'_1 & \textbf{e}_3 \cdot \textbf{e}'_2 & \textbf{e}_3 \cdot \textbf{e}'_3 
\end{bmatrix} = 
\begin{bmatrix}
0 & -1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{bmatrix}
\label{eqn:rotate903ax}
\end{align}
For a general rank-$2$ tensor $F$, by (\ref{eqn:rank2aF}) its transformation will then be
\begin{align}
\begin{bmatrix}
F'_{11} & F'_{12} & F'_{13} \\
F'_{21} & F'_{22} & F'_{23} \\
F'_{31} & F'_{32} & F'_{33}  
\end{bmatrix} &= 
\begin{bmatrix}
0 & 1 & 0 \\
-1 & 0 & 0 \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
F_{11} & F_{12} & F_{13} \\
F_{21} & F_{22} & F_{23} \\
F_{31} & F_{32} & F_{33}  
\end{bmatrix} 
\begin{bmatrix}
0 & -1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{bmatrix} \nonumber \\
&= \begin{bmatrix}
F_{22} & -F_{21} & F_{23} \\
-F_{12} & F_{11} & -F_{13} \\
F_{32} & -F_{31} & F_{33}
\end{bmatrix}
\end{align}
If the tensor is isotropic such that $F'_{ij} = F_{ij}$, then by comparing the entries on both sides, we must have
\begin{subequations}
\begin{align}
F_{11} &= F_{22} \\
F_{13} = F_{23} = -F_{13} \implies F_{13} &= F_{23} = 0 \\
F_{31} = F_{32} = -F_{31} \implies F_{31} &= F_{32} = 0 
\end{align}
\end{subequations}
By the same technique, we can carry out a \SI{90}{\degree} positive rotation about the $2$-axis to obtain
\begin{subequations}
\begin{align}
F_{11} &= F_{33} \\
F_{12} &= F_{32} = 0 \\
F_{21} &= F_{23} = 0 
\end{align}
\end{subequations}
Therefore, $F_{11} = F_{22} = F_{33}$ and $F_{12} = F_{21} = F_{13} = F_{31} = F_{23} = F_{32} = 0$, the diagonal entries have the same value and all off-diagonal entries are zero, and any isotropic rank-$2$ tensor must be in the form of $\lambda \delta_{ij}$, a constant multiple of the Kronecker delta symbol. We can show the same for the rank-$3$ isotropic tensor to be found. First, do a rotation about the $(1,1,1)^T$ axis so that the effect is a subscript permutation $1 \to 2 \to 3 \to 1$, then we have
\begin{subequations}
\begin{align}
T_{111} &= T_{222} = T_{333} \\
T_{112} &= T_{223} = T_{331} \\
T_{122} &= T_{233} = T_{311} \\
T_{212} &= T_{323} = T_{131} \\
T_{211} &= T_{322} = T_{133} \\
T_{121} &= T_{232} = T_{313} \\
T_{221} &= T_{332} = T_{113} \\
T_{123} &= T_{231} = T_{312} \\
T_{132} &= T_{321} = T_{213}
\end{align}
\end{subequations}
Next, apply a \SI{90}{\degree} positive rotation about the $3$-axis as suggested by (\ref{eqn:rotate903ax}) before. After some tedious algebra, we can obtain the following relations:
\begin{subequations}
\begin{align}
T_{111} = -T_{222} &\implies T_{111} = T_{222} = T_{333} = 0 \\
T_{112} = T_{221} \quad \text{and}\quad T_{221} = -T_{112} &\implies T_{112} = T_{221} = 0 \\
T_{122} = T_{211} \quad \text{and}\quad  T_{211} = -T_{122} &\implies T_{122} = T_{211} = 0 \\
T_{121} = T_{212} \quad \text{and}\quad  T_{212} = -T_{121} &\implies T_{121} = T_{212} = 0 \\
T_{123} = -T_{213} &\implies \begin{aligned}
T_{123} &= T_{231} = T_{312}  \\
= -T_{213} &= -T_{321} = -T_{132}    
\end{aligned}
\end{align}
The first four of them show that the rank-$3$ isotropic tensor must take a value of zero when there are repeated subscripts. The last relation implies that it will have an opposite sign when the subscripts belong to an odd/even permutation respectively. This coincides with Definition \ref{defn:epsilon} for the epsilon symbol and thus any isotropic rank-$3$ tensor must be in the form of $\lambda\epsilon_{ijk}$, a constant multiple of the epsilon symbol.
\end{subequations}

\chapter{Cyclic Subspaces and Jordan Normal Forms}

Note that this chapter is basically a rewritten summary of Friedberg's Section 5.4 and Chapter 7. 

\section{From Cyclic Subspaces to Cayley-Hamilton Theorem}
\label{section:B1}

A \index{Cyclic Subspace}\keywordhl{cyclic subspace} of a vector space is generated by repeatedly applying the same linear operator on a fixed initial vector and taking the span of all the resulting vectors.
\begin{defn}[Cyclic Subspace]
Given a vector space $\mathcal{V}$ with a linear operator $T: \mathcal{V} \to \mathcal{V}$ on it, and choose a non-zero vector $\vec{v} \in \mathcal{V}$, then the subspace
\begin{align}
\mathcal{W} = \text{span}(\{\vec{v}, T(\vec{v}), T^2(\vec{v}), \ldots\})
\end{align}
is known as the \textit{$T$-cyclic subspace of $\mathcal{V}$ generated by $\vec{v}$}. 
\end{defn}
It is obviously $T$-invariant (Definition \ref{defn:invariantsub}) because for any vector $\vec{w} \in \mathcal{W}$, if we write $\vec{w} = c_0\vec{v} + c_1T(\vec{v}) + c_2T^2(\vec{v}) + \cdots$ as a span (Definition \ref{defn:span}, extended to an infinite one if needed), then $T(\vec{w}) = c_0T(\vec{v}) + c_1T^2(\vec{v}) + c_2T^3(\vec{v}) + \cdots \in \mathcal{W}$ lies in the $T$-cyclic subspace too. It is also the smallest $T$-invariant subspace that contains $\vec{v}$: For other $T$-invariant subspaces $\mathcal{W}'$ that also contain $\vec{v}$, by definition they contain $T(\vec{v}), T^2(\vec{v}), \ldots$ as well, and by Properties \ref{proper:WcontainsspanS} they have to completely contain the corresponding $T$-cyclic subspace $\mathcal{W} = \text{span}(\{\vec{v}, T(\vec{v}), T^2\vec{v}, \ldots\})$ too, i.e.\ $\mathcal{W} \subseteq \mathcal{W}'$.\par

To proceed further, we need some more intermediate results. First,
\begin{proper}
\label{proper:B12}
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a \textit{finite-dimensional} vector space $\mathcal{V}$, and $\mathcal{W}$ be a $T$-invariant subspace of $\mathcal{V}$. Then, the characteristic polynomial (Definition \ref{defn:charactereqn}) of the restriction of $T$ to $\mathcal{W}$, that is, $T|_W$ (Definition \ref{defn:restrictionTW}), divides the characteristic polynomial of $T$.
\end{proper}
We refer the readers to Friedberg for the proof. It means that if the characteristic polynomial of $T|_W$ is $g(t)$ and that of $T$ is $f(t)$, then there exists some other polynomial $q(t)$ such that $f(t) = q(t)g(t)$. Next,
\begin{proper}
\label{proper:B13}
Given $T: \mathcal{V} \to \mathcal{V}$ as a linear operator on a \textit{finite-dimensional} vector space $\mathcal{V}$ and let $\mathcal{W}$ be a $T$-cyclic subspace of $\mathcal{V}$ generated by a non-zero vector $\vec{v} \in \mathcal{V}$. If $\dim(\mathcal{W}) = k$, then
\begin{enumerate}[label=(\alph*)]
    \item $\{\vec{v}, T(\vec{v}), T^2(\vec{v}), \ldots, T^{k-1}(\vec{v})\}$ is a basis for $\mathcal{W}$, and
    \item If $a_0 \vec{v} + a_1 T(\vec{v}) + \cdots + a_{k-1}T^{k-1}(\vec{v}) + T^k(\vec{v}) = \textbf{0}$, then the characteristic polynomial of $T|_W$ is $g(t) = (-1)^k (a_0 + a_1t + \cdots + a_{k-1}t^{k-1} + t^k)$.
\end{enumerate}
\end{proper}
Again, we refer the readers to Friedberg for the proof of part (a) as it is explained in detail there. Meanwhile, to show (b), denote the basis in part (a) as $\mathcal{\beta} = \{\vec{v}, T(\vec{v}), T^2(\vec{v}), \ldots, T^{k-1}(\vec{v})\}$. Then we have $T(T^j(\vec{v})) = T^{j+1}(\vec{v})$ for $0 \leq j < k-1$, and $T(T^{k-1}(\vec{v})) = T^k(\vec{v}) = -a_0 \vec{v} - a_1 T(\vec{v}) - \cdots - a_{k-1}T^{k-1}(\vec{v})$. Subsequently, the linear operator $T|_W$ will have a matrix representation of
\begin{align}
[T|_W]_\beta = 
\begin{bmatrix}
0 & 0 & 0 & \cdots & 0 & -a_0 \\
1 & 0 & 0 & \cdots & 0 & -a_1 \\
0 & 1 & 0 & \cdots & 0 & -a_2 \\
\vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 1 & -a_{k-1}
\end{bmatrix}
\end{align}
This is called a \index{Companion Matrix}\textit{companion matrix}. Its characteristic polynomial (in $t$) is then directly obtained by repeated cofactor expansion (Properties \ref{proper:cofactorex}) along the rows from top to bottom:
\begin{align}
&\quad \det([T|_W]_\beta - t I) \nonumber \\
&= \begin{vmatrix}
-t & 0 & 0 & \cdots & 0 & -a_0 \\
1 & -t & 0 & \cdots & 0 & -a_1 \\
0 & 1 & -t & \cdots & 0 & -a_2 \\
\vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 1 & -a_{k-1}-t
\end{vmatrix} \nonumber \\
&=
-t
\begin{vmatrix}
-t & 0 & \cdots & 0 & -a_1 \\
1 & -t & \cdots & 0 & -a_2 \\
\vdots & \vdots  & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & -a_{k-1}-t
\end{vmatrix}
+ (-a_0) (-1)^{k+1}
\begin{vmatrix}
1 & -t & 0 & \cdots & 0 \\
0 & 1 & -t & \cdots & 0 \\
\vdots & \vdots & \vdots & & \vdots \\
0 & 0 & 0 & \cdots & 1 
\end{vmatrix} \nonumber \\
&= -t
\begin{vmatrix}
-t & 0 & \cdots & 0 & -a_1 \\
1 & -t & \cdots & 0 & -a_2 \\
\vdots & \vdots  & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & -a_{k-1}-t
\end{vmatrix}
+ a_0 (-1)^k (1) \nonumber \\
&\quad 
\begin{aligned}
&\text{(The determinant of an upper-triangular block is simply} \\    
&\text{the product of its diagonal entries, which are all $1$ here)}
\end{aligned} \nonumber \\
&= (-t)(-t)
\begin{vmatrix}
-t & \cdots & 0 & -a_2 \\
\vdots & & \vdots & \vdots \\
0 & \cdots & 1 & -a_{k-1}-t
\end{vmatrix} + (-t)(-a_1)(-1)^{(k+1)-1} 
\begin{vmatrix}
1 & -t & \cdots & 0 \\
\vdots & \vdots  & & \vdots \\
0 & 0 & \cdots & 1 
\end{vmatrix} \nonumber \\
&\quad + a_0 (-1)^k \nonumber \\
&= t^2 \begin{vmatrix}
-t & \cdots & 0 & -a_2 \\
\vdots & & \vdots & \vdots \\
0 & \cdots & 1 & -a_{k-1}-t
\end{vmatrix}
+ (a_0 + a_1t)(-1)^k \nonumber \\
& \quad \vdots \nonumber\\
&= (-t)^{k-1}(-a_{k-1}-t) + (a_0 + a_1t + \cdots + a_{k-2}t^{k-2})(-1)^k \nonumber \\
&= (a_0 + a_1t + \cdots + a_{k-2}t^{k-2} + a_{k-1}t^{k-1} + t^k)(-1)^k
\end{align}

The Cayley-Hamilton Theorem then easily follows.
\begin{thm}[Cayley-Hamilton Theorem for Linear Operators]
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a finite-dimensional vector space $\mathcal{V}$ and $f(t)$ be its characteristic polynomial. Then $f(T) = T_0$ yields the zero transformation.
\end{thm}
\begin{proof}
We have to show that for all $\vec{v} \in \mathcal{V}$, $[f(T)](\vec{v}) = \textbf{0}$. Let $\mathcal{W}$ be the $T$-cyclic subspace generated by $\vec{v}$ and $\dim(\mathcal{W}) = k$. By part (a) of Properties \ref{proper:B13}, $T^k(\vec{v})$ can be written as a linear combination of $\{\vec{v}, T(\vec{v}), T^2(\vec{v}), \ldots, T^{k-1}(\vec{v})\}$, and thus there exist scalars $a_j$ such that $a_0 \vec{v} + a_1 T(\vec{v}) + \cdots + a_{k-1}T^{k-1}(\vec{v}) + T^k(\vec{v}) = \textbf{0}$. By part (b) of Properties \ref{proper:B13}, the characteristic polynomial of $T|_W$ is then $g(t) = (-1)^k(a_0 + a_1t + \cdots + a_{k-1}t^{k-1} + t^k)$. These two statements together imply that
\begin{align*}
[g(T)](\vec{v}) &= (-1)^k(a_0\vec{v} + a_1T(\vec{v}) + \cdots + a_{k-1}T^{k-1}(\vec{v}) + T^k(\vec{v})) \\
&= (-1)^k[a_0I + a_1T + \cdots + a_{k-1}T^{k-1} + T^k](\vec{v}) = \textbf{0}
\end{align*}
By Properties \ref{proper:B12}, $g(t)$ divides the characteristic polynomial $f(t)$ of $T$ so that $f(t) = q(t)g(t)$, and therefore
\begin{align*}
[f(T)](\vec{v}) = [q(T)g(T)](\vec{v}) = [q(T)][g(T)(\vec{v})] = [q(T)](\textbf{0}) = \textbf{0}
\end{align*}
so we are done. The Cayley-Hamilton Theorem \ref{thm:CHthmmat} for matrices is then an immediate corollary of this.
\end{proof}

\section{Jordan Normal Forms}
\label{sec:JNF}

\subsection{Generalized Eigenvectors, Existence of Jordan Normal Forms}

We have seen that the matrix in Example \ref{exmp:dummyjordan}
\begin{align*}
A &=
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}
\end{align*}
has only one eigenvalue $\lambda = 1$ and its geometric multiplicity $= 1$, is less than the algebraic multiplicity $= 2$: there is only one eigenvector $\vec{v}^{(1)} = (1,0)^T$, so that $A\vec{v}^{(1)} = \lambda\vec{v}^{(1)}$, or equivalently $(A - \lambda I)\vec{v}^{(1)} = \textbf{0}$. Now, let's iterate the procedure and consider $(A - \lambda I)\vec{v}^{(2)} = \vec{v}^{(1)}$, then we have
\begin{align*}
\begin{bmatrix}
1-\lambda & -1 \\
0 & 1-\lambda
\end{bmatrix} 
\vec{v}^{(2)}
=
\begin{bmatrix}
1-1 & -1 \\
0 & 1-1
\end{bmatrix}
\vec{v}^{(2)} 
=
\begin{bmatrix}
0 & -1 \\
0 & 0
\end{bmatrix} 
\vec{v}^{(2)}
=
\begin{bmatrix}
1 \\
0
\end{bmatrix} 
\end{align*}
The solution to $\vec{v}^{(2)}$ is clearly $(0,-1)^T$. $\{\vec{v}^{(1)}, \vec{v}^{(2)}\}$ then forms a basis for $\mathbb{R}^2$ in this case. We have successfully salvaged the problem where an eigenvalue is deficient by finding a \index{Generalized Eigenvector}\keywordhl{generalized eigenvector} $\vec{v}^{(2)}$ that fulfills
\begin{align*}
(A - \lambda I)[(A - \lambda I)\vec{v}^{(2)}] &= (A - \lambda I)\vec{v}^{(1)} \\
(A - \lambda I)^2\vec{v}^{(2)} &= \textbf{0}
\end{align*}
This motivates us to make the following definition.
\begin{defn}
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a vector space $\mathcal{V}$. A non-zero vector $\vec{v} \in \mathcal{V}$ is known as a generalized eigenvector of $T$ corresponding to an eigenvalue of $\lambda$ if
\begin{align}
(T - \lambda I)^p(\vec{v}) &= \textbf{0}
\end{align}
for some integer $p$.
\end{defn}
Clearly, generalized eigenvectors of $p=1$ are just the ordinary eigenvectors. In parallel, we also have
\begin{defn}
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a vector space $\mathcal{V}$. The \index{Generalized Eigenspace}\keywordhl{generalized eigenspace} of $T$ corresponding to an eigenvalue of $\lambda$ is denoted by $\mathcal{K}_\lambda$, which is a subspace of $\mathcal{V}$ such that
\begin{align}
\mathcal{K}_\lambda = \{\vec{v} \in \mathcal{V} \mid (T - \lambda I)^p(\vec{v}) &= \textbf{0} \text{ for some integer p.}\}
\end{align}
\end{defn}
It is not hard to show that $\mathcal{K}_\lambda$ is $T$-invariant and contains the ordinary eigenspace $\mathcal{E}_\lambda$. We note two of the most relevant results of generalized eigenvectors/eigenspaces without proof.
\begin{proper}
The generalized eigenvectors/eigenspaces corresponding to two distinct eigenvalues $\lambda \neq \mu$ are linearly independent.
\end{proper}
\begin{proper}
The generalized eigenspace $\mathcal{K}_\lambda$ for the linear operator $T$, corresponding to some fixed eigenvalue $\lambda$ with an algebraic multiplicity of $m$, satisfies the following properties:
\begin{enumerate}
    \item $\mathcal{K}_\lambda = \mathcal{N}((T-\lambda I)^m)$,
    \item $\dim(\mathcal{K}_\lambda) = m$.
\end{enumerate}
\end{proper}
As a result, if $T$ has a characteristic polynomial
\begin{align}
p(\lambda) = (\lambda_1-\lambda)^{m_1}(\lambda_2-\lambda)^{m_2}\cdots(\lambda_k-\lambda)^{m_k} \label{eqn:lambdammm}
\end{align}
then the respective generalized eigenspaces of all distinct eigenvalues form a direct sum that reproduces the full vector space:
\begin{align}
\mathcal{V} = \mathcal{K}_{\lambda_1} \oplus \mathcal{K}_{\lambda_2} \oplus \cdots \oplus \mathcal{K}_{\lambda_k} \label{eqn:VKKK}
\end{align} 
On the other hand, extending the above example and definitions of generalized eigenvectors/eigenspaces, we can make a \index{Cycle of Generalized Eigenvectors}\keywordhl{cycle of generalized eigenvectors}.
\begin{defn}
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a vector space $\mathcal{V}$ and $\vec{v} \in \mathcal{V}$ be a generalized eigenvector of $T$ corresponding to the eigenvalue $\lambda$. If $p$ is the smallest positive integer such that $(T-\lambda I)^p (\vec{v}) = \textbf{0}$, then
\begin{align}
&\quad \{(T-\lambda I)^{p-1}(\vec{v}), (T-\lambda I)^{p-2}(\vec{v}), \ldots, (T-\lambda I)(\vec{v}),\vec{v}\} \label{eqn:B6}\\
&\coloneqq \{\vec{v}^{(1)}, \vec{v}^{(2)}, \ldots, \vec{v}^{(p-1)}, \vec{v}^{(p)}\} \nonumber
\end{align}
is called a \textit{cycle of generalized eigenvectors} of $T$ corresponding to $\lambda$ with a length of $p$. The vectors $\vec{v}^{(1)} \coloneqq (T-\lambda I)^{p-1}(\vec{v})$ and $\vec{v}^{(p)} \coloneqq \vec{v}$ are called the \textit{initial vector} and \textit{end vector} of that cycle respectively.
\end{defn}
The initial vector of any cycle of generalized eigenvectors is then the only eigenvector in that cycle.\par

Now we are going to deduce the form of a \textit{Jordan block} for a cycle of generalized eigenvectors of some linear operator $T$ corresponding to an eigenvalue of $\lambda$. Consider the linear operator $T-\lambda I$ and the cycle described in (\ref{eqn:B6}). It is not hard to see that (\ref{eqn:B6}) is $(T-\lambda I)$-cyclic (as well as $T$-cyclic) and thus it forms a basis $\gamma$ by part (a) of Properties \ref{proper:B13}. Then the matrix representation of $T-\lambda I$ with respect to $\gamma$ is readily seen to be
\begin{align}
[T-\lambda I]_\gamma = 
\begin{bmatrix}
0 & 1 & 0 & \cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & & \vdots \\
0 & 0 & 0 & \cdots & 1 \\
0 & 0 & 0 & \cdots & 0
\end{bmatrix}
\end{align}
and thus the matrix representation of $T$, also with respect to $\gamma$, will be
\begin{align}
[T]_\gamma = 
\begin{bmatrix}
\lambda & 1 & 0 & \cdots & 0 \\
0 & \lambda  & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & & \vdots \\
0 & 0 & 0 & \cdots & 1 \\
0 & 0 & 0 & \cdots & \lambda 
\end{bmatrix} 
\label{eqn:B8}
\end{align}
This $[T]_\gamma$ is then known as the \index{Jordan Block}\keywordhl{Jordan block} of size $p$ corresponding to $\lambda$, which is "almost diagonal" with $1$s filling above the main diagonal. Note that it has a characteristic polynomial of $g_\gamma(t) = (\lambda-t)^p$. Moreover, it is not hard to check that $([T-\lambda I]_\gamma)^p = [\textbf{0}] = [T_0]_\gamma$ is the zero transformation and $T-\lambda I$ is known as \index{Nilpotent}\keywordhl{nilpotent} (in the generalized eigenspace $\mathcal{K}_\lambda$). To proceed, we first note some intermediate observations.
\begin{thm}
Given $T$ as a linear operator on a vector space $\mathcal{V}$, if $\gamma_1, \gamma_2, \ldots, \gamma_q$ are cycles of generalized eigenvectors of $T$ corresponding to the same eigenvalue $\lambda$ such that the initial vectors of these $\gamma_j$ are linearly independent, then these $\gamma_j$ will be disjoint and their union $\gamma = \bigcup_{j=1}^{q} \gamma_j$ will also be linearly independent.
\end{thm}
\begin{thm}
\label{thm:B27}
Given $T$ as a linear operator on a \textit{finite-dimensional} vector space $\mathcal{V}$ and $\lambda$ as an eigenvalue of $T$, then the associated generalized eigenspace $\mathcal{K}_\lambda$ has a basis consisting of a union of disjoint cycles of generalized eigenvectors corresponding to $\lambda$.
\end{thm}
These two results are derived in Friedberg and we refer the readers to it for the proof. Using our previous example, the matrix
\begin{align*}
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}    
\end{align*}
has a basis made up of two generalized eigenvectors $\gamma = \{(1,0)^T, (0,-1)^T\}$ for the generalized eigenspace $\mathcal{K}_{\lambda=1}$, and thus we can express the matrix as a Jordan block using this \index{Jordan Canonical Basis}\keywordhl{Jordan canonical basis} following (\ref{eqn:B8}) by a coordinate transformation to the $\gamma$ system (refer to Properties \ref{proper:endomorph}):
\begin{align*}
J &= Q^{-1}AQ \\
\begin{bmatrix}
1 & 1 \\
0 & 1
\end{bmatrix} 
&=
\begin{bmatrix}
1 & 0 \\
0 & -1
\end{bmatrix}^{-1} 
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix} 
\begin{bmatrix}
1 & 0 \\
0 & -1
\end{bmatrix}
\end{align*}\par
The final result is that we can derive the \index{Jordan Normal Form}\keywordhl{Jordan Normal Form} of a matrix/linear operator by composing Jordan blocks as a direct sum.
\begin{thm}
Let $T: \mathcal{V} \to \mathcal{V}$ be a linear operator on a finite-dimensional vector space $\mathcal{V}$ and the characteristic polynomial of $T$ \textit{splits} over the chosen field (scalar: $\mathbb{R}$ or $\mathbb{C}$). Then $T$ can be transformed into its Jordan Normal Form as the direct sum of smaller Jordan blocks. Using the matrix representation $A$, we have
\begin{align}
J = Q^{-1}AQ
\end{align}
where $Q$ consists of the union of Jordan canonical bases for all the Jordan blocks in columns.
\end{thm}
By Theorem \ref{thm:B27}, $\mathcal{K}_{\lambda_i}$ has a basis $\gamma_{\lambda_i} = \bigcup_{j=1}^{q} (\gamma_j)_{\lambda_i}$ formed by smaller cycles to the same eigenvalue $\lambda_i$ and hence $[T|_{\mathcal{K}_{\lambda_i}}]_\gamma = {(J_1)}_{\lambda_i} \oplus {(J_2)}_{\lambda_i} \oplus \cdots {(J_q)}_{\lambda_i} = J_{\lambda_i}$, and by (\ref{eqn:VKKK}), we have
\begin{align}
[T]_\beta = J_{\lambda_1} \oplus J_{\lambda_2} \oplus \cdots \oplus J_{\lambda_k}
\end{align}
as its Jordan Normal Form, where $\beta = \bigcup_{i=1}^{k} \gamma_{\lambda_i}$. The Jordan Normal Form of the last example is simply its only Jordan block:
\begin{align*}
J = \begin{bmatrix}
1 & 1 \\
0 & 1
\end{bmatrix} 
\end{align*}
To conclude, Jordan Normal Forms are the generalization and remedy of diagonalized forms where the constraint is relaxed to accommodate the case of deficient eigenvalues, so that it exists whenever the characteristic polynomial splits. But the price is that there will be some extra "1" above the main diagonal in the Jordan blocks of size $\geq 2$. Conversely, if all the Jordan blocks are of size $1$, then it will coincide with the usual diagonalization. For example,
\begin{align*}
J = 
\begin{bmatrix}
\mathcolor{red}{2} & \mathcolor{red}{1} & 0 & 0 & 0 & 0 & 0 \\
\mathcolor{red}{0} & \mathcolor{red}{2} & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & \mathcolor{blue}{3} & \mathcolor{blue}{1} & \mathcolor{blue}{0} & 0 & 0 \\
0 & 0 & \mathcolor{blue}{0} & \mathcolor{blue}{3} & \mathcolor{blue}{1} & 0 & 0 \\
0 & 0 & \mathcolor{blue}{0} & \mathcolor{blue}{0} & \mathcolor{blue}{3} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \mathcolor{Green}{3} & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & \mathcolor{purple}{3}
\end{bmatrix}
\end{align*}
consists of one Jordan block corresponding to $\lambda = 2$ of size $2$ (red), one Jordan block corresponding to $\lambda = 3$ of size $3$ (blue), and two Jordan blocks corresponding to $\lambda = 3$ of size $1$ (green and purple). The positions of blocks can be swapped with each other. However, note that a characteristic polynomial (\ref{eqn:lambdammm}) does not uniquely determine the Jordan Normal Form. In the above example, the characteristic polynomial is $p(\lambda) = (2-\lambda)^2(3-\lambda)^5$. But a different Jordan Normal Form like
\begin{align*}
J' =
\begin{bmatrix}
2 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 2 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 3 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 3 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 3 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 3 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 3
\end{bmatrix}
\end{align*}
also shares the same characteristic polynomial. We note in passing that there is a closely related concept called a \index{Minimal Polynomial}\textit{minimal polynomial}, which the readers can search about if they are interested in exploring the deeper theoretical aspect.

\subsection{Methods to Compute Jordan Normal Forms}

In this subsection, we will investigate the techniques to compute the Jordan Normal Form of any given matrix practically.

\begin{exmp}
Derive the Jordan Normal Form of
\begin{align*}
A = 
\begin{bmatrix}
3&3&-2\\ 
-1&0&1\\ 
0&1&1
\end{bmatrix}
\end{align*}
\end{exmp}
\begin{solution}
We leave it to the readers to check that the characteristic polynomial is $p(\lambda) = (1-\lambda)^2(2-\lambda)$. Hence, both the algebraic and geometric multiplicity of $\lambda = 2$ is $1$, which is easy to deal with. The problem will be $\lambda = 1$ which has an algebraic multiplicity of $2$, and there are two possibilities: if the geometric multiplicity is also $2$, then the situation is the same as the usual diagonalization where we have two ordinary eigenvectors. Otherwise, if the geometric multiplicity is $1$, then we will have a single Jordan block of length $2$. Now consider the system $A-I = \textbf{0}$:
\begin{align*}
\left[\begin{array}{@{\,}wc{10pt}wc{10pt}wc{10pt}|wc{10pt}@{\,}}
2 & 3 & -2 & 0  \\ 
-1 & -1 & 1 & 0\\ 
0 & 1 & 0 & 0
\end{array}\right] 
& \to 
\left[\begin{array}{@{\,}wc{10pt}wc{10pt}wc{10pt}|wc{10pt}@{\,}}
1 & \frac{3}{2} & -1 & 0 \\ 
-1 & -1 & 1 & 0 \\ 
0 & 1 & 0 & 0
\end{array}\right] 
& \frac{1}{2}R_1 \to R_1 \\
& \to 
\left[\begin{array}{@{\,}wc{10pt}wc{10pt}wc{10pt}|wc{10pt}@{\,}}
1 & \frac{3}{2} & -1 & 0 \\ 
0 & \frac{1}{2} & 0 & 0 \\ 
0 & 1 & 0 & 0
\end{array}\right] 
& R_2 + R_1 \to R_2 \\
& \to 
\left[\begin{array}{@{\,}wc{10pt}wc{10pt}wc{10pt}|wc{10pt}@{\,}}
1 & \frac{3}{2} & -1 & 0 \\ 
0 & 1 & 0 & 0 \\ 
0 & 0 & 0 & 0
\end{array}\right]
&
\begin{aligned}
2R_2 \to R_2 \\
R_3 - R_2 \to R_3
\end{aligned}
\end{align*}
The nullity is $1$ and there is only one eigenvector. By letting the third entry be the free variable, it is not hard to see that the eigenvector will be $\vec{v}^{(1)} = (1,0,1)^T$. So there will be only one Jordan block, and due to this, we can obtain the next generalized eigenvector in a relatively simple reverse fashion by solving
\begin{align*}
(A - I)\vec{v}^{(2)} &= \vec{v}^{(1)} \\
\begin{bmatrix}
2 & 3 & -2\\ 
-1 & -1 & 1\\ 
0 & 1 & 0
\end{bmatrix}
\vec{v}^{(2)} &=
\begin{bmatrix}
1 \\
0 \\
1
\end{bmatrix}
\end{align*}
We can pick any particular solution for this, and again, we will leave it to the readers to confirm that $\vec{v}^{(2)} = (-1,1,0)^T$ is a valid choice. Finally, we simply note that the eigenvector for the other eigenvalue $\lambda = 2$ is $(-1,1,1)^T$. Thus, the required Jordan Normal Form is
\begin{align*}
Q^{-1}AQ &= 
\begin{bmatrix}
1 & -1 & -1\\
0 & 1 & 1\\
1 & 0 & 1
\end{bmatrix}^{-1}
\begin{bmatrix}
3&3&-2\\ 
-1&0&1\\ 
0&1&1
\end{bmatrix}
\begin{bmatrix}
1 & -1 & -1\\
0 & 1 & 1\\
1 & 0 & 1
\end{bmatrix} \\
J &=
\begin{bmatrix}
1 & 1 & 0 \\
0 & 1 & 0 \\
0 & 0 & 2
\end{bmatrix}
\end{align*}
\end{solution}

Next, we will tackle a much harder case where there are multiple Jordan blocks corresponding to the same eigenvalue.
\begin{exmp}
Convert the matrix
\begin{align*}
A = \begin{bmatrix}
5&-1&-1&-2&0&-1\\ 
-1&19&9&24&6&17\\ 
-1&3&5&4&2&1\\ 
2&-10&-6&-12&-4&-10\\ 
1&-5&-3&-8&2&-5\\ 
-1&1&1&2&0&5
\end{bmatrix}    
\end{align*}
to its Jordan Normal Form.
\end{exmp}
\begin{solution}
We will leave it to the readers to verify that the characteristic polynomial is $p(\lambda) = (4-\lambda)^6$, so the only eigenvalue is $\lambda = 4$. There can be many possible scenarios: the geometric multiplicity may be $6$, so there will be $6$ eigenvectors and $6$ singleton Jordan blocks, which reduces to the usual diagonalization. The other extreme will have a single Jordan block of length $6$. In between these two extremes, we may have $2$ to $5$ Jordan blocks of varying sizes. To begin, we compute the eigenvectors that are also the generalized eigenvectors of order $1$, by solving the system $A-4I = \textbf{0}$:
\begin{align*}
& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0\\ 
-1 & 15 & 9 & 24 & 6 & 17 & 0\\ 
-1 & 3 & 1 & 4 & 2 & 1 & 0\\  
2 & -10 & -6 & -16 & -4 & -10 & 0\\  
1 & -5 & -3 & -8 & -2 & -5 & 0\\  
-1 & 1 & 1 & 2 & 0 & 1 & 0
\end{array}\right] \\
\to& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0\\ 
0 & 14 & 8 & 22 & 6 & 16 & 0\\ 
0 & 2 & 0 & 2 & 2 & 0 & 0\\  
0 & -8 & -4 & -12 & -4 & -8 & 0\\  
0 & -4 & -2 & -6 & -2 & -4 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] &
\begin{aligned}
R_2 + R_1 \to R_2\\
R_3 + R_1 \to R_3\\
R_4 - 2R_1 \to R_4\\
R_5 - R_1 \to R_5\\
R_6 + R_1 \to R_6
\end{aligned} \\
\to& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0\\ 
0 & 1 & 0 & 1 & 1 & 0 & 0\\  
0 & 14 & 8 & 22 & 6 & 16 & 0\\ 
0 & -8 & -4 & -12 & -4 & -8 & 0\\  
0 & -4 & -2 & -6 & -2 & -4 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] &
\begin{aligned}
R_2 \leftrightarrow R_3 \\
\frac{1}{2} R_2 \to R_2
\end{aligned} \\
\to& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0\\ 
0 & 1 & 0 & 1 & 1 & 0 & 0\\  
0 & 0 & 8 & 8 & -8 & 16 & 0\\ 
0 & 0 & -4 & -4 & 4 & -8 & 0\\  
0 & 0 & -2 & -2 & 2 & -4 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] &
\begin{aligned}
R_3 - 14R_2 \to R_3 \\
R_4 + 8R_2 \to R_4 \\
R_5 + 4R_2 \to R_5
\end{aligned} \\
\to& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0\\ 
0 & 1 & 0 & 1 & 1 & 0 & 0\\  
0 & 0 & 1 & 1 & -1 & 2 & 0\\ 
0 & 0 & 0 & 0 & 0 & 0 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] &
\begin{aligned}
\frac{1}{8} R_3 \to R_3 \\
R_4 + 4R_3 \to R_4 \\
R_5 + 2R_3 \to R_5
\end{aligned} \\
\to& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}@{\,}}
1 & 0 & 0 & 0 & 0 & 1 & 0\\ 
0 & 1 & 0 & 1 & 1 & 0 & 0\\  
0 & 0 & 1 & 1 & -1 & 2 & 0\\ 
0 & 0 & 0 & 0 & 0 & 0 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0\\  
0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_1 + R_3 \to R_1 \\
R_1 + R_2 \to R_1
\end{aligned}
\end{align*}
The nullity is $3$. We can assign the last three components as the free variables to obtain three eigenvectors: $(0,-1,-1,1,0,0)^T$, $(0,-1,1,0,1,0)^T$, $(-1,0,-2,0,0,1)^T$ (the readers should do the checking), and there will be three Jordan blocks of unknown lengths. Now here is the trick. We create an augmented system where the left is $A-\lambda I$ and the right consists of the three eigenvectors we just obtained. We will derive its RREF but we can reuse the previous reduction step:
\begin{align*}
& \left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0 & 0 & -1\\ 
-1 & 15 & 9 & 24 & 6 & 17 & -1 & -1 & 0 \\ 
-1 & 3 & 1 & 4 & 2 & 1 & -1 & 1 & -2\\  
2 & -10 & -6 & -16 & -4 & -10 & 1 & 0 & 0\\  
1 & -5 & -3 & -8 & -2 & -5 & 0 & 1 & 0 \\  
-1 & 1 & 1 & 2 & 0 & 1 & 0 & 0 & 1
\end{array}\right] \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0 & 0 & -1\\  
0 & 14 & 8 & 22 & 6 & 16 & -1 & -1 & -1 \\ 
0 & 2 & 0 & 2 & 2 & 0 & -1 & 1 & -3\\  
0 & -8 & -4 & -12 & -4 & -8 & 1 & 0 & 2\\ 
0 & -4 & -2 & -6 & -2 & -4 & 0 & 1 & 1 \\    
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_2 + R_1 \to R_2\\
R_3 + R_1 \to R_3\\
R_4 - 2R_1 \to R_4\\
R_5 - R_1 \to R_5\\
R_6 + R_1 \to R_6
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0 & 0 & -1\\  
0 & 1 & 0 & 1 & 1 & 0 & -\frac{1}{2} & \frac{1}{2} & -\frac{3}{2}\\
0 & 14 & 8 & 22 & 6 & 16 & -1 & -1 & -1 \\ 
0 & -8 & -4 & -12 & -4 & -8 & 1 & 0 & 2\\ 
0 & -4 & -2 & -6 & -2 & -4 & 0 & 1 & 1 \\    
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_2 \leftrightarrow R_3 \\
\frac{1}{2} R_2 \to R_2
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0 & 0 & -1\\  
0 & 1 & 0 & 1 & 1 & 0 & -\frac{1}{2} & \frac{1}{2} & -\frac{3}{2}\\ 
0 & 0 & 8 & 8 & -8 & 16 & 6 & -8 & 20 \\ 
0 & 0 & -4 & -4 & 4 & -8 & -3 & 4 & -10\\   
0 & 0 & -2 & -2 & 2 & -4 & -2 & 3 & -5 \\     
0 & 0 & 0 & 0 & 0 & 0 &  0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_3 - 14R_2 \to R_3 \\
R_4 + 8R_2 \to R_4 \\
R_5 + 4R_2 \to R_5
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & -1 & -1 & -2 & 0 & -1 & 0 & 0 & -1\\  
0 & 1 & 0 & 1 & 1 & 0 & -\frac{1}{2} & \frac{1}{2} & -\frac{3}{2}\\ 
0 & 0 & 1 & 1 & -1 & 2 & \frac{3}{4} & -1 & \frac{5}{2} \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\   
0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 1 & 0 \\     
0 & 0 & 0 & 0 & 0 & 0 &  0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
\frac{1}{8} R_3 \to R_3 \\
R_4 + 4R_3 \to R_4 \\
R_5 + 2R_3 \to R_5
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & 0 & 0 & 0 & 0 & 1 & \frac{1}{4} & -\frac{1}{2} & 0\\  
0 & 1 & 0 & 1 & 1 & 0 & -\frac{1}{2} & \frac{1}{2} & -\frac{3}{2}\\ 
0 & 0 & 1 & 1 & -1 & 2 & \frac{3}{4} & -1 & \frac{5}{2} \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\   
0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{2} & 1 & 0 \\     
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_1 + R_3 \to R_1 \\
R_1 + R_2 \to R_1
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & 0 & 0 & 0 & 0 & 1 & \frac{1}{4} & -\frac{1}{2} & 0\\  
0 & 1 & 0 & 1 & 1 & 0 & -\frac{1}{2} & \frac{1}{2} & -\frac{3}{2}\\ 
0 & 0 & 1 & 1 & -1 & 2 & \frac{3}{4} & -1 & \frac{5}{2} \\ 
0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 0 \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\   
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_4 \leftrightarrow R_5 \\
-2R_4 \to R_4
\end{aligned} \\
\to&
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\  
0 & 1 & 0 & 1 & 1 & 0 & 0 & -\frac{1}{2} & -\frac{3}{2}\\ 
0 & 0 & 1 & 1 & -1 & 2 & 0 & \frac{1}{2} & \frac{5}{2} \\ 
0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 0 \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\   
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right] & 
\begin{aligned}
R_1 - \frac{1}{4}R_4 \to R_1 \\
R_2 + \frac{1}{2}R_4 \to R_1 \\
R_3 - \frac{3}{4}R_4 \to R_1
\end{aligned}
\end{align*}
We only need to look at the top three rows where pivots exist in the left portion and discard the bottom three rows. The first column to the right is trivially all zeros, so it is a dead end. Meanwhile, by considering the second/third column to the right, we can simply take the particular solutions $(0,-\frac{1}{2},\frac{1}{2},0,0,0)^T$ and $(0,-\frac{3}{2},\frac{5}{2},0,0,0)^T$ as candidates for the generalized eigenvectors of order $2$. Now there are $5$ generalized eigenvectors, and the last generalized eigenvector will be found with an order of $3$. We repeat the same procedure and the result is simply shown below:
\begin{align*}
\left[\begin{array}{@{\,}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}wc{13pt}|wc{13pt}wc{13pt}wc{13pt}@{\,}}
1 & 0 & 0 & 0 & 0 & 1 & -\frac{1}{2} & 0 \\  
0 & 1 & 0 & 1 & 1 & 0 & \frac{1}{2} & 0\\ 
0 & 0 & 1 & 1 & -1 & 2 & -1 & 0 \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\   
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{array}\right]
\end{align*}
So the end vector of the Jordan block with length $3$ is $(-\frac{1}{2}, \frac{1}{2}, -1,0,0,0)^T$. There will be two other Jordan blocks with lengths $1$ and $2$ respectively. We will compute the cycle of generalized eigenvectors for all three Jordan blocks in the opposite direction. We begin from the largest Jordan block of length $3$, where the end vector has been found to be $(-\frac{1}{2}, \frac{1}{2}, -1,0,0,0)^T$, and then compute:
\begin{align*}
(A-4I)\left(-\frac{1}{2}, \frac{1}{2}, -1,0,0,0\right)^T &= (0,-1,1,0,0,0)^T \\
(A-4I)(0,-1,1,0,0,0)^T &= (0,-6,-2,4,2,0)^T
\end{align*}
For the second Jordan block of length $2$, we note that at the third stage, the inconsistent $1$ to the right appears in the second column, and hence we will choose the second candidate generalized eigenvector $(0,-\frac{3}{2},\frac{5}{2},0,0,0)^T$ as the new end vector at the second stage, and calculate:
\begin{align*}
(A-4I)\left(0,-\frac{3}{2},\frac{5}{2},0,0,0\right)^T &= (-1,0,-2,0,0,1)^T    
\end{align*}
By the same logic, we choose the first eigenvector $(0,-1,-1,1,0,0)^T$ for the last singleton Jordan block. The final answer will then be
\begin{align*}
Q^{-1}AQ =&
\small
\begin{bmatrix}
0 & 0 & -\frac{1}{2} & -1 & 0 & 0 \\
-6 & -1 & \frac{1}{2} & 0 & -\frac{3}{2} & -1\\
-2 & 1 & -1 & -2 & \frac{5}{2} & -1\\
4 & 0 & 0 & 0 & 0 & 1\\
2 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0
\end{bmatrix}^{-1}
\begin{bmatrix}
5&-1&-1&-2&0&-1\\ 
-1&19&9&24&6&17\\ 
-1&3&5&4&2&1\\ 
2&-10&-6&-12&-4&-10\\ 
1&-5&-3&-8&2&-5\\ 
-1&1&1&2&0&5
\end{bmatrix} \\ 
& \small\begin{bmatrix}
0 & 0 & -\frac{1}{2} & -1 & 0 & 0 \\
-6 & -1 & \frac{1}{2} & 0 & -\frac{3}{2} & -1\\
-2 & 1 & -1 & -2 & \frac{5}{2} & -1\\
4 & 0 & 0 & 0 & 0 & 1\\
2 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0
\end{bmatrix} \\
J =&
\begin{bmatrix}
4&1&0&0&0&0\\ 
0&4&1&0&0&0\\ 
0&0&4&0&0&0\\ 
0&0&0&4&1&0\\ 
0&0&0&0&4&0\\ 
0&0&0&0&0&4
\end{bmatrix}
\end{align*}
\end{solution}

\section{Matrix Exponentials}
\label{sec:matexp}

In Section \ref{section:sysode}, we have seen how to solve a system of first-order ODEs $\textbf{y}' = A\textbf{y}$ by diagonalization: the new system $\textbf{z}' = D\textbf{z}$ where $D$ is a diagonal matrix can be solved according to (\ref{eqn:ODEDz}) so that $z_{i} = c_ie^{D_{ii}x}$. Since in the one-dimensional case, we simply have $z = ce^{\alpha x}$, we may like to ask if we can also write something along the lines of $\textbf{z} = e^{Dx}\textbf{c}$ (or even one step further, $\textbf{y} = e^{Ax}\textbf{c}$) without resorting to dealing with each variable separately. This motivates the definition of a \index{Matrix Exponential}\keywordhl{matrix exponential}. Recall that from elementary Calculus, the usual exponential function can be expressed in terms of a Taylor series:
\begin{align*}
e^{x} = 1 + x + \frac{x^2}{2!} + \cdots = \sum_{k=0}^{\infty} \frac{x^k}{k!}
\end{align*}
and thus in parallel, we demand that
\begin{defn}
The exponential of a matrix $X$ is computed by
\begin{align}
e^{X} = I + X + \frac{X^2}{2!} + \cdots = \sum_{k=0}^{\infty} \frac{X^k}{k!}
\end{align}
\end{defn}
When we set $X = Ax$ where $A$ is a matrix and $x$ is an independent variable, it becomes
\begin{align}
e^{Ax} = I + Ax + \frac{A^2}{2!}x^2 + \cdots = \sum_{k=0}^{\infty} \frac{A^k}{k!}x^k  
\end{align}
Now, if we are given a diagonal matrix $D$ such that
\begin{align*}
D = \begin{bmatrix}
D_{11} & 0 & \cdots & 0 \\
0 & D_{22} & & 0 \\
\vdots & & \ddots \\
0 & 0 & & D_{nn}
\end{bmatrix}
\end{align*}
Then it is obvious that
\begin{align*}
D^k =  \begin{bmatrix}
D_{11}^k & 0 & \cdots & 0 \\
0 & D_{22}^k & & 0 \\
\vdots & & \ddots \\
0 & 0 & & D_{nn}^k
\end{bmatrix}
\end{align*}
and thus
\begin{align*}
e^{Dx} &= I + Dx + \frac{D^2}{2!}x^2 + \cdots = \sum_{k=0}^{\infty} \frac{D^k}{k!}x^k \\
&= \begin{bmatrix}
1 & 0 & \cdots & 0 \\
0 & 1 & & 0 \\
\vdots & & \ddots \\
0 & 0 & & 1
\end{bmatrix}
+
\begin{bmatrix}
D_{11}x & 0 & \cdots & 0 \\
0 & D_{22}x & & 0 \\
\vdots & & \ddots \\
0 & 0 & & D_{nn}x
\end{bmatrix} \\
&\quad + 
\begin{bmatrix}
\frac{D_{11}^2}{2!}x^2 & 0 & \cdots & 0 \\
0 & \frac{D_{22}^2}{2!}x^2 & & 0 \\
\vdots & & \ddots \\
0 & 0 & & \frac{D_{nn}^2}{2!}x^2
\end{bmatrix}
+ \cdots \\
&=
\footnotesize \begin{bmatrix}
1 + D_{11}x + \frac{D_{11}^2}{2!}x^2 + \cdots & 0 & \cdots & 0 \\
0 & 1 + D_{22}x + \frac{D_{22}^2}{2!}x^2 + \cdots & & 0 \\
\vdots & & \ddots \\
0 & 0 & & 1 + D_{nn}x + \frac{D_{nn}^2}{2!}x^2 + \cdots
\end{bmatrix} \\
&=
\begin{bmatrix}
e^{D_{11}x} & 0 & \cdots & 0 \\
0 & e^{D_{22}x} & & 0 \\
\vdots & & \ddots \\
0 & 0 & & e^{D_{nn}x}
\end{bmatrix}
\end{align*}
so the proposed solution form $\textbf{z} = e^{Dx}\textbf{c}$ will be consistent with what we have established in Section \ref{section:sysode}.\par

Since diagonalization is not always possible but we can compensate with a Jordan Normal Form, the next step is to derive $e^{Jx}$ where $J$ is a Jordan block with an eigenvalue of $\lambda$, and we will start from the easier case of length $3$:
\begin{align*}
J = 
\begin{bmatrix}
\lambda & 1 & 0\\ 
0 & \lambda & 1 \\
0 & 0 & \lambda
\end{bmatrix}
\end{align*}
then
\begin{align*}
J^2 &= 
\begin{bmatrix}
\lambda & 1 & 0\\ 
0 & \lambda & 1 \\
0 & 0 & \lambda
\end{bmatrix}
\begin{bmatrix}
\lambda & 1 & 0\\ 
0 & \lambda & 1 \\
0 & 0 & \lambda
\end{bmatrix}
=
\begin{bmatrix}
\lambda^2 & 2\lambda & 1 \\ 
0 & \lambda^2 & 2\lambda \\
0 & 0 & \lambda^2
\end{bmatrix} \\
J^3 &=
\begin{bmatrix}
\lambda^2 & 2\lambda & 1 \\ 
0 & \lambda^2 & 2\lambda \\
0 & 0 & \lambda^2
\end{bmatrix}
\begin{bmatrix}
\lambda & 1 & 0\\ 
0 & \lambda & 1 \\
0 & 0 & \lambda
\end{bmatrix}
=
\begin{bmatrix}
\lambda^3 & 3\lambda^2 & 3\lambda \\
0 & \lambda^3 & 3\lambda^2 \\
0 & 0 & \lambda^3
\end{bmatrix}
\end{align*}
Going on, we have
\begin{align*}
J^k =
\begin{bmatrix}
\lambda^k & k\lambda^{k-1} & \frac{k(k-1)}{2}\lambda^{k-2}\\
0 & \lambda^k & k\lambda^{k-1} \\
0 & 0 & \lambda^k 
\end{bmatrix}
\end{align*}
and thus
\begin{align*}
e^{Jx} &= I + Jx + \frac{J^2}{2!}x^2 + \frac{J^3}{3!}x^3 + \cdots = \sum_{k=0}^{\infty} \frac{J^k}{k!}x^k  \\
&= 
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}
+ 
\begin{bmatrix}
\lambda x & x & 0\\ 
0 & \lambda x & x \\
0 & 0 & \lambda x
\end{bmatrix}
+ \frac{1}{2!}
\begin{bmatrix}
\lambda^2x^2 & 2\lambda x^2 & x^2 \\ 
0 & \lambda^2x^2 & 2\lambda x^2 \\
0 & 0 & \lambda^2x^2
\end{bmatrix} \\
&\quad + \frac{1}{3!}
\begin{bmatrix}
\lambda^3x^3 & 3\lambda^2x^3 & 3\lambda x^3\\
0 & \lambda^3x^3 & 3\lambda^2x^3 \\
0 & 0 & \lambda^3x^3
\end{bmatrix} + \cdots \\
&= 
\footnotesize
\begin{bmatrix}
1 + \lambda x + \frac{\lambda^2}{2!}x^2 + \frac{\lambda^3}{3!}x^3 + \cdots & x(1 + \lambda x + \frac{\lambda^2}{2!}x + \cdots) & \frac{x^2}{2!}(1 + \lambda x + \cdots) \\
0 & 1 + \lambda x + \frac{\lambda^2}{2!}x^2 + \frac{\lambda^3}{3!}x^3 + \cdots & x(1 + \lambda x + \frac{\lambda^2}{2!}x + \cdots) \\
0 & 0 & 1 + \lambda x + \frac{\lambda^2}{2!}x^2 + \frac{\lambda^3}{3!}x^3 + \cdots
\end{bmatrix} \\
&=
\begin{bmatrix}
e^{\lambda x} & xe^{\lambda x} & \frac{x^2}{2!}e^{\lambda x} \\
0 & e^{\lambda x} & xe^{\lambda x} \\
0 & 0 & e^{\lambda x}
\end{bmatrix}
\end{align*}
In general, we have the following formula.
\begin{proper}
The matrix exponential $e^{Jx}$ where $J$ is a Jordan block of size $n$ corresponding to an eigenvalue of $\lambda$ is
\begin{align}
e^{Jx} =
\begin{bmatrix}
e^{\lambda x} & xe^{\lambda x} & \frac{x^2}{2!}e^{\lambda x} & \cdots & \frac{x^{n-1}}{(n-1)!} e^{\lambda x} \\
0 & e^{\lambda x} & xe^{\lambda x} & & \frac{x^{n-2}}{(n-2)!} e^{\lambda x} \\
0 & 0 &e^{\lambda x} & & \frac{x^{n-3}}{(n-3)!} e^{\lambda x} \\
\vdots & & & \ddots & \vdots \\
0 & 0 & 0 & \cdots & e^{\lambda x}
\end{bmatrix}
\end{align}
\end{proper}
Given a Jordan Normal Form $J = J_1 \oplus J_2 \oplus \cdots$ as a direct sum composed of smaller Jordan blocks, its matrix exponential $e^{Jx}$ will also be a direct sum of the respective smaller matrix exponential $e^{J_ix}$, i.e.\ $e^{Jx} = e^{J_1x} \oplus e^{J_2x} \oplus \cdots$. For example, if $J = J_1 \oplus J_2$ with
\begin{align*}
J_1 &= \begin{bmatrix}
2 & 1 \\
0 & 2
\end{bmatrix}
& \text{and} & & J_2 &=
\begin{bmatrix}
5 & 1 & 0 & 0 \\
0 & 5 & 1 & 0 \\
0 & 0 & 5 & 1 \\
0 & 0 & 0 & 5
\end{bmatrix}
\end{align*}
then
\begin{align*}
e^{Jx} =
\begin{bmatrix}
e^{2x} & xe^{2x} & 0 & 0 & 0 & 0 \\
0 & e^{2x} & 0 & 0 & 0 & 0 \\
0 & 0 & e^{5x} & xe^{5x} & \frac{x^2}{2}e^{5x} & \frac{x^3}{6}e^{5x} \\
0 & 0 & 0 & e^{5x} & xe^{5x} & \frac{x^2}{2}e^{5x} \\
0 & 0 & 0 & 0 & e^{5x} & xe^{5x} \\
0 & 0 & 0 & 0 & 0 & e^{5x} 
\end{bmatrix}
\end{align*}
Extending the above results, if we have a square matrix $A$ that has a Jordan Normal Form of $J$ such that $J = Q^{-1}AQ$, where $Q$ forms its Jordan canonical basis, then
\begin{align}
e^{Ax} &= I + Ax + \frac{A^2}{2!}x^2 + \cdots = \sum_{k=0}^{\infty} \frac{A^k}{k!}x^k  \nonumber \\
&= QIQ^{-1} + QJQ^{-1}x + \frac{1}{2!}(QJQ^{-1})(QJQ^{-1})x^2 + \cdots \nonumber \\
&= QIQ^{-1} + QJQ^{-1}x + \frac{1}{2!}QJ^2Q^{-1}x^2 + \cdots \nonumber \\
&= Q(I + Jx + \frac{J^2}{2!}x^2 + \cdots)Q^{-1} \nonumber \\
&= Qe^{Jx}Q^{-1}
\end{align}
Now let's apply these to solve Example \ref{exmp:chemtrace} when $k_1 = k_2 = k$ and diagonalization is no longer possible. The system becomes
\begin{align*}
\textbf{y}'
= A\textbf{y} =
\begin{bmatrix}
-k & 0 & 0 \\
k & -k & 0 \\
0 & k & 0 
\end{bmatrix}
\textbf{y}
\end{align*}
It remains the same that the eigenvector for $\lambda = 0$ is $(0,0,1)^T$. For $\lambda = -k$, it will form a Jordan block of length $2$: the initial eigenvector can be seen to be $(0,-1,1)^T$ and the end generalized eigenvector is found by solving
\begin{align*}
\left[\begin{array}{ccc|c}
0 & 0 & 0 & 0 \\
k & 0 & 0 & -1 \\
0 & k & k & 1
\end{array}\right]
\end{align*}
where one of the possible choices is $(-\frac{1}{k},0,\frac{1}{k})^T$. The solution of the system is then 
\begin{align*}
\textbf{y} = e^{At} \textbf{c} = Qe^{Jt}Q^{-1}\textbf{c}
\end{align*}
with
\begin{align*}
Q &= \begin{bmatrix}
0 & 0 & -\frac{1}{k}\\
0 & -1 & 0 \\
1 & 1 & \frac{1}{k}
\end{bmatrix} &
J &=
\begin{bmatrix}
0 & 0 & 0 \\
0 & -k & 1 \\
0 & 0 & -k
\end{bmatrix}
\end{align*}
and 
\begin{align*}
e^{Jt} =
\begin{bmatrix}
1 & 0 & 0 \\
0 & e^{-kt} & te^{-kt} \\
0 & 0 & e^{-kt}
\end{bmatrix}
\end{align*}
Therefore,
\begin{align*}
\textbf{y} &= Qe^{Jt}Q^{-1}\textbf{c} \\
&= \begin{bmatrix}
0 & 0 & -\frac{1}{k}\\
0 & -1 & 0 \\
1 & 1 & \frac{1}{k}
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 \\
0 & e^{-kt} & te^{-kt} \\
0 & 0 & e^{-kt}
\end{bmatrix}
\begin{bmatrix}
0 & 0 & -\frac{1}{k}\\
0 & -1 & 0 \\
1 & 1 & \frac{1}{k}
\end{bmatrix}^{\mathrlap{-1}}
\textbf{c} \\
&= \begin{bmatrix}
0 & 0 & -\frac{1}{k}\\
0 & -1 & 0 \\
1 & 1 & \frac{1}{k}
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 \\
0 & e^{-kt} & te^{-kt} \\
0 & 0 & e^{-kt}
\end{bmatrix}
\begin{bmatrix}
1 & 1 & 1\\ 
0 & -1 & 0\\ 
-k & 0 & 0
\end{bmatrix}
\textbf{c}\\
&= 
\begin{bmatrix}
e^{-kt} & 0 & 0 \\ 
ke^{-kt}t & e^{-kt} & 0 \\ 
-ke^{-kt}t-e^{-kt}+1 & -e^{-kt}+1 & 1
\end{bmatrix}
\begin{bmatrix}
c_1 \\
c_2 \\
c_3
\end{bmatrix}
\end{align*}
Plugging in the initial condition gives $c_1 = 1, c_2 = 0, c_3 = 0$, and thus the new solution will be
\begin{align*}
\left\{\begin{alignedat}{2}
&[P]& &= e^{-kt} \\
&[Q]& &= ke^{-kt}t \\
&[R]& &= 1-ke^{-kt}t-e^{-kt} 
\end{alignedat}\right.
\end{align*}