# MATLAB TUTORIAL, part 2.1: Polynomial Interpolation

Polynomial Interpolation

Originally, spectral decomposition was developed for symmetric or self-adjoint matrices. Following tradition, we present this method for symmetric/self-adjoint matrices, and later expand it for arbitrary matrices.

A matrix A is said to be unitary diagonalizable if there is a unitary matrix U such that $${\bf U}^{\ast} {\bf A}\,{\bf U} = {\bf \Lambda} ,$$ where Λ is a diagonal matrix and $${\bf U}^{\ast} = {\bf U}^{-1} .$$ A matrix A is said to be orthogonally diagonalizable if there is an orthogonal matrix P such that $${\bf P}^{\mathrm T} {\bf A}\,{\bf P} = {\bf \Lambda} ,$$ where Λ is a diagonal matrix and $${\bf P}^{\mathrm T} = {\bf P}^{-1} .$$

Theorem: The matrix A is orthogonally diaginalisable if and only if A is symmetric ($${\bf A} = {\bf A}^{\mathrm T}$$ ). ■

Theorem: The matrix A is unitary diaginalisable if and only if A is normal ($${\bf A}\, {\bf A}^{\ast} = {\bf A}^{\ast} {\bf A}$$ ). ■

Example: The matrix

${\bf A} = \begin{bmatrix} 1&{\bf j}&0 \\ {\bf j}&1&0 \\ 0&0&1 \end{bmatrix}$
is symmetric, normal, but not self-adjoint. Another matrix
${\bf B} = \begin{bmatrix} 2&1 \\ -1&2 \end{bmatrix}$
is normal, but not self-adjoint. Therefore, both matrices are unitary diagonalizable but not orthogonally diagonalizable.

Example. Consider a symmetric matrix

${\bf A} = \begin{bmatrix} 2&1&1 \\ 1&2&1 \\ 1&1&2 \end{bmatrix}$
that has the characteristic polynomial $$\chi_{A} (\lambda ) = \det \left( \lambda {\bf I} - {\bf A} \right) = \left( \lambda -1 \right)^2 \left( \lambda -4 \right) .$$ Thus, the distinct eigenvalues of A are $$\lambda_1 =1,$$ which has geometrical multiplicity 2, and $$\lambda_3 =4.$$ The corresponding eigenvectors are
${\bf u}_1 = \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix} , \quad {\bf u}_2 = \begin{bmatrix} 1 \\ -1 \\ 0 \end{bmatrix} \qquad \mbox{and} \qquad {\bf u}_3 = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} .$
The vectors $${\bf u}_1 , \ {\bf u}_2$$ form the basis for the two-dimensional eigenspace corresponding $$\lambda_1 =1 ,$$ while $${\bf u}_3$$ is the eigenvectors corresponding to $$\lambda_3 =4 .$$ Applying the Gram--Schmidt process to $${\bf u}_1 , \ {\bf u}_2$$ yields the following orthogonal basis:
${\bf v}_1 = {\bf u}_1 = \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix} \quad \mbox{and} \quad {\bf v}_2 = {\bf u}_2 - \frac{\langle {\bf u}_2 , {\bf v}_1 \rangle}{\| {\bf v}_1 \|^2} \, {\bf v}_1 = \frac{1}{2} \begin{bmatrix} 1 \\ -2 \\ 1 \end{bmatrix}$
because $$\langle {\bf u}_2 , {\bf u}_1 \rangle = -1$$ and $$\| {\bf v}_1 \|^2 =2 .$$ Normalizing these vectors, we obtain orthonormal basis:
${\bf q}_1 = \frac{{\bf v}_1}{\| {\bf v}_1 \|} = \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix} , \quad {\bf q}_2 = \frac{{\bf v}_2}{\| {\bf v}_2 \|} = \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\ -2 \\ 1 \end{bmatrix} , \quad {\bf q}_3 = \frac{{\bf v}_3}{\| {\bf v}_3 \|} = \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} .$
Finally, using $${\bf q}_1 ,\ {\bf q}_2 , \ {\bf q}_3$$ as column vectors, we obtain the unitary matrix
${\bf U} = \begin{bmatrix} \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} \\ 0& \frac{-2}{\sqrt{6}} & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} \end{bmatrix} ,$
which orthogonally diagonalizes A. As a check, we confirm
${\bf U}^{\mathrm T} {\bf A} \,{\bf U} = \begin{bmatrix} \frac{-1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & \frac{-2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \\ \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} \end{bmatrix} \, \begin{bmatrix} 2&1&1 \\ 1&2&1 \\ 1&1&2 \end{bmatrix} \, \begin{bmatrix} \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} \\ 0& \frac{-2}{\sqrt{6}} & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} \end{bmatrix} = \begin{bmatrix} 1 &0&0 \\ 0 &1&0 \\ 0 &0&4 \end{bmatrix} ,$

