Skip to main content

Scientific Computing

Essential Maths

Linear algebra 1 []

Scientific Computing

Essential Maths

Systems of di... []

This material has been adapted from material by Fergus Cooper from the "Essential Mathematics" module of the SABS R³ Center for Doctoral Training.

This material has been adapted from material by Fergus Cooper from the "Essential Mathematics" module of the SABS R³ Center for Doctoral Training.

Creative Commons License
This course material was developed as part of UNIVERSE-HPC, which is funded through the SPF ExCALIBUR programme under grant number EP/W035731/1

This course material was developed as part of UNIVERSE-HPC, which is funded through the SPF ExCALIBUR programme under grant number EP/W035731/1

Creative Commons License

Linear algebra 2

Eigenvalues and Eigenvectors


YouTube lecture recording from October 2020

The following YouTube video was recorded for the 2020 iteration of the course. The material is still very similar:

Youtube lecture thumbnail


Using Python to calculate the inverse

To find A1A^{-1} for

A=(302303011)A = \begin{pmatrix}3&0&2\\ 3&0&-3\\ 0&1&1\end{pmatrix}

NumPy

import numpy as np A = np.array([[3, 0, 2], [3, 0, -3], [0, 1, 1]]) np.linalg.inv(A)
array([[ 0.2 , 0.13333333, 0. ], [-0.2 , 0.2 , 1. ], [ 0.2 , -0.2 , -0. ]])

SymPy

import sympy as sp A = sp.Matrix([[3, 0, 2], [3, 0, -3], [0, 1, 1]]) A.inv()

[1521501515115150]\displaystyle \left[\begin{matrix}\frac{1}{5} & \frac{2}{15} & 0\\- \frac{1}{5} & \frac{1}{5} & 1\\\frac{1}{5} & - \frac{1}{5} & 0\end{matrix}\right]

It doesn't always work

Consider the following system

eq1:x+y+z=aeq2:2x+5y+2z=beq3:7x+10y+7z=c\begin{array}{cccccccc} eq1: & x & + & y & + & z & = & a \\ eq2: & 2x & + & 5y & + & 2z & = & b \\ eq3: & 7x & +& 10y & + & 7z & = & c \end{array}
A = np.array([[1, 1, 1], [2, 5, 2], [7, 10, 7]]) np.linalg.inv(A)
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) Cell In[5], line 2 1 A = np.array([[1, 1, 1], [2, 5, 2], [7, 10, 7]]) ----> 2 np.linalg.inv(A) File ~/GitRepos/gutenberg-material/DtcMathsStats2018/venv2/lib/python3.10/site-packages/numpy/linalg/linalg.py:561, in inv(a) 559 signature = 'D->D' if isComplexType(t) else 'd->d' 560 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) --> 561 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) 562 return wrap(ainv.astype(result_t, copy=False)) File ~/GitRepos/gutenberg-material/DtcMathsStats2018/venv2/lib/python3.10/site-packages/numpy/linalg/linalg.py:112, in _raise_linalgerror_singular(err, flag) 111 def _raise_linalgerror_singular(err, flag): --> 112 raise LinAlgError("Singular matrix") LinAlgError: Singular matrix

Singular matrices

The rank of an   n×n  \;n\,\times\,n\; matrix   A  \;A\; is the number of linearly independent rows in   A  \;A\; (rows not combinations of other rows).

When   rank(A)<n  \;\text{rank}(A) < n\; then

  • The system   Ax=b  \;A\textbf{x} = \textbf{b}\; has fewer equations than unknowns

  • The matrix is said to be singular

  • A  A\; has no inverse

  • The determinant of   A  \;A\; is 0

  • The equation   Au=0  \;A\textbf{u} = \textbf{0}\; has non-trivial solutions (solutions where u0\textbf{u} \neq \textbf{0})

Singular matrix example

An under-determined system (fewer equations than unknowns) may mean that there are many solutions or that there are no solutions.

An example with many solutions is

