6. Linear Algebra

Linear Algebra is a discipline where vector spaces and linear mapping between them are studied. In physics and astronomy, several phenomena can be readily written in terms of linear variables, what makes Computational Linear Algebra a very important topic to be covered throughout this course. We shall cover linear systems of equations, techniques for calculating inverses and determinants and factorization methods.

An interesting fact of Computational Linear Algebra is that it does not comprises numerical approaches as most of the methods are exact. The usage of a computer is then necessary because of the large number of calculations rather than the non-soluble nature of the problems. Numerical errors come then from round-off approximations.

See: http://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture14.pdf

---

Contents

• Linear Algebra

  • Matrices in Python

  • Gaussian Elimination

  • Linear Systems of Equations

  • Matrix Inversion

  • Determinant of a Matrix

  • Matrix diagonalization

  • LU Factorization

---

from IPython.display import HTML, Markdown, YouTubeVideo,Latex
%pylab inline

import numpy as np
import matplotlib.pyplot as plt
## JSAnimation import available at https://github.com/jakevdp/JSAnimation
#from JSAnimation import IPython_display
from matplotlib import animation
#Interpolation add-on
import scipy.interpolate as interp
xrange=range

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

https://stackoverflow.com/a/50748163

try:
    from google.colab.output._publish import javascript
    from IPython.display import Math as math
    from IPython.display import Latex as latex
    url = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/latest.js?config=default"
    def Math(*args,**kwargs):
        javascript(url=url)
        return math(*args,**kwargs)
    def Latex(*args,**kwargs):
        #print(args[0].replace('$',''))
        javascript(url=url)
        return math(args[0],**kwargs)
except:
    from IPython.display import Math, Latex

---

6.1. Matrices in Python

One of the most useful advantages of high level languages like Python, is the manipulation of complex objects like matrices and vectors. For this part we are going to use advanced capabilities for handling matrices, provided by the library NumPy.

NumPy, besides the extreme useful NumPy array objects, which are like tensors of any rank, also provides the Matrix objects that are overloaded with proper matrix operations.

6.1.1. Numpy arrays

• A matrix can be represented by an array of two indices

#NumPy Arrays
M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )
M2 = np.array( a[[3,-2,1],[-1,5,4],[1,-2,3]] )
print( f'M1=\n{M1}, with shape={M1.shape},\n\nM2=\n{M2}' )

M1=
[[ 5 -4  0]
 [-4  7 -3]
 [ 0 -3  5]], with shape=(3, 3),

M2=
[[ 3 -2  1]
 [-1  5  4]
 [ 1 -2  3]]

Note tha the previous definitions correspond to integer arrays.

WARNING: Be careful with integer arrays. In the next example, a single float entry convert the full matrix from an integer array to a float array:

M1[0,0]

5

float

np.array([[ 5, -4,  0],
          [-4,  7, -3.],
          [ 0, -3,  5]])#[0,0]

array([[ 5., -4.,  0.],
       [-4.,  7., -3.],
       [ 0., -3.,  5.]])

M1/3

array([[ 1.66666667, -1.33333333,  0.        ],
       [-1.33333333,  2.33333333, -1.        ],
       [ 0.        , -1.        ,  1.66666667]])

(M1/3).astype(int)

array([[ 1, -1,  0],
       [-1,  2, -1],
       [ 0, -1,  1]])

6.1.1.1. Special arrays:

Let $n$ the range of the matrix

• zero matrix

n=3
np.zeros( (n,n) )

array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

• Ones matrix

np.ones( (n,n) )

array([[1., 1., 1.],
       [1., 1., 1.],
       [1., 1., 1.]])

np.ones( (n,n) )*1j

array([[0.+1.j, 0.+1.j, 0.+1.j],
       [0.+1.j, 0.+1.j, 0.+1.j],
       [0.+1.j, 0.+1.j, 0.+1.j]])

• Identity matrix

np.identity(n)

array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

np.diag([1,2,4])

array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 4]])

• Random matrix with entries between 0 and 1

np.random.random()

0.9460452544468382

#np.random.seed(986554575)
np.random.random((n,n))

array([[0.05016015, 0.81396332, 0.17984221],
       [0.46070592, 0.46019579, 0.9506249 ],
       [0.28245888, 0.83429461, 0.47399034]])

np.random.uniform(0,10,(n,n))

array([[5.16683546, 5.73486742, 4.15495224],
       [0.75064986, 0.18470119, 3.9795813 ],
       [0.59178719, 6.74641986, 7.40144976]])

• Integer Random matrix

np.random.randint(0,10,(n,n)) #0 to 9

array([[9, 1, 9],
       [1, 8, 4],
       [0, 4, 2]])

np.random.randint(-10,10,(n,n))

array([[  4,  -3,  -5],
       [  7,   2,  -9],
       [-10,   5,  -8]])

6.1.1.2. Analytical matrices

Some times it is necessary to work with analytical Matrices. In such a case we can use the symbolic module Sympy

import sympy
x =sympy.Symbol('x') ## declare analytical varibles
sympy.init_printing() ## Use LaTeX to print sympy obejects

import numpy as np

#NumPy Arrays
M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )
M2 = np.array( [[3,-2,1],[-1,5,4],[1,-2,3]] )

sympy.Matrix(M2)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}3& -2& 1\\ -1& 5& 4\\ 1& -2& 3\end{array}\right]}\end{array}$$

In the following we will focus only in numerical matrices as numpy arrays, but also use sympy to print the matrices in a more readable way

6.1.1.3. Matrix operations (numpy)

Both as functions or array’s methods

• Transpose

$M_2^{\mathrm{T}}=$

from IPython import display

#Latex(r'M_2^{\rm T}=')
sympy.Matrix( M2.transpose())

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}3& -1& 1\\ -2& 5& -2\\ 1& 4& 3\end{array}\right]}\end{array}$$

• Matrix addition

sympy.Matrix(M1)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}5& -4& 0\\ -4& 7& -3\\ 0& -3& 5\end{array}\right]}\end{array}$$

sympy.Matrix(M2)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}3& -2& 1\\ -1& 5& 4\\ 1& -2& 3\end{array}\right]}\end{array}$$

sympy.Matrix( M1+M2 )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}8& -6& 1\\ -5& 12& 1\\ 1& -5& 8\end{array}\right]}\end{array}$$

Complex arrays are allowed

Mc=M1+M2*1j
sympy.Matrix( M1+M2*1j )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}5.0+3.0i& -4.0-2.0i& 1.0i\\ -4.0-1.0i& 7.0+5.0i& -3.0+4.0i\\ 1.0i& -3.0-2.0i& 5.0+3.0i\end{array}\right]}\end{array}$$

with the corresponding complex matrix operations

display ( sympy.Matrix(  Mc.conjugate() ) )
display ( sympy.Matrix(  Mc.conjugate().transpose() ) ) ## hermitian-conjugate

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}5.0-3.0i& -4.0+2.0i& -1.0i\\ -4.0+1.0i& 7.0-5.0i& -3.0-4.0i\\ -1.0i& -3.0+2.0i& 5.0-3.0i\end{array}\right]}\end{array}$$

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}5.0-3.0i& -4.0+1.0i& -1.0i\\ -4.0+2.0i& 7.0-5.0i& -3.0+2.0i\\ -1.0i& -3.0-4.0i& 5.0-3.0i\end{array}\right]}\end{array}$$

