The Linear Quadratic Regulator

Optimal Control and Linear-Quadratic-Regulator (LQR)

Today, I will not write an introductory passage to write off my blog. Because, writing an introduction to Optimal Control in itself will required a blog. However, I will add in small tidbits as and when needed.
To understand the topic, we need some basic definitions with us.

---

1.
A control system can be represented in terms of State Space, as follows,
$$\boxed{\dot{x}(t) = A(t)x(t)+B(t)u(t) \\ y(t)=C(t)x(t)+D(t)u(t)}$$
In the above formulation,
$x(.)$ is the state vector; $x(t) \in \mathcal{R}^{n}$.
$y(.)$ is the output vector; $y(t) \in \mathcal{R}^{q}$.
$u(.)$ is the input vector; $u(t) \in \mathcal{R}^{p}$.
$A(.)$ is the System Matrix; $\mathrm{dim}[A(.)]=n\times n$.
$B(.)$ is the Input Matrix; $\mathrm{dim}[B(.)]=n\times p$.
$C(.)$ is the Output Matrix; $\mathrm{dim}[C(.)]=q\times n$.
$A(.)$ is the Feed-forward Matrix; $\mathrm{dim}[D(.)]=q\times p$.
Now, for a system to be controllable, we first define a matrix $\mathcal{Q_c}$, called the controllability matrix, such that,
$$\boxed{\mathcal{Q}_c=\begin{bmatrix} B & AB & A^{2}B & \ldots & A^{n-1}B \\ \end{bmatrix}}$$
The system is controllable if $\mathcal{Q}_c$ has full row rank (i.e. rank($\mathcal{Q}_c$) $=n$).

---

We will assume that we deal with Controllable systems only.

---

Usually, a single input system’s state feedback controller is designed using the Eigen-value method, or Pole Placement method.

---

2.
Pole placement method is the methodology of finding the control vector $\mathcal{U}$ in the form $-\mathcal{k}x$
So, the state space representation changes as,
$$\dot{x}=(A-Bk)x$$
$k$ is found as,
$$|sI-(A-Bk)|=(s-\mu_1)(s-\mu_2)...(s-\mu_n)$$
Here, $\mu$ are the desired pole locations. Note that $k$ is defined as $k=\begin{bmatrix} k_1 & k_2 & k_3 & \ldots & k_n\ \end{bmatrix}$

---

However, for a multi-input system the feedback gain i.e. $k$ is not unique.
Linear Quadratic Control strategy is used to deal with this issue.

---

Now, we dive into the Linear Quadratic Regulator (LQR) formulation, for an $m$-input and $n$-state system with $x \in \mathcal{R}^n$,$u \in \mathcal{R}^m$. Consider a system,
$$\dot{x}=A(t)x(t)+B(t)u(t) \text{ provided }x(0)=x_0 \tag*{(1)}$$
Our aim is to find an open loop control $u(\tau)$, for $\tau \in [t_0, t_f]$ such that we minimize:
$$\boxed{J(u, x_0, t_0, t_f) = \int_{t_0}^{t_f}[x^{T}(t)Q(t)x(t)+u^{T}(t)R(t)u(t)]dt+x(t_f)^{T}Sx(t_f) } \tag*{(2)}$$
where $Q(t)$ and $S$ are symmetric positive semi-definite $n \times n$ matrices.
$R(t)$ is a symmetric positive definite $m \times m$ matrix. Note that $x_0$, $t_0$ and $t_f$ are fixed and given data.
The controller aim is to basically keep $x(t)$ close to 0 especially at $t_f$, which is the final time.
In $(2)$,

• $x^T(t)Q(t)x(t)$ works against the transient response.
• $x^T(t_f)Sx(t_f)$ works against the finite state.
• $u^T(t)R(t)u(t)$ works against the control effort.

The above formulation can regulate the output $y(t) = C(t)x(t)$ near $0$.
Note that, we can define, $S$ and $Q(t)$ as $C^T(t)W(t)C(t)$ where, $W(t) \in R^{r \times r}$
We can now have a theorem as follows,

For a system with fixed initial and final conditions, $\dot{x}=f(x,u,t)$; and clearly $x(t_0)=x_0$. We define our time horizon as $[t_0,t_f]$ such that $t \in [t_0,t_f]$. We find $u(t)$ such that our cost function, $J$ is minimized. $J$ is defined as,
$$J(u(.),x_0)= \phi(x(t_f))+\int_{t_0}^{t_f}L(x(t), u(t), t)dt$$
Here, the first term of $J$ is the final cost and the second term is the recurring cost.

---

Now, we will formulate some important functions that will convert the $J$ which is a constrained optimal control problem to a unconstrained optimal control problem. [THIS MAY NOT MAKE SENSE TO YOU, WHICH IS NATURAL. HOLD ON].
$$\boxed{\dot{\lambda}=-H_x=-\frac{\partial{L}}{\partial{x}}-\lambda^T\frac{\partial{f}}{\partial{x}}} \tag*{(3)}$$
Note that, $\lambda(t)$ ( $\in R^n$ ) is called the Lagrangian.
$H$ is the Hamiltonian operator. Defined in terms of $L$ and $\lambda$ as in $(3)$. Or it can be defined as,
$$\boxed{H(x,u,t) := L(x,u,t)+\lambda^T(t)f(x,u,t)} \tag*{(4)}$$
The above definition is in terms of $L$ as defined in the theorem. So, we define $\lambda$ in the same lines. Just for convenience of computation.
$$\boxed{\lambda^T(t_f)=\frac{\partial \phi}{\partial x}(x(t_f))} \tag*{(5)}$$
$(5)$ can be written as $\lambda^T(t_f)=\phi_x(x(t_f))$
Equation $(3)$,$(4)$ and $(5)$ together form a set of $2n$ differential equations (in $x$ and $\lambda$, obviously) with split boundary conditions at $t_0$ and $t_f$. Now, we can easily define $u(t)$ in terms of $x$ or/and $\lambda$.
As mentioned earlier, the solution is found by converting $J$ from a constrained optimal problem to a constrained optimal problem using a Lagrange multiplier function $\lambda(t)$:
$$\boxed{\bar{J}(u,x_0)=J(u(.),x_0)+\int_{t_0}^{t_f} \lambda^T [f(x,u,t)-\dot{x}]dt} \tag*{(6)}$$
Notice that,
$$\frac{d}{dt}(\lambda^T(t)\dot{x}(t))=\dot{\lambda}^T(t)x(t)+\lambda^T(t)\dot{x} \tag*{(7)}$$
Therefore,
$$\int_{t_0}^{t_f}\lambda^T\dot{x}\ dt=\lambda^T(t_f)\dot{x}(t_f)-\lambda^T(t_0)\dot{x}(t_0)-\int_{t_0}^{t_f}\dot{\lambda}^Tx\ dt \tag*{(8)}$$
As the Hamiltonian Function is defined in $(4)$, thus,
$$\boxed{\bar{J}=\phi(x(t_f))-\lambda^T(t_f)x(t_f) + \lambda^T(t_0)x(t_0)+\int_{t_0}^{t_f}[H(x(t),u(t),t)+\dot{\lambda}(t)x(t)] \ dt} \tag*{(9)}$$
The necessary condition for an optimal solution is $\delta \bar{J}$ of the modified cost with respect to all variations of the system be minimal at all times from $t_0$ to $t_f$.
We will define $\delta \bar{J}$ analytically in the next post and formulate the Riccati Equation that will lay the foundation to some amazing control strategies.
Cheers!

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.