x+y=1,2x+2y=2,\begin{align*} x+y &=& 1,\\ 2x+2y &=& 2, \end{align*}

has infinitely many solutions (x=0, y=1; x=-89.3, y=90.3...)

An example with no solutions is

x+y=1,2x+2y=0,\begin{align*} x+y &=& 1,\\ 2x+2y &=& 0, \end{align*}

where the second equation is inconsistent with the first.

These examples use the matrix

(1122),\begin{pmatrix}1&1\\ 2&2\end{pmatrix},

which has rank 1.

Null space

When a matrix is singular we can find non-trivial solutions to   Au=0\;A\textbf{u} = \textbf{0}.

These are vectors which form the null space for   A\;A.

Any vector in the null space makes no difference to the effect that AA is having:

A(x+u)=Ax+Au=Ax+0=Ax.\displaystyle A(\textbf{x} + \textbf{u}) = A\textbf{x} + A\textbf{u} = A\textbf{x} + \textbf{0} = A\textbf{x}.

Note that any combination or scaling of vectors in the null space is also in the null space.

That is, if   Au=0  \;A\textbf{u} = \textbf{0}\; and   Av=0  \;A\textbf{v} = \textbf{0}\; then

A(αu+βv)=0\displaystyle A(\alpha\textbf{u} + \beta\textbf{v}) = \textbf{0}

The number of linearly independent vectors in the null space is denoted  null(A) ~\text{null}(A)~ and

null(A)+rank(A)=order(A).\text{null}(A) + \text{rank}(A) = \text{order}(A).

Null space example

Previous example of a singular system:

A=(1112527107)\displaystyle A = \left(\begin{matrix} 1& 1& 1\\ 2& 5& 2\\ 7& 10&7 \end{matrix}\right)

A = np.array([[1, 1, 1], [2, 5, 2], [7, 10, 7]]) np.linalg.matrix_rank(A)
2
import scipy.linalg scipy.linalg.null_space(A)
array([[-7.07106781e-01], [-1.11022302e-16], [ 7.07106781e-01]])

remember, scaled vectors in the null space are also in the null space, for example,   x=1,y=0,z=1  \;x=1, y=0, z=-1\; is in the null space.

Try it:

(1112527117)(100001000)=?\left(\begin{matrix} 1& 1& 1\\ 2& 5& 2\\ 7& 11&7 \end{matrix}\right) \left(\begin{matrix} -1000\\ 0 \\ 1000 \end{matrix}\right) = \quad ?

np.matmul(A,np.array([-1000,0,1000]))
array([0, 0, 0])

Eigenvectors: motivation

The eigenvalues and eigenvectors give an indication of how much effect the matrix has, and in what direction.

A=(cos(45)sin(45)sin(45)cos(45))has no scaling effect.\displaystyle A=\left(\begin{matrix} \cos(45)&-\sin(45)\\ \sin(45)&\cos(45)\\\end{matrix}\right)\qquad\text{has no scaling effect.} > B=(20012)doubles in x, but halves in y.\displaystyle B=\left(\begin{matrix} 2& 0 \\ 0&\frac{1}{2}\\\end{matrix}\right)\qquad\qquad\text{doubles in }x\text{, but halves in }y\text{.}

Repeated applications of   A  \;A\; stay the same distance from the origin, but repeated applications of   B  \;B\; move towards   (,0).\;(\infty, 0).

  • Transitions with probability

  • Markov chains

  • Designing bridges

  • Solution of systems of linear ODEs

  • Stability of systems of nonlinear ODEs

Eigenvectors and Eigenvalues

A  A\; is a matrix,   v  \;\textbf{v}\; is a non-zero vector,   λ  \;\lambda\; is a scalar,

If:

Av=λv\displaystyle A \textbf{v} = \lambda \textbf{v}

then   v  \;\textbf{v}\; is called an eigenvector and   λ  \;\lambda\; is the corresponding eigenvalue.

Note that if   v  \;\textbf{v}\; is a solution, then so is a scaling   av\;a\textbf{v}:

