<jupyter_start><jupyter_text>Week 3 MA544---**Objectives and Plan**1. Pseudo-inverse (Perron-Frobenius inverse) Properties1. Linear Systems of Equations and Gaussian Elimination with pivoting1. LU Decomposition of A1. QR Decomposition of a matrix1. Iterative solution of Linear Systems<jupyter_code>#IMPORT
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

## Set a seed for the random number generator
np.random.seed(100)<jupyter_output><empty_output><jupyter_text>Linear System of Equations---Consider the following system of $m$ linear equations in $n$ variables.\begin{align}a_{11} x_1 + a_{12} x_2  + \cdots + a_{1n} x_n  &= b_1 \\a_{21} x_1 + a_{22} x_2  + \cdots + a_{2n} x_n  &= b_2 \\ \vdots \qquad \qquad   & \ \\a_{m1} x_1 + a_{m2} x_2  + \cdots + a_{mn} x_n  &= b_ m,\end{align}-  The solution of a linear system represents the **point of intersection of hyperplanes**$$ x_1 \begin{pmatrix}a_{11}\\a_{21}\\ \vdots \\a_{m1}\end{pmatrix} + x_2 \begin{pmatrix}a_{12}\\a_{22}\\ \vdots \\a_{m2}\end{pmatrix} + \cdots + x_n \begin{pmatrix}a_{1n}\\a_{2n}\\ \vdots \\a_{mn}\end{pmatrix} = \begin{pmatrix}b_1\\b_2\\ \vdots \\b_m\end{pmatrix} $$  -  The solution also represents the **linear coding of the columns** of a matrix $A$ to get a vector $b$ in the column space ($\mathcal{R}(A)$). - The system could be represented in a compact form as $Ax = b$, where  $$ A=\begin{pmatrix}a_{11} & a_{12} & \cdots & a_{1n} \\a_{21} & a_{22} & \cdots & a_{2n} \\\vdots & \vdots & \ddots & \vdots \\a_{m1} & a_{m2} & \cdots & a_{mn}\end{pmatrix},\quad\mathbf{x}=\begin{pmatrix}x_1 \\x_2 \\\vdots \\x_n\end{pmatrix},\quad\mathbf{b}=\begin{pmatrix}b_1 \\b_2 \\\vdots \\b_m\end{pmatrix}$$ >- The system is called...>>- overdetermined when $m > n$.>>- underdetermined when $ m < n$.  Gaussian elimination with scaled-row partial pivoting--->- What is the need for pivoting?>- How can you find a factorization of the matrix by using the followign code.>- How can you modify the code to find the determinant of a matrix?>- How can you modify the following code to solve a linear system?>- How can you modify the code to find the inverse of a matrix?<jupyter_code>## Gaussian Elimination: Scaled Row Pivoting
## This function is based on the pseudo-code on page-148 in the Text by <NAME> Cheney
def GE_rsp(A):
    '''
    This function returns the P'LU factorization of a square matrix A
    by scaled row partial pivoting. 
    In place of returning L and U, elements of modified A are used to hold values of L and U.
    '''
    m,n = A.shape
    
    L = np.eye(n) # Not being used
    U = np.zeros_like(A) # Not being used
    if m !=n:
        sys.exit("This function needs a square matrix as an input.")
        
    # The initial ordering of rows
    p = list(range(n))
    
    # Scaling vector: absolute maximum elements of each row
    s = np.max(np.abs(A), axis=1)
    
    print("Scaling Vector: ",s)
    
    # Start the k-1 passes of Guassian Elimination on A
    for k in range(n-1):
        
        print("\n PASS {}: \n".format(k+1), A)
        # Find the pivot element and interchange the rows
        pivot_index = k + np.argmax(np.abs(A[p[k:], k])/s[p[k:]])        
        
        # Interchange element in the permutation vector
        if pivot_index !=k:
            temp = p[k]
            p[k]=p[pivot_index]
            p[pivot_index] = temp
            print("permutation vector: ",p)
            
        print("\n Pivot Element: {0:.2f} \n".format(A[p[k],k]))
        if np.abs(A[p[k],k]) < 10**(-20):
             sys.exit("ERROR!! Provided matrix is singular.")
        
        # For the k-th pivot row Perform the Gaussian elimination on the following rows
        for i in range(k+1, n):
            # Find the multiplier
            z = A[p[i],k]/A[p[k],k]
            
            #Save z in A itself. You can save this in L also
            A[p[i],k] = z
            
            #Elimination operation: Changes all elements in a row simultaneously
            ##
            A[p[i],k+1:] -= z*A[p[k],k+1:]
    return A, p