• Matrix multiplication

dot(a, b, out=None)

Dot product of two arrays. Specifically,

• If both a and b are 1-D arrays, it is inner product of vectors (without complex conjugation).

• If both a and b are 2-D arrays, it is matrix multiplication, but using :func:matmul or a @ b is preferred.

• If either a or b is 0-D (scalar), it is equivalent to :func:multiply and using numpy.multiply(a, b) or a * b is preferred.

• …

M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )
M2 = np.array( [[3,-2,1],
                [-1,5,4],
                [1,-2,3]] )

#Multiplication
sympy.Matrix( np.dot( M1, M2 ) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}19& -30& -11\\ -22& 49& 15\\ 8& -25& 3\end{array}\right]}\end{array}$$

sympy.Matrix( M1.dot(M2) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}19& -30& -11\\ -22& 49& 15\\ 8& -25& 3\end{array}\right]}\end{array}$$

Recommended way

sympy.Matrix(M1@M2)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}19& -30& -11\\ -22& 49& 15\\ 8& -25& 3\end{array}\right]}\end{array}$$

np.array( [1,2])@np.array( [[1],[2]])

array([5])

#Multiplication is not commutative
sympy.Matrix ( M2@M1 )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}23& -29& 11\\ -25& 27& 5\\ 13& -27& 21\end{array}\right]}\end{array}$$

• Complex matrices

σ_2=np.array(
    [[0,-1j],
     [1j,0]])
sympy.Matrix(σ_2)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cc}0& -1.0i\\ 1.0i& 0\end{array}\right]}\end{array}$$

• Trace

np.trace(σ_2)

0j

• Determinant

np.linalg.det(σ_2)

(-1+0j)

Raise a square matrix to the (integer) power n.

sympy.Matrix ( np.linalg.matrix_power(σ_2,2) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cc}1.0& 0\\ 0& 1.0\end{array}\right]}\end{array}$$

(σ_2@σ_2).astype(float)

/tmp/ipykernel_12325/501469815.py:1: ComplexWarning: Casting complex values to real discards the imaginary part
  (σ_2@σ_2).astype(float)

array([[1., 0.],
       [0., 1.]])

odd power

sympy.Matrix ( np.linalg.matrix_power(σ_2,7) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cc}0& -1.0i\\ 1.0i& 0\end{array}\right]}\end{array}$$

even power

sympy.Matrix ( np.linalg.matrix_power(σ_2,14) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cc}1.0& 0\\ 0& 1.0\end{array}\right]}\end{array}$$

Inverse of a matrix

sympy.Matrix( np.linalg.inv(M2) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}0.442307692307692& 0.0769230769230769& -0.25\\ 0.134615384615385& 0.153846153846154& -0.25\\ -0.0576923076923077& 0.0769230769230769& 0.25\end{array}\right]}\end{array}$$

sympy.Matrix ( ( M2@np.linalg.inv(M2) ).round(15) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}1.0& 0& 0\\ 0& 1.0& 0\\ 0& 0& 1.0\end{array}\right]}\end{array}$$

( M2@np.linalg.inv(M2) ).round(15).astype(int)

array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])

sympy.Matrix ( ( M2@np.linalg.inv(M2) ).round(15).astype(int) )

6.1.1.4. Element by element operations

print(M1)
M2

[[ 5 -4  0]
 [-4  7 -3]
 [ 0 -3  5]]

array([[ 3, -2,  1],
       [-1,  5,  4],
       [ 1, -2,  3]])

M1*M2

array([[ 15,   8,   0],
       [  4,  35, -12],
       [  0,   6,  15]])

#Division
print( M1/M2 )

[[ 1.66666667  2.          0.        ]
 [ 4.          1.4        -0.75      ]
 [ 0.          1.5         1.66666667]]

sympy.Matrix( M1@np.linalg.inv(M2).round(3) )

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}1.67& -0.231& -0.25\\ -0.649& 0.539& -1.5\\ -0.695& -0.077& 2.0\end{array}\right]}\end{array}$$

6.1.1.5. Submatrix from matrix:

M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )

Elements of the matrix

M1[1,2]

-3

6.1.1.6. Extract rows

M1[1,:]

array([-4,  7, -3])

or

M1[1]

array([-4,  7, -3])

In general: M1[list of rows to extract]

M1[[1]]

array([[-4,  7, -3]])

M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )

M1[[0,2]]

array([[ 5, -4,  0],
       [ 0, -3,  5]])

6.1.1.7. Extract columns

M1[:,1]

array([-4,  7, -3])

and convert back into a column

np.c_[M1[:,1]]

array([[-4],
       [ 7],
       [-3]])

6.1.1.8. Reordering

M1 = np.array( [[ 5 ,-4, 0],
                [-4 , 7,-3],
                [ 0 ,-3, 5]] )

Reverse row order to 2,1,0

M1[[2,1,0]]

array([[ 0, -3,  5],
       [-4,  7, -3],
       [ 5, -4,  0]])

or

np.r_[[M1[2]],[M1[1]],[M1[0]]]

array([[ 0, -3,  5],
       [-4,  7, -3],
       [ 5, -4,  0]])

Reverse column order to 2,1,0

np.c_[ M1[:,2],M1[:,1],M1[:,0]  ]

array([[ 0, -4,  5],
       [-3,  7, -4],
       [ 5, -3,  0]])

6.1.1.9. It is possible to reasign a full set of rows to a matrix

M1 = np.array( [[ 5,-4, 0],
                [-4, 7,-3],
                [ 0,-3, 5]] )

M1[[0,2]]=[[0  ,0,0.3],
           [0.3,0,0  ]]
M1

array([[ 0,  0,  0],
       [-4,  7, -3],
       [ 0,  0,  0]])

Note that because M1 has only integer numbers, it is assumed to be an integer matrix. In order to allow float assignments, it is first necessary to convert to a float matrix:

M1=M1.astype(float)
M1[[0,2]]=[[0  ,0,0.3],
           [0.3,0,0  ]]
M1

array([[ 0. ,  0. ,  0.3],
       [-4. ,  7. , -3. ],
       [ 0.3,  0. ,  0. ]])

WARNING
In some cases a reasignment of a list or array points to the same space in memory. To really copy an array to other variable it is always recommended to use the .copy() method.

A=np.array([1,2])
A

array([1, 2])

Reassign memory space

C=A
C

array([1, 2])

A change is reflected in both variables

C[0]=5
print(f'A={A} → C={C}')

A=[5 2] → C=[5 2]

A[1]=8
print(f'A={A} → C={C}')

A=[5 8] → C=[5 8]

To keep the first memory space

B=A.copy() ## two different memory spaces
B[0]=4
print('A =',A)
print('B =',B)

A = [5 2]
B = [4 2]

6.1.1.10. Matrix object

numpy also has a Matrix object with simplified operations. However it is recommended to use the general array object even for specific matrix operations. This helps to avoid incompatibilities with usual array operations. For example, as shown below the multiplication symbol, *, behaves different for arrays that for matrices, and could induce to errors.

#NumPy Matrix
M1 = np.matrix( [[5,-4,0],[-4,7,-3],[0,-3,5]] )
M2 = np.matrix( [[3,-2,1],[-1,5,4],[1,-2,3]] )