A(av)=λ(av).\displaystyle A (a \textbf{v}) = \lambda (a \textbf{v}).

Finding Eigenvalues

Another way to write previous equation:

Av=λv,AvλIv=0,(AλI)v=0.\begin{align*} A \textbf{v} &=& \lambda \textbf{v},\\ A \textbf{v} - \lambda I \textbf{v}&=& \textbf{0},\\ (A - \lambda I) \textbf{v}&=& \textbf{0}. \end{align*}

Consider the equation:

Bx=0\displaystyle B\textbf{x}=\textbf{0}

where   B  \;B\; is a matrix and   x  \;\textbf{x}\; is a non-zero column vector.

Assume   B  \;B\; is nonsingular:

B1(Bx)=B10=0\displaystyle B^{-1}(B\textbf{x})=B^{-1}\textbf{0}=\textbf{0}

But:

B1(Bx)=(B1B)x=Ix=x\displaystyle B^{-1}(B\textbf{x})=(B^{-1}B)\textbf{x}=I\textbf{x}=\textbf{x}

This means that:

B1(Bx)=0=x\displaystyle B^{-1}(B\textbf{x})=\textbf{0}=\textbf{x}

but we stated above that   x0  \;\textbf{x}\neq\textbf{0}\; so our assumption that   B  \;B\; is nonsingular must be false:   B  \;B\; is singular.

We can state that:

(AλI)v=0withv0,\displaystyle (A-\lambda I)\textbf{v} = \textbf{0} \quad \text{with} \quad \textbf{v}\neq \textbf{0},

so (AλI)\displaystyle (A-\lambda I) must be singular:

AλI=0.\displaystyle |A-\lambda I|=0.

Example

What are the eigenvalues for this matrix?

A=(2215)\displaystyle A=\left(\begin{matrix}-2&-2\\ 1&-5\\\end{matrix}\right) > AλI=2λ215λ=(2λ)(5λ)(2)\displaystyle |A-\lambda I|=\left\vert\begin{matrix}-2-\lambda&-2\\ 1&-5-\lambda\end{matrix}\right\vert=(-2-\lambda)(-5-\lambda)-(-2) > =10+5λ+λ2+2λ+2=λ2+7λ+12=(λ+3)(λ+4)=0\displaystyle =10+5\lambda+\lambda^2+2\lambda+2=\lambda^2+7\lambda+12=(\lambda+3)(\lambda+4)=0

So the eigenvalues are λ1=3\lambda_1=-3 and λ2=4\lambda_2=-4.

Eigenvalues using Python

Numpy:

A = np.array([[-2, -2], [1, -5]]) np.linalg.eig(A)[0]
array([-3., -4.])

SymPy:

A2 = sp.Matrix([[-2, -2], [1, -5]]) A2.eigenvals()

{4:1, 3:1}\displaystyle \left\{ -4 : 1, \ -3 : 1\right\}

Finding Eigenvectors

For an eigenvalue, the corresponding vector comes from substitution into   Av=λv\;A \textbf{v} = \lambda \textbf{v}:

Example

What are the eigenvectors for

A=(2215)?\displaystyle A=\left(\begin{matrix}-2&-2\\ 1&-5\\\end{matrix}\right)?

The eigenvalues are   λ1=3  \;\lambda_1=-3\; and   λ2=4.  \;\lambda_2=-4.\; The eigenvectors are   v1  \;\textbf{v}_1\; and   v2  \;\textbf{v}_2\; where:

Av1=λ1v1.\displaystyle A\textbf{v}_1=\lambda_1 \textbf{v}_1. > (2215)(u1v1)=(3u13v1)\displaystyle \left(\begin{matrix}-2&-2\\ 1&-5\\\end{matrix}\right) \left(\begin{matrix}u_1\\ v_1\\\end{matrix}\right) = \left(\begin{matrix}-3u_1\\ -3v_1\\\end{matrix}\right) > u1=2v1. (from the top or bottom equation)\displaystyle u_1 = 2v_1. \text{ (from the top or bottom equation)} (u1v1)=(21),(10.5),(4.42.2),(2αα)\displaystyle \left(\begin{matrix}u_1\\ v_1\\\end{matrix}\right) = \left(\begin{matrix}2 \\ 1\\\end{matrix}\right), \left(\begin{matrix}1 \\ 0.5\\\end{matrix}\right), \left(\begin{matrix}-4.4 \\ -2.2\\\end{matrix}\right), \left(\begin{matrix}2\alpha \\ \alpha\\\end{matrix}\right)\ldots

Eigenvectors in Python

Numpy:

A = np.array([[-2, -2], [1, -5]]) np.linalg.eig(A)[1]
array([[0.89442719, 0.70710678], [0.4472136 , 0.70710678]])

SymPy:

c = sp.symbols('c') A2 = sp.Matrix([[-2, c], [1, -5]]) A2.eigenvects()

[(4c+9272, 1, [[324c+921]]), (4c+9272, 1, [[4c+92+321]])]\displaystyle \left[ \left( - \frac{\sqrt{4 c + 9}}{2} - \frac{7}{2}, \ 1, \ \left[ \left[\begin{matrix}\frac{3}{2} - \frac{\sqrt{4 c + 9}}{2}\\1\end{matrix}\right]\right]\right), \ \left( \frac{\sqrt{4 c + 9}}{2} - \frac{7}{2}, \ 1, \ \left[ \left[\begin{matrix}\frac{\sqrt{4 c + 9}}{2} + \frac{3}{2}\\1\end{matrix}\right]\right]\right)\right]

Diagonalising matrices

Any nonsingular matrix AA can be rewritten as a product of eigenvectors and eigenvalues.

If   A  \;A\; has eigenvalues   λ1  \;\lambda_1\; and   λ2  \;\lambda_2\; with corresponding eigenvectors   (u1v1)  \;\left(\begin{matrix}u_1\\ v_1\\\end{matrix}\right)\; and

(u2v2)\displaystyle \left(\begin{matrix}u_2\\ v_2\\\end{matrix}\right)

then

A=(u1u2v1v2)(λ100λ2)(u1u2v1v2)1.\begin{align*} A = \left(\begin{matrix}u_1 & u_2\\v_1 & v_2\\\end{matrix}\right) \left(\begin{matrix}\lambda_1 & 0\\0 & \lambda_2\\\end{matrix}\right) \left(\begin{matrix}u_1 & u_2\\v_1 & v_2\\\end{matrix}\right)^{-1}. \end{align*}

This is like a scaling surrounded by rotations and separates how much effect the matrix has   (λi)  \;(\lambda_i)\; from the directions   (vi).\;(\textbf{v}_i).

For example

A=(2215)\displaystyle A=\left(\begin{matrix}-2&-2\\ 1&-5\\\end{matrix}\right)

A = np.array([[-2, -2], [1, -5]]) w, v = np.linalg.eig(A) inv_v = np.linalg.inv(v) np.matmul( np.matmul(v, np.diag(w)) , inv_v )
array([[-2., -2.], [ 1., -5.]])

Orthogonal eigenvectors

If   x1  \;\textbf{x}_1\; and   x2  \;\textbf{x}_2\; are perpendicular or orthogonal vectors, then the scalar/dot product is zero.

x1.x2=0\displaystyle \textbf{x}_1.\textbf{x}_2=0

e.g.

x1.x2=(213).(271)=2×2+(1×7)+(3×1)=47+3=0.\displaystyle \textbf{x}_1.\textbf{x}_2=\left(\begin{matrix}2\\ -1\\ 3\end{matrix}\right).\left(\begin{matrix}2\\ 7\\ 1\end{matrix}\right)=2\times 2+(-1\times 7)+(3\times1)=4-7+3=0.

Since this dot product is zero, the vectors x1\textbf{x}_1 and x2\textbf{x}_2 are orthogonal.