Theorem: Let A be a symmetric or normal $$n \times n$$ matrix with eigenvalues $$\lambda_1 , \ \lambda_2 , \ \ldots , \ \lambda_n$$ and corresponding eigenvectors $${\bf u}_1 , \ {\bf u}_2 , \ \ldots , \ {\bf u}_n .$$ Then

${\bf A} = \begin{bmatrix} \uparrow & \uparrow & \cdots & \uparrow \\ {\bf u}_1 & {\bf u}_2 & \cdots & {\bf u}_n \\ \downarrow & \downarrow & \cdots & \downarrow \end{bmatrix} \, \begin{bmatrix} \lambda_1 &&&0 \\ &\lambda_2 && \\ &&\ddots & \\ 0&&& \lambda_n \end{bmatrix} \, \begin{bmatrix} \longleftarrow & {\bf u}_1 & \longrightarrow \\ \longleftarrow & {\bf u}_2 & \longrightarrow & \\ \vdots & \\ \longleftarrow & {\bf u}_n & \longrightarrow \end{bmatrix} = \sum_{i=1}^n \lambda_i {\bf u}_i {\bf u}_i^{\ast} ,$
which is called spectral decomposition for a symmetric/ normal matrix A. ■
The term was cointed around 1905 by a German mathematician David Hilbert (1862--1942). Denoting the rank one projection matrix $${\bf u}_i {\bf u}_i^{\ast}$$ by $${\bf E}_i = {\bf u}_i {\bf u}_i^{\ast} ,$$ we obtain spectral decomposition of A:
${\bf A} = \lambda_1 {\bf E}_1 + \lambda_2 {\bf E}_2 + \cdots + \lambda_n {\bf E}_n .$
This formular allows one to define a function of a square symmetric/self-adjoint matrix:
$f\left( {\bf A} \right) = f(\lambda_1 )\, {\bf E}_1 + f(\lambda_2 )\, {\bf E}_2 + \cdots + f(\lambda_n )\,{\bf E}_n$
because the matrices $${\bf E}_k , \quad k=1,2,\ldots n,$$ are projection matrices:
${\bf E}_i {\bf E}_j = \delta_{i,j} {\bf E}_i = \begin{cases} {\bf E}_i , & \mbox{ if } i=j, \\ {\bf 0} , & \mbox{ if } i \ne j , \end{cases} \qquad i,j =1,2,\ldots n.$

Example. Consider a self-adjoint 2-by-2 matrix

${\bf A} = \begin{bmatrix} 1 &2+{\bf j} \\ 2- {\bf j} & 5\end{bmatrix} ,$
where $${\bf j}^2 =-1 .$$ We have $${\bf S}^{\ast} {\bf A} {\bf S} = {\bf S}^{-1} {\bf A} {\bf S} = {\bf \Lambda}$$ for the matrices
${\bf S} = \begin{bmatrix} \left( 2+{\bf j} \right) / \sqrt{6} & \left( 2+{\bf j} \right) / \sqrt{30} \\ - 1/\sqrt{6} & 5\sqrt{30} \end{bmatrix} \quad \mbox{and} \quad {\bf \Lambda} = \begin{bmatrix} 0&0 \\ 0& 6 \end{bmatrix} .$
Here column vectors in matrix S are normalized eigenvectors corresponding eigenvalues $$\lambda_1 =0 \quad \mbox{and} \quad \lambda_2 =6 .$$

Therefore, the spectral decomposition of A becomes $${\bf A} = 0\,{\bf E}_1 + 6\,{\bf E}_2 ,$$ which is clearly matrix A itself.

In our case, projection matrices are

${\bf E}_1 = \frac{1}{6} \begin{bmatrix} 5 & -2 - {\bf j} \\ -2+{\bf j} & 1 \end{bmatrix} , \qquad {\bf E}_2 = \frac{1}{6} \begin{bmatrix} 1 &2+{\bf j} \\ 2- {\bf j} & 5\end{bmatrix} = \frac{1}{6}\, {\bf A} .$
It is easy to check that
${\bf E}_1^2 = {\bf E}_1 , \qquad {\bf E}_2^2 = {\bf E}_2 , \qquad \mbox{and} \qquad {\bf E}_1 {\bf E}_2 = {\bf 0} .$
The exponential matrix-function is
$e^{{\bf A}\,t} = {\bf E}_1 + e^{6t} \,{\bf E}_2 .$

