Preparation of exam.

Preparation Notes - Hell lot of Math and Algorithms; very helpful.

Note : These notes are made during my preparation for my exams. They might be crude and may not be in their best form of presentation. Read at your own risk.

---

Subject : IC 021 - Computational Techniques in Control Engineering

Professor : Dr. A Ramakalyan.

Difficulty : Runs to Heavens and goes around it $\infty$ times.

Books:
1. B N Datta’s Numerical Methods for Linear Control Systems. - referenced as [1]
2. Gilbert Strang’s Linear ALgebra and it’s applications. - referenced as [2]
3. S Kumaresen’s Linear Algebra - A Geometric Approach. - referenced as [3]

---

1. Computing $e^\mathcal{A}$[1]

1. There are many ways of doing this, one of the very popular one is the Taylor Series method.
   We simply expand $e^\mathcal{A}$ as per the Taylor series.
   Hence, $e^\mathcal{A}$ can be written as,
   $$\boxed{e^\mathcal{A} = I + \mathcal{A} + \frac{\mathcal{A}^{2}}{2!} + \frac{\mathcal{A^3}}{3!} + ...}$$

   • The drawbacks to this method are that a large number of terms is needed for convergence, and even when convergence occurs, the answer can be totally wrong.

2. Hence, I will now get into the Pade’s Scaling and Squaring algorithm for this task.
   Suppose that, $f(x)$ is a power series presented as,
   $$f(x) = f_0 + f_1 x+ f_2 x^2+...$$
   Then the Pade’s Scaling and Squaring Algorithm, given an approximation of $f(x)$ as ,
   $$\boxed{f(x) \approx \frac{c(x)}{d(x)}= \frac{\sum_{k=0}^p c_{k} x^{k}}{\sum_{k=0}^q d_{k} x^{k}}}$$
   The $c_k$ and $d_k$ are chosen such that the terms containing $x^0, x^1, x^2,...,x^{p+q}$ are cancelled out in $f(x)d(x) - c(x)$. The order of Pade’s approximation is $p+q$. The ratio $(\frac{c(x)}{d(x)})$ is unique and the coefficients $c_0, c_1,...,c_p$ and $d_0, d_1,...,d_q$ always exist.
   
   Similarly, Pade’s approximation for $e^\mathcal{A}$ is given by,
   $$\boxed{R_{pq}(\mathcal{A}) = [D_{pq}(\mathcal{A})]^{-1}N_{pq}(\mathcal{A})}$$
   Where,
   $$\boxed{D_{pq}(\mathcal{A})=\sum_{k=0}^q\frac{(p+q-j)!q!}{(p+q)!j!(q-j)!}(-\mathcal{A})^{j}}$$
   and
   $$\boxed{N_{pq}(\mathcal{A})=\sum_{k=0}^p\frac{(p+q-j)!p!}{(p+q)!j!(p-j)!}(\mathcal{A})^{j}}$$

   • It can be shown (Exercise 5.16) that if p and q are large, or if A is a matrix having all its eigenvalues with negative real parts, then $D_{pq}(\mathcal{A})$ is nonsingular.
   • The difficulty of Round off errors in Pade’s method can be controlled by Scaling and squaring.
   • Since $e^\mathcal{A}$ is same as $(e^\frac{\mathcal{A}}{m})^m$, the idea is it to find out $m$ in terms of powers of $2$ so that $e^\frac{\mathcal{A}}{m}$ can be computed accurately and then find $(e^\frac{\mathcal{A}}{m})^m$ with repetitive squaring.
   • This method has favorable numerical properties when $p=q$.
     
     Let $m = 2^j$ be chosen such that $||A||_{\infty} \leq 2^{j-1}$, then Moler and Van Loan(1978) have shown that there exists an $E$ such that ,
     $$\boxed{[R_{pp}(A/2^j)]^{2^j} = e^{\mathcal{A} + E}},$$
     Where $||E||_{\infty} \leq \epsilon ||\mathcal{A}||_{\infty}$, with
     $$\epsilon = 2^{3-2p} (\frac{(p!)^2}{(2p)!(2p+1)!)}) \\ \boxed{\epsilon = \frac{8}{2^{2p}} (\frac{(p!)^2}{(2p)!(2p+1)!)})}$$
     Given an error tolerance of $\delta$ we should choose $p$ such that $||E||_{\infty} \leq \delta ||\mathcal{A}||_{\infty}$
     
     The Algorithm is as follows:
     Input : $A \in \mathcal{R}^{n \times n}$, $\delta > 0$
     Output : $F = e^{\mathcal{A} + E}$ with $||E||_{\infty} \leq \delta ||\mathcal{A}||_{\infty}$
     
     Step 1: Choose $j$ such that $||\mathcal{A}||_{\infty} \leq 2^{j-1}$, set $\mathcal{A} \equiv \mathcal{A}/2^j$.
     Step 2: Choose $p>0$ and smallest non-negative integer satisfying,
     $$\boxed{(\frac{8}{2^{2p}}) (\frac{(p!)^2}{(2p)!(2p+1)!)}) \leq \delta}.$$
     Step 3: Set $\boxed{D \equiv I, N \equiv I, Y \equiv I, c = 1}$
     Step 4: For $k = 1, 2, ..., p$ do
     $$\boxed{c \equiv c(p - k + 1)/[(2p - k + 1)k] \\ Y \equiv AY,\\ N \equiv N + cY , \\D = D+(-1)^{k}cY.}$$
     Step 5: Solve for $F$ : $\boxed{DF = N}$.
     Step 6: For $k = 1, 2, ..., j$ do,
     $$\boxed{F \equiv F^{2}}$$
   • The algorithm requires about $2(p + j + (1/3)n^3$ flops.
   • The algorithm is numerically stable when A is normal. When A is non-normal, an analysis of the stability property becomes difficult, because, in this case e a may grow before it decays during the squaring process; which is known as “hump” phenomenon.
     Example at Page 143 of [1]

3. We have another effective method called the Matrix decomposition method uses the Real Schur form to do so.
   
   Let $\mathcal{A}$ be transformed to a real Schur matric $R$ using an orthogonal similarity transformation:
   $$\boxed{P^{T}\mathcal{A}P = R}$$
   Then,
   $$\boxed{e^{\mathcal{A}}=Pe^{R}P^{T}}$$
   Note that $R$ is upper-triangular.
   
   The algorithm is:
   Step 1: Transform $\mathcal{A}$ to $R$ using the $QR$ iteration algorithm.
   $$\boxed{P^{T}\mathcal{A}P = R}$$
   Step 2: Compute $e^R = G = g_{ij}$:
   for $i = 1,...,n$ do,
   $$g_{ii} = e^{r_{ii}}$$
   end.
   
   for $k = 1, 2, ..., n-1$ do
   for $i = 1, 2, ..., n-k$ do
   set $j = i+k$
   $$\boxed{g_{ij}=\frac{1}{(r_{ii}-r_{jj})}\Big[r_{ij}(g_{ii}-g_{jj})+\sum_{p=i+1}^{j-1}(g_{ip}r_{pj}-r_{ip}g_{pj})\Big]}$$
   end
   end
   Step 3: Now we have $e^R$ then we compute $e^\mathcal{A}$

---

2.Lyapunov Stability theorem: [1]

Lyapunov Stability theorem,
The system :
$$\boxed{\dot{x}(t)=\mathcal{A}x(t)}$$
is asymptotically stable $iff$, for any symmetric positive definite matrix $M$, there exists a unique symmetric positive definite matrix $X$ satisfying the equation:
$$\boxed{X\mathcal{A}+\mathcal{A}^{T}X=-M}$$

Kronecker Product and Lyapunov:
$$\boxed{\mathcal{A} \bigotimes B = \big[ a_{ij} [B] \big]}$$
For example,
if
$$\mathcal{A} = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{bmatrix}$$
and,
$$B = \begin{bmatrix} b_{11} &b_{12} &b_{13}\\ b_{21} &b_{22} &b_{23}\\ \end{bmatrix}$$
Then, the Kronecker product will look like,
$$\mathcal{A} \bigotimes B = \begin{bmatrix} a_{11}b_{11} &a_{11}b_{12} &a_{11}b_{13} &a_{12}b_{11} &a_{12}b_{12} &a_{12}b_{13} \\ a_{21}b_{21} &a_{21}b_{22} &a_{21}b_{23}&a_{22}b_{11} &a_{22}b_{12} &a_{22}b_{13} \\ \end{bmatrix}$$
It is evident that,
$$\boxed{\mathcal{A} \bigotimes B \neq B\bigotimes\mathcal{A}}$$
Now, we formulate the $X\mathcal{A}+\mathcal{A}^{T}X=M$ using the Kronecker product. Note that $M$ is negative semi-definite. $X$ is positive definite and symmetric.

Vectorisation of a matrix $\large \triangleq$ Stack all elements of
the matrix into one column.

Going directly to the final answer,
$$\boxed{X\mathcal{A}+\mathcal{A}^{T}X=M = [(A^{T} \bigotimes I )+( I\bigotimes A^{T})]\overrightarrow{X} = \overrightarrow{M}}$$
Verify it once.

If $\mathcal{A}$ is stable then $X$ has unique solution.


Schur method of Lyapunov:
Method proposed by Bartels and Stewart. It is based on the reduction of $A^{T}$ to $RSF$.

Step 1: Reduce the problem.
For,
$$X\mathcal{A}+\mathcal{A}^{T}X = C$$
is reduced to,
$$\boxed{YR^{T}+RY=\hat{C}}$$
Where,
$R = U^{T}A^{T}U$, $\hat{C}=U^{T}CU$ and $Y=U^{T}XU$.

Step 2: Solution of the reduced form.
Let,
$$Y = (y_1, y_2,...,y_n),$$
$$\hat{C}=(c_1,...,c_n)$$
and $$R=(r_{ij})$$
Assume that the columns $y_{k+1}$ through $y_{n}$ have been computed, and consider the following two cases,

Case 1: $r_{kk-1} = 0$. Then $y_k$ is determined by solving the quasi-triangular system,
$$\boxed{(R+r_{kk}I)y_k=c_k-\sum_{j=k+1}^nr_{kj}y_j}$$
If, in particular, R is upper triangular, that is there are no “Schur bumps” on the diagonal, then each $y_i$, $i=n, n - 1 ..... 2, 1$ can be obtained
by solving an $n \times n$ upper triangular system as follows:
$$\boxed{(R+r_{nn}I)y_n=c_n} \\ \boxed{(R+r_{n-1,n-1}I)y_{n-1}=c_{n-1}-r_{n-1,n}y_n}$$

Case 2: $r_{kk-1} \neq 0$ this shows a Schur Bump.
$$\boxed{R(y_{k-1}, y_k)+(y_{k-1}, y_k) \begin{bmatrix} r_{k-1,k-1} &r_{k,k-1}\\ r_{k-1,k} &r_{k,k}\\ \end{bmatrix}\\ =(c_{k-1},c_k)-\sum_{j=k+1}^n(r_{k-1,j}y_j, r_{kj}y_j)=(d_{k-1},d_k)}$$
Example at Page 274


If $\mathcal{A}$ doesn’t have Eigen Values in the Limit circle , we modify $\mathcal{A}$ as,
$$A_{1} = (I-A)^{-1}(I+A) \text{ and } Q_{1} = (I+A^{T}_{1})(I+A_{1})\\ \text { solve for A } \implies (I-A)A_{1}=(I+A)\\ \implies A_{1}-AA_{1}=I+A\\ \implies A_{1} - I = AA_{1}+A\\ \implies A(A_{1}+I) = A_{1}-I\\ \implies A=(A_{1}-I)(A_{1}+I)^{-1}$$
Use this $A$ in the Discrete Lyapunov which is given by,
$$\mathcal{L}_{D} : A^{T}PA-P=Q$$
The resulting equation is,
$$A^{T}_1P+PA_{1}=Q_{1} \text{ assume, } Q = -2I$$
You will get a similar equation to that of the normal Lyapunov equation.

---

3. Modified QR algorithms:[4]

1. Modified Gram-Schmidt method:
   for
   $$k = 1:n\\ q^{T}_k = a^{T}_k$$
   end
   
   for
   $$k = 1:n\\ r_{kk} = (q^{T}_{k}q_k)^{1/2}\\ q^{T}_k= \frac{q^{T}_k}{r_{kk}}$$
   -for
   $$j = k+1:n\\ r_{kj}=a^{T}_{j}q_{k}\\ q^{T}_{j}=q^{T}_{j}-r_{kj}q^{T}_{k}$$
   -end
   end
   It requires $\boxed{O(\frac{8}{3}n^{3})}$ operations.
   Practice example from the above link

2. Original algorithm of QR decomposition:
   Aim is to write matrix $\mathcal{A}$ in terms of $\boxed{Q\times R}$.
   Here, $Q$ is an orthogonal matrix with orthonormal columns; and $R$ is an upper-triangular matrix.
   
   Now, we take an example,
   $$\mathcal{A}=\begin{bmatrix} 1 & 2 & 3\\ -1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1 \end{bmatrix}$$
   We will first find $q_1$, $q_2$ and $q_3$ as the three orthonormal vectors using the GS-strategy.
   Now, we will define $a$, $b$ and $c$ as,
   $$a = \begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, b=\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix} \text{and}, c = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}$$
   We define a partial basis as $\mathcal{X}$ which will be updated at each step.
   
   For the first step,
   we find $v_1$ as,
   $$\boxed{v_1 = \frac{a}{1} = \begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}} \text{now the partial basis is,} \mathcal{X} = \big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix} \big]$$
   similarly, we find $v_2$ as follows,
   $$v_2 = b - proj_{\mathcal{X}}\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}\\ = \begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix} = \begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{3}{4}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}\\ \text{Now,} \boxed{4v_2=\begin{bmatrix} 8\\ 4\\ 4\\ 4 \end{bmatrix}-\begin{bmatrix} 3\\ -3\\ 3\\ 3 \end{bmatrix}\\ =\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}$$
   Now, we update $\mathcal{X}$,
   $$\mathcal{X} = \Big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix} \Big]$$
   Now, we find $v_3$ as follows,
   $$v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix} - proj_{\mathcal{X}}\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}\\ v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}\\ v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{4}{4}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}-\frac{24}{76}\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}\\ \text{now},\\ \boxed{19v_3 = \begin{bmatrix} 8\\ -4\\ -6\\ -6 \end{bmatrix}}$$
   Hence, the partial basis is,
   $$\boxed{\mathcal{X} = \Big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 8\\ -4\\ -6\\ -6 \end{bmatrix}\Big]}$$