print (M1, "\n")
print (M2)

[[ 5 -4  0]
 [-4  7 -3]
 [ 0 -3  5]] 

[[ 3 -2  1]
 [-1  5  4]
 [ 1 -2  3]]

#Addition
sympy.Matrix( M1+M2 )

#Multiplication
sympy.Matrix( M1*M2 )

6.1.2. Basic operations with matrices

In order to simplify the following methods, we introduce here 3 basic operations over linear systems:

1. The $i$-th row $E_i$ can be multiplied by a non-zero constant $\lambda$, and then a new row used in place of $E_i$, i.e. $(E_i)\to (\lambda E_i)$. We denote this operation as $\mathrm{Lamb}(E_i,\lambda )$.

2. The $j$-th row $E_j$ can be multiplied by a non-zero constant $\lambda$ and added to some row $E_i$. The resulting value take the place of $E_i$, i.e. $(E_i)\to (E_i+\lambda E_j)$. We denote this operation as $\mathrm{Comb}(E_i,E_j,\lambda )$.

3. Finally, we define a swapping of two rows as $(E_i)\leftrightarrow (E_j)$, denoted as $\mathrm{Swap}(E_i,E_j)$

6.1.2.1. Extract rows from an array

6.1.3. Activity

Write three routines that perform, given a matrix $A$, the previous operations over matrices:

• row_lamb(i, λ, A): i is the row to be changed, λ the multiplicative factor and A the matrix. This function should return the new matrix with the performed operation $(E_i)\to (\lambda E_i)$.

• row_comb(i, j, λ, A): i is the row to be changed, j the row to be added, λ the multiplicative factor and A the matrix. This function should return the new matrix with the performed operation $(E_i)\to (E_i+\lambda E_j)$.

• row_swap(i, j, A): i and j are the rows to be swapped. This function should return the new matrix with the performed operation $(E_i)\leftrightarrow (E_j)$.

import numbers

isinstance(np.array([2])[0],numbers.Integral)

True

M1=M1.astype(float)
M1

array([[ 2., -4.,  0.],
       [-4.,  7., -3.],
       [ 0., -3.,  5.]])

M1.dtype==np.dtype('int64')

False

def row_sawp(i,j,A):
    '''
    i and j are the rows to be swapped
    '''
    B=A.copy()
    B[[i,j]]=B[[j,i]]
    return B
    
def row_comb(i,j,λ,A):
    '''
    `i` is the row to be changed, 
    `j` the row to be added
    `λ` the multiplicative factor 
    `A` the matrix.
    operation
    ---------
    E_i→ E_i+λ*E_j 
    '''
    B=A.copy()    
    if B.dtype==np.dtype('int64'):
        B=B.astype(float)
    B[i]=B[i]+λ*B[j]
    return B
    
def row_lamb(i,λ,A):
    '''
    `i` is the row to be changed, 
    `λ` the multiplicative factor 
    `A` the matrix.
    operation
    ---------
    λ → λ*E_i 
    '''
    B=A.copy()    
    if B.dtype==np.dtype('int64'):
        B=B.astype(float)
    B[i]=λ*B[i]
    return B

Check with

import numpy as np
import sympy

M1=np.array( [[5,-4,0],[-4,7,-3],[0,-3,5]] )
M1

array([[ 5, -4,  0],
       [-4,  7, -3],
       [ 0, -3,  5]])

row_lamb(1,3.4,M1)

array([[  5. ,  -4. ,   0. ],
       [-13.6,  23.8, -10.2],
       [  0. ,  -3. ,   5. ]])

row_comb(0,1,5/4,M1)

array([[ 0.  ,  4.75, -3.75],
       [-4.  ,  7.  , -3.  ],
       [ 0.  , -3.  ,  5.  ]])

row_sawp(0,1,M1)

array([[-4,  7, -3],
       [ 5, -4,  0],
       [ 0, -3,  5]])

6.2. Gaussian Elimination

6.2.1. Example 1

Using Kirchhoff’s circuit laws, it is possible to solve the next system:

obtaining the next equations:

(6.1)$$\begin{array}{rlrl}A:& & I_A(R+4R)-I_B(4R)=& E\\ C:& & -I_A(4R)+I_B(4R+3R)-I_C(3R)=& 0\\ B:& & -I_B(3R)+I_C(3R+2R)=& -2E\end{array}$$

or

(6.2)$$\begin{array}{rl}5\frac{I_{A}R}{E}-4\frac{I_{B}R}{E}=& 1\\ -4\frac{I_{A}R}{E}+7\frac{I_{B}R}{E}-3\frac{I_{C}R}{E}=& 0\\ -3\frac{I_{B}R}{E}+5\frac{I_{C}R}{E}=& -2.\end{array}$$

Defining the variables $x_1=R{I}_A/E$, $x_2=R{I}_B/E$ and $x_3=R{I}_C/E$ we have

(6.3)$$\begin{array}{rl}5x_1-4x_2=& 1\\ -4x_1+7x_2-3x_3=& 0\\ -3x_2+5x_3=& -2.\end{array}$$

A first method to solve linear systems of equations is the Gaussian elimination. This procedure consists of a set of recursive steps performed in order to diagonalise the matrix of the problem. A suitable way to introduce this method is applying it to some basic problem. To do so, let’s take the result of the Example 1:

$$\begin{array}{r}\left[\begin{array}{ccc}5& -4& 0\\ -4& 7& -3\\ 0& -3& 5\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ -2\end{array}\right]\end{array}$$

Constructing the associated augmented matrix, we obtain

$$\begin{array}{r}\left[\begin{array}{ccccc}5& -4& 0& ⋮& 1\\ -4& 7& -3& ⋮& 0\\ 0& -3& 5& ⋮& -2\end{array}\right]\end{array}$$

6.2.1.1. Activity

Use numpy functions to build the augmented matrix M1

M1 = np.array( [[5,-4,0],[-4,7,-3],[0,-3,5]] )
sympy.Matrix(M1)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{ccc}5& -4& 0\\ -4& 7& -3\\ 0& -3& 5\end{array}\right]}\end{array}$$

np.array( [[1],[0],[-2]] )

array([[ 1],
       [ 0],
       [-2]])

M1_aug=np.c_[ M1, [[1],[0],[-2]] ]
sympy.Matrix(M1_aug)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cccc}5& -4& 0& 1\\ -4& 7& -3& 0\\ 0& -3& 5& -2\end{array}\right]}\end{array}$$

or

M1_aug=np.hstack( (M1,[[1],[0],[-2]]))
sympy.Matrix(M1_aug)

$$\begin{array}{r}{\displaystyle \left[\begin{array}{cccc}5& -4& 0& 1\\ -4& 7& -3& 0\\ 0& -3& 5& -2\end{array}\right]}\end{array}$$

At this point, we can eliminate the coefficients of the variable $x_1$ in the second and third equations. For this, we apply the operations $\mathrm{Comb}(E_2,E_1,4/5)$.

$E_2$ → [-4,7,-3⋮0]+[4,-16/5,0⋮4/5]=[0,19/5,-3⋮4/5]

The coefficient in the third equation is already null. We then obtain:

$$\begin{array}{r}\left[\begin{array}{ccccc}5& -4& 0& ⋮& 1\\ 0& 19/5& -3& ⋮& 4/5\\ 0& -3& 5& ⋮& -2\end{array}\right]\end{array}$$