Symmetric matricies

Symmetric matricies have orthogonal eigenvectors, e.g.

A=(1920162013416431)\displaystyle A=\left(\begin{matrix}19&20&-16\\ 20&13&4 \\ -16&4&31\\\end{matrix}\right)

A = np.array([[19, 20, -16], [20, 13, 4], [-16, 4, 31]]) w, v = np.linalg.eig(A) print(v) print('\ndot products of eigenvectors:\n') print(np.dot(v[:,0],v[:,1])) print(np.dot(v[:,0],v[:,2])) print(np.dot(v[:,1],v[:,2]))
[[ 0.66666667 -0.66666667 0.33333333] [-0.66666667 -0.33333333 0.66666667] [ 0.33333333 0.66666667 0.66666667]] dot products of eigenvectors: -1.942890293094024e-16 1.3877787807814457e-16 1.1102230246251565e-16

Normalised eigenvectors

If

(xyz),\displaystyle \left(\begin{matrix}x\\ y\\ z\\\end{matrix}\right),

is an eigenvector, then

(xx2+y2+z2yx2+y2+z2zx2+y2+z2)\displaystyle \left(\begin{matrix}\frac{x}{\sqrt{x^2+y^2+z^2}}\\ \frac{y}{\sqrt{x^2+y^2+z^2} }\\ \frac{z}{\sqrt{x^2+y^2+z^2} }\end{matrix}\right)

is the corresponding normalised vector: a vector of unit length (magnitude).

Orthogonal matrices

A matrix is orthogonal if its columns are normalised orthogonal vectors. One can prove that if   M  \;M\; is orthogonal then:

MTM=IMT=M1\displaystyle M^TM=I\qquad M^T=M^{-1}

Note that if the eigenvectors are written in orthogonal form then the diagonalising equation is simplified:

A=(u1u2v1v2)(λ100λ2)(u1u2v1v2)1=(u1u2v1v2)(λ100λ2)(u1v1u2v2).\begin{align*} A &= \left(\begin{matrix}u_1 & u_2\\v_1 & v_2\\\end{matrix}\right) \left(\begin{matrix}\lambda_1 & 0\\0 & \lambda_2\\\end{matrix}\right) \left(\begin{matrix}u_1 & u_2\\v_1 & v_2\\\end{matrix}\right)^{-1} \\ &= \left(\begin{matrix}u_1 & u_2\\v_1 & v_2\\\end{matrix}\right) \left(\begin{matrix}\lambda_1 & 0\\0 & \lambda_2\\\end{matrix}\right) \left(\begin{matrix}u_1 & v_1\\u_2 & v_2\\\end{matrix}\right). \end{align*}

Summary

  • Matrix representation of simultaneous equations

  • Matrix-vector and matrix-matrix multiplication

  • Determinant, inverse and transpose

  • Null space of singular matrices

  • Finding eigenvalues and eigenvectors

  • Python for solving systems, finding inverse, null space and eigenvalues/vectors

  • Diagonalising matrices (we will use this for systems of differential equations)

Introductory problems

Introductory problems 1

Given

A=(100i);\displaystyle A = \left(\begin{array}{cc} 1 & 0 \\ 0 & i \end{array}\right); > B=(0ii0);\displaystyle B = \left(\begin{array}{cc} 0 & i \\ i & 0 \end{array}\right); > C=(1212(1i)12(1+i)12);\displaystyle C = \left(\begin{array}{cc} \frac{1}{\sqrt{2}} & \frac{1}{2}\left(1-i\right) \\\frac{1}{2}\left(1+i\right) & -\frac{1}{\sqrt{2}} \end{array}\right);

verify by hand, and using the numpy.linalg module, that

  1. AA1=A1A=IAA^{-1}=A^{-1}A =I;

  2. BB1=B1B=IBB^{-1}=B^{-1}B =I;

  3. CC1=C1C=ICC^{-1}=C^{-1}C =I;

