# A = QR avec Q qui est orthogonale et R triangulaire supérieure
def decompositionQR(A):
  m, n = A.shape
  Q = np.zeros((m,n))
  R = np.zeros((n,n))

  Q[:,0] = A[:,0]/np.linalg.norm(A[:,0])

  for k in range(1,n):
    Q[:,k] = A[:,k]
    #je dois orthogonaliser cette colonne Q[:,k]
    #par rapoort à toutes les précédentes
    for i in range(k):
      projection = (Q[:,k]@Q[:,i])*Q[:,i] #projection sur le vecteur i
      Q[:,k] = Q[:,k] - projection
    Q[:,k] = Q[:,k]/np.linalg.norm(Q[:,k])

  R = Q.T@A

  return Q, R

from scipy.linalg import hilbert
import matplotlib.pyplot as plt

erreur = np.zeros(200)
for n in range(2,200):
  A = hilbert(n)
  Q, R = decompositionQR(A)
  erreur[n] = np.linalg.norm(Q.T@Q-np.eye(n))

plt.plot(erreur)

def algoQR(A,maxiter=10**3):
  n = A.shape[0]
  B = A.copy()
  Q = np.eye(n) #on initialise la matrice Q à l'identité
  for _ in range(maxiter):
    Qk, R = decompositionQR(B)
    B     = R@Qk #-> converger vers une matrice triangulaire où les valeurs propres seront sur la diagonale
    Q     = Q@Qk #-> converge vers une matrice contenant tous les vecteurs propres
  return B, Q


A = np.array([[3,1,-1],
              [1,3,-1],
              [-1,-1,5]])
B, Q = algoQR(A,maxiter=500)
print('B=')
print(B)
print('Q=')
print(Q)