M1_LU=row_comb(1,0,4/5,M1_aug)
M1_LU

array([[ 5. , -4. ,  0. ,  1. ],
       [ 0. ,  3.8, -3. ,  0.8],
       [ 0. , -3. ,  5. , -2. ]])

Now, we proceed to eliminate the coefficient of $x_2$ in the third equation. For this, we apply again $\mathrm{Comb}(E_3,E_2,3\cdot 5/19)$

$E_3$ → [0,-3,5⋮-2]+[0,(3*5/19)*(19/5),-(3*5/19)*3⋮(3*5/19)*4/5]→[0,0,50/19⋮-26/19]

$$\begin{array}{r}\left[\begin{array}{ccccc}5& -4& 0& ⋮& 1\\ 0& 16/5& -3& ⋮& 4/5\\ 0& 0& 50/19& ⋮& -26/19\end{array}\right]\end{array}$$

M1_LU=row_comb(2,1,3*5/19,M1_LU)
M1_LU

array([[ 5.        , -4.        ,  0.        ,  1.        ],
       [ 0.        ,  3.8       , -3.        ,  0.8       ],
       [ 0.        ,  0.        ,  2.63157895, -1.36842105]])

M1_LU

array([[ 5.        , -4.        ,  0.        ,  1.        ],
       [ 0.        ,  3.8       , -3.        ,  0.8       ],
       [ 0.        ,  0.        ,  2.63157895, -1.36842105]])

The final step is to solve for $x_3$ in the last equation, doing the operation $\mathrm{Lamb}(E_3,19/50)$, yielding:

$$\begin{array}{r}\left[\begin{array}{ccccc}5& -4& 0& ⋮& 1\\ 0& 3.8& -3& ⋮& 0.8\\ 0& 0& 1& ⋮& -26/50\end{array}\right]\end{array}$$

From this, it is direct to conclude that $x_3=-26/50=-0.52$, for $x_2$ and $x_1$ it is only necessary to replace the found value of $x_3$.

row_lamb(2,19/50,M1_LU)

array([[ 5.  , -4.  ,  0.  ,  1.  ],
       [ 0.  ,  3.8 , -3.  ,  0.8 ],
       [ 0.  ,  0.  ,  1.  , -0.52]])

The full implementation in numpy is with np.linalg.solve

np.linalg.solve?

np.linalg.solve(M1,[1,0,-2])

array([ 0.04, -0.2 , -0.52])

np.linalg.solve(M1,[[1],[0],[-2]])

array([[ 0.04],
       [-0.2 ],
       [-0.52]])

6.3. Linear Systems of Equations

A linear system is a set of equations of $n$ variables that can be written in a general form:

$$a_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1$$

$$⋮$$

$$a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mn}{x}_n=b_n$$

where $n$ is, again, the number of variables, and $m$ the number of equations.

A linear system has solution if and only if $m\ge n$. This leads us to the main objetive of a linear system, find the set $\{x_i{\}}_{i=0}^n$ that fulfills all the equations.

Although there is an intuitive way to solve this type of systems, by just adding and subtracting equations until reaching the desire result, the large number of variables of some systems found in physics and astronomy makes necessary to develop iterative and general approaches. Next, we shall introduce matrix and vector notation and some basic operations that will be the basis of the methods to be developed in this section.

Quadratic eqution…

6.3.1. Matrices and vectors

A $m\times n$ matrix can be defined as a set of numbers arranged in columns and rows such as:

$$\begin{array}{r}A=[a_{ij}]=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]\end{array}$$

In the same way, it is possible to define a $n$-dimensional column vector as

$$\begin{array}{r}\mathit{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]\end{array}$$

and a column vector of constant coefficients

$$\begin{array}{r}\mathit{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]\end{array}$$

The system of equations

$$a_{11}{x}_1+a_{12}{x}_2+\cdots a_{1n}{x}_n=b_1$$

$$⋮$$

$$a_{m1}{x}_1+a_{m2}{x}_2+\cdots a_{mn}{x}_n=b_n$$

can be then written in a more convenient way as

$$\begin{array}{r}A\mathbf{\text{x}}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\mathit{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]\end{array}$$

We can also introducing the $n\times (n+1)$ augmented matrix as

$$\begin{array}{r}[A:\mathbf{\text{b}}]=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& ⋮& b_1\\ a_{21}& a_{22}& \cdots & a_{2n}& ⋮& b_2\\ ⋮& ⋮& & ⋮& ⋮& ⋮\\ a_{m1}& a_{n2}& \cdots & a_{nn}& ⋮& b_n\end{array}\right]\end{array}$$

6.3.1.1. Implementation in Scipy

import scipy.linalg as linalg

linalg.lu?

M1_aug

array([[ 5, -4,  0,  1],
       [-4,  7, -3,  0],
       [ 0, -3,  5, -2]])

P,L,U=linalg.lu(M1_aug)

U

array([[ 5.        , -4.        ,  0.        ,  1.        ],
       [ 0.        ,  3.8       , -3.        ,  0.8       ],
       [ 0.        ,  0.        ,  2.63157895, -1.36842105]])

$x_3$ is now

x3=U[2,3]/U[2,2]
'x3=-26/50={:.2f}'.format(x3)

'x3=-26/50=-0.52'

x2=(U[1,3]-U[1,2]*x3)/U[1,1]
round(x2,2)

x1=( U[0,3]-U[0,1]*x2-U[0,2]*x3 )/U[0,0]
round(x1,2)

6.3.2. General Gaussian elimination

Now, we shall describe the general procedure for Gaussian elimination:

1. Give an augmented matrix $\hat{A}=[A:\mathbf{\text{b}}]$.

2. Find the first non-zero coefficient $a_{i1}$ associated to $x_1$. This element is called pivot.

3. Apply the operation $\mathrm{Swap}(E_1,E_i)$. This guarantee the first row has a non-zero coefficient $a_{11}$.

4. Apply the operation $\mathrm{Comb}(E_j,E_1,-a_{j1}/a_{11})$. This eliminates the coefficients associated to $x_1$ in all the rows but in the first one.

5. Repeat steps 2. to 4. but for the coefficients of $x_2$ and then $x_3$ and so up to $x_n$. When iterating the coefficients of $x_k$, do not take into account the first $k$-th rows as they are already sorted.

6. Once you obtain a diagonal form of the matrix, apply the operation $\mathrm{Lamb}(E_n,1/a_{nn})$. This will make the coefficient of $x_n$ in the last equation equal to 1.

7. The final result should be an augmented matrix of the form: $$\left[\begin{array}{ccccccc}{a}_{11}& a_{12}& \cdots & a_{1(n-1)}& a_{1n}& ⋮& {\hat{b}}_1\\ 0& a_{22}& \cdots & a_{2(n-1)}& a_{2n}& ⋮& {\hat{b}}_2\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & a_{(n-1)(n-1)}& a_{(n-1)n}& ⋮& {\hat{b}}_{n-1}\\ 0& 0& \cdots & 0& 1& ⋮& {\hat{b}}_n\end{array}\right]$$

