Reduction to Hessenberg form

Reduction to Hessenberg form#

The upper-triangular Hessenberg form has zeros below the first sub-diagonal

\[\begin{split} \newcommand{\star}{\color{red}*} \begin{bmatrix} \star & * & * & * & * & * \\ * & \star & * & * & * & * \\ 0 & * & \star & * & * & * \\ 0 & 0 & * & \star & * & * \\ 0 & 0 & 0 & * & \star & * \\ 0 & 0 & 0 & 0 & * & \star \\ \end{bmatrix} \end{split}\]
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))

[[ 3.60e-01 -7.11e-01  7.19e-01  5.49e-01 -5.56e-01  6.03e-02 -5.06e-01]
 [-1.01e+00  2.13e+00 -1.20e+00 -5.06e-02  1.88e-01  1.89e-01  1.72e-01]
 [ 1.11e-16 -9.62e-01  5.47e-01 -3.44e-01 -1.53e-01  1.09e-01 -6.03e-01]
 [ 0.00e+00  0.00e+00 -8.04e-01 -2.18e-01 -4.06e-01 -1.63e-01  1.46e-01]
 [ 0.00e+00  0.00e+00 -6.94e-18  5.17e-01 -1.39e-01  1.75e-01  2.09e-01]
 [ 2.78e-17  0.00e+00 -1.11e-16  0.00e+00 -5.89e-01  1.98e-01 -5.40e-01]
 [ 5.55e-17  0.00e+00 -2.78e-17  0.00e+00 -1.11e-16 -8.54e-02  2.20e-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]
 [ 1.11e-16  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 -6.94e-18  0.00e+00  0.00e+00  0.00e+00  0.00e+00]
 [ 2.78e-17  0.00e+00 -1.11e-16  0.00e+00  0.00e+00  0.00e+00  0.00e+00]
 [ 5.55e-17  0.00e+00 -2.78e-17  0.00e+00 -1.11e-16  0.00e+00  0.00e+00]]

Compute the maximum over the zero elements.

print(np.abs(L).max())

1.1102230246251565e-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())

8.881784197001252e-16

Hermitian case#

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.37 -1.3   0.    0.    0.    0.    0.  ]
 [-1.3   2.21  0.83  0.    0.   -0.   -0.  ]
 [ 0.    0.83 -0.23 -0.55  0.    0.   -0.  ]
 [ 0.    0.   -0.55  0.38  0.2   0.   -0.  ]
 [ 0.    0.    0.    0.2   0.5   0.1   0.  ]
 [ 0.   -0.    0.   -0.    0.1   0.26 -0.29]
 [ 0.   -0.    0.   -0.    0.   -0.29 -0.03]]

Remark#

The above function may not work correctly if the input matrix is complex. Write a complex version and test it.

Contents