QR Solver

qr_solver.py

import scipy
import numpy as np

#Ref: https://inst.eecs.berkeley.edu/~ee127/sp21/livebook/l_lineqs_solving.html
def qr_solve(A,y,silent=False):
    A = np.array(A)
    y = np.array(y)
    m,n = A.shape
    
    Q,R,P= scipy.linalg.qr(A,pivoting=True)
    rank = r = np.linalg.matrix_rank(A)
    P = np.eye(len(P))[:,P]
    
    Q_1 = Q[:,:rank]
    R_1 = R[:rank,:rank]
    R_2 = R[:,:-(rank-m)]
    
    z_1 = np.linalg.inv(R_1)@Q_1.T@y
    padding = np.zeros(n-r) #zero padding
    x_0 = P@np.hstack([z_1,padding]) #Add padding
    return x_0

qr_solver_documented.py

import scipy
import numpy as np
import logging,unittest

#Ref: https://inst.eecs.berkeley.edu/~ee127/sp21/livebook/l_lineqs_solving.html
def qr_solve(A,y,silent=False):
    A = np.array(A)
    y = np.array(y)
    
    #Solving y=Ax for an x_0
    #A has m rows and n columns -> y has m rows and x has n rows
    m,n = A.shape
    
    #QR decomposition
    Q,R,P= scipy.linalg.qr(A,pivoting=True)
    #Q is an m by m orthogonal matrix
    #R is an m by n upper right triangle matrix
    
    #P is a permuntation matrix with n rows and n columns
    #P can order A by rank for *rank revealing*
    P = np.eye(len(P))[:,P]
    
    #Let r be the number of linearly independent columns & rows in A (the rank)
    rank = r = np.linalg.matrix_rank(A)
    
    #Q is a m by m square orthogonal matrix
    #Q_1 has m rows and r columns
    Q_1 = Q[:,:rank]
    
    #R_1 is an r by r upper right triangular matrix
    R_1 = R[:rank,:rank]
    
    #R_2 is an r by m-r matrix
    R_2 = R[:,:-(rank-m)]
    
    #z_1 is a column vector with r elements
    z_1 = np.linalg.inv(R_1)@Q_1.T@y
    
    #We pad z_1 with r-m zeros at the bottom for a solution vector
    padding = np.zeros(n-r) #zero padding
    x_0 = P@np.hstack([z_1,padding]) #Add padding
    
    if not silent:
        #Log if there was a solution
        is_solution = np.allclose(y,A@x_0)
        logging.warning("Solution Found!") if is_solution else print("No Solution!")
    return x_0


class TestQRSolver(unittest.TestCase):
    def test_is_no_soultion(self):
        m,n = (100,150)
        A = np.random.rand( m, n )
        y = np.random.rand( m )
        x_0 = qr_solve(A, y, silent=True)
        return np.allclose(y, A@x_0)
    
    def test_is_solution(self):
        m,n = (100,150)
        A = np.random.rand( m, n )
        x = np.random.rand( n )
        y = A@x
        x_0 = qr_solve(A, y, silent=True)
        return np.allclose(y, A@x_0)
if __name__=="__main__":
    unittest.main()