Householder comparison

householder.py

import perfplot
import numpy as np


def householder(A):
    m, n = A.shape
    R = A.copy()
    Q = np.identity(m)

    for j in range(n):
        hr = HouseholderReflection(R[j:, j])
        hr.apply(R[j:])
        hr.apply(Q[j:])

    Q = Q[:n].T
    R = np.triu(R[:n])
    return Q, R


class HouseholderReflection:
    def __init__(self, x):
        v = x.copy()
        v = v.reshape(-1)

        # Householder only works with the Euclidean inner product
        self.dot = lambda x, y: np.dot(x.conj(), y)

        gamma = v[0].copy()
        v[0] = 1

        sigma2 = self.dot(v[1:], v[1:])
        xnorm = np.sqrt(np.abs(gamma) ** 2 + sigma2)

        # is x a multiple of first unit vector?
        if sigma2 == 0:
            beta = 0
            xnorm = np.abs(gamma)
            alpha = 1 if gamma == 0 else gamma / xnorm
        else:
            beta = 2
            if gamma == 0:
                v[0] = -np.sqrt(sigma2)
                alpha = 1
            else:
                v[0] = gamma + gamma / np.abs(gamma) * xnorm
                alpha = -gamma / np.abs(gamma)

        self.xnorm = xnorm
        self.v = v / np.sqrt(np.abs(v[0]) ** 2 + sigma2)
        self.alpha = alpha
        self.beta = beta

    def apply(self, x: np.ndarray) -> None:
        if self.beta == 0:
            return
        x -= self.beta * np.outer(self.v, self.dot(self.v, x))


b = perfplot.bench(
    setup=lambda n: np.random.rand(n, n),
    kernels=[householder, np.linalg.qr],
    n_range=[2 ** k for k in range(12)]
)
b.save("out.png")
b.show()

numpy.linalg.qr is much faster (about factor 50):