qr_decomposition.py

# -*- coding: utf-8 -*-
"""QR-Algorithm.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1Og_eOm5rgWFhTmlR2Cp5qGl-t2TqXpYm
"""

import numpy as np

"""### Utility"""

def assembleQ(q):
  m,n=q.shape
  Q=np.eye(m)
  for i in range(n):
    reflectMult(Q,q[:,n-1-i])
  return Q

def sgn(x):
  if x==0:
    return 1
  else:
    return np.sign(x)

"""### Householder Reflection"""

def reflectMult(A, w):
  u = w/np.linalg.norm(w) # normalized vector
  Q = np.eye(len(w), dtype=float) - 2*np.outer(u, u)
  A[:,:] = Q@A

"""### QR Algorithm by HouseHolder Reflection"""

def myQR(A):
    nRows, nCols = A.shape
    R = np.copy(A)
    Q = np.eye(nRows)
    q = np.zeros((nRows, nCols))
    for i in range(nCols - (nRows ==nCols)):
        Qi = np.eye(nRows) # householder transform matrix at step i
        v = R[i:, i]
        # construct e
        e = np.zeros(nRows-i)
        e[0] = 1
        # constuct the reflection vecotr w (unnormalized)
        w = -sgn(v[0])*e*np.linalg.norm(v) - v
        q[i:, i] = w
        # construct the reflection matrix
        Qi[i:, i:] = np.eye(nRows-i) - 2*np.outer(w, w)/np.dot(w, w)
        # apply changes onto the input matrix to get R = Q...A
        R = Qi@R # update the transformed matrix
    return q, R[:nCols, :nCols]

"""### Testing"""

# QR Algorithm
for i in range(10):
  A=np.random.rand(100,50)
  q,R=myQR(A)
  Q=assembleQ(q)
  print(np.allclose(Q[:,:50]@R,A))