QR Analysis

output.txt

Matrix Rank: 9
QR is Solution? True
These a solution in the span of Q_2.
Diagonal of R (9); Rank of A (9)
Is non-zero diag of R equal to rank? (True)
Is x_1 a solution? True
Is x_1 == x_0? True
Is our new solution x_0? True
Is our new solution x_0? True
Is our new solution x_0? True
Is our new solution x_0? True

qr_analysis.ipynb

qr_analysis.py

#Setup
import numpy as np
import scipy,random
import matplotlib.pyplot as plt
np.set_printoptions(precision=2)

#n>m (over-determined)
m,n = (10,15)

#Make exmaple ternary matrix
A = (2*np.random.rand(m,n)-1).round()

#Reduce Rank
for i in range(n):
    ids = set(range(n))-set([i]) #Remove current column
    n_cols = 2#Number of columns in convex combination
    col_ids = random.choices(list(ids),k=n_cols) #List of random col ids
    cols = A[:,col_ids] #Get a list of columns
    weights = np.random.rand(n_cols)
    weights /= sum(weights) #Calculate convex combinations of weights
    convex_combination = np.add.reduce([w*v for w,v in zip(weights.T,cols.T)])
    A[:,i] = convex_combination

print("Matrix Rank:",np.linalg.matrix_rank(A))


#Generate solution
x = (2*np.random.rand( n ) - 1).round() #Ternary soultion
y = A@x

# QR pivoting
rank = np.linalg.matrix_rank(A)
Q,R,P= scipy.linalg.qr(A,pivoting=True)
rank = r = sum(np.around(np.diag(R),12)!=0)
P = np.eye(len(P))[:,P]
Q_1 = Q[:,:rank]
Q_2 = Q[:,:-rank]
R_1 = R[:rank,:rank]
R_2 = R[:,:-(rank-m)]
# Get a single solution
z_1 = scipy.linalg.solve(R_1,Q_1.T@y)
padding = np.zeros((n-r)) #zero padding
x_0 = P@np.hstack([z_1,padding]) #One solution
print("QR is Solution?",np.allclose(y,A@x_0)) 
#>>> QR is Solution? True

if np.allclose(Q_2.T@y,0):
    print("There exists a solution to this system")
else:
    print("These a solution in the span of Q_2.")

print(f"Diagonal of R ({r}); Rank of A ({np.linalg.matrix_rank(A)})")
print(f"Is non-zero diag of R equal to rank? ({r==np.linalg.matrix_rank(A)})")
#>>> Is non-zero diag of R equal to rank? (True)


# Trying to get soultion set
z_2 = scipy.linalg.solve(R_1,Q_1.T@R_2)
zeros = np.zeros((n-r,m-r)) #padding
L = -P@np.vstack([z_2,zeros])

x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()
while not np.allclose(y,A@x_1):
    x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()
print("Is x_1 a solution?", np.allclose(y,A@x_1))
print("Is x_1 == x_0?", np.allclose(x_0,x_1))

#Show other solutions
for i,col in enumerate(Q_2):
    x_n = x_0 + L@col
    if np.allclose(y,A@x_n):
        print("Is our new solution x_0?",np.allclose(x_n,x_0))