[](https://colab.research.google.com/drive/1TnyJDxbMlMu01WZiBKEBAH-CpSBkfPMi?usp=sharing)


# ECHELON FORM AND ROW REDUCED FORM

12/02/2022

Leonardo Saccotelli

License: Creative Commons non-commercial not derivative works - 2.5 Italy License http://creativecommons.org/licenses/by-nc-nd/2.5/it/

##  Row Reduced Echelon Form from Sympy library

We use the symbolic package 'sympy' that compute the row reduce echelon form with different type of matrices.

### Import of library


```python
import sympy as sym
```

### Definition of matrices and computation of the row reduce echelon form

Definition of the matrix A (with 3 rows and 7 columns) and transformation of it in the row reduce echelon form.


```python
# Definition of the matrix A
A = sym.Matrix([[1, 2, 2, 3, 1, 0, 0], 
                [2, 4, 1, 3, 0,1,0], 
                [3, 6, 1, 4,0,0,1]])  
print("Matrix A: {} ".format(A)) 

# Use sympy.rref() method to compute the row reduce echelon form. 
# The method rref return a matrix in the row reduce echelon form and
# a vector with the index of pivots/ basic columns in the matrix.
A_rref = A.rref()   
      
# print the row reduce echelon form of the matrix A and the indexes of the basic columns
print("\nThe Row echelon form of matrix A and the pivot columns are: {}"
      .format(A_rref))
```

    Matrix A: Matrix([[1, 2, 2, 3, 1, 0, 0], [2, 4, 1, 3, 0, 1, 0], [3, 6, 1, 4, 0, 0, 1]]) 
    
    The Row echelon form of matrix A and the pivot columns are: (Matrix([
    [1, 2, 0, 1, 0, -1,  1],
    [0, 0, 1, 1, 0,  3, -2],
    [0, 0, 0, 0, 1, -5,  3]]), (0, 2, 4))


Definition of the matrix B (with 4 rows and 6 columns) and transformation of it in the row reduce echelon form.


```python
# Definition of the matrix B
B = sym.Matrix([[1, 1, 2, 2, 1, 1], [2, 2, 4, 4, 3, 1], 
                [2, 2, 4, 4, 2, 2], [3, 5, 8, 6, 5, 3]]) 
print("Matrix B: {} ".format(B))

B_rref = B.rref()

# print the row reduce echelon form of the matrix B 
# and the indexes of the basic columns
print("\nThe Row echelon form of matrix B and the pivot columns are: {}"
      .format(B_rref))
```

    Matrix B: Matrix([[1, 1, 2, 2, 1, 1], [2, 2, 4, 4, 3, 1], [2, 2, 4, 4, 2, 2], [3, 5, 8, 6, 5, 3]]) 
    
    The Row echelon form of matrix B and the pivot columns are: (Matrix([
    [1, 0, 1, 2, 0,  1],
    [0, 1, 1, 0, 0,  1],
    [0, 0, 0, 0, 1, -1],
    [0, 0, 0, 0, 0,  0]]), (0, 1, 4))


## Implementation of the Echelon Form Algorithm with partial pivoting

We define a function which takes in Input a Matrix of size m x n and that gives in output the echelon form of it. The permutation is performed in order to have a stable algorithm using the technique of partial pivoting.

### Import of library


```python
import numpy as np
```

### Function to compute the Echelon form with partial pivoting 


```python
def compute_echelon_form_with_partial_pivot(A): 
    # The number of column is greater or equal the number of rows
    # Copying the matrix A in U 
    U = np.copy(A) 
    
    # Retrive the dimension of the matrix A, n_rows, n_cols
    (n_rows, n_cols)=A.shape
    
    # Index related to the column
    j = 0
    for i in range(0,n_rows): 
        ech=1
        while (ech == 1) & (j < n_cols):
            # find the pivotal index, maximum element in the column j
            # and return the index of the row in 
            # which we have the pivotal element
            indm=np.argmax(abs(U[i:n_rows,j])) 
            indm=indm+i
            # perform the permutation to work well, 
            # we should change values != 0 
            # with an absolute value less then a constant very small
            if (indm != i) & (abs(U[indm,j]) > 1e-15):
                # perform permutation between row i and row in 
                # which we have the pivotal element
                U[ [i, indm],:]=U[[indm,i],:]
                
            # If the pivotal element is different from zero
            # we can compute the M value, which is the constant 
            # used for the elimination
            if (abs(U[i,j]) > 1e-15):
                
                # we put in vector all the value to compute the elimination
                # of all the row below the pivotal element
                M = U[i+1:n_rows,j]/U[i,j] 
                # in numpy array there is no 
                # difference row vectors or column vectors
                # np.outer to perform the outer product between two vectors
                # We perform all the row linear combination with one statement
                U[i+1:n_rows,j+1:n_cols]=U[i+1:n_rows,j+1:n_cols]-np.outer(M,U[i,j+1:n_cols]) 
                
                # We put 0 below the pivotal element
                U[i+1:n_rows,j]=0
                
                # Move on the next column
                j=j+1
                ech=0
            else:
                j=j+1
                ech=1      
    return(U)    


```

### Test with a matrix of size 3x4


```python
C = np.array([[1., 2., 2., 3.], [2, 4, 1, 3], [3, 6, 1, 4]]) 
print("Matrix C: \n{}".format(C))

U=compute_echelon_form_with_partial_pivot(C)
print("\nThe Row echelon form of matrix C: \n{}".format(U))
```

    Matrix C: 
    [[1. 2. 2. 3.]
     [2. 4. 1. 3.]
     [3. 6. 1. 4.]]
    
    The Row echelon form of matrix C: 
    [[3.00000000e+00 6.00000000e+00 1.00000000e+00 4.00000000e+00]
     [0.00000000e+00 0.00000000e+00 1.66666667e+00 1.66666667e+00]
     [0.00000000e+00 0.00000000e+00 0.00000000e+00 1.11022302e-16]]


## Implementation of the Row Reduced Echelon Form Algorithm with partial pivoting

We define a function which takes in Input a Matrix of size m x n and that gives in output the row reduced echelon form echelon form of it and the related permutation matrix. The permutation is performed in order to have a stable algorithm using the technique of partial pivoting.

### Function to compute the Row Reduced Echelon Form with partial pivoting 


```python
def compute_row_reduced_form_with_partial_pivot(A):
    # Copy the original matrix
    U = np.copy(A)
    
    # Find number of rows and columns
    (m,n)=A.shape
    
    # Create an identity matrix of size mxm
    # to use as permutation matrix
    P = np.eye(m)
    
    # Set the first column of the matrix
    j = 0
    for i in range(0,m): 
        ech=1
        while (ech == 1) & (j < n):
            # find the pivotal index, maximum element in the column j
            # and return the index of the row 
            # in which we have the pivotal element
            indm=np.argmax(abs(U[i:m,j]))
            indm=indm+i
            
            # perform the permutation to work well,
            # we should change values != 0 
            # with an absolute value less then a constant very small
            if (indm != i) & (abs(U[indm,j]) > 1e-15):
                
                # perform permutation between row i and row in 
                # which we have the pivotal element
                U[ [i, indm],:]=U[[indm,i],:] 
                
                # Update the permutation matrix to save all the 
                # permutation that we will perform
                P[ [i, indm],:]=P[[indm,i],:] 
                
            # If the pivotal element is different from zero
            # we can compute the M value, which is the constant 
            # used for the elimination
            if ( abs(U[i,j]) > 1e-15):
                
                # we put in vector all the value to compute the elimination
                # of all the row below the pivotal element
                M=U[i+1:m,j]/U[i,j]
                
                # in numpy array there is no difference 
                # row vectors or column vectors
                # np.outer to perform the outer product between two vectors
                # We perform all the row linear combination with one statement
                U[i+1:m,j+1:n]=U[i+1:m,j+1:n]-np.outer(M,U[i,j+1:n])
                
                # We put 0 below the pivotal element
                U[i+1:m,j]=0 
                
                # All the pivotal element must be equal to 1
                # so we compute the division of all the pivotal
                # element at the same time
                U[i,j:n] = U[i,j:n]/U[i,j]
                
                # If the actual row is not the first row, we have to
                # compute linear combination of row in order to put 
                # 0 above the pivotal element
                if i>0:
                    # We compute the constant to be used in the 
                    # elimination on the element above the pivotal 
                    # element
                    M= U[0:i,j]/U[i,j]
                    
                    # Computing the elimination on the upper part
                    U[0:(i),j:n]=U[0:i,j:n]-np.outer(M,U[i,j:n])     
                j=j+1
                ech=0
            else:
                j=j+1
                ech=1
            
    return(U,P)
```

### Test with a matrix of size 3x7


```python
D = np.array([[1., 2., 2., 3., 1., 0, 0], [2, 4, 1, 3, 0,1,0], [3, 6, 1, 4,0,0,1]]) 
print("Matrix D: \n{}".format(D))

(U,P)=compute_row_reduced_form_with_partial_pivot(D)
print("\nThe Row Reduced Echelon form of matrix D: \n{}".format(U))
print("\nThe Permutation matrix used to perform the partial pivoting on D: \n{}".format(P))
```

    Matrix D: 
    [[1. 2. 2. 3. 1. 0. 0.]
     [2. 4. 1. 3. 0. 1. 0.]
     [3. 6. 1. 4. 0. 0. 1.]]
    
    The Row Reduced Echelon form of matrix D: 
    [[ 1.00000000e+00  2.00000000e+00  0.00000000e+00  1.00000000e+00
       0.00000000e+00 -1.00000000e+00  1.00000000e+00]
     [ 0.00000000e+00  0.00000000e+00  1.00000000e+00  1.00000000e+00
       0.00000000e+00  3.00000000e+00 -2.00000000e+00]
     [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.11022302e-16
       1.00000000e+00 -5.00000000e+00  3.00000000e+00]]
    
    The Permutation matrix used to perform the partial pivoting on D: 
    [[0. 0. 1.]
     [1. 0. 0.]
     [0. 1. 0.]]


## Comparing the Sympy function sys.rref() with our function to compute the row reduced form

Exercise:  use the function sp.rref and the function compute_row_reduced_form_with_partial_pivot to compute the row reduced form of the matrix in the exercise 4.2.1  on page 178 of the book of Carl Meyer and find the  spanning sets for all the four fundamental subspaces


```python
E = sym.Matrix([[1, 2, 1, 1, 5],[-2, -4, 0, 4, -2],[1, 2, 2, 4, 9]])
print("Matrix E: {}".format(E))

E_rref = E.rref()
print("\n1) Using sym.rref function")
print("----------------------------------------------------------")
print("The Row Reduced form of matrix E: \n{}".format(E_rref))
print("\nThe spanning sets of all the four fundamental subspaces: ")

(U,P)=compute_row_reduced_form_with_partial_pivot(E)
print("\n2) Using compute_row_reduced_form_with_partial_pivot function")
print("----------------------------------------------------------")
print("The Row Reduced form of matrix E: \n{}".format(U))
print("\nThe spanning sets of all the four fundamental subspaces: ")
```

    Matrix E: Matrix([[1, 2, 1, 1, 5], [-2, -4, 0, 4, -2], [1, 2, 2, 4, 9]])
    
    1) Using sym.rref function
    ----------------------------------------------------------
    The Row Reduced form of matrix E: 
    (Matrix([
    [1, 2, 0, -2, 1],
    [0, 0, 1,  3, 4],
    [0, 0, 0,  0, 0]]), (0, 2))
    
    The spanning sets of all the four fundamental subspaces: 
    
    2) Using compute_row_reduced_form_with_partial_pivot function
    ----------------------------------------------------------
    The Row Reduced form of matrix E: 
    [[1 2 0 -2 1]
     [0 0 1 3 4]
     [0 0 0 0 0]]
    
    The spanning sets of all the four fundamental subspaces: 