8. The solution to the problem is then obtained through backward substitutions, i.e. $$x_n={\hat{b}}_n$$$$x_{n-1}=\frac{{\hat{b}}_{n-1}-a_{(n-1)n}{x}_n}{a_{(n-1)(n-1)}}$$$$⋮$$$$x_i=\frac{{\hat{b}}_i-\sum _{j=i+1}^n{a}_{ij}{x}_j}{a_{ii}}\text{for}i=n-1,n-2,\cdots ,1$$

See custom implementation of general Gaussian elimination

#Gaussian Elimination
import scipy
def Gaussian_Elimination( A0 ):
    #Making local copy of matrix
    P,L,A=scipy.linalg.lu(A0)
    n = len(A)    
    #Finding solution
    x = np.zeros( n )
    x[n-1] = A[n-1,n]
    for i in range( n-1, -1, -1 ):
        x[i] = ( A[i,n] - sum(A[i,i+1:n]*x[i+1:n]) )/A[i,i]
    
    #Upper diagonal matrix and solutions x
    return A, x

np.random.seed(3)
M =  np.random.random( (4,5) )
M

array([[0.5507979 , 0.70814782, 0.29090474, 0.51082761, 0.89294695],
       [0.89629309, 0.12558531, 0.20724288, 0.0514672 , 0.44080984],
       [0.02987621, 0.45683322, 0.64914405, 0.27848728, 0.6762549 ],
       [0.59086282, 0.02398188, 0.55885409, 0.25925245, 0.4151012 ]])

U,x=Gaussian_Elimination(M)

Upper triangular

U

array([[ 0.89629309,  0.12558531,  0.20724288,  0.0514672 ,  0.44080984],
       [ 0.        ,  0.63097203,  0.16354802,  0.47919953,  0.62205662],
       [ 0.        ,  0.        ,  0.52490983, -0.06699671,  0.21531002],
       [ 0.        ,  0.        ,  0.        ,  0.32582312,  0.00303691]])

solutions

x

array([0.27395594, 0.8721633 , 0.41137442, 0.00932074])

M1_aug

array([[ 5, -4,  0,  1],
       [-4,  7, -3,  0],
       [ 0, -3,  5, -2]])

Gaussian_Elimination(M1_aug)

(array([[ 5.        , -4.        ,  0.        ,  1.        ],
        [ 0.        ,  3.8       , -3.        ,  0.8       ],
        [ 0.        ,  0.        ,  2.63157895, -1.36842105]]),
 array([ 0.04, -0.2 , -0.52]))

6.3.3. Pivoting Strategies

The previous method of Gaussian Elimination for finding solutions of linear systems is mathematically exact, however round-off errors that appear in computational arithmetics can represent a real problem when high accuracy is required.

In order to illustrate this, consider the next situation:

$$E_1:0.00300x_1+59.14x_2=59.17$$

$$E_2:5.291x_1-6.130x_2=46.78$$

Using four-digit arithmetic we obtain:

1. Constructing the augmented matrix:

$$\begin{array}{r}\left[\begin{array}{cccc}0.003& 59.14& ⋮& 59.17\\ 5.291& -6.130& ⋮& 46.78\end{array}\right]\end{array}$$

2. Applying the reduction with the first pivot, we obtain:

$$(E_1+m{E}_2)\to (E_1)$$

where:

$$m=-\frac{a_{21}}{a_{11}}=-\frac{5.291}{0.003}=1763.666\cdots \approx 1764$$

In this step, we have taken the first four digits. This leads us to:

$$\begin{array}{r}\left[\begin{array}{cccc}0.003& 59.14& ⋮& 59.17\\ 0& -104300& ⋮& -104400\end{array}\right]\end{array}$$

The exact system is instead

$$\begin{array}{r}\left[\begin{array}{cccc}0.003& 59.14& ⋮& 59.17\\ 0& -104309.37\overline{6}& ⋮& -104309.37\overline{6}\end{array}\right]\end{array}$$

Using the solution $x_2\approx 1.001$, we obtain

$$x_1\approx \frac{59.17-(59.14)(1.001)}{0.00300}=-10$$

The exact solution is however:

$$x_1=10.00$$

The source of such a large error is that the factor $59.14/0.00300\approx 20000$. This quantity is propagated through the combination steps of the Gaussian Elimination, yielding a complete wrong result.

6.3.4. Computing time

As mentioned before, Algebra Linear methods do not invole numerical approximations excepting round-off errors. This implies the computing time depends on the number of involved operations (multiplications, divisions, additions and subtractions). It is possible to demonstrate that the numbers of required divisions/multiplications is given by:

$$\frac{n^3}{3}+n^2-\frac{n}{3}$$

and the required additions/subtractions:

$$\frac{n^3}{3}+\frac{n^2}{2}-\frac{5n}{6}$$

These numbers are calculated separately as the computing time per operation for divisions and multiplications are similar and much larger than computing times for additions and subtractions. According to this, the overall computing time, for large $n$, scales as $n^3$/3.

6.3.5. Example 2

Using the library datetime of python, evaluate the computing time required for Gaussian_Elimination to find the solution of a system of $n=20$ equations and $20$ unknowns. For this purpose you can generate random systems. For a more robust result, repeat $500$ times, give a mean value and plot an histrogram of the computing times.

#Importing datatime
from datetime import datetime

#Number of repeats
Nrep = 500
#Size of matrix
n = 20

def Gaussian_Time( n, Nrep ):
    #Arrays of times
    Times = []
    #Cicle for number of repeats
    for i in xrange(Nrep):
        #Generating random matrix
        M = np.matrix( np.random.random( (n,n+1) ) )
        
        #Starting time counter
        tstart = datetime.now()
        #Invoking Gaussian Elimination routine
        Gaussian_Elimination(M)
        #Killing time counter
        tend = datetime.now()
        
        #Saving computing time
        Times.append( (tend-tstart).microseconds )
        
    #Numpy Array
    Times = np.array(Times)
    
    print "The mean computing time for a %dx%d matrix is: %lf microseconds"%(n,n,Times.mean())
    
    #Histrogram
    plt.figure( figsize=(8,5) )
    histo = plt.hist( Times, bins = 30 )
    plt.xlabel( "Computing Time [microseconds]" )
    plt.ylabel( "Ocurrences" )
    plt.grid()
    
    return Times.mean()
    
    
Gaussian_Time( n, Nrep )

The mean computing time for a 20x20 matrix is: 2762.628000 microseconds

2762.6280000000002

6.3.6. Activity

Using the previous code, estimate the computing time for random matrices of size $n=5,10,50,100,500,1000$. For each size, compute $500$ times in order to reduce spurious errors. Plot the results in a figure of $n$ vs computing time. Is it verified the some of scaling laws (Multiplication/Division or Addition/Sustraction). Note that for large values of $n$, both scaling laws converge to the same.

---

6.3.7. Partial pivoting

A suitable method to reduce round-off errors is to choose a pivot more conveniently. As we saw before, a small pivot generally implies larger propagated errors as they appear usually as dividends. The partial pivoting method consists then of a choosing of the largest absolute coefficient associated to $x_i$ instead of the first non-null one, i.e.

$$|a_{ii}|=\underset{i\le j\le n}{max}|a_{ji}|$$

This way, propagated multiplication errors would be minimized.

6.3.8. Activity

Create a new routine Gaussian_Elimination_Pivoting from Gaussian_Elimination in order to include the partial pivoting method. Compare both routines with some random matrix and with the exact solution.

