GitHub Repository: restrepo/ComputationalMethods
Path: blob/master/material/linear-algebra.ipynb
⁹³⁴ views

Kernel: Python 3 (ipykernel)

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

:depth: 2

In [1]:

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

Out[1]:

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

https://stackoverflow.com/a/50748163

In [2]:

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

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.

Numpy arrays

A matrix can be represented by an array of two indices

In [3]:

#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}' )

Out[3]:

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:

In [4]:

M1[0,0]

Out[4]:

5

float

In [5]:

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

Out[5]:

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

In [6]:

M1/3

Out[6]:

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

In [7]:

(M1/3).astype(int)

Out[7]:

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

Special arrays:

Let $n$ the range of the matrix

zero matrix

In [8]:

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

Out[8]:

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

Ones matrix

In [9]:

np.ones( (n,n) )

Out[9]:

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

In [10]:

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

Out[10]:

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

In [11]:

np.identity(n)

Out[11]:

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

In [12]:

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

Out[12]:

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

In [ ]:

Random matrix with entries between 0 and 1

In [13]:

np.random.random()

Out[13]:

0.9460452544468382

In [14]:

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

Out[14]:

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

In [15]:

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

Out[15]:

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

Integer Random matrix

In [16]:

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

Out[16]:

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

In [17]:

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

Out[17]:

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

Analytical matrices

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

In [1]:

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

In [2]:

import numpy as np

In [3]:

#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]] )

In [4]:

sympy.Matrix(M2)

Out[4]:

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

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

Matrix operations (numpy)

Both as functions or array's methods

Transpose

$M_2^{\rm T}=$

In [21]:

from IPython import display

In [7]:

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

Out[7]:

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

Matrix addition

In [8]:

sympy.Matrix(M1)

Out[8]:

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

In [9]:

sympy.Matrix(M2)

Out[9]:

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

In [10]:

sympy.Matrix( M1+M2 )

Out[10]:

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

Complex arrays are allowed

In [11]:

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

Out[11]:

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

with the corresponding complex matrix operations

In [12]:

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

Out[12]:

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

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

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.
...

In [13]:

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

In [14]:

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

Out[14]:

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

In [15]:

sympy.Matrix( M1.dot(M2) )

Out[15]:

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

Recommended way

In [16]:

sympy.Matrix(M1@M2)

Out[16]:

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

In [22]:

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

Out[22]:

array([5])

In [23]:

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

Out[23]:

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

Complex matrices

In [24]:

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

Out[24]:

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

Trace

In [25]:

np.trace(σ_2)

Out[25]:

0j

Determinant

In [26]:

np.linalg.det(σ_2)

Out[26]:

(-1+0j)

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

In [27]:

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

Out[27]:

$\displaystyle \left[\begin{matrix}1.0 & 0\\0 & 1.0\end{matrix}\right]$

In [29]:

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

Out[29]:

/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

In [30]:

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

Out[30]:

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

even power

In [31]:

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

Out[31]:

$\displaystyle \left[\begin{matrix}1.0 & 0\\0 & 1.0\end{matrix}\right]$

Inverse of a matrix

In [32]:

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

Out[32]:

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

In [33]:

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

Out[33]:

$\displaystyle \left[\begin{matrix}1.0 & 0 & 0\\0 & 1.0 & 0\\0 & 0 & 1.0\end{matrix}\right]$

In [35]:

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

Out[35]:

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

In [52]:

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

Out[52]:

$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

Element by element operations

In [36]:

print(M1)
M2

Out[36]:

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

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

In [37]:

M1*M2

Out[37]:

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

In [38]:

#Division
print( M1/M2 )

Out[38]:

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

In [39]:

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

Out[39]:

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

Submatrix from matrix:

In [57]:

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

Elements of the matrix

In [58]:

M1[1,2]

Out[58]:

-3

Extract rows

In [59]:

M1[1,:]

Out[59]:

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

In [60]:

M1[1]

Out[60]:

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

In general: M1[list of rows to extract]

In [61]:

M1[[1]]

Out[61]:

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

In [62]:

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

In [63]:

M1[[0,2]]

Out[63]:

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

Extract columns

In [64]:

M1[:,1]

Out[64]:

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

and convert back into a column

In [65]:

np.c_[M1[:,1]]

Out[65]:

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

Reordering

In [53]:

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

Reverse row order to 2,1,0

In [66]:

M1[[2,1,0]]

Out[66]:

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

In [67]:

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

Out[67]:

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

Reverse column order to 2,1,0

In [68]:

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

Out[68]:

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

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

In [69]:

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

In [70]:

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

Out[70]:

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:

In [71]:

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

Out[71]:

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.

In [72]:

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

Out[72]:

array([1, 2])

Reassign memory space

In [73]:

C=A
C

Out[73]:

array([1, 2])

A change is reflected in both variables

In [74]:

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

Out[74]:

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

In [64]:

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

Out[64]:

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

