In [1]:
import numpy as np
from scipy import linalg
import pylab as pl
In [2]:
def house(x):
    """
    Given a vetor x, computes vectors v with v[0] = 1
    and scalar beta such that P = I - beta v v^T
    is orthogonal and P x = ||x|| e_1

    Parameters
    ----------
    x : array, shape (n,) or (n, 1)

    Returns
    -------
    beta : scalar
    v : array, shape (n, 1)
    """
    x = np.asarray(x)
    if x.ndim == 1:
        x = x[:, np.newaxis]
    sigma = linalg.norm(x[1:, 0]) ** 2
    v = np.vstack((1, x[1:]))
    if sigma == 0:
        beta = 0
    else:
        mu = np.sqrt(x[0, 0] ** 2 + sigma)
        if x[0, 0] <= 0:
            v[0, 0] = x[0, 0] - mu
        else:
            v[0, 0] = - sigma / (x[0, 0] + mu)
        beta = 2 * (v[0, 0] ** 2) / (sigma + v[0, 0] ** 2)
        v /= v[0, 0]
    return beta, v
In [3]:
n = 5
np.random.seed(0)
x = np.random.randn(n)
beta, v = house(x)
P = np.eye(n) - beta * np.dot(v, v.T)
print np.round(P.dot(x) / linalg.norm(x), decimals=15) 
[ 1. -0. -0. -0. -0.]

QR decomposition using Householder reflections. To generate the animated gif from the saved pictures I run the command-line convert -delay 30 -loop 0 house_*.png house.gif

In [10]:
n = 20
X = np.random.randn(n, n)
pl.matshow(np.abs(X), cmap=pl.cm.Paired, vmax=5.)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(8.5,6.0)
pl.axis('off')
pl.colorbar()
pl.savefig('house_random.png')

X_0 = X.copy()
title = 'X'
for i in range(X.shape[1] -1):
    title = ('P_{%s} ' % i) + title
    beta, v = house(X_0[i:, i])
    P = np.eye(X.shape[1])
    P[i:, i:] = np.eye(n - i) - beta * np.dot(v, v.T)
    X_0 = P.dot(X_0)
    pl.matshow(np.abs(X_0), cmap=pl.cm.Paired, vmax=5.)
    fig = matplotlib.pyplot.gcf()
    fig.set_size_inches(8.5,6.0)
    pl.title('$' + title + '$', fontsize='xx-large',x =1., ha='right')
    pl.axis('off')
    pl.colorbar()
    pl.savefig('house_%03d.png' % i)
In [5]:
X = (X + X.T) / 2.

X_0 = X.copy()
title = 'X'
for i in range(X.shape[1] - 2):
    title = ('P_{%s} ' % i) + title + ('P_{%s} ' % i)
    beta, v = house(X_0[i+1:, i])
    P = np.eye(X.shape[1])
    P[i+1:, i+1:] = np.eye(n - i - 1) - beta * np.dot(v, v.T)
    X_0 = P.T.dot(X_0).dot(P)
    pl.matshow(np.abs(X_0), cmap=pl.cm.Paired, vmax=5.)
    fig = matplotlib.pyplot.gcf()
    fig.set_size_inches(8.5,6.0)
    pl.title('$' + title + '$', fontsize='xx-large',x =.5, ha='center')
    pl.axis('off')
    pl.colorbar()
    pl.savefig('house_tridiag_%03d.png' % i)