---

6.4. Matrix Inversion

Asumming a nonsingular matrix $A$, if a matrix $A^{-1}$ exists, with $A{A}^{-1}=I$ and $A^{-1}A=I$, where $I$ is the identity matrix, then $A^{-1}$ is called the inverse matrix of $A$. If such a matrix does not exist, $A$ is said to be a singular matrix.

A corollary of this definition is that $A$ is also the inverse matrix of $A^{-1}$.

Once defined the Gaussian Elimination method, it is possible to extend it in order to find the inverse of any nonsingular matrix. Let’s consider the next equation:

$$\begin{array}{r}A{A}^{-1}=AB=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{m1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1n}\\ b_{21}& b_{22}& \cdots & b_{2n}\\ ⋮& ⋮& & ⋮\\ b_{n1}& b_{n2}& \cdots & b_{nn}\end{array}\right]=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right]\end{array}$$

This can be rewritten as a set of $n$ systems of equations, i.e.

$$\begin{array}{r}\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{c}{b}_{11}\\ b_{21}\\ ⋮\\ b_{n1}\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right],\end{array}$$

$$\begin{array}{r}\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{c}{b}_{12}\\ b_{22}\\ ⋮\\ b_{n2}\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ ⋮\\ 0\end{array}\right]\end{array}$$

$$⋮$$

$$\begin{array}{r}\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]\left[\begin{array}{c}{b}_{1n}\\ b_{2n}\\ ⋮\\ b_{nn}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 1\end{array}\right]\end{array}$$

These systems can be solved individually by using Gaussian Elimination, however we can mix all the problems, obtaining the augmented matrix:

$$\begin{array}{r}\left[\begin{array}{ccccccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& ⋮& 1& 0& \cdots & 0\\ a_{21}& a_{22}& \cdots & a_{2n}& ⋮& 0& 1& \cdots & 0\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}& ⋮& 0& 0& \cdots & 1\end{array}\right]\end{array}$$

Now, applying Gaussian Elimination we can obtain a upper diagonal form for the first matrix. Completing the steps using forwards elimination we can convert the first matrix into the identity matrix, obtaining

$$\begin{array}{r}\left[\begin{array}{ccccccccc}1& 0& \cdots & 0& ⋮& b_{11}& b_{12}& \cdots & b_{1n}\\ 0& 1& \cdots & 0& ⋮& b_{21}& b_{22}& \cdots & b_{2n}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 1& ⋮& b_{n1}& b_{n2}& \cdots & b_{nn}\end{array}\right]\end{array}$$

Where the second matrix is then the inverse $B=A^{-1}$.

6.4.1. Activity

Using the previous routine Gaussian_Elimination_Pivoting, create a new routine Inverse that calculates the inverse of any given squared matrix.

---

6.5. Determinant of a Matrix

The determinant of a matrix is a scalar quantity calculated for square matrix. This provides important information about the matrix of coefficients of a system of linear of equations. For example, any system of $n$ equations and $n$ unknowns has an unique solution if the associated determinant is nonzero. This also implies the determinant allows to evaluate whether a matrix is singular or nonsingular.

6.5.1. Calculating determinants

Next, we shall define some properties of determinants that will allow us to calculate determinants by using a recursive code:

1. If $A=[a]$ is a $1\times 1$ matrix, its determinant is then $detA=a$.

2. If $A$ is a $n\times n$ matrix, the minor matrix $M_{ij}$ is the determinant of the $(n-1)\times (n-1)$ matrix obtained by deleting the $i$th row and the $j$th column.

3. The cofactor $A_{ij}$ associated with $M_{ij}$ is defined by $A_{ij}=(-1{)}^{i+j}{M}_{ij}$.

4. The determinant of a $n\times n$ matrix $A$ is given by:

$$detA=\sum _{j=1}^n{a}_{ij}{A}_{ij}$$

or

$$detA=\sum _{i=1}^n{a}_{ij}{A}_{ij}$$

This is, it is possible to use both, a row or a column for calculating the determinant.

6.5.2. Computing time of determinants

Using the previous recurrence, we can calculate the computing time of the previous algorithm. First, let’s consider the number of required operations for a $2\times 2$ matrix: let $A$ be a $2\times 2$ matrix given by:

$$\begin{array}{r}A=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\end{array}$$

The determinant is then given by:

$$det(A)=a_{11}{a}_{22}-a_{12}{a}_{21}$$

the number of required multiplications was $2$ and subtractions is $1$.

Now, using the previous formula for the determinant

$$detA=\sum _{j=1}^n{a}_{ij}{A}_{ij}$$

For a $3\times 3$ matrix, it is necessary to calculate $3$ times $2\times 2$ determinants. Besides, it is necessary to multiply the cofactor $A_{ij}$ with the coefficient $a_{ij}$, that leads us with $t_{n=3}=3\times 2+3$ multiplications. Additions are not important as their computing time is far less than multiplications.

For a $4\times 4$ matrix, we need four deteminants of $3\times 3$ submatrices, leading $t_{n=4}=4\times (3\times 2+3)+4=4!+\frac{4!}{2!}+\frac{4!}{3!}$. In general, for a $n\times n$ matrix, we have then:

$$t_n=\frac{n!}{(n-1)!}+\frac{n!}{(n-2)!}+\cdots +\frac{n!}{1!}=n!\left(\sum _{i=1}^{n-1}\frac{1}{i!}\right)$$

If $n$ is large enough, we can approximate $t_n\approx n!$.

In computers, this is a prohibitive computing time so other schemes have to be proposed.

6.5.3. Activity

Evaluate the computing time of the Determinant routine for matrix sizes of $n=1,2,3,\cdots ,10$ and doing several repeats. Plot your result ($n$ vs $t_n$). What can you conclude about the behaviour of the computing time?

6.5.4. Properties of determinants

Determinants have a set of properties that can reduce considerably computing times. Suppose $A$ is a $n\times n$ matrix:

1. If any row or column of $A$ has only zero entries, then $detA=0$.

2. If two rows or columns of $A$ are the same, then $detA=0$.

3. If $\hat{A}$ is obtained from $A$ by using a swap operation $(E_i)\leftrightarrow (E_j)$, then $det\hat{A}=-detA$.

4. If $\hat{A}$ is obtained from $A$ by using a escalation operation $(\lambda E_i)\leftrightarrow (E_i)$, then $det\hat{A}=\lambda detA$.

5. If $\hat{A}$ is obtained from $A$ by using a combination operation $(E_i+\lambda E_j)\leftrightarrow (E_i)$, then $det\hat{A}=detA$.

6. If $B$ is also a $n\times n$ matrix, then $det(AB)=(detA)(detB).$

7. $det{A}^t=detA.$

8. $det{A}^{-1}=(detA{)}^{-1}$

9. Finally and most importantly: if $A$ is an upper, lower or diagonal matrix, then:

$$detA=\prod _{i=1}^n{a}_{ii}$$

As we analysed before, Gaussian Elimination takes a computing time scaling like $n^3$ for large matrix sizes. According to the previous properties, the determinant of a upper diagonal matrix just takes $n-1$ multiplications, far less than a nondiagonal matrix. Combining these properties, we can track back and relate the determinant of the resulting upper diagonal matrix and the original one. Leading us to a computing time scaling like $n^3$, much better than the original $n!$.

