QR decomposition by Householder projection for tridiagonal matrices in Julia and Python.

td_qr.jl

function householder!(x)
    x[1] = x[1] + sign(x[1]) .* norm(x)
    x ./= norm(x);
end

function tridiag_qr(T)
    Q = eye(size(T)...)
    R = copy(T)

    for i in 1:(size(R, 1) - 1)
        u = householder!(R[i:i+1, i])
        M = u * u'

        for j in 1:size(R, 2)
            # 2 optimized matrix multiplications
            # equivalent to R[i:i+1, j] -= 2 .* M * R[i:i+1, j] 
            tmp = 2 .* (M[1, 1] .* R[i, j] + M[1, 2] .* R[i+1, j])
            R[i+1, j] -= 2 .* (M[2, 1] .* R[i, j] + M[2, 2] .* R[i+1, j])
            R[i,j] -= tmp

            # similar to Q[i:i+1, j] -= 2 .* M * R[i:i+1, j], except all transposed
            tmp = 2 .* (M[1, 1] .* Q[j, i] + M[1, 2] .* Q[j, i+1])
            Q[j, i+1] -= 2 .* (M[2, 1] .* Q[j, i] + M[2, 2] .* Q[j, i+1])
            Q[j, i] -= tmp
        end
    end
    Q, R
end

function rand_tridiag(size)
    full(Tridiagonal(rand(size-1), rand(size), rand(size-1)))
end


function main()
    const SIZE = 1500
    const TRIALS = 20

    subsup = rand(SIZE - 1)
    diagonal = rand(SIZE)

    tridiag = Tridiagonal(subsup, diagonal, subsup)
    T = full(tridiag)

    for i=1:TRIALS
        println("$(@elapsed tridiag_qr(T))")
        #println("$(@elapsed qr(T))")
    end
end

function test()
    T = rand_tridiag(5)
    show(qr(T))
    println()
    show(tridiag_qr(T))
end

main()

td_qr.py

import numpy as np

def householder(x):
    x[0] = x[0] + np.sign(x[0]) * np.linalg.norm(x)
    x /= np.linalg.norm(x)
    return x

def tridiag_qr(T):
    R = T.copy()
    Qt = np.eye(T.shape[0])

    for i in xrange(T.shape[0] - 1):
        u = householder(T[i:i+2, i])

        R[i:i+2, :] = R[i:i+2, :] - 2 * np.outer(u, np.dot(u, R[i:i+2, :]))
        Qt[i:i+2, :] = Qt[i:i+2, :] - 2 * np.outer(u, np.dot(u, Qt[i:i+2, :]))

    return Qt.T, R


def main():
    SIZE = 1500
    TRIALS = 20

    subsup = np.random.rand(SIZE - 1)
    diagonal = np.random.rand(SIZE)
    T = np.diag(subsup, 1) + np.diag(subsup, -1) + np.diag(diagonal)

    import time

    for _ in xrange(TRIALS):
        tic = time.clock()
        tridiag_qr(T)
        #np.linalg.qr(T)
        print time.clock() - tic

if __name__ == '__main__':
    main()

You should be able to speed this up a bit by using this definition of householder!:

function householder!(x)
    n = norm(x)
    x[1] = x[1] + sign(x[1]) * n
    for i in 1:length(x)
        x[i] /= n
    end
end

A little bit of devectorization of these two lines would also go a long way:

    R[i:i+1, :] = R[i:i+1, :] - 2 .* u * u' * R[i:i+1, :]
    Qt[i:i+1, :] = Qt[i:i+1, :] - 2 .* u * u' * Qt[i:i+1, :]