$$proj_{\overrightarrow{a}}{\overrightarrow{b}} = \frac{{\overrightarrow{a}}\cdot {\overrightarrow{b}}}{{\overrightarrow{a}}^2}{\overrightarrow{a}}$$

This partial basis is the orthogonal set, hence we now make it orthonoamal,
$$Q_{orthogonal} = \begin{bmatrix} 1 & 5 & 8\\ -1 & 7 & -2\\ 1 & 1 & -3\\ 1 & 1 & -3\\ \end{bmatrix}$$
Hence,
norm of first column, $$||n_1|| = \sqrt{1^2+(-1)^2+1^2+1^2} = 2$$
Similarly,
$$||n_2|| = \sqrt{76}$$
and
$$||n_2|| = \sqrt{38}$$
Hence, the orthonormal $Q$ is,
$$Q_{orthonormal}=\begin{bmatrix} 1/2 & 5/\sqrt{76} & 8/\sqrt{38} \\ -1/2 & 7/\sqrt{76} & -2/\sqrt{38}\\ 1/2 & 1/\sqrt{76} & -3/\sqrt{38}\\ 1/2 & 1/\sqrt{76} & -3/\sqrt{38}\\ \end{bmatrix}$$
Now, we have, $\mathcal{A}$ and $Q$; we now find $R$; which by theory has to be upper-triangular,
$$\mathcal{A}= Q \cdot R \\ Q^{T} \cdot \mathcal{A} = Q^{T} \cdot Q \cdot R\\ \text{NOTE:}\text{because } Q \text{ is orthonormal } \boxed{Q^{T} \cdot Q = I}.\\ \therefore \boxed{R = Q^{T} \cdot \mathcal{A}}$$
Make up the algorithm from the above steps.

