from numpy import *

A = floor(random.rand(4, 4) * 20 - 10)  # random matrix

Q, R = linalg.qr(A)  # qr decomposition of A

b = floor(random.rand(4, 1) * 20 - 10)
# solve Ax = b using the standard numpy function
x = linalg.solve(A, b)

# solve Ax = b using Q and R
y = dot(Q.T, b)
xQR = linalg.solve(R, y)

print("\nSolution compared")
print(x.T, "Ax=b")
print(xQR.T, "Rx=y")