Householder Similarity Transforms

In [2]:
import numpy as np
import numpy.linalg as la

np.set_printoptions(precision=2, linewidth=150)
In [3]:
n = 8

e1 = np.zeros(n); e1[0] = 1
e2 = np.zeros(n); e2[1] = 1

A = np.random.randn(n, n)
A
Out[3]:
array([[ 1.54,  1.4 , -0.32,  0.35,  1.16,  0.86, -0.07,  1.67],
       [ 0.34,  0.56, -1.68, -0.9 , -0.24, -1.35,  1.57,  0.07],
       [-0.69, -0.81, -0.99, -0.56,  0.45,  1.55, -0.98,  0.82],
       [-0.85,  0.57, -0.51,  1.42,  0.31,  1.41,  0.18, -0.43],
       [-0.46, -0.72, -1.13,  0.91,  1.56,  1.17,  0.58,  0.  ],
       [-0.68, -1.53,  0.84,  1.05,  0.36,  0.  , -1.01,  2.12],
       [ 0.02,  0.87,  1.88, -0.91,  0.3 , -1.26,  0.08, -0.36],
       [-0.12,  1.21, -0.75,  0.18, -1.83, -0.47,  0.11, -1.83]])

Now try to zero the first column with a similarity transform.

Starting with the first row

Lets first try to proceed as in Householder QR, creating a transformation $$H=I-2\frac{vv^T}{v^Tv}$$ where $v = a_1 - ||a_1||_2e_1$ with $a_1$ being the first column of $A$.

In [4]:
a = A[:, 0].copy()
v = a-la.norm(a)*e1

H1 = np.eye(n) - 2*np.outer(v, v)/(v@v)

We can apply the transformation from the left as in QR to reduce the first column to a multiple of the first elementary vector.

In [5]:
H1@A
Out[5]:
array([[  3.66e+00,   1.95e+00,  -1.71e-01,  -5.52e-01,   1.24e+00,   1.29e-01,   2.61e-02,  -2.09e+00],
       [  2.24e-16,   9.07e-01,  -1.63e+00,   4.23e-01,   3.42e-01,   7.02e-01,   6.64e-01,   2.02e+00],
       [ -1.14e-16,   4.67e-02,  -4.70e-01,   1.46e+00,  -4.03e-01,  -8.46e-01,   2.05e-01,   5.92e-01],
       [  3.85e-16,  -1.10e+00,   3.98e-01,   5.12e-01,  -8.57e-01,  -8.08e-01,  -1.83e-01,   1.17e+00],
       [  3.34e-17,  -7.27e-01,   1.25e-01,   9.67e-01,  -7.66e-01,  -2.40e-01,   2.24e+00,   1.55e+00],
       [ -1.88e-16,   1.02e+00,   1.60e+00,  -1.42e+00,  -1.82e-01,   2.04e-02,  -3.10e-02,   4.73e-01],
       [  3.06e-17,  -7.74e-01,  -1.29e+00,  -2.87e-01,  -5.34e-01,  -9.70e-01,  -2.01e+00,  -5.13e-01],
       [ -1.46e-16,  -1.06e+00,  -9.55e-01,  -1.74e+00,   6.56e-01,  -1.42e+00,   2.06e+00,  -6.85e-01]])

However, to ensure we do not perturb the eigenvalues of $A$, we must also apply the matrix from the right, resulting in a similarity transformation.

In [5]:
H1@A@H1.T
Out[5]:
array([[ 2.38,  1.58,  0.42, -0.4 ,  0.43, -0.56,  0.64,  0.62],
       [ 0.73, -0.11, -1.01,  0.96,  0.96,  0.64,  1.09,  0.92],
       [ 1.37, -1.21,  1.19,  0.04,  0.36,  1.01, -0.12, -0.28],
       [ 1.15,  0.39,  1.51,  1.36, -0.23,  0.1 ,  1.25, -1.89],
       [ 0.23, -0.57, -0.53,  0.27,  0.95, -0.01,  1.17, -0.87],
       [ 0.41, -1.35,  1.82,  0.21, -0.5 , -1.68, -0.14,  0.83],
       [ 0.21,  0.73,  2.13, -0.54,  0.51, -0.93,  0.05, -0.27],
       [ 1.32,  0.47,  0.98,  1.94, -0.96,  0.76,  0.22, -1.78]])