# hint import numpy as np # In Python the imaginary unit is "1j" A = np.array([[1, 0], [0, 1j]]) print(A * np.linalg.inv(A))

Introductory problems 2

Let A be an n×nn \times n invertible matrix. Let II be the n×nn\times n identity matrix and let BB be an n×nn\times n matrix.

Suppose that ABA1=IABA^{-1}=I. Determine the matrix BB in terms of the matrix AA.

Main problems

Main problems 1

Let AA be the coefficient matrix of the system of linear equations

x12x2=1,2x1+3x2=1.\begin{aligned} -x_1 - 2 x_2 &= 1,\\ 2 x_1 + 3 x_2 &= -1. \end{aligned}
  1. Solve the system by finding the inverse matrix A1A^{-1}.

  2. Let x=(x1x2)\displaystyle \mathbf{x} = \left(\begin{array}{cc} x_1 \\ x_2 \end{array}\right) be the solution of the system obtained in part 1.

Calculate and simplify A2017xA^{2017} \mathbf{x}.

Main problems 2

For each of the following matrices

A=(2314);\displaystyle A = \left(\begin{array}{cc} 2 & 3 \\ 1 & 4 \end{array}\right); > B=(4268);\displaystyle B = \left(\begin{array}{cc} 4 & 2 \\ 6 & 8 \end{array}\right); > C=(1411);\displaystyle C = \left(\begin{array}{cc} 1 & 4 \\ 1 & 1 \end{array}\right); > D=(x00y),\displaystyle D = \left(\begin{array}{cc} x & 0 \\ 0 & y \end{array}\right),

compute the determinant, eigenvalues and eigenvectors by hand.

Check your results by verifying that Qx=λixQ\mathbf{x} = \lambda_i \mathbf{x}, where Q=AQ=A, BB, CC or DD, and by using the numpy.linalg module.

# hint import numpy as np A = np.array([[2, 3], [1, 4]]) e_vals, e_vecs = np.linalg.eig(A) print(e_vals) print(e_vecs)

Main problems 3

Orthogonal vectors.

Two vectors x1{\bf x_1} and x2{\bf x_2} are said to be perpendicular or orthogonal if their dot/scalar product is zero:

x1.x2=(213).(271)=(2×2)+(1×7)+(3×1)=47+3=0,\mathbf{x_1}.\mathbf{x_2}=\left(\begin{array}{c} 2 \\ -1 \\ 3 \\ \end{array} \right).\left(\begin{array}{c} 2 \\ 7 \\ 1 \\ \end{array} \right)=(2{\times}2)+(-1{\times}7)+(3{\times}1)=4-7+3=0,

thus x1\mathbf{x_1} and x2\mathbf{x_2} are orthogonal.

Find vectors that are orthogonal to (12)and(121).\displaystyle \left({\begin{array}{c} 1 \\ 2 \\ \end{array} } \right) \text{and} \left({\begin{array}{c} 1 \\ 2 \\ -1 \\ \end{array} } \right).

How many such vectors are there?

Main problems 4

Diagonalize the 2×22 \times 2 matrix A=(2112)\displaystyle A=\left(\begin{array}{cc} 2 & -1 \\ -1 & 2 \\ \end{array} \right) by finding a non-singular matrix SS and a diagonal matrix DD such that A=SDS1\displaystyle A = SDS^{-1}.

Main problems 5

Find a 2×22{\times}2 matrix AA such that A(21)=(14)andA(53)=(02)\displaystyle \quad A\left(\begin{array}{c} 2 \\ 1 \end{array} \right) = \left(\begin{array}{c} -1 \\ 4 \end{array} \right)\quad\text{and}\quad A \left(\begin{array}{c} 5 \\ 3\end{array} \right)=\left(\begin{array}{c} 0 \\ 2\end{array} \right).

Extension problems

Extension problems 1

If there exists a matrix MM whose columns are those of normalised (unit length) and orthogonal vectors, prove that MTM=IM^TM=I which implies that MT=M1M^T=M^{-1}.