#!/usr/bin/env python
# coding: utf-8

# # QR algorithm

# In[71]:


import numpy as np


# ## Classical Gram-Schmidt algorithm

# In[72]:


def qr(A):
    n = A.shape[0]
    Q = np.empty_like(A)
    R = np.zeros((n,n))
    for j in range(n):
        Q[:,j] = A[:,j]
        for i in range(j):
            R[i,j] = Q[:,i].dot(A[:,j])
            Q[:,j] -= R[i,j]*Q[:,i]
        R[j,j] = np.linalg.norm(Q[:,j])
        Q[:,j] /= R[j,j]
    return Q,R


# We generate a random matrix and test this.

# In[73]:


n = 4
A = np.random.rand(n,n)
print('A =\n',A)
Q,R = qr(A)
print('Q =\n',Q)
print('R =\n',R)
print('A-QR =\n',A-Q.dot(R))


# ## Modified Gram-Schmidt

# In[74]:


def mqr(A):
    n = A.shape[0]
    Q = A.copy()
    R = np.zeros((n,n))
    for i in range(n):
        R[i,i] = np.linalg.norm(Q[:,i])
        Q[:,i] /= R[i,i]
        for j in range(i+1,n):
            R[i,j] = Q[:,i].dot(Q[:,j])
            Q[:,j] -= R[i,j]*Q[:,i]
    return Q,R


# In[75]:


Q,R = mqr(A)
print('Q =\n',Q)
print('R =\n',R)
print('A-QR =\n',A-Q.dot(R))