## Example on page number 146 (Kincaid Cheney).
## Example solved manually in class
A = np.array([[2, 3, -6], [1,-6,8], [3, -2, 1]], dtype=float)
print("\n Given A: \n ",A)
A,p =GE_rsp(A)
print("\n After Gaussian Elimination with RSPP: \n", A)
print("\n The permutation Vector is: \n", p)<jupyter_output>Given A: 
  [[ 2.  3. -6.]
 [ 1. -6.  8.]
 [ 3. -2.  1.]]
Scaling Vector:  [6. 8. 3.]

 PASS 1: 
 [[ 2.  3. -6.]
 [ 1. -6.  8.]
 [ 3. -2.  1.]]
permutation vector:  [2, 1, 0]

 Pivot Element: 3.00 


 PASS 2: 
 [[ 0.66666667  4.33333333 -6.66666667]
 [ 0.33333333 -5.33333333  7.66666667]
 [ 3.         -2.          1.        ]]
permutation vector:  [2, 0, 1]

 Pivot Element: 4.33 


 After Gaussian Elimination with RSPP: 
 [[ 0.66666667  4.33333333 -6.66666667]
 [ 0.33333333 -1.23076923 -0.53846154]
 [ 3.         -2.          1.        ]]

 The permutation Vector is: 
 [2, 0, 1]<jupyter_text>LU Factorization ---One can use the Gaussian elimination to decompose the matrix A as follows$$P A = L U,\text{ or }\ A = P^T L U;$$where $P$ is a permutation matrix, $L$ is a unit lower-triangular matrix and $U$ is an upper triangular matrix.<jupyter_code>print("\n Upper triangular, U:\n ", np.triu(A[p,:]))
print("\n Lower triangular, L:\n", np.tril(A[p,:], -1)+np.eye(3))
print("The Permutation Matrix, P:\n",np.eye(3)[p,:])<jupyter_output>Upper triangular, U:
  [[ 3.         -2.          1.        ]
 [ 0.          4.33333333 -6.66666667]
 [ 0.          0.         -0.53846154]]

 Lower triangular, L:
 [[ 1.          0.          0.        ]
 [ 0.66666667  1.          0.        ]
 [ 0.33333333 -1.23076923  1.        ]]
The Permutation Matrix, P:
 [[ 0.  0.  1.]
 [ 1.  0.  0.]
 [ 0.  1.  0.]]<jupyter_text>Iterative Methods for Linear Systems--- Jacobi Method<jupyter_code># You can modify this code to answer the following
'''
Jacobi's iteration method for solving the system of equations Ax=b.
p0 is the initialization for the iteration.
'''
def jacobi(A, b, p0, tol, maxIter=100):
    n=len(A)
    p = p0

    for k in range(maxIter):
        p_old = p.copy() # In python assignment is not the same as copy
        
        # Update every component of iterant p
        for i in range(n):
            sumi = b[i];
            for j in range(n):
                if i==j: # Diagonal elements are not included in Jacobi
                    continue;
                sumi = sumi - A[i,j] * p_old[j]
            p[i] = sumi/A[i,i]
                
        rel_error = np.linalg.norm(p-p_old)/n
        # print("Relative error in iteration", k+1,":",rel_error)
        if rel_error<tol:
            print("TOLERANCE MET BEFORE MAX-ITERATION")
            break
    return p;
# Example System
A = np.array([[10, -1, 2, 0],
              [-1, 11, -1, 3],
              [2, -1, 10, -1],
              [0, 3, -1, 8]],dtype=float)
b = np.array([6, 25, -11, 15],dtype=float)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14504783342712438, 'fraction_numerical': 0.049943725379853685, 'mean_word_length': 2.380408939610081, 'pattern_counts': {'":': 1, '<': 19, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 14, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '657309', 'license': 'MIT', 'n_tokens_mistral': 2735, 'n_tokens_neox': 2458, 'n_words': 818, 'path': 'course_notes/MA544 Share/NB3 MA544.ipynb', 'repo_name': 'jschmidtnj/ma544-final-project', 'size': '7113'}