---

4. Condition Number:[5]

As a rule of thumb, if the condition number $\kappa (A)=10^{k}$, then you may lose up to $k$ digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods.

A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned. The condition number is a property of the problem.

Condition number is mathematically defined as,
$$\kappa(A) = ||A||_{1}||A^{-1}||_{1}$$
We define the $1$-norm as follows,
For a matrix $M$
$$M = \begin{bmatrix} 2 & 3\\ 1 & -1 \end{bmatrix}$$
The $1$-norm is,
$$\text{Maximum of }(2+1,3+1) = 4.\\ \text{Hence, }||M||_{1} = 4.$$

So, for a given matrix,
$$P = \begin{bmatrix} 1 & 10^{4}\\ -1 & 2 \end{bmatrix}$$
We get,
$$P^{-1} = \frac{1}{2+10^{4}}\begin{bmatrix} 2 & -10^{4}\\ 1 & 1 \end{bmatrix}$$
Hence,
$$||P||_{1} = \text{maximum of }(1+1, 2+10^{2}) = 2+10^{4}$$
and,
$$||P^{-1}||_{1} = \frac{1}{2+10^{4}}\text{maximum of }(2+1, 10^{4}+1) = \\ \frac{1}{2+10^{4}}(1+10^{4})$$
Therefore,
$$\text{Condition number},\\ \kappa(P) = ||P||_{1}||P^{-1}||_{1} = 1+10^{4}$$
Thus, the matrix (or system) is ill-conditioned.

---

The question paper that came for this exam was pretty sadistic.
Give it a try, if you can and tell me how awesome it feels. I will score around 30/50.