6.5.5. Activity

Using the Gaussian_Elimination routine and tracking back the performed operations, construct a new routine called Gaussian_Determinant. Make the same analysis of the computing time as the previous activity. Compare both results.

6.5.6. Existence of inverse

A matrix $A$ has an inverse if $detA\ne 0$. See for example here

If the matrix $A$ has an inverse, then

(6.4)$$\begin{array}{rl}A\mathit{x}=& \mathit{b}\\ A^{-1}A\mathit{x}=& A^{-1}\mathit{b}\\ \mathit{x}=& A^{-1}\mathit{b},\end{array}$$

Example From the previous example

np.random.seed(3)
M =  np.random.random( (4,5) )
M

array([[0.5507979 , 0.70814782, 0.29090474, 0.51082761, 0.89294695],
       [0.89629309, 0.12558531, 0.20724288, 0.0514672 , 0.44080984],
       [0.02987621, 0.45683322, 0.64914405, 0.27848728, 0.6762549 ],
       [0.59086282, 0.02398188, 0.55885409, 0.25925245, 0.4151012 ]])

M[:,0]

array([0.5507979 , 0.89629309, 0.02987621, 0.59086282])

np.c_[M[:,0]]

array([[0.5507979 ],
       [0.89629309],
       [0.02987621],
       [0.59086282]])

A=np.c_[ tuple( [ np.c_[M[:,i]]    for i in range(4) ] ) ]
A

array([[0.5507979 , 0.70814782, 0.29090474, 0.51082761],
       [0.89629309, 0.12558531, 0.20724288, 0.0514672 ],
       [0.02987621, 0.45683322, 0.64914405, 0.27848728],
       [0.59086282, 0.02398188, 0.55885409, 0.25925245]])

b=np.c_[ M[:,4] ]
b

array([[0.89294695],
       [0.44080984],
       [0.6762549 ],
       [0.4151012 ]])

such that

(6.5)$$\begin{array}{rl}A\mathit{x}=& \mathit{b},\end{array}$$

Check that $A$ has an inverse and calculate $\mathit{x}$

np.linalg.det(A)

$\mathit{x}=$

sympy.Matrix( np.dot( np.linalg.inv(A) , b).round(4) ) 

$$\begin{array}{r}{\displaystyle \left[\begin{array}{c}0.274\\ 0.8722\\ 0.4114\\ 0.0093\end{array}\right]}\end{array}$$

6.5.6.1. Numpy implementation

np.linalg.solve(A,b)

array([[0.27395594],
       [0.8721633 ],
       [0.41137442],
       [0.00932074]])

 np.linalg.solve(A,M[:,4])

array([0.27395594, 0.8721633 , 0.41137442, 0.00932074])

6.6. Matrix diagonalization

See linear-algebra-diagonalization

---

6.7. LU Factorization

As we saw before, the Gaussian Elimination algorithm takes a computing time scaling as $\mathcal{O}(n^3/3)$ in order to solve a system of $n$ equations and $n$ unknowns. Let’s assume a system of equations $\mathbf{A}\mathbf{x}=\mathbf{b}$ where $\mathbf{b}$ is already in upper diagonal form.

$$\begin{array}{r}\left[\begin{array}{ccccccc}{a}_{11}& a_{12}& \cdots & a_{1(n-1)}& a_{1n}& ⋮& {\hat{b}}_1\\ 0& a_{22}& \cdots & a_{2(n-1)}& a_{2n}& ⋮& {\hat{b}}_2\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & a_{(n-1)(n-1)}& a_{(n-1)n}& ⋮& {\hat{b}}_{n-1}\\ 0& 0& \cdots & 0& a_{nn}& ⋮& {\hat{b}}_n\end{array}\right]\end{array}$$

The Gauss-Jordan algorithm can reduce even more this problem in order to solve it directly, yielding:

$$\begin{array}{r}\left[\begin{array}{ccccccc}1& 0& \cdots & 0& 0& ⋮& x_1\\ 0& 1& \cdots & 0& 0& ⋮& x_2\\ ⋮& ⋮& & ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 1& 0& ⋮& x_{n-1}\\ 0& 0& \cdots & 0& 1& ⋮& x_n\end{array}\right]\end{array}$$

From the upper diagonal form to the completely reduced one, it is necessary to perform $n+(n-1)+(n-2)+\cdots \propto n(n-1)$ backwards substitutions. The computing time for solving a upper diagonal system is then $\mathcal{O}(n^2)$.

Now, let $\mathbf{A}\mathbf{x}=\mathbf{b}$ be a general system of equations of $n$ dimensions. Let’s assume $\mathbf{A}$ can be written as a multiplication of two matrices, one lower diagonal $\mathbf{L}$ and other upper diagonal $\mathbf{U}$, such that $\mathbf{A}=\mathbf{L}\mathbf{U}$. Defining a vector $\mathbf{y}=\mathbf{U}\mathbf{x}$, it is obtained for the original system

$$\mathbf{A}\mathbf{x}=\mathbf{L}(\mathbf{U}\mathbf{x})=\mathbf{L}\mathbf{y}=\mathbf{b}$$

For solving this system we can then:

1. Solve the equivalent system $\mathbf{L}\mathbf{y}=\mathbf{b}$, what takes a computing time of $\mathcal{O}(n^2)$.

2. Once we know $\mathbf{y}$, we can solve the system $\mathbf{U}\mathbf{x}=\mathbf{y}$, with a computing time of $\mathcal{O}(n^2)$.

The global computing time is then $\mathcal{O}(2n^2)$

6.7.1. Activity

In order to compare the computing time that Gaussian Elimination takes and the previous time for the LU factorization, make a plot of both computing times. What can you conclude when $n$ becomes large enough?

6.7.2. Derivation of LU factorization

Although the LU factorization seems to be a far better method for solving linear systems as compared with say Gaussian Elimination, we was assuming we already knew the matrices $\mathbf{L}$ and $\mathbf{U}$. Now we are going to see the algorithm for perfoming this reduction takes a computing time of $\mathcal{O}(n^3/3)$.

You may wonder then, what advantage implies the use of this factorization? Well, matrices $\mathbf{L}$ and $\mathbf{U}$ do not depend on the specific system to be solved, i.e. there is not dependence on the $\mathbf{b}$ vector, so once we have both matrices, we can use them to solve any system we want, just taking a $\mathcal{O}(2n^2)$ computing time.

First, let’s assume a matrix $\mathbf{A}$ with all its pivots are nonzero, so there is not need to swap rows. Now, when we want to eliminate all the coefficients associated to $x_1$, we perform the next operations:

$$(E_j-m_{j1}{E}_1)\to (E_j),\text{where}{m}_{j1}=\frac{a_{j1}^{(1)}}{a_{11}^{(1)}}$$

henceforth, $a_{ij}^{(1)}$ denotes the components of the original matrix $\mathbf{A}={\mathbf{A}}^{(1)}$, $a_{ij}^{(2)}$ the components of the matrix after eliminating the coefficients of $x_1$, and generally, $a_{ij}^{(k)}$ the components of the matrix after eliminating the coefficients of $x_{k-1}$.

The previous operation over the matrix $\mathbf{A}$ can be also reproduced defining the matrix ${\mathbf{M}}^{(1)}$

