The upper-triangular Hessenberg form has zeros below the first sub-diagonal
$$ \begin{bmatrix} * & * & * & * & * & * \\ * & * & * & * & * & * \\ 0 & * & * & * & * & * \\ 0 & 0 & * & * & * & * \\ 0 & 0 & 0 & * & * & * \\ 0 & 0 & 0 & 0 & * & * \\ \end{bmatrix} $$import numpy as np
# We need sign function with sign(0) = 1
def sign(x):
if x >= 0.0:
return 1.0
else:
return -1.0
# Algorithm 26.1
def hessenberg(A):
m = A.shape[0]
for k in range(m-2):
x = A[k+1:m, k]
e1 = np.zeros_like(x)
e1[0] = 1.0
v = sign(x[0]) * np.linalg.norm(x) * e1 + x
v = v / np.linalg.norm(v)
A[k+1:m,k:m] = A[k+1:m,k:m] - 2.0 * np.outer(v, np.dot(v, A[k+1:m,k:m]))
A[0:m,k+1:m] = A[0:m,k+1:m] - 2.0 * np.outer(A[0:m,k+1:m] @ v, v)
return A
Test on a random matrix.
m = 7
A = np.random.rand(m,m)
H = hessenberg(A)
print(np.array_str(H, precision=2))
[[ 1.46e-02 -9.73e-01 3.66e-01 -2.42e-01 4.67e-01 -6.16e-01 6.15e-02] [-1.24e+00 2.15e+00 -8.66e-01 4.98e-03 -2.90e-01 -2.23e-01 -1.51e-01] [ 0.00e+00 -1.76e+00 1.74e-01 -1.22e-01 5.24e-01 -8.45e-03 1.60e-01] [-2.78e-17 2.22e-16 5.69e-01 -1.04e-02 1.19e-01 2.12e-01 -3.13e-01] [ 0.00e+00 2.22e-16 0.00e+00 2.86e-01 1.78e-01 1.56e-01 -5.46e-01] [ 0.00e+00 -2.22e-16 0.00e+00 0.00e+00 5.48e-02 1.30e-01 -1.65e-01] [ 0.00e+00 -8.67e-19 0.00e+00 0.00e+00 -8.67e-19 1.04e-01 -1.25e-01]]
Extract the lower triangular part which is supposed to be zero.
L = np.tril(H,k=-2)
print(np.array_str(L, precision=2))
[[ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [ 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [-2.78e-17 2.22e-16 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [ 0.00e+00 2.22e-16 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [ 0.00e+00 -2.22e-16 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00] [ 0.00e+00 -8.67e-19 0.00e+00 0.00e+00 -8.67e-19 0.00e+00 0.00e+00]]
Compute the maximum over the zero elements.
print(np.abs(L).max())
2.220446049250313e-16
Test on a larger matrix without printing.
m = 100
A = np.random.rand(m, m)
H = hessenberg(A)
L = np.tril(H,k=-2)
print(np.abs(L).max())
1.1102230246251565e-15
In this case, the Hessenberg reduction should give a tridiagonal matrix.
m = 7
B = np.random.rand(m,m)
A = 0.5 * (B + B.T)
H = hessenberg(A)
print(np.array_str(H, precision=2, suppress_small=True))
[[ 0.91 -1.34 0. 0. -0. 0. -0. ] [-1.34 2.56 -1.57 -0. -0. 0. -0. ] [ 0. -1.57 0.79 0.31 0. -0. -0. ] [ 0. -0. 0.31 -0.09 0.26 -0. 0. ] [ 0. 0. 0. 0.26 -0.21 0.23 0. ] [ 0. 0. 0. 0. 0.23 0.41 -0.51] [-0. 0. 0. 0. -0. -0.51 0.03]]
The above function may not work correctly if the input matrix is complex. Write a complex version and test it.