Note that applying the Householder transformation from the right filled in the elements annihilated by applying it from the left.

Starting in the second row

To avoid this, we define the Householder transformation to annihilate elements below the first subdiagonal. That way, the first transformation does not affect the first row when applied from the left, and consequently does not affect the first column when applied for the right, preserving the zeros we've annihilated.

In [8]:
a = A[:, 0].copy()
a[0] = 0
v = a-la.norm(a)*e2

H2 = np.eye(n) - 2*np.outer(v, v)/(v@v)
In [9]:
H2 @ A
Out[9]:
array([[  1.54e+00,   1.40e+00,  -3.17e-01,   3.51e-01,   1.16e+00,   8.59e-01,  -6.80e-02,   1.67e+00],
       [  1.42e+00,   1.07e+00,   4.34e-01,  -1.62e+00,  -9.81e-01,  -2.29e+00,   1.04e+00,  -9.96e-01],
       [  1.11e-16,  -4.76e-01,   3.73e-01,  -1.03e+00,  -2.16e-02,   9.48e-01,  -1.31e+00,   1.34e-01],
       [  1.73e-18,   9.70e-01,   1.16e+00,   8.44e-01,  -2.72e-01,   6.72e-01,  -2.33e-01,  -1.27e+00],
       [  3.56e-17,  -5.06e-01,  -2.30e-01,   6.00e-01,   1.25e+00,   7.73e-01,   3.61e-01,  -4.46e-01],
       [  5.81e-17,  -1.21e+00,   2.18e+00,   5.93e-01,  -1.03e-01,  -5.88e-01,  -1.35e+00,   1.44e+00],
       [ -4.20e-18,   8.63e-01,   1.84e+00,  -9.00e-01,   3.13e-01,  -1.25e+00,   9.43e-02,  -3.43e-01],
       [  1.39e-17,   1.27e+00,  -5.09e-01,   9.98e-02,  -1.92e+00,  -5.76e-01,   5.30e-02,  -1.95e+00]])
In [10]:
H2 @ A @ H2.T
Out[10]:
array([[  1.54e+00,  -6.49e-01,  -1.64e+00,  -1.27e+00,   2.94e-01,  -4.46e-01,  -2.60e-02,   1.44e+00],
       [  1.42e+00,   2.54e+00,   1.38e+00,  -4.60e-01,  -3.59e-01,  -1.35e+00,   1.01e+00,  -8.29e-01],
       [  1.11e-16,  -1.63e-01,   5.75e-01,  -7.82e-01,   1.11e-01,   1.15e+00,  -1.32e+00,   1.70e-01],
       [  1.73e-18,  -9.71e-01,  -8.72e-02,  -6.90e-01,  -1.09e+00,  -5.63e-01,  -1.93e-01,  -1.49e+00],
       [  3.56e-17,  -1.10e+00,  -6.13e-01,   1.31e-01,   9.97e-01,   3.96e-01,   3.73e-01,  -5.14e-01],
       [  5.81e-17,  -1.54e+00,   1.96e+00,   3.26e-01,  -2.46e-01,  -8.03e-01,  -1.34e+00,   1.40e+00],
       [ -4.20e-18,   3.79e-01,   1.53e+00,  -1.28e+00,   1.08e-01,  -1.55e+00,   1.04e-01,  -3.98e-01],
       [  1.39e-17,   1.56e+00,  -3.22e-01,   3.29e-01,  -1.79e+00,  -3.91e-01,   4.71e-02,  -1.91e+00]])

To generalize this process, we continue to eliminate everything below the subdiagonal in the next column and applying the two-sided transformations, finally resulting in an upper-Hessenberg matrix.