$$\begin{array}{r}{\mathbf{M}}^{(1)}=\left[\begin{array}{ccccc}1& 0& \cdots & 0& 0\\ -m_{21}& 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ -m_{(n-1)1}& 0& \cdots & 1& 0\\ -m_{n1}& 0& \cdots & 0& 1\end{array}\right]\end{array}$$

This is called the first Gaussian transformation matrix. From this, we have

$${\mathbf{A}}^{(2)}\mathbf{x}={\mathbf{M}}^{(1)}{\mathbf{A}}^{(1)}\mathbf{x}={\mathbf{M}}^{(1)}{\mathbf{b}}^{(1)}={\mathbf{b}}^{(2)}$$

where ${\mathbf{A}}^{(2)}$ is matrix with null coefficients associated to $x_1$ but the first one.

Repeating the same procedure for the next pivots, we obtain then

$${\mathbf{A}}^{(n)}={\mathbf{M}}^{(n-1)}{\mathbf{M}}^{(n-2)}\cdots {\mathbf{M}}^{(1)}{\mathbf{A}}^{(1)}$$

where the $k$th Gaussian transformation matrix is defined as

$$\begin{array}{r}{\mathbf{M}}_{ij}^{(k)}=\{\begin{array}{cc}1& \text{if}i=j\\ -m_{ij}& \text{if}j=k\text{and}k+1\le i\le n\\ 0& \text{otherwise}\end{array}\end{array}$$

and

$$m_{ij}=\frac{a_{ij}^{(j)}}{a_{jj}^{(j)}}$$

Note ${\mathbf{A}}^{(n)}$ is a upper diagonal matrix given by

$$\begin{array}{r}{\mathbf{A}}^{(n)}=\left[\begin{array}{ccccc}{a}_{11}^{(n)}& a_{12}^{(n)}& \cdots & a_{1(n-1)}^{(n)}& a_{1n}^{(n)}\\ 0& a_{22}^{(n)}& \cdots & a_{2(n-1)}^{(n)}& a_{2n}^{(n)}\\ ⋮& ⋮& & ⋮& ⋮& \\ 0& 0& \cdots & a_{(n-1)(n-1)}^{(n)}& a_{(n-1)n}^{(n)}\\ 0& 0& \cdots & 0& a_{nn}^{(n)}\end{array}\right]\end{array}$$

so we can define $\mathbf{U}\equiv {\mathbf{A}}^{(n)}$.

Now, taking the equation

$${\mathbf{A}}^{(n)}={\mathbf{M}}^{(n-1)}{\mathbf{M}}^{(n-2)}\cdots {\mathbf{M}}^{(1)}{\mathbf{A}}^{(1)}$$

and defining the inverse of ${\mathbf{M}}^{(k)}$ as

$$\begin{array}{r}{\mathbf{L}}_{ij}^{(k)}={\left({\mathbf{M}}^{(k)}\right)}_{ij}^{-1}=\{\begin{array}{cc}1& \text{if}i=j\\ m_{ij}& \text{if}j=k\text{and}k+1\le i\le n\\ 0& \text{otherwise}\end{array}\end{array}$$

we obtain

$${\mathbf{L}}^{(1)}\cdots {\mathbf{L}}^{(n-2)}{\mathbf{L}}^{(n-1)}{\mathbf{A}}^{(n)}={\mathbf{L}}^{(1)}\cdots {\mathbf{L}}^{(n-2)}{\mathbf{L}}^{(n-1)}{\mathbf{M}}^{(n-1)}{\mathbf{M}}^{(n-2)}\cdots {\mathbf{M}}^{(1)}{\mathbf{A}}^{(1)}=\mathbf{L}\mathbf{U}$$

where the lower diagonal matrix $\mathbf{L}$ is given by:

$$\mathbf{L}={\mathbf{L}}^{(1)}\cdots {\mathbf{L}}^{(n-2)}{\mathbf{L}}^{(n-1)}$$.

6.7.3. Algorithm for LU factorization

The algorithm is then given by:

1. Give a square matrix $\mathbf{A}$ where the pivots are nonzero.

2. Apply the operation $Comb(E_j,E_1,-a_{j1}^{(1)}/a{a}_{11}^{(1)})$. This eliminates the coefficients associated to $x_1$ in all the rows but in the first one.

3. Construct the matrix ${\mathbf{L}}^{(1)}$ given by

$$\begin{array}{r}{\mathbf{L}}_{ij}^{(k)}=\{\begin{array}{cc}1& \text{if}i=j\\ m_{ij}=\frac{a_{ij}^{(j)}}{a_{jj}^{(j)}}& \text{if}j=k\text{and}k+1\le i\le n\\ 0& \text{otherwise}\end{array}\end{array}$$

with $k=1$.

4. Repeat the steps 2 and 3 for the next column until reaching the last one.

5. Return the matrices $\mathbf{U}={\mathbf{A}}^{(n)}$ and $\mathbf{L}={\mathbf{L}}^{(1)}\cdots {\mathbf{L}}^{(n-2)}{\mathbf{L}}^{(n-1)}$.

6.7.4. Activity

Create a routine called LU_Factorization that, given a matrix $\mathbf{A}$ and the previous algorithm, calculate the LU factorization of the matrix. Test your routine with a random square matrix, verify that $\mathbf{A}=\mathbf{L}\mathbf{U}$.

---

M1

array([[ 5, -4,  0],
       [-4,  7, -3],
       [ 0, -3,  5]])

import scipy

P,L,U=scipy.linalg.lu(M1)

U

array([[ 5.        , -4.        ,  0.        ],
       [ 0.        ,  3.8       , -3.        ],
       [ 0.        ,  0.        ,  2.63157895]])

The same obtained before

scipy.linalg.lu?

6.7.5. Eigenvalues and Eigenvectors activity

6.7.5.1. Electron interacting with a magnetic field

An electron is placed to interact with an uniform magnetic field. To give account of the possible allowed levels of the electron in the presence of the magnetic field it is necessary to solve the next equation

(6.6)$$\hat{H}|\mathrm{\Psi }⟩=E|\mathrm{\Psi }⟩$$

where the hamiltonian is equal to $H=-\mu \cdot B=-\gamma B\cdot S$, with the gyromagnetic ratio $\gamma$, $\mathbf{\text{B}}$ the magnetic field and $\mathbf{\text{S}}$ the spin. It can be shown that the hamiltonian expression is transformed in

(6.7)$$\widehat{\mathbf{\text{H}}}=-\frac{\gamma \hslash }{2}\left(\begin{array}{cc}{B}_z& B_x-i{B}_y\\ B_x+i{B}_y& -B_z\end{array}\right)$$

Then, by solving the problem $|H-EI|=0$ is found the allowed energy levels, i.e., finding the determinant of the matrix $H-EI$ allows to get the values $E_1$ and $E_2$.

And solving the problem $\hat{H}\mathrm{\Psi }$ - E$\mathrm{\Psi }=0$ gives the autofunctions $\mathrm{\Psi }$, i.e., the column vector $\mathrm{\Psi }=\{{\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2\}$.

The function scipy.optimize.root can be used to get roots of a given equation. The experimental value of $\gamma$ for the electron is 2. The order of magnitude of the magnetic field is $1g$.

1. Find the allowed energy levels.

2. Find the autofunctions and normalize them.

Hint: An imaginary number in python can be written as 1j.