## Block Matrix Definition and Multiplication

Definition of a block matrix using the numpy function. First, we define single matrix and define we create a block matrix merging them. After that, we will perform the row reduced form on this matrix and at the end we will multiply the result (the inverse of the original matrix) with the original matrix in order to check if we are able to find the identity. The experiment show us that we are not able to find again the perfect identity.


```python
Random_Matrix = np.random.rand(3,3)
Identity_Matrix = np.eye(3)

# Defining a matrix using its blocks
Block_Matrix = np.block([Random_Matrix,Identity_Matrix]) 

print("Random Matrix (first block):\n{}".format(Random_Matrix))
print("\nIdentity Matrix (second block):\n{}".format(Identity_Matrix))
print("\nBlock Matrix: \n{}".format(Block_Matrix))
```

    Random Matrix (first block):
    [[0.32854889 0.98950404 0.39290949]
     [0.25902266 0.13538146 0.37090808]
     [0.37575205 0.96043677 0.77677363]]
    
    Identity Matrix (second block):
    [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 1.]]
    
    Block Matrix: 
    [[0.32854889 0.98950404 0.39290949 1.         0.         0.        ]
     [0.25902266 0.13538146 0.37090808 0.         1.         0.        ]
     [0.37575205 0.96043677 0.77677363 0.         0.         1.        ]]