Example. Consider a symmetric 3-by-3 matrix from the previous example

${\bf A} = \begin{bmatrix} 2&1&1 \\ 1&2&1 \\ 1&1&2 \end{bmatrix} .$
Its spectral decomposition is $${\bf A} = 1\,{\bf E}_1 + 1\,{\bf E}_2 + 4\,{\bf E}_3 ,$$ where the projection matrices $${\bf E}_i = {\bf q}_i {\bf q}_i^{\ast}$$ are obtained from the orthonormal eigenvectors:
\begin{align*} {\bf E}_1 &= \frac{1}{6} \begin{bmatrix} -1 \\ -1 \\ 2 \end{bmatrix} \left[ -1 \ -1 \ 2 \right] = \frac{1}{6} \begin{bmatrix} 1&1& -2 \\ 1&1& -2 \\ -2&-2& 4 \end{bmatrix} , \\ {\bf E}_2 &= \frac{1}{2} \begin{bmatrix} -1 \\ 1 \\ 0 \end{bmatrix} \left[ -1 \ 1 \ 0 \right] = \frac{1}{2} \begin{bmatrix} 1&-1&0 \\ -1&1&0 \\ 0&0&0 \end{bmatrix} , \\ {\bf E}_3 &= \frac{1}{3} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \left[ 1 \ 1 \ 1 \right] = \frac{1}{3} \begin{bmatrix} 1&1& 1 \\ 1&1& 1 \\ 1&1& 1 \end{bmatrix} . \end{align*}
Indeed, these diagonalizable matrices satisfy the following relations:
${\bf E}_1^2 = {\bf E}_1 , \quad {\bf E}_2^2 = {\bf E}_2 , \quad {\bf E}_3^2 = {\bf E}_3 , \quad {\bf E}_1 {\bf E}_2 = {\bf 0} , \quad {\bf E}_1 {\bf E}_3 = {\bf 0} , \quad {\bf E}_3 {\bf E}_2 = {\bf 0} ,$
and they all have eigenvalues $$\lambda = 1, 0, 0 .$$ Using this spectral decomposition, we define two matrix-functions corresponding to $${\Phi}(\lambda ) = \cos \left( \sqrt{\lambda} \,t \right)$$ and $${\Psi}(\lambda ) = \frac{1}{\sqrt{\lambda}} \,\sin \left( \sqrt{\lambda} \,t \right)$$ that do not depend on the branch of the choisen square root:
\begin{align*} {\bf \Phi} (t) &= \cos \left( \sqrt{\bf A} \,t \right) = \cos t\, {\bf E}_1 + \cos t\, {\bf E}_2 + \cos (2t) \,{\bf E}_3 = \frac{\cos t}{3} \, \begin{bmatrix} 2&-1&-1 \\ -1&2&-1 \\ -1&-1& 2 \end{bmatrix} + \frac{\cos 2t}{3} \begin{bmatrix} 1&1& 1 \\ 1&1& 1 \\ 1&1& 1 \end{bmatrix} , \\ {\bf \Psi} (t) &= \frac{1}{\sqrt{\bf A}} \,\sin \left( \sqrt{\bf A} \,t \right) = \sin t\, {\bf E}_1 + \sin t\, {\bf E}_2 + \frac{\sin (2t)}{2} \,{\bf E}_3 = \frac{\sin t}{3} \, \begin{bmatrix} 2&-1&-1 \\ -1&2&-1 \\ -1&-1& 2 \end{bmatrix} + \frac{\sin 2t}{6} \begin{bmatrix} 1&1& 1 \\ 1&1& 1 \\ 1&1& 1 \end{bmatrix} . \end{align*}
These matrix-functions are solutions of the following initial value problems for the second order matrix differential equations:
\begin{align*} & \ddot{\bf \Phi}(t) + {\bf A}\,{\bf \Phi} (t) ={\bf 0} , \qquad {\bf \Phi}(0) = {\bf I}, \quad \dot{\bf \Phi}(0) = {\bf 0}, \\ &\ddot{\bf \Psi}(t) + {\bf A}\,{\bf \Phi} (t) ={\bf 0} , \qquad {\bf \Phi}(0) = {\bf 0}, \quad \dot{\bf \Psi}(0) = {\bf I} . \end{align*}
Since the given matrix A is positive definite, we can define four square roots:
\begin{align*} {\bf R}_1 &= {\bf E}_1 + {\bf E}_2 + 2\,{\bf E}_3 = \frac{1}{3} \begin{bmatrix} 4&1&1 \\ 1&4&1 \\ 1&1&4 \end{bmatrix} , \\ {\bf R}_2 &= {\bf E}_1 + {\bf E}_2 - 2\,{\bf E}_3 = \begin{bmatrix} 0&-1&-1 \\ -1&0&-1 \\ -1&-1&0 \end{bmatrix} , \\ {\bf R}_3 &= {\bf E}_1 - {\bf E}_2 + 2\,{\bf E}_3 = \frac{1}{3} \begin{bmatrix} 1&4&1 \\ 4&1&1 \\ 1&1&4 \end{bmatrix} , \\ {\bf R}_4 &= -{\bf E}_1 + {\bf E}_2 + 2\,{\bf E}_3 = \begin{bmatrix} 1&0&1 \\ 0&1&1 \\ 1&1&0 \end{bmatrix} , \end{align*}

