# Eigenvalue Decomposition - Definitions and Facts

## Prerequisites

The reader should be familiar with basic linear algebra concepts.

## Competences

The reader should be able to understand and check the facts about eigenvalue decomposition.

## Selected references

There are many excellent books on the subject. Here we list a few:

J. W. Demmel, Applied Numerical Linear Algebra

G. H. Golub and C. F. Van Loan, Matrix Computations

N. Higham, Accuracy and Stability of Numerical Algorithms

L. Hogben, ed., Handbook of Linear Algebra

B. N. Parlett, The Symmetric Eigenvalue Problem

G. W. Stewart, Matrix Algorithms, Vol. II: Eigensystems

L. N. Trefethen and D. Bau, III, Numerical Linear Algebra

J. H. Wilkinson, The Algebraic Eigenvalue Problem

## General matrices

For more details and the proofs of the Facts below, see L. M. DeAlba, Determinants and Eigenvalues and the references therein.

### Definitions

We state the basic definitions:

Let $\mathbb{F}=\mathbb{R}$ or $F=\mathbb{C}$ and let $A\in \mathbb{F}^{n\times n}$ with elements $a_{ij}\in \mathbb{F}$.

An element $\lambda \in \mathbb{F}$ is an eigenvalue of $A$ if $\exists x\in \mathbb{F}^n$, $x\neq 0$ such that

$$Ax=\lambda x,$$

and $x$ is an eigenvector of $\lambda$.

Characteristic polynomial of $A$ is $p_A(x)=\det(A-xI)$.

Algebraic multiplicty, $\alpha(\lambda)$, is the multiplicity of $\lambda$ as a root of $p_A(x)$.

Spectrum of $A$, $\sigma(A)$, is the multiset of all eigenvalues of $A$, with each eigenvalue appearing $\alpha(A)$ times.

Spectral radius of $A$ is

$$\rho(A)=\max\{|\lambda|, \lambda \in \sigma(A)\}.$$

Eigenspace of $\lambda$ is

$$E_{\lambda}(A)=\ker(A-\lambda I).$$

Geometric multiplicity of $\lambda$ is

$$\gamma(\lambda)=\dim(E_{\lambda}(A)).$$

$\lambda$ is simple if $\alpha(\lambda)=1$.

$\lambda$ is semisimple if $\alpha(\lambda)=\gamma(\lambda)$.

$A$ is nonderogatory if $\gamma(\lambda)=1$ for all $\lambda$.

$A$ is nondefective if every $\lambda$ is semisimple.

$A$ is diagonalizable if there exists nonsingular $B$ matrix and diagonal matrix $D$ such that $A=BDB^{-1}$.

Trace of $A$ is

$$\mathop{\mathrm{tr}}(A)=\sum_i a_{ii}.$$

$Q\in\mathbb{C}^{n\times n}$ is unitary if $Q^*Q=QQ^*=I$, where $Q^*=(\bar Q)^T$.

Schur decomposition of $A$ is $A=QTQ^*$, where $Q$ is unitary and $T$ is upper triangular.

$A$ and $B$ are similar if $B=QAQ^{-1}$ for some nonsingular matrix $Q$.

$A$ is normal if $AA^*=A^*A$.

### Facts

There are many facts related to the eigenvalue problem for general matrices. We state some basic ones:

1. $\lambda\in\sigma(A) \Leftrightarrow p_A(\lambda)=0$.

2. Cayley-Hamilton Theorem. $p_A(A)=0$.

3. A simple eigenvalue is semisimple.

4. $\mathop{\mathrm{tr}}(A)=\sum_{i=1}^n \lambda_i$.

5. $\det(A)=\prod_{i=1}^n \lambda_i$.

6. $A$ is singular $\Leftrightarrow$ $\det(A)=0$ $\Leftrightarrow$ $0\in\sigma(A)$.

7. If $A$ is triangular, $\sigma(A)=\{a_{11},a_{22},\ldots,a_{nn}\}$.

8. For $A\in\mathbb{C}^{n\times n}$, $\lambda\in\sigma(A)$ $\Leftrightarrow$ $\bar\lambda\in\sigma(A^*)$.

9. Corollary of the Fundamental theorem of algebra. For $A\in\mathbb{R}^{n\times n}$, $\lambda\in\sigma(A)$ $\Leftrightarrow$ $\bar\lambda\in\sigma(A)$.

10. If $(\lambda,x)$ is an eigenpair of a nonsingular $A$, then $(1/\lambda,x)$ is an eigenpair of $A^{-1}$.

11. Eigenvectors corresponding to distinct eigenvalues are linearly independent.

12. $A$ is diagonalizable $\Leftrightarrow$ $A$ is nondefective $\Leftrightarrow$ $A$ has $n$ linearly independent eigenvectors.

13. Every $A$ has Schur decomposition. Moreover, $T_{ii}=\lambda_i$.

14. If $A$ is normal, matrix $T$ from its Schur decomposition is normal. Consequently:

• $T$ is diagonal, and has eigenvalues of $A$ on diagonal,
• matrix $Q$ of the Schur decomposition is the unitary matrix of eigenvectors,
• all eigenvalues of $A$ are semisimple and $A$ is nondefective.
15. If $A$ and $B$ are similar, $\sigma(A)=\sigma(B)$. Consequently, $\mathop{\mathrm{tr}}(A)=\mathop{\mathrm{tr}}(B)$ and $\det(A)=\det(B)$.

16. Eigenvalues and eigenvectors are continous and differentiable: if $\lambda$ is a simple eigenvalue of $A$ and $A(\varepsilon)=A+\varepsilon E$ for some $E\in F^{n\times n}$, for small $\varepsilon$ there exist differentiable functions $\lambda(\varepsilon)$ and $x(\varepsilon)$ such that

$$A(\varepsilon) x(\varepsilon) = \lambda(\varepsilon) x(\varepsilon).$$
1. Classical motivation for the eigenvalue problem is the following: consider the system of linear differential equations with constant coefficients, $$\dot y(t)=Ay(t).$$ If the solution is $y=e^{\lambda t} x$ for some constant vector $x$, then $$\lambda e^{\lambda t} x=Ae^{\lambda t} x \quad \textrm{or} \quad Ax=\lambda x.$$

### Examples

We shall illustrate above Definitions and Facts on several small examples, using symbolic computation.



In [1]:

using SymPy




In [2]:

A=[-3 7 -1; 6 8 -2; 72 -28 19]




Out[2]:

3×3 Array{Int64,2}:
-3    7  -1
6    8  -2
72  -28  19




In [3]:

@vars x




Out[3]:

(x,)




In [4]:

using LinearAlgebra
eye(n)=Matrix{Int}(I,n,n)
A-x*eye(3)




Out[4]:

$\left[ \begin{array}{rrr}- x - 3&7&-1\\6&8 - x&-2\\72&-28&19 - x\end{array}\right]$




In [5]:

# Characteristic polynomial p_A(λ)
p(x)=det(A-x*eye(3))
p(x)




Out[5]:

\begin{equation*}26 x + \left(8 - x\right) \left(19 - x\right) \left(- x - 3\right) - 894\end{equation*}




In [6]:

# Characteristic polynomial in nicer form
p(x)=factor(det(A-x*eye(3)))
p(x)




Out[6]:

\begin{equation*}- \left(x - 15\right)^{2} \left(x + 6\right)\end{equation*}




In [7]:

λ=solve(p(x),x)




Out[7]:

$\left[ \begin{array}{r}-6\\15\end{array} \right]$



The eigenvalues are $\lambda_1=-6$ and $\lambda_2=15$ with algebraic multiplicities $\alpha(\lambda_1)=1$ and $\alpha(\lambda_2)=2$.



In [8]:

g=nullspace(map(Float64,A-λ[1]*eye(3)))




Out[8]:

3×1 Array{Float64,2}:
0.23570226039551587
-0.23570226039551587
-0.9428090415820634




In [9]:

h=nullspace(map(Float64,A-λ[2]*eye(3)))




Out[9]:

3×1 Array{Float64,2}:
0.21821789023599256
0.43643578047198495
-0.8728715609439692



The geometric multiplicities are $\gamma(\lambda_1)=1$ and $\gamma(\lambda_2)=1$. Thus, $\lambda_2$ is not semisimple, therefore $A$ is defective and not diagonalizable.



In [10]:

# Trace and determinant
tr(A), λ[1]+λ[2]+λ[2]




Out[10]:

(24, 24)




In [11]:

det(A), λ[1]*λ[2]*λ[2]




Out[11]:

(-1350.0, -1350)




In [12]:

# Schur decomposition
T,Q=schur(A)
T




Out[12]:

3×3 Array{Float64,2}:
-6.0  25.4662       -72.2009
0.0  15.0          -12.0208
0.0   1.48587e-15   15.0




In [13]:

Q




Out[13]:

3×3 Array{Float64,2}:
-0.235702  -0.0571662  -0.970143
0.235702  -0.971825   -5.90663e-16
0.942809   0.228665   -0.242536




In [14]:

# ?schur




In [15]:

F=schur(A)
fieldnames(typeof(F))




Out[15]:

(:T, :Z, :values)




In [16]:

F.Z




Out[16]:

3×3 Array{Float64,2}:
-0.235702  -0.0571662  -0.970143
0.235702  -0.971825   -5.90663e-16
0.942809   0.228665   -0.242536




In [17]:

println(diag(T))




[-6.000000000000002, 14.999999999999988, 14.999999999999988]




In [18]:

Q'*Q




Out[18]:

3×3 Array{Float64,2}:
1.0          1.11022e-16  2.77556e-17
1.11022e-16  1.0          1.52656e-16
2.77556e-17  1.52656e-16  1.0




In [19]:

Q*Q'




Out[19]:

3×3 Array{Float64,2}:
1.0          1.35877e-16  0.0
1.35877e-16  1.0          3.22346e-17
0.0          3.22346e-17  1.0




In [20]:

# Similar matrices
M=rand(-5:5,3,3)
B=M*A*inv(M)
eigvals(B), tr(B), det(B)




Out[20]:

(Complex{Float64}[-6.000000000000157 + 0.0im, 15.00000000000006 - 6.732659804061512e-7im, 15.00000000000006 + 6.732659804061512e-7im], 24.000000000000014, -1349.9999999999982)



### Example

This matrix is nondefective and diagonalizable.



In [21]:

A=[57 -21 21; -14 22 -7; -140 70 -55]




Out[21]:

3×3 Array{Int64,2}:
57  -21   21
-14   22   -7
-140   70  -55




In [22]:

p(x)=factor(det(A-x*eye(3)))
p(x)




Out[22]:

\begin{equation*}- \left(x - 15\right)^{2} \left(x + 6\right)\end{equation*}




In [23]:

λ=solve(p(x),x)




Out[23]:

$\left[ \begin{array}{r}-6\\15\end{array} \right]$




In [24]:

h=nullspace(map(Float64,A-λ[2]*eye(3)))




Out[24]:

3×2 Array{Float64,2}:
0.0       -0.57735
0.707107  -0.57735
0.707107   0.57735




In [25]:

F=schur(A)
F.T




Out[25]:

3×3 Array{Float64,2}:
-6.0  -174.816        -36.5844
0.0    15.0           -8.88178e-16
0.0     7.99361e-14   15.0



### Example

Let us try some random examples of dimension $n=4$ (the largest $n$ for which we can compute eigevalues symbolically).



In [26]:

using Random
Random.seed!(123)
A=rand(-4:4,4,4)




Out[26]:

4×4 Array{Int64,2}:
0  3  -3   3
-2  3   3   1
-3  3  -1  -2
4  0   1   0




In [27]:

p(x)=factor(det(A-x*eye(4)))
p(x)




Out[27]:

\begin{equation*}x^{4} - 2 x^{3} - 25 x^{2} + 30 x + 324\end{equation*}




In [28]:

λ=solve(p(x),x)




Out[28]:

$\left[ \begin{array}{r}\frac{1}{2} - \frac{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}{2} - \frac{\sqrt{\frac{106}{3} - 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}} + \frac{8}{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}} - \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}}{2}\\\frac{1}{2} + \frac{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}{2} - \frac{\sqrt{\frac{106}{3} - 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}} - \frac{8}{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}} - \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}}{2}\\\frac{1}{2} + \frac{\sqrt{\frac{106}{3} - 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}} - \frac{8}{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}} - \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}}{2} + \frac{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}{2}\\\frac{1}{2} + \frac{\sqrt{\frac{106}{3} - 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}} + \frac{8}{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}} - \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}}{2} - \frac{\sqrt{\frac{53}{3} + \frac{4693}{18 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}} + 2 \sqrt[3]{\frac{298871}{216} + \frac{\sqrt{43320759} i}{12}}}}{2}\end{array} \right]$




In [29]:

length(λ)




Out[29]:

4



Since all eigenvalues are distinct, they are all simple and the matrix is diagonalizable. With high probability, all eigenvalues of a random matrix are simple.

Do not try to use nullspace() here.



In [30]:

A=rand(4,4)
p(x)=factor(det(A-x*eye(4)))
p(x)




Out[30]:

\begin{equation*}\frac{3.38186862787728 \left(0.11230475307369 x^{8} - 0.577482168557925 x^{7} + 1.0 x^{6} - 0.650050059179624 x^{5} + 0.0234706828935998 x^{4} + 0.137035547738385 x^{3} - 0.0431687459613977 x^{2} + 0.00168390430737259 x + 6.01501218583027 \cdot 10^{-5}\right)}{0.379799921181417 x^{4} - 1.0 x^{3} + 0.836624012915075 x^{2} - 0.247969313739554 x + 0.0143453492392233}\end{equation*}



In this case, symbolic computation does not work well with floating-point numbers - the degree of $p_A(x)$ is 8 instead of 4.

Let us try Rational numbers:



In [31]:

A=map(Rational,A)




Out[31]:

4×4 Array{Rational{Int64},2}:
2766357556280301//4503599627370496  …  2798440949260299//4503599627370496
3614880418241535//4503599627370496      196003842078327//562949953421312
1251253776444781//2251799813685248     2569813346554741//4503599627370496
529613335220473//562949953421312       918719114247981//4503599627370496




In [32]:

p(x)=factor(det(A-x*eye(4)))
p(x)




Out[32]:

\begin{equation*}\frac{205688069665150755269371147819668813122841983204197482918576128 x^{4} - 516098893015258486191992613103377934939181547116050259862618112 x^{3} + 19556503958273307413339192236881133420393877152071708923920384 x^{2} + 132070114958308006066088439196161677871309827249193117647110144 x + 2916696370953252746150072496481024006168149914567049629396175}{205688069665150755269371147819668813122841983204197482918576128}\end{equation*}




In [33]:

λ=solve(p(x),x)




Out[33]:

$\left[ \begin{array}{r}\frac{11300134159719059}{18014398509481984} - \frac{\sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}{2} - \frac{\sqrt{\frac{367651734667466909083179859605059}{121694457621910022543683507716096} - 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}} - \frac{886575462337898459234562458973361265771256866127}{365375409332725729550921208179070754913983135744 \sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}} - \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}}{2}\\\frac{11300134159719059}{18014398509481984} + \frac{\sqrt{\frac{367651734667466909083179859605059}{121694457621910022543683507716096} - 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}} - \frac{886575462337898459234562458973361265771256866127}{365375409332725729550921208179070754913983135744 \sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}} - \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}}{2} - \frac{\sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}{2}\\\frac{11300134159719059}{18014398509481984} - \frac{\sqrt{\frac{367651734667466909083179859605059}{121694457621910022543683507716096} - 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}} + \frac{886575462337898459234562458973361265771256866127}{365375409332725729550921208179070754913983135744 \sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}} - \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}}{2} + \frac{\sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}{2}\\\frac{11300134159719059}{18014398509481984} + \frac{\sqrt{\frac{367651734667466909083179859605059}{121694457621910022543683507716096} - 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}} + \frac{886575462337898459234562458973361265771256866127}{365375409332725729550921208179070754913983135744 \sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}} - \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}}{2} + \frac{\sqrt{\frac{367651734667466909083179859605059}{243388915243820045087367015432192} + \frac{1031004536818940592052623113611754218686054429798863358692187005}{3702385253972713594848680660754038636211155697675554692534370304 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}} + 2 \sqrt[3]{\frac{241595952437880192075309903440465131856058991805221929547017039828558083201874506011805762843}{7039996334211983914143748299142810924593874768509475682198501392968741771982091656569821331456} + \frac{17 \sqrt{52780695677775735142518817340259915367992921746115475579682933189934552104701294351811171096509941338520401207545834934618082773875552670825799311238498173448668923623452644749787115626} i}{100124392308792660112266642476697755372001774485468098591267575366666549645967525782326347825152}}}}{2}\end{array} \right]$




In [34]:

length(λ)




Out[34]:

4



### Example - Circulant matrix

For more details, see A. Böttcher and I. Spitkovsky, Special Types of Matrices and the references therein.

Given $a_0,a_1,\ldots,a_{n-1}\in \mathbb{C}$, the circulant matrix is

$$C(a_0,a_1,\ldots,a_{n-1})=\begin{bmatrix} a_0 & a_{n-1} & a_{n-2} & \cdots & a_{1} \\ a_1 & a_0 & a_{n-1} & \cdots & a_{2} \\ a_2 & a_1 & a_{0} & \cdots & a_{3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n-1} & a_{n-2} & a_{n-3} & \cdots & a_{0} \end{bmatrix}.$$

Let $a(z)=a_0+a_1z+a_2z^2+\cdots +a_{n-1}z^{n-1}$ be the associated complex polynomial.

Let $\omega_n=\exp\big(\displaystyle\frac{2\pi i}{n}\big)$. The Fourier matrix of order $n$ is

$$F_n=\frac{1}{\sqrt{n}} \bigg[ \omega_n^{(j-1)(k-1)}\bigg]_{j,k=1}^n= \frac{1}{\sqrt{n}} \begin{bmatrix} 1 & 1 & \cdots & 1 \\ 1& \omega_n & \omega_n^2 & \cdots \omega_n^{n-1} \\ 1& \omega_n^2 & \omega_n^4 & \cdots \omega_n^{2(n-1)} \\ \vdots & \vdots & \ddots & \vdots \\ 1& \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots \omega_n^{(n-1)(n-1)} \end{bmatrix}.$$

Fourier matrix is unitary. Fourier matrix is a Vandermonde matrix,

$$F_n=\displaystyle\frac{1}{\sqrt{n}} V(1,\omega_n,\omega_n^2,\ldots, \omega_n^{n-1}).$$

Circulant matrix is normal and, thus, unitarily diagonalizable, with the eigenvalue decomposition