```python
m = 3
n = 3
(Row_Reduced_Block_Matrix,P)=compute_row_reduced_form_with_partial_pivot(Block_Matrix)

print("First Block of the Row reduced form of the block matrix: \n{}"
      .format(Row_Reduced_Block_Matrix[0:m,0:n]))

print("\nSecond Block of the Row reduced form of the block matrix: \n{}"
      .format(Row_Reduced_Block_Matrix[0:m,n:m+n]))

X=Row_Reduced_Block_Matrix[0:m,n:m+n]
```

    First Block of the Row reduced form of the block matrix: 
    [[1. 0. 0.]
     [0. 1. 0.]
     [0. 0. 1.]]
    
    Second Block of the Row reduced form of the block matrix: 
    [[ 3.80907055  5.935809   -4.76105175]
     [ 0.93807126 -1.63198622  0.30477318]
     [-3.00245028 -0.85349552  3.21362554]]



```python
np.dot(Random_Matrix,X)
```




    array([[ 1.00000000e+00, -1.73255491e-16,  2.13847669e-16],
           [-1.11823876e-16,  1.00000000e+00,  5.36994453e-18],
           [-2.45580527e-16, -8.08890910e-17,  1.00000000e+00]])




```python
np.dot(X,Random_Matrix)
```




    array([[ 1.00000000e+00, -3.50121445e-16,  2.74577539e-16],
           [ 1.19842123e-17,  1.00000000e+00,  6.57938721e-17],
           [-8.31437201e-17, -2.81521386e-16,  1.00000000e+00]])



1ccf3eb6e3f26b0b4706705a5ca262911ba1072d
ipynb
Jupyter Notebook
Numerical-Methods-Exercises/1_Echelon-Row-Reduced-Form/echelon_form_and_row_reduced_form.ipynb
LeonardoSaccotelli/Numerical-Methods-for-Computer-Science
00802e4397545ea65b4cf946390dbe198d8eba34
MIT
Numerical-Methods-Exercises/1_Echelon-Row-Reduced-Form/echelon_form_and_row_reduced_form.ipynb
LeonardoSaccotelli/Numerical-Methods-for-Computer-Science
00802e4397545ea65b4cf946390dbe198d8eba34
MIT
Numerical-Methods-Exercises/1_Echelon-Row-Reduced-Form/echelon_form_and_row_reduced_form.ipynb
LeonardoSaccotelli/Numerical-Methods-for-Computer-Science
00802e4397545ea65b4cf946390dbe198d8eba34
MIT
Qwen/Qwen-72B
1. YES
2. YES