Orthogonal-Iteration.py

# -*- coding: utf-8 -*-
"""Lab4_Math170A_P3.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1BUC8PsCggw7UbKex-SKcCdFwtlyWPRFO

Problem 3
"""

import numpy as np

"""### Utility"""

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

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

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 myQR(A):
    m,n=A.shape
    R=np.copy(A)
    q=np.zeros((m,n))
    for i in range(n):
        tau=np.linalg.norm(R[i:,i])
        q[i:,i]=R[i:,i]
        q[i,i]+=sgn(R[i,i])*tau
        reflectMult(R[i:,i:],q[i:,i])
    return q, R[:n,:n]

"""### Orthogonal Iteration to Calculate EigenValues"""

def myOrth(A):
  n=len(A)
  Y=np.random.rand(n,n)
  i = 0
  while True:
    q, R=myQR(Y)
    Q=assembleQ(q)[:, :n]
    Y=A@Q # cols of Y converages 
    if np.allclose(np.zeros((n,n)),np.tril(np.transpose(Q)@A@Q,-1))==True:
      break
    i+=1
  return np.diag(np.transpose(Q)@A@Q) # n-dimensional vector of all eigenvalues

"""### Testing"""

# Personal Little Check 
P=np.random.rand(4,4)
d=np.array([345, 2, 3, 4])
A=P@np.diag(d)@np.linalg.inv(P)
eig=myOrth(A)
d[::-1].sort()
print(eig)
print(d)
np.allclose(eig,d)

for i in range(10):
  P=np.random.rand(10,10)
  d=np.random.rand(10)
  d[::-1].sort()
  D=np.diag(d)
  A=P@D@np.linalg.inv(P)
  eig=myOrth(A)
  print(np.allclose(eig,d))