$$C(a_0,a_1,\ldots,a_{n-1})=F_n^*\mathop{\mathrm{diag}}[(a(1),a(\omega_n),a(\omega_n^2),\ldots, a(\omega_n^{n-1})]F_n.$$

We shall use the package SpecialMatrices.jl.



In [35]:




In [37]:

using SpecialMatrices
using Polynomials




In [38]:

varinfo(SpecialMatrices)




Out[38]:

name
size
summary

Cauchy
40 bytes
UnionAll

Circulant
40 bytes
UnionAll

Companion
40 bytes
UnionAll

Frobenius
40 bytes
UnionAll

Hankel
40 bytes
UnionAll

Hilbert
40 bytes
UnionAll

InverseHilbert
40 bytes
UnionAll

Kahan
80 bytes
UnionAll

Riemann
40 bytes
UnionAll

SpecialMatrices
185.473 KiB
Module

Strang
40 bytes
UnionAll

Toeplitz
40 bytes
UnionAll

Vandermonde
40 bytes
UnionAll

embed
0 bytes
typeof(embed)




In [39]:

import Random
Random.seed!(127)
n=6
a=rand(-9:9,n)




Out[39]:

6-element Array{Int64,1}:
-1
3
2
-2
9
7




In [40]:

C=Circulant(a)




Out[40]:

6×6 Circulant{Int64}:
-1   7   9  -2   2   3
3  -1   7   9  -2   2
2   3  -1   7   9  -2
-2   2   3  -1   7   9
9  -2   2   3  -1   7
7   9  -2   2   3  -1




In [41]:

# Check for normality
C*C'-C'*C




Out[41]:

6×6 Array{Int64,2}:
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0




In [42]:

Cm=Matrix(C)




Out[42]:

6×6 Array{Int64,2}:
-1   7   9  -2   2   3
3  -1   7   9  -2   2
2   3  -1   7   9  -2
-2   2   3  -1   7   9
9  -2   2   3  -1   7
7   9  -2   2   3  -1




In [43]:

Cm*Cm'-Cm'*Cm




Out[43]:

6×6 Array{Int64,2}:
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0
0  0  0  0  0  0




In [48]:

p1=Poly(a)




Out[48]:

-1 + 3∙x + 2∙x2 - 2∙x3 + 9∙x4 + 7∙x5




In [49]:

# Polynomials are easy to manipulate
p2=p1*p1




Out[49]:

1 - 6∙x + 5∙x2 + 16∙x3 - 26∙x4 + 32∙x5 + 82∙x6 - 8∙x7 + 53∙x8 + 126∙x9 + 49∙x10




In [50]:

# n-th complex root of unity
ω=exp(2*π*im/n)




Out[50]:

0.5000000000000001 + 0.8660254037844386im




In [51]:

v=[ω^k for k=0:n-1]
F=Vandermonde(v)




Out[51]:

6×6 Vandermonde{Complex{Float64}}:
1.0+0.0im   1.0+0.0im           1.0+0.0im          …   1.0+0.0im
1.0+0.0im   0.5+0.866025im     -0.5+0.866025im         0.5-0.866025im
1.0+0.0im  -0.5+0.866025im     -0.5-0.866025im        -0.5-0.866025im
1.0+0.0im  -1.0+3.88578e-16im   1.0-7.77156e-16im     -1.0+1.94289e-15im
1.0+0.0im  -0.5-0.866025im     -0.5+0.866025im        -0.5+0.866025im
1.0+0.0im   0.5-0.866025im     -0.5-0.866025im     …   0.5+0.866025im




In [52]:

Fn=Matrix(F)/sqrt(n)
Λ=Fn*C*Fn'




Out[52]:

6×6 Array{Complex{Float64},2}:
18.0+0.0im          …   1.01955e-14-9.84146e-15im
-2.35916e-15+1.74243e-15im     -4.41314e-15-3.55271e-15im
-2.33153e-15+1.74029e-15im      8.76935e-16+4.87854e-15im
-2.81676e-16+3.4972e-15im      -2.61444e-15+5.86306e-16im
1.71426e-15+5.79502e-15im      3.79963e-15+9.37263e-15im
1.06966e-14+8.37778e-15im  …           0.5+9.52628im




In [54]:

# Check
[diag(Λ) p1.(v) eigvals(Matrix(C))]




Out[54]:

6×3 Array{Complex{Float64},2}:
18.0+0.0im       18.0+0.0im          -13.5-2.59808im
0.5-9.52628im    0.5-9.52628im      -13.5+2.59808im
-13.5+2.59808im  -13.5+2.59808im        0.5-9.52628im
2.0+0.0im        2.0-3.10862e-15im    0.5+9.52628im
-13.5-2.59808im  -13.5-2.59808im        2.0+0.0im
0.5+9.52628im    0.5+9.52628im       18.0+0.0im



## Hermitian and real symmetric matrices

For more details and the proofs of the Facts below, see W. Barrett, Hermitian and Positive Definite Matrices and the references therein.

### Definitions

Matrix $A\in \mathbb{C}^{n\times n}$ is Hermitian or self-adjoint if $A^*=A$, or element-wise, $\bar a_{ij}=a_{ji}$. We say $A\in\mathcal{H}_n$.

Matrix $A\in \mathbb{R}^{n\times n}$ is symmetric if $A^T=A$, or element-wise, $a_{ij}=a_{ji}$. We say $A\in\mathcal{S}_n$.

Rayleigh qoutient of $A\in\mathcal{H}_n$ and nonzero vector $x\in\mathbb{C}^n$ is

$$R_A(x)=\frac{x^*Ax}{x^*x}.$$

Matrices $A,B \in\mathcal{H}_n$ are congruent if there exists nonsingular matrix $C$ such that $B=C^*AC$.

Inertia of $A\in\mathcal{H}_n$ is the ordered triple $$\mathop{\mathrm{in}}(A)=(\pi(A),\nu(A),\zeta(A)),$$

where $\pi(A)$ is the number of positive eigenvalues of $A$, $\nu(A)$ is the number of negative eigenvalues of $A$, and $\zeta(A)$ is the number of zero eigenvalues of $A$.

Gram matrix of a set of vectors $x_1,x_2,\ldots,x_k\in\mathbb{C}^{n}$ is the matrix $G$ with entries $G_{ij}=x_i^*x_j$.

### Facts

Assume $A$ is Hermitian and $x\in\mathbb{C}^n$ is nonzero. Then

1. Real symmetric matrix is Hermitian, and real Hermitian matrix is symmetric.

2. Hermitian and real symmetric matrices are normal.

3. $A+A^*$, $A^*A$, and $AA^*$ are Hermitian.

4. If $A$ is nonsingular, $A^{-1}$ is Hermitian.

5. Main diagonal entries of $A$ are real.

6. Matrix $T$ from the Schur decomposition of $A$ is Hermitian. Consequently:

• $T$ is diagonal and real, and has eigenvalues of $A$ on diagonal,
• matrix $Q$ of the Schur decomposition is the unitary matrix of eigenvectors,
• all eigenvalues of $A$ are semisimple and $A$ is nondefective,
• eigenvectors corresponding to distinct eigenvalues are orthogonal.
7. Spectral Theorem. To summarize:

• if $A\in\mathcal{H}_n$, there is a unitary matrix $U$ and real diagonal matrix $\Lambda$ such that $A=U\Lambda U^*$. The diagonal entries of $\Lambda$ are the eigenvalues of $A$, and the columns of $U$ are the corresponding eigenvectors.
• if $A\in\mathcal{S}_n$, the same holds with orthogonal matrix $U$, $A=U\Lambda U^T$.
• if $A\in\mathcal{H}_n$ with eigenpairs $(\lambda_i,u_i)$, then

$$A=\lambda_1u_1u_1^*+\lambda_2 u_2u_2^* +\cdots +\lambda_n u_nu_n^*.$$

• similarly, if $A\in\mathcal{S}_n$, then

$$A=\lambda_1u_1u_1^T+\lambda_2 u_2u_2^T +\cdots +\lambda_n u_nu_n^T.$$

1. Since all eigenvalues of $A$ are real, they can be ordered:
$$\lambda_1\geq \lambda_2\geq \cdots \geq \lambda_n.$$
1. Rayleigh-Ritz Theorem. It holds:
\begin{align*} \lambda_n &\leq \frac{x^*Ax}{x^*x} \leq \lambda_1, \\ \lambda_1&=\max_x\frac{x^*Ax}{x^*x} =\max_{\|x\|_2=1} x^*Ax, \\ \lambda_n&=\min_x\frac{x^*Ax}{x^*x} =\min_{\|x\|_2=1} x^*Ax. \end{align*}
1. Courant-Fischer Theorem. It holds:
\begin{align*} \lambda_k& =\max_{\dim(W)=k}\min_{x\in W} \frac{x^*Ax}{x^*x}\\ & =\min_{\dim(W)=n-k+1}\max_{x\in W} \frac{x^*Ax}{x^*x}. \end{align*}

1. Cauchy Interlace Theorem. Let $B$ be the principal submatrix of $A$ obtained by deleting the $i$-th row and the $i$-th column of $A$. Let $\mu_1\geq \mu_2\geq \cdots \geq \mu_{n-1}$ be the eignvalues of $B$. Then
$$\lambda_1\geq \mu_1\geq \lambda_2\geq \mu_2\geq \lambda_3\geq\cdots\geq \lambda_{n-1}\geq\mu_{n-1}\geq\lambda_n.$$
1. Weyl Inequalities. For $A,B\in\mathcal{H}_n$, it holds:
\begin{align*} \lambda_{j+k-1}(A+B)& \leq \lambda_j(A)+\lambda_k(B), & \textrm{for} \ j+k\leq n+1,\\ \lambda_{j+k-n}(A+B)& \geq \lambda_j(A)+\lambda_k(B), & \textrm{for} \ j+k\geq n+1, \end{align*}

and, in particular,

$$\lambda_j(A)+\lambda_n(B) \leq \lambda_j(A+B) \leq \lambda_j(A)+\lambda_1(B),\quad \textrm{for} \ j=1,2,\ldots,n.$$
1. $\pi(A)+\nu(A)+\zeta(A)=n$.

2. $\mathop{\mathrm{rank}}(A)=\pi(A)+\nu(A)$.

3. If $A$ is nonsingular, $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(A^{-1})$.

4. If $A,B \in\mathcal{H}_n$ are similar, $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(B)$.

5. Sylvester's Law of Inertia. $A,B \in\mathcal{H}_n$ are congruent if and only if $\mathop{\mathrm{in}}(A)=\mathop{\mathrm{in}}(B)$.

6. Subadditivity of Inertia. For $A,B \in\mathcal{H}_n$,

$$\pi(A+B)\leq \pi(A)+\pi(B), \qquad \nu(A+B)\leq \nu(A)+\nu(B).$$
1. Gram matrix is Hermitian.

### Example - Hermitian matrix



In [55]:

# Generating Hermitian matrix
Random.seed!(432)
n=4
A=rand(ComplexF64,n,n)




Out[55]:

4×4 Array{Complex{Float64},2}:
1.84102+0.0im        1.06576+0.33703im   …  0.323092+0.0339997im
1.06576-0.33703im    1.54757+0.0im           1.60398+0.359329im
1.4439-0.305891im   1.37017+0.022169im     0.570542+0.0799815im
0.323092-0.0339997im  1.60398-0.359329im      1.71878+0.0im




In [56]:

ishermitian(A)




Out[56]:

true




In [57]:

# Diagonal entries
diag(A)




Out[57]:

4-element Array{Complex{Float64},1}:
1.8410195749746157 + 0.0im
1.5475690889240252 + 0.0im
0.9901494279719274 + 0.0im
1.718783744539766 + 0.0im




In [58]:

# Schur decomposition
F=schur(A)




Out[58]:

Schur{Complex{Float64},Array{Complex{Float64},2}}
T factor:
4×4 Array{Complex{Float64},2}:
4.82827-3.02154e-17im  …  4.71551e-16+2.51989e-16im
0.0+0.0im             -3.6118e-17+4.659e-17im
0.0+0.0im             -2.3491e-17+8.73366e-17im
0.0+0.0im               -0.436509+3.71066e-17im
Z factor:
4×4 Array{Complex{Float64},2}:
0.49437+1.88747e-17im  -0.626715+1.54319e-16im  …  -0.167364+9.09494e-12im
0.568825-0.134676im      0.236761-0.0057747im       -0.596825-0.289187im
0.454329-0.101044im     -0.270866+0.0654529im        0.602523+0.198302im
0.409987-0.170835im      0.605081-0.327664im         0.350606+0.0829377im
eigenvalues:
4-element Array{Complex{Float64},1}:
4.8282732952424965 - 3.021541012371012e-17im
1.761576281473247 - 1.3354010200483882e-16im
-0.05581915748485714 + 2.9070704761235285e-17im
-0.4365085828205522 + 3.710661184361823e-17im




In [59]:

F.T




Out[59]:

4×4 Array{Complex{Float64},2}:
4.82827-3.02154e-17im  …  4.71551e-16+2.51989e-16im
0.0+0.0im             -3.6118e-17+4.659e-17im
0.0+0.0im             -2.3491e-17+8.73366e-17im
0.0+0.0im               -0.436509+3.71066e-17im




In [60]:

λ,U=eigen(A)




Out[60]:

Eigen{Complex{Float64},Float64,Array{Complex{Float64},2},Array{Float64,1}}
values:
4-element Array{Float64,1}:
-0.43650858282055194
-0.05581915748485722
1.7615762814732465
4.828273295242489
vectors:
4×4 Array{Complex{Float64},2}:
-0.162869+0.0385275im  -0.526681-0.239632im   …  -0.456339-0.190149im
-0.647367-0.14403im     0.376493+0.143651im      -0.576866-0.0944707im
0.63199+0.0542742im   0.550588+0.0147649im     -0.458243-0.0814773im
0.360282-0.0im        -0.446582-0.0im           -0.444155-0.0im




In [61]:

λ




Out[61]:

4-element Array{Float64,1}:
-0.43650858282055194
-0.05581915748485722
1.7615762814732465
4.828273295242489




In [62]:

# Spectral theorem
norm(A-U*Diagonal(λ)*U')




Out[62]:

6.702874538353148e-15




In [63]:

# Spectral theorem
A-sum([λ[i]*U[:,i]*U[:,i]' for i=1:n])




Out[63]:

4×4 Array{Complex{Float64},2}:
2.22045e-15+0.0im          …  2.22045e-16+3.33067e-16im
2.44249e-15-7.77156e-16im     2.22045e-16-1.11022e-16im
1.9984e-15-6.66134e-16im     6.66134e-16+2.77556e-17im
3.33067e-16-3.33067e-16im     6.66134e-16+0.0im




In [64]:

# Cauchy Interlace Theorem (repeat several times)
@show i=rand(1:n)
μ=eigvals(A[[1:i-1;i+1:n],[1:i-1;i+1:n]])




i = rand(1:n) = 4

Out[64]:

3-element Array{Float64,1}:
-0.31223338776364146
0.5740391221926207
4.116932357441587




In [65]:

λ




Out[65]:

4-element Array{Float64,1}:
-0.43650858282055194
-0.05581915748485722
1.7615762814732465
4.828273295242489




In [66]:

# Inertia
inertia(A)=[sum(eigvals(A).>0), sum(eigvals(A).<0), sum(eigvals(A).==0)]




Out[66]:

inertia (generic function with 1 method)




In [67]:

inertia(A)




Out[67]:

3-element Array{Int64,1}:
2
2
0




In [68]:

# Similar matrices
C=rand(ComplexF64,n,n)
B=C*A*inv(C)
eigvals(B)




Out[68]:

4-element Array{Complex{Float64},1}:
-0.43650858282055194 - 4.760622708291153e-16im
-0.055819157484860334 - 1.1732254580750932e-15im
1.7615762814732525 - 1.5950102848760153e-16im
4.828273295242495 + 2.113110034174248e-15im




In [69]:

B




Out[69]:

4×4 Array{Complex{Float64},2}:
1.85492-0.696005im  -0.566638-1.70282im  …   0.657662+1.85676im
1.50219-1.52789im     -1.6778-3.07909im      0.351122+4.20872im
1.58685-1.61565im    -1.02112-3.37356im     -0.736184+4.17921im
0.917302-0.848219im   -1.91749-1.04715im         2.067+2.07822im



This did not work numerically due to rounding errors!



In [70]:

# Congruent matrices - this does not work either, without some preparation
B=C'*A*C
inertia(A)==inertia(B)




MethodError: no method matching isless(::Int64, ::Complex{Float64})
Closest candidates are:
isless(!Matched::Missing, ::Any) at missing.jl:87
isless(!Matched::PyCall.PyObject, ::Any) at C:\Users\Ivan_Slapnicar\.julia\packages\PyCall\zqDXB\src\pyoperators.jl:75
isless(!Matched::Sym, ::Number) at C:\Users\Ivan_Slapnicar\.julia\packages\SymPy\JaxDZ\src\generic.jl:14
...

Stacktrace:
[1] <(::Int64, ::Complex{Float64}) at .\operators.jl:268
[2] >(::Complex{Float64}, ::Int64) at .\operators.jl:294
[8] inertia(::Array{Complex{Float64},2}) at .\In[66]:2
[9] top-level scope at In[70]:3




In [71]:

# However,
inertia(A)==inertia(Hermitian(B))




Out[71]:

true




In [72]:

# or,
inertia((B+B')/2)




Out[72]:

3-element Array{Int64,1}:
2
2
0




In [73]:

# Weyl's Inequalities
B=rand(ComplexF64,n,n)
B=(B+B')/10
@show λ
μ=eigvals(B)
γ=eigvals(A+B)
μ,γ




λ = [-0.43650858282055194, -0.05581915748485722, 1.7615762814732465, 4.828273295242489]

Out[73]:

([-0.0988735390914506, -0.01182765702077481, 0.07296810055246222, 0.3287937288227195], [-0.4545242135758396, -0.07251608580244176, 1.7740633000178827, 5.141559469033688])




In [74]:

# Theorem uses different order!
j=rand(1:n)
k=rand(1:n)
@show j,k
if j+k<=n+1
@show sort(γ,rev=true)[j+k-1], sort(λ,rev=true)[j]+sort(μ,rev=true)[k]
end
if j+k>=n+1
sort(γ,rev=true)[j+k-n], sort(λ,rev=true)[j]+sort(μ,rev=true)[k]
end




(j, k) = (1, 3)
((sort(γ, rev = true))[(j + k) - 1], (sort(λ, rev = true))[j] + (sort(μ, rev = true))[k]) = (-0.07251608580244176, 4.8164456382217145)




In [75]:

sort(λ,rev=true)




Out[75]:

4-element Array{Float64,1}:
4.828273295242489
1.7615762814732465
-0.05581915748485722
-0.43650858282055194



### Example - Real symmetric matrix



In [76]:

# Generating real symmetric matrix
Random.seed!(531)
n=6
A=rand(-9:9,n,n)
A=A+A'




Out[76]:

6×6 Array{Int64,2}:
8   17   -6  -11  -2    9
17  -14   12    3  10    1
-6   12    0  -10   3   -1
-11    3  -10   12   9  -10
-2   10    3    9  10    7
9    1   -1  -10   7   10




In [77]:

issymmetric(A)




Out[77]:

true




In [78]:

T,Q=schur(A)




Out[78]:

Schur{Float64,Array{Float64,2}}
T factor:
6×6 Array{Float64,2}:
-33.0112    1.89247e-16   1.62178e-15  …  -5.61774e-15  -1.8153e-15
0.0     -11.1569        2.6645e-15      -1.71972e-15   4.61092e-16
0.0       0.0          32.326            3.6099e-15    4.48108e-16
0.0       0.0           0.0             -2.85135e-15  -1.00712e-15
0.0       0.0           0.0             10.2386       -1.1609e-15
0.0       0.0           0.0          …   0.0          21.1844
Z factor:
6×6 Array{Float64,2}:
-0.44514     0.0429978  -0.571264    0.506595   0.451608  -0.114327
0.755485   -0.226136   -0.212836    0.45081   -0.096524  -0.346784
-0.408122   -0.397543   -0.141593    0.19712   -0.78213   -0.0690242
-0.21589    -0.531446    0.60184     0.133565   0.344137  -0.415303
-0.128703    0.544883    0.0240726  -0.149403  -0.181727  -0.794111
0.0368434  -0.457866   -0.495499   -0.679054   0.153464  -0.242522
eigenvalues:
6-element Array{Float64,1}:
-33.01123470268622
-11.156924850235118
32.32596330116704
6.41917104929743
10.238647969944779
21.184377232512148




In [79]:

T




Out[79]:

6×6 Array{Float64,2}:
-33.0112    1.89247e-16   1.62178e-15  …  -5.61774e-15  -1.8153e-15
0.0     -11.1569        2.6645e-15      -1.71972e-15   4.61092e-16
0.0       0.0          32.326            3.6099e-15    4.48108e-16
0.0       0.0           0.0             -2.85135e-15  -1.00712e-15
0.0       0.0           0.0             10.2386       -1.1609e-15
0.0       0.0           0.0          …   0.0          21.1844




In [80]:

Q




Out[80]:

6×6 Array{Float64,2}:
-0.44514     0.0429978  -0.571264    0.506595   0.451608  -0.114327
0.755485   -0.226136   -0.212836    0.45081   -0.096524  -0.346784
-0.408122   -0.397543   -0.141593    0.19712   -0.78213   -0.0690242
-0.21589    -0.531446    0.60184     0.133565   0.344137  -0.415303
-0.128703    0.544883    0.0240726  -0.149403  -0.181727  -0.794111
0.0368434  -0.457866   -0.495499   -0.679054   0.153464  -0.242522




In [81]:

λ,U=eigen(A)
λ




Out[81]:

6-element Array{Float64,1}:
-33.01123470268616
-11.15692485023515
6.4191710492974146
10.238647969944783
21.184377232512134
32.325963301167015




In [82]:

cond(A)




Out[82]:

5.142600882445594




In [83]:

U




Out[83]:

6×6 Array{Float64,2}:
-0.44514    -0.0429978   0.506595   0.451608  -0.114327    0.571264
0.755485    0.226136    0.45081   -0.096524  -0.346784    0.212836
-0.408122    0.397543    0.19712   -0.78213   -0.0690242   0.141593
-0.21589     0.531446    0.133565   0.344137  -0.415303   -0.60184
-0.128703   -0.544883   -0.149403  -0.181727  -0.794111   -0.0240726
0.0368434   0.457866   -0.679054   0.153464  -0.242522    0.495499




In [84]:

A-sum([λ[i]*U[:,i]*U[:,i]' for i=1:n])




Out[84]:

6×6 Array{Float64,2}:
-3.19744e-14   2.4869e-14   -3.01981e-14  …  -9.32587e-15  -1.24345e-14
2.84217e-14  -4.79616e-14   3.90799e-14      5.32907e-15   4.44089e-16
-3.01981e-14   3.90799e-14  -3.18634e-14     -1.02141e-14   3.10862e-15
-2.66454e-14   2.4869e-14   -2.4869e-14       1.77636e-15   1.59872e-14
-9.99201e-15   5.32907e-15  -1.06581e-14      1.95399e-14  -1.15463e-14
-1.24345e-14   0.0           3.10862e-15  …  -1.24345e-14   6.57252e-14




In [85]:

inertia(A)




Out[85]:

3-element Array{Int64,1}:
4
2
0




In [86]:

C=rand(n,n)
inertia(C'*A*C)




Out[86]:

3-element Array{Int64,1}:
4
2
0



## Positive definite matrices

These matrices are an important subset of Hermitian or real symmteric matrices.

### Definitions

Matrix $A\in\mathcal{H}_n$ is positive definite (PD) if $x^*Ax>0$ for all nonzero $x\in\mathbb{C}^n$.

Matrix $A\in\mathcal{H}_n$ is positive semidefinite (PSD) if $x^*Ax\geq 0$ for all nonzero $x\in\mathbb{C}^n$.

### Facts

1. $A\in\mathcal{S}_n$ is PD if $x^TAx>0$ for all nonzero $x\in \mathbb{R}^n$, and is PSD if $x^TAx\geq 0$ for all $x\in \mathbb{R}^n$.

2. If $A,B\in \mathrm{PSD}_n$, then $A+B\in \mathrm{PSD}_n$. If, in addition, $A\in \mathrm{PD}_n$, then $A+B\in \mathrm{PD}_n$.

3. If $A\in \mathrm{PD}_n$, then $\mathop{\mathrm{tr}}(A)>0$ and $\det(A)>0$.

4. If $A\in \mathrm{PSD}_n$, then $\mathop{\mathrm{tr}}(A)\geq 0$ and $\det(A)\geq 0$.

5. Any principal submatrix of a PD matrix is PD. Any principal submatrix of a PSD matrix is PSD. Consequently, all minors are positive or nonnegative, respectively.

6. $A\in\mathcal{H}_n$ is PD iff every leading principal minor of $A$ is positive. $A\in\mathcal{H}_n$ is PSD iff every principal minor is nonnegative.

7. For $A\in \mathrm{PSD}_n$, there exists unique PSD $k$-th root, $A^{1/k}=U\Lambda^{1/k} U^*$.

8. Cholesky Factorization. $A\in\mathcal{H}_n$ if PD iff there is an invertible lower triangular matrix $L$ with positive diagonal entries such that $A=LL^*$.

9. Gram matrix is PSD. If the vectors are linearly independent, Gram matrix is PD.

### Example - Positive definite matrix



In [87]:

# Generating positive definite matrix as a Gram matrix
n=5
A=rand(n,n)+im*rand(n,n)
A=A*A'




Out[87]:

5×5 Array{Complex{Float64},2}:
4.78771+0.0im       2.49491-1.08523im   …  2.78758-0.586195im
2.49491+1.08523im   1.99217+0.0im           1.3238+0.166899im
2.70349+0.882787im   1.4953-0.238075im     1.89543+0.154661im
3.23928+1.3186im    1.78581-0.505257im     2.20113+0.809557im
2.78758+0.586195im   1.3238-0.166899im     2.69021+0.0im




In [88]:

ishermitian(A)




Out[88]:

true




In [89]:

eigvals(A)




Out[89]:

5-element Array{Float64,1}:
0.017539281952608544
0.3247567476856929
0.5585201067285008
1.607345246729099
12.73056538580112




In [90]:

# Positivity of principal leading minors
[det(A[1:i,1:i]) for i=1:n]




Out[90]:

5-element Array{Complex{Float64},1}:
4.787713836028915 + 0.0im
2.135651400648032 + 0.0im
0.5577068310530778 + 6.668576318172403e-17im
0.22986275838388015 - 1.249000902703301e-16im
0.06509770344359479 - 2.2854981795994433e-16im




In [91]:

# Square root
λ,U=eigen(A)
Ar=U*Diagonal(λ)*U'
A-Ar




Out[91]:

5×5 Array{Complex{Float64},2}:
0.0+0.0im          …  -4.44089e-16+2.22045e-16im
4.44089e-16+2.22045e-16im     -2.22045e-16+1.94289e-16im
0.0-1.11022e-16im     -6.66134e-16+0.0im
-8.88178e-16-1.33227e-15im     -8.88178e-16-2.22045e-16im
-4.44089e-16-2.22045e-16im      1.77636e-15+0.0im




In [92]:

# Cholesky factorization - the upper triangular factor is returned
L=cholesky(A).U




Out[92]:

5×5 UpperTriangular{Complex{Float64},Array{Complex{Float64},2}}:
2.18808+0.0im   1.14022-0.495971im  …   1.27398-0.267903im
⋅      0.667884+0.0im          -0.39183-0.238799im
⋅               ⋅              0.346592-0.232844im
⋅               ⋅              0.570903-0.0370585im
⋅               ⋅              0.532168+0.0im




In [93]:

norm(A-L'*L)




Out[93]:

7.889616961715915e-16



### Example - Positive semidefinite matrix



In [94]:

# Generating positive semidefinite matrix as a Gram matrix, try it several times
n=6
m=4
A=rand(-9:9,n,m)




Out[94]:

6×4 Array{Int64,2}:
7   1   1   1
5   5   1  -7
6   7  -6  -6
9  -3  -1   4
0  -8   5   8
5  -2  -3   1




In [95]:

A=A*A'




Out[95]:

6×6 Array{Int64,2}:
52   34    37   63     5  31
34  100   101    1   -91   5
37  101   157   15  -134  28
63    1    15  107    51  58
5  -91  -134   51   153   9
31    5    28   58     9  39




In [96]:

# There are rounding errors!
eigvals(A)




Out[96]:

6-element Array{Float64,1}:
-3.7549491324113813e-14
4.0880160397987686e-14
4.449546314098275
36.197204039470954
201.65495386016767
365.6982957862633




In [97]:

s=svdvals(A)




Out[97]:

6-element Array{Float64,1}:
365.6982957862633
201.65495386016755
36.19720403947091
4.449546314098264
1.0124582481354459e-14
1.3670083042179748e-15




In [98]:

s[1]*eps()




Out[98]:

8.120133360961808e-14




In [99]:

rank(A)




Out[99]:

4




In [100]:

# Cholesky factorization - this can fail
L=cholesky(A).U




PosDefException: matrix is not positive definite; Cholesky factorization failed.

Stacktrace:
[1] checkpositivedefinite at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\factorization.jl:18 [inlined]
[2] cholesky!(::Hermitian{Float64,Array{Float64,2}}, ::Val{false}; check::Bool) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\cholesky.jl:226
[3] cholesky!(::Array{Float64,2}, ::Val{false}; check::Bool) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\cholesky.jl:258
[4] #cholesky#136 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\cholesky.jl:348 [inlined]
[5] cholesky at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\LinearAlgebra\src\cholesky.jl:348 [inlined] (repeats 2 times)
[6] top-level scope at In[100]:1



### Example - Covariance and correlation matrices

Correlation matrix is diagonally scaled covariance matrix.



In [101]:

Random.seed!(651)
x=rand(10,5)




Out[101]:

10×5 Array{Float64,2}:
0.828124   0.525628  0.664775   0.734664   0.964204
0.760442   0.778445  0.0804647  0.0480193  0.731189
0.718784   0.331317  0.413675   0.96678    0.0723234
0.226571   0.363896  0.513771   0.642912   0.414969
0.655436   0.127381  0.529349   0.9076     0.548456
0.926612   0.973135  0.76572    0.941035   0.902337
0.948178   0.796415  0.232122   0.46858    0.647936
0.182908   0.977579  0.179693   0.286669   0.872155
0.0559694  0.857692  0.514281   0.227538   0.773259
0.317514   0.58723   0.832853   0.249614   0.603709




In [102]:

using Statistics
A=cov(x)




Out[102]:

5×5 Array{Float64,2}:
0.11085     -0.00555573  -0.00317136   0.0516332   0.00382669
-0.00555573   0.0850514   -0.0163683   -0.0504624   0.0509808
-0.00317136  -0.0163683    0.0622391    0.03573     0.00499945
0.0516332   -0.0504624    0.03573      0.113258   -0.029674
0.00382669   0.0509808    0.00499945  -0.029674    0.0705389




In [103]:

# Covariance matrix is a Gram matrix
(x.-mean(x,dims=1))'*(x.-mean(x,dims=1))/(size(x,1)-1)-A




Out[103]:

5×5 Array{Float64,2}:
0.0          -8.67362e-19  0.0           6.93889e-18  0.0
-8.67362e-19   1.38778e-17  0.0          -6.93889e-18  0.0
0.0           0.0          0.0           0.0          0.0
6.93889e-18  -6.93889e-18  0.0           0.0          0.0
3.03577e-18  -6.93889e-18  8.67362e-19   0.0          0.0




In [104]:

B=cor(x)




Out[104]:

5×5 Array{Float64,2}:
1.0        -0.0572179  -0.0381809   0.460814   0.0432754
-0.0572179   1.0        -0.224973   -0.514153   0.65819
-0.0381809  -0.224973    1.0         0.425565   0.075453
0.460814   -0.514153    0.425565    1.0       -0.331991
0.0432754   0.65819     0.075453   -0.331991   1.0




In [105]:

# Diagonal scaling
D=1 ./sqrt.(diag(A))




Out[105]:

5-element Array{Float64,1}:
3.0035294964691293
3.4289352335301344
4.008376039395491
2.971425838234638
3.765179232014274




In [106]:

Diagonal(D)*A*Diagonal(D)




Out[106]:

5×5 Array{Float64,2}:
1.0        -0.0572179  -0.0381809   0.460814   0.0432754
-0.0572179   1.0        -0.224973   -0.514153   0.65819
-0.0381809  -0.224973    1.0         0.425565   0.075453
0.460814   -0.514153    0.425565    1.0       -0.331991
0.0432754   0.65819     0.075453   -0.331991   1.0




In [107]:

eigvals(A)




Out[107]:

5-element Array{Float64,1}:
0.018962686335371806
0.026747954810784666
0.0754000977625162
0.1179434475533185
0.20288364944175208




In [108]:

eigvals(B)




Out[108]:

5-element Array{Float64,1}:
0.23570502689844522
0.2900555550445134
1.0485883146491362
1.253733349516822
2.1719177538910843




In [109]:

C=cov(x')




Out[109]:

10×10 Array{Float64,2}:
0.0274039    0.0123544   -0.00911327  …  -0.0134788    -0.00888699
0.0123544    0.144262    -0.0658908       0.0407097    -0.00597273
-0.00911327  -0.0658908    0.12114        -0.10068      -0.0611635
-0.00365153  -0.0514985    0.0163041       0.00148382    0.00103769
0.0247812   -0.0607921    0.0659274      -0.0739685    -0.0382523
-0.00073939   0.0164507    0.00580677  …   0.000502438  -0.0133413
0.00684369   0.0930205   -0.00177452     -0.00984111   -0.0331679
-0.00541184   0.0918037   -0.0952559       0.116714      0.0212344
-0.0134788    0.0407097   -0.10068         0.1183        0.05371
-0.00888699  -0.00597273  -0.0611635       0.05371       0.0558742




In [110]:

eigvals(C)




Out[110]:

10-element Array{Float64,1}:
-2.8678282441648895e-17
-2.0606155317260784e-17
-1.142644713428881e-17
3.2024131844992747e-19
8.633927427481286e-18
4.0501282659044764e-17
0.05465143722525408
0.06765627053231633
0.2131797129097473
0.4740510413291936




In [111]:

inertia(C)




Out[111]:

3-element Array{Int64,1}:
7
3
0




In [112]:

rank(C)




Out[112]:

4



Explain the function rank().



In [113]:

@which rank(C)




Out[113]:

rank(A::AbstractArray{T,2} where T; atol, rtol) in LinearAlgebra at C:\Users\Ivan_Slapnicar\AppData\Local\Programs\Julia\Julia-1.4.1\share\julia\stdlib\v1.4\LinearAlgebra\src\generic.jl:984