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