and four others are just negative of these four; so total number of square roots is 8. Note that we cannot obtain $${\bf R}_3$$ and $${\bf R}_4$$ using neither Sylvester's method nor the Resolvent method because they are based on the minimal polynomial $$\psi (\lambda ) = (\lambda -1)(\lambda -4) .$$

We expand spectral decomposition for arbitary square matrices. Let $$f (\lambda )$$ be an analytic function in a neighborhood of the origin and A be a square $$n \times n$$ matrix. We choose the origin as an example; application of the spectral decomposition requirs functions to be expressed as convergent power series in neighborhoods of every eigenvalue. Using a Maclaurin series

$f(\lambda ) = f_0 + f_1 \lambda + f_2 \lambda^2 + \cdots + f_k \lambda^k + \cdots = \sum_{k\ge 0} f_k \lambda^k ,$
we can define the matrix-valued function $$f ({\bf A} )$$ as
$f({\bf A} ) = \sum_{k\ge 0} f_k {\bf A}^k .$
Let $$\psi (\lambda )$$ be a minimal polynomial of degree m for the matrix A. Then every power $${\bf A}^p$$ of matrix A can be expressed as a polynomial of degree not higher than $$m-1.$$ Therefore,
$f({\bf A} ) = \sum_{j= 0}^{m-1} b_j {\bf A}^j ,$
where coefficients $$b_j , \quad j=0,1,\ldots , m-1;$$ should satisfy the following equations
$f(\lambda_k ) = \sum_{j= 0}^{m-1} b_k \,\lambda_k^j ,$
for each eigenvalue $$\lambda_k , \quad k=1,2,\ldots , s ,$$ where s is the number of distinct eigenvalues. If the eigenvalue $$\lambda_k$$ is of multiplicity $$m_k$$ in the minimal polynomial $$\psi (\lambda ) ,$$ then we need to add $$m_k -1$$ auxiliary algebraic equations
$\left. \frac{{\text d}^p f(\lambda )}{{\text d} \lambda^p} \right\vert_{\lambda = \lambda_k} = \left. \frac{{\text d}^p}{{\text d} \lambda^p} \, \sum_{j= 0}^{m-1} b_k \,\lambda_k^j \right\vert_{\lambda = \lambda_k} , \quad p=1,2,\ldots , m_k -1 .$

Example: Consider a diagonalizable 4-by-4 matrix

${\bf A} = \begin{bmatrix} -4&7&1&4 \\ 6&-16&-3&-9 \\ 12&-27&-4&-15 \\ -18&43&7&24 \end{bmatrix} .$
Its characteristic polynomial is
$\chi (\lambda ) = \det \left( \lambda {\bf I} - {\bf A} \right) = \lambda \left( \lambda -2 \right) \left( \lambda +1 \right)^2 ,$
but its minimal polynomial is a product of simple terms:
$\psi (\lambda ) = \det \left( \lambda {\bf I} - {\bf A} \right) = \lambda \left( \lambda -2 \right) \left( \lambda +1 \right) .$
Therefore, matrix A is diagonalizable because its minimal polynomail is a product of simple terms, and we don't need to find eigenspaces. ■


Complete

Applications