To keep the first memory space

In [75]:

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

Out[75]:

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

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.

In [76]:

#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)

Out[76]:

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

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

In [77]:

#Addition
sympy.Matrix( M1+M2 )

Out[77]:

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

In [78]:

#Multiplication
sympy.Matrix( M1*M2 )

Out[78]:

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

Basic operations with matrices

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

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)\rightarrow (\lambda E_i)$ . We denote this operation as $\operatorname{Lamb}(E_i,\lambda)$ .
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) \rightarrow (E_i + \lambda E_j)$ . We denote this operation as $\operatorname{Comb}(E_i,E_j,\lambda)$ .
Finally, we define a swapping of two rows as $(E_i)\leftrightarrow (E_j)$ , denoted as $\operatorname{Swap}(E_i,E_j)$

Extract rows from an array

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)\rightarrow (\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) \rightarrow (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)$ .

In [30]:

import numbers

In [35]:

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

Out[35]:

True

In [43]:

M1=M1.astype(float)
M1

Out[43]:

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

In [44]:

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

Out[44]:

False

In [71]:

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

In [24]:

import numpy as np
import sympy

In [74]:

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

Out[74]:

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

In [59]:

row_lamb(1,3.4,M1)

Out[59]:

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

In [60]:

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

Out[60]:

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

In [75]:

row_sawp(0,1,M1)

Out[75]:

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

Gaussian Elimination

Example 1

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

obtaining the next equations:

\begin{align} 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{align}

or $\begin{align} 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{align}$

Defining the variables $x_1 = R I_A/E$ , $x_2 = R I_B/E$ and $x_3 = R I_C/E$ we have $\begin{align} 5x_1 - 4x_2 =& 1 \\ -4 x_1 + 7x_2 - 3 x_3 =& 0 \\ -3x_2 + 5 x_3 =& -2\,. \end{align}$

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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 5 & -4 & 0 \\…

Constructing the associated augmented matrix, we obtain

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 5 & -4 & 0 & …

Activity

Use numpy functions to build the augmented matrix M1

In [76]:

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

Out[76]:

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

In [77]:

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

Out[77]:

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

In [78]:

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

Out[78]:

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

In [35]:

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

Out[35]:

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

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 $\operatorname{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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 5 & -4 & 0 & …

In [79]:

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

Out[79]:

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 $\operatorname{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]

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 5 & -4 & 0 & …

In [80]:

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

Out[80]:

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

In [42]:

M1_LU

Out[42]:

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 $\operatorname{Lamb}(E_3,19/50)$ , yielding:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 5 & -4 & 0 & …

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$ .

In [81]:

row_lamb(2,19/50,M1_LU)

Out[81]:

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

In [41]:

np.linalg.solve?

In [82]:

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

Out[82]:

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

In [83]:

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

Out[83]:

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

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

\vdots

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\geq 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...

Matrices and vectors

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 23: …{ij}] = \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 26: …l{x} = \left[ \̲m̲a̲t̲r̲i̲x̲{ x_{1} \\ x_{2…

and a column vector of constant coefficients

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 25: …ol{b} = \left[ \̲m̲a̲t̲r̲i̲x̲{ b_{1} \\ b_{2…

The system of equations

a_{11}x_1 + a_{12}x_2 + \cdots a_{1n}x_n = b_1

\vdots

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 24: …f{x} = \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 26: …{b}] = \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

Implementation in Scipy

In [91]:

import scipy.linalg as linalg

In [ ]:

linalg.lu?

In [87]:

M1_aug

Out[87]:

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

In [93]:

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

In [94]:

Out[94]:

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

$x_3$ is now

In [95]:

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

Out[95]:

'x3=-26/50=-0.52'

In [67]:

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

Out[67]:

$\displaystyle -0.2$

In [68]:

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

Out[68]:

$\displaystyle 0.04$

General Gaussian elimination

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

Give an augmented matrix $\hat{A}=[A:\textbf{b}]$ .
Find the first non-zero coefficient $a_{i1}$ associated to $x_1$ . This element is called pivot.
Apply the operation $\operatorname{Swap}(E_1,E_i)$ . This guarantee the first row has a non-zero coefficient $a_{11}$ .
Apply the operation $\operatorname{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.
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.
Once you obtain a diagonal form of the matrix, apply the operation $\operatorname{Lamb}(E_n,1/a_{nn})$ . This will make the coefficient of $x_n$ in the last equation equal to 1.
The final result should be an augmented matrix of the form: ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…
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)}}$ $\vdots$ $x_i = \frac{\hat b_i - \sum_{j=i+1}^n a_{ij}x_j}{a_{ii}}\ \ \ \ \mbox{for}\ \ \ i=n-1, n-2, \cdots, 1$

See custom implementation of general Gaussian elimination

In [105]:

#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

In [106]:

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

Out[106]:

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 ]])

In [107]:

U,x=Gaussian_Elimination(M)

Upper triangular

In [108]:

Out[108]:

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

In [109]:

Out[109]:

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

In [62]:

M1_aug

Out[62]:

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

In [63]:

Gaussian_Elimination(M1_aug)

Out[63]:

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

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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 0.003 & 59.14…

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

(E_1 + m E_2)\rightarrow (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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 0.003 & 59.14…

The exact system is instead

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 0.003 & 59.14…

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.

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.

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.

In [ ]:

#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

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.

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}| = \max_{i\leq j\leq n}|a_{ji}|

This way, propagated multiplication errors would be minimized.

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.

Matrix Inversion

Asumming a nonsingular matrix $A$ , if a matrix $A^{-1}$ exists, with $AA^{-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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 22: …} = AB= \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

\vdots

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 1 & 0 & \cdot…

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

Activity

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

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.

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 $\det A = 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:

\det A = \sum_{j=1}^n a_{ij}A_{ij}

\det A = \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.

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:

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 12: A = \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

\det A = \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.

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?

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 $\det A = 0$ .

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

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

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

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=\det A$ .

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

7. $\det A^t=\det A.$

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

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

\det A = \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!$ .

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.

Existence of inverse

A matrix $A$ has an inverse if $\det{A}\ne 0$ . See for example here

If the matrix $A$ has an inverse, then $\begin{align} A \boldsymbol{x}=&\boldsymbol{b}\\ A^{-1}A\boldsymbol{x}=&A^{-1}\boldsymbol{b}\\ \boldsymbol{x}=&A^{-1}\boldsymbol{b}\,, \end{align}$

Example From the previous example

In [74]:

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

Out[74]:

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 ]])

In [75]:

M[:,0]

Out[75]:

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

In [76]:

np.c_[M[:,0]]

Out[76]:

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

In [90]:

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

Out[90]:

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]])

In [78]:

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

Out[78]:

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

such that $\begin{align} A \boldsymbol{x}=&\boldsymbol{b}\,, \end{align}$

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

In [82]:

np.linalg.det(A)

Out[82]:

$\displaystyle -0.0967223339297434$

$\boldsymbol{x}=$

In [83]:

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

Out[83]:

$\displaystyle \left[\begin{matrix}0.274\\0.8722\\0.4114\\0.0093\end{matrix}\right]$

Numpy implementation

In [79]:

np.linalg.solve(A,b)

Out[79]:

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

In [91]:

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

Out[91]:

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

Matrix diagonalization

See linear-algebra-diagonalization

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.

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11} & a_{1…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 8: \left[ \̲m̲a̲t̲r̲i̲x̲{ 1 & 0 & \cdot…

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)$

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?

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)\rightarrow (E_j), \ \ \ \ \mbox{where}\ \ \ \ m_{j1} = \frac{a^{(1)}_{j1}}{a^{(1)}_{11}}

henceforth, $a^{(1)}_{ij}$ denotes the components of the original matrix $\mathbf{A}=\mathbf{A}^{(1)}$ , $a^{(2)}_{ij}$ the components of the matrix after eliminating the coefficients of $x_1$ , and generally, $a^{(k)}_{ij}$ 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)}$

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 27: …{(1)} = \left[ \̲m̲a̲t̲r̲i̲x̲{ 1 & 0 & \cdot…

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 33: …{ij} = \left\{ \̲m̲a̲t̲r̲i̲x̲{ 1 & \mbox{if}…

and

m_{ij} = \frac{a^{(j)}_{ij}}{a^{(j)}_{jj}}

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 27: …{(n)} = \left[ \̲m̲a̲t̲r̲i̲x̲{ a_{11}^{(n)} …

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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 76: …ij} = \left\{ \̲m̲a̲t̲r̲i̲x̲{ 1 & \mbox{if}…

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)}

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^{(1)}_{j1}/aa^{(1)}_{11})$ . 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

ParseError: KaTeX parse error: Undefined control sequence: \matrix at position 34: …ij} = \left\{ \̲m̲a̲t̲r̲i̲x̲{ 1 & \mbox{if}…

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)}$ .

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}$ .

In [86]:

M1

Out[86]:

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

In [87]:

import scipy

In [88]:

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

In [89]:

Out[89]:

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

In [ ]:

The same obtained before

In [ ]:

scipy.linalg.lu?

Eigenvalues and Eigenvectors activity

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

\begin{equation} \hat{H}|\Psi\rangle = E|\Psi\rangle \end{equation}

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

\begin{equation} \hat{\textbf{H}} = - \frac{\gamma \hbar}{2} \left( \begin{array}{cc} B_z & B_x -i B_y \\ B_x + i B_y & -B_z \end{array} \right) \end{equation}

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}\Psi$ - E $\Psi = 0$ gives the autofunctions $\Psi$ , i.e., the column vector $\Psi= \{\Psi_1, \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$ .

Find the allowed energy levels.
Find the autofunctions and normalize them.

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

In [ ]: