Solves the matrix equation Ax = b for x using QR decomposition and Gram-Schmidt orthogonalisation.

qr_decomp.py

"""
Solves the matrix equation Ax = b for x using QR decomposition and 
Gram-Schmidt orthogonalisation.

Author: spdskatr
"""
import random
import numpy as np

inner_product = np.dot

def qr_decompose(A : np.ndarray) -> (np.ndarray, np.ndarray):
  """
  Decomposes square matrix A into QR, where Q is an orthogonal matrix
  and R is an upper triangular matrix.
  """
  # U is Q's transpose (to make calculations easier)
  U = np.array(A.T, dtype=np.float64)

  # Make orthogonal basis via Gram-Schmidt Orthogonalisation
  for i in range(1, U.shape[0]):
    result = np.array(U[i], dtype=np.float64)
    for j in range(0, i):
      result -= (inner_product(U[j], U[i])/inner_product(U[j], U[j])) * U[j]
    U[i] = result

  # Make orthonormal basis
  for i in range(U.shape[0]):
    U[i] = U[i] / np.linalg.norm(U[i])
  
  # Re-transpose
  Q = U.T

  # Create upper triangular matrix
  R = np.ndarray((A.shape[0], A.shape[0]))
  for i in range(R.shape[0]):
    for j in range(i, R.shape[1]):
      R[i, j] = inner_product(U[i], A.T[j])
  
  return Q, R

def qr_solve(A : np.ndarray, b : np.ndarray) -> np.ndarray:
  """
  Finds solution x for Ax = b
  """
  # Find QR Decomposition
  Q, R = qr_decompose(A)

  # Remember that the inverse of Q is Q's transpose since Q is orthogonal
  # i.e. Q^-1 = Q^T
  QT_dot_b = Q.T.dot(b)

  # Rx = Q.T.dot(b)
  # Find x where R is upper triangular matrix
  # See https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution
  x = np.ndarray((Q.shape[0], 1))
  for i in range(Q.shape[0] - 1, -1, -1):
    s = np.array(QT_dot_b[i])
    for j in range(i+1, Q.shape[0]):
      s -= R[i, j] * x[j]
    x[i] = s / R[i, i]
  return x

#print(qr_solve(np.array([[2., 3.],
#                         [1., 8.]]),
#               np.array([[12.],
#                         [19.]])))

def build_test_matrix(x : np.ndarray) -> (np.ndarray, np.ndarray):
  x_shaped = x.reshape((-1, 1))
  m = x_shaped.shape[0]
  mat = np.random.random((m, m)) * 100
  while np.linalg.det(mat) == 0:
    mat = np.random.random((m, m)) * 100
  b = mat.dot(x_shaped)
  return mat, b

def do_solve_test():
  x = np.random.random((25, 1)) * 100
  A, b = build_test_matrix(x)
  print(x)
  print(A)
  print(b)
  x_sol = qr_solve(A, b)
  print("Solution: ")
  print(x_sol)

if __name__ == "__main__":
  do_solve_test()

You should implement modified Gram-Schmidt. Also, if you normalize after every single step, instead of at the end, you get slightly better performance (you can skip the denominator in (Uj.Ui / Ui.Ui)) and, if I'm not mistaken, increase numerical stability.