Numerical Linear Algebra: QR, Hessenberg, and Schur

import numpy as np
np.set_printoptions( precision=3, linewidth = 132, suppress = True )

# Set dimensions and create identity matrices
m,n = (8,6)
I_m = np.identity(m)
I_n = np.identity(n)

# Seed the RNG.
rng = np.random.default_rng(111)

# A: a "tall" m x n matrix
data = rng.negative_binomial(10, 0.25, size=m*n)    
A0 = np.array(data, dtype=np.float64).reshape((m,n))

# B: a square n x n matrix
data = rng.negative_binomial(10, 0.25, size=n*n)
B0 = np.array(data, dtype=np.float64).reshape((n,n))

def house(x):
  """
  Golub and Van Loan, 3rd Ed.  
  Algorithm 5.1.1
  """
  n = len(x)
  sigma = x[1:] @ x[1:]
  v = np.concatenate([[1], x[1:]])

  if sigma == 0:
    beta = 0
  else:
    mu = np.sqrt(x[0]**2 + sigma)
    if x[0] <= 0:
      v[0] = x[0] - mu
    else:
      v[0] = -sigma / (x[0] + mu)
    beta = 2 * v[0]**2 / (sigma + v[0]**2)
    v = v/v[0]

  return(v, beta)

# A householder matrix to zero out elements of the first column of A

# Let `x0` be the first column of A
x0 = A0[:,0].copy()

# Compute the householder matrix
v, beta = house(x0)
P0 = (I_m - beta * np.outer(v,v))

# Its effect: zeroes!
P0 @ x0

P1 = np.identity(m)

# Let `x1` be the next column of A. 
x1 = A0[:,1].copy()

# A Householder matrix to introduce zeros below the diagonal in the second column.  At the
# second iteration, we'll be working with P0 @ x:
v1, b1 = house((P0 @ x1)[1:m])
P1[1:m, 1:m] = np.identity(m-1) - b1 * np.outer(v1, v1)

# P1 @ P0 @ A0
P1 @ P0 @ A0

"""
Golub and Van Loan, 3rd Ed.  
Algorithm 5.2.1: Householder QR
"""

# The "R" in a "QR" factorization
A = A0.copy()
R = A
V = np.zeros_like(A)
beta = np.zeros(n)

for j in range(n):
  x = A[j:m,j]
  v,beta[j] = house(x) 

  R[j:m, j:n] = R[j:m, j:n] - beta[j] * np.outer(v, v @ A[j:m, j:n])

  if j < m:
    R[(j+1):m,j] = 0
    V[(j+1):m,j] = v[1:(m-j+1)]

# The "Q" in "QR" factorization
Q = np.identity(m)
for j in reversed(range(n)):
  v = V[j:m,j].copy()
  v[0] = 1
  Q[j:m,j:m] = Q[j:m,j:m] - beta[j] * np.outer(v, v @ Q[j:m,j:m])

# Q is unitary
np.allclose(I_m, Q.T @ Q)

# R is upper diagonal
R

# Q @ R = A
np.allclose(Q @ R, A0)

"""
Golub and Van Loan, 3rd Ed.  
Algorithm 7.4.2: Householder Reduction to Hessenberg Form
"""

H = B0.copy()

V = np.zeros_like(H)
beta = np.zeros(n-2)

for k in range(n-2):
  x = H[(k+1):n,k].copy()
  v, beta[k] = house(x)

  V[(k+1):n,k] = v

  H[(k+1):n,k:n] = H[(k+1):n,k:n] - beta[k] * np.outer(v, v @ H[(k+1):n,k:n])
  H[0:n,(k+1):n] = H[0:n,(k+1):n] - beta[k] * np.outer( H[0:n,(k+1):n] @ v, v)

# The "U" in U' @ B0 @ U = H
U = np.identity(n)
for j in reversed(range(n-2)):
  v = V[j:n,j]
  U[j:n,j:n] = U[j:n,j:n] - beta[j] * np.outer(v, v @ U[j:n,j:n])

H

# U' @ B0 @ U = H 
np.allclose(U.T @ B0 @ U, H)

# U' @ U = I
np.allclose(U.T @ U, np.identity(n))

def givens(a,b):
  """
  Golub and Van Loan, 3rd Ed.  
  Algorithm 5.1.3
  """
  if b == 0:
    c = 1
    s = 0
  else:
    if np.abs(b) > np.abs(a):
      tau = -a/b
      s = 1 / np.sqrt(1 + tau**2)
      c = s * tau
    else:
      tau = -b/a
      c = 1/np.sqrt(1 + tau**2)
      s = c*tau

  return (c,s)

"""
Golub and Van Loan, 3rd Ed.  
Algorithm 7.4.1
"""

def qrstep(H):
  G = np.zeros((n-1,2,2))
  for k in range(n-1):
    c,s = givens(H[k,k], H[k+1,k])
    G[k] = np.array([[c,s], [-s,c]])
    H[k:(k+2), k:n] = G[k].T @ H[k:(k+2),k:n]
  
  for k in range(n-1):
    H[0:(k+2), k:(k+2)] = H[0:(k+2), k:(k+2)] @ G[k]

  return H

# Set the starting point and apply one step
H = U.T @ B0 @ U
qrstep(H)

# The goal after one step: check
H = U.T @ B0 @ U
Q, R = np.linalg.qr(H)
R @ Q

# Apply 400 steps to obtain a Schur decomposition.
S = U.T @ B0 @ U
for i in range(400):
  S = qrstep(S)
S

# E.g., the eigenvalues from the second diagonal block:
np.linalg.eigvals(S[1:3, 1:3])

# All eigenvalues; these include the ones from the second diagonal block, as expected
np.linalg.eigvals(S)

# similar matrices have the same eigenvalues: compare the spectra
eb = np.linalg.eigvals(B0)
es = np.linalg.eigvals(S)

import matplotlib as mpl
import matplotlib.pyplot as plt

colors = [c['color'] for c in mpl.rcParams['axes.prop_cycle']]

fig = plt.figure(figsize = (8,8))
ax = fig.gca()
ax.plot(np.real(eb), np.imag(eb), ls="none", marker="o", ms = 8, color = colors[0])
ax.plot(np.real(es), np.imag(es), ls="none", marker="o", ms = 16, color = "none", markeredgecolor = colors[1])
ax.set_xlabel("real")
ax.set_ylabel("imaginary")
ax.set_title("spectra of 'S' and 'B0'")
ax;