Skeleton for Benchmarking QRupdate

gistfile1.txt

## This skeleton provides a framework to develop a benchmark for various routines using QRupdate

## The following example tests the efficiency of obtaining the subset of linearly independent columns from
## various matrices. An application of these procedures is for instrumental least squares estimators (e.g., 2SLS).

using QRupdate
srand(0)
m = Int(1e6)
X = rand(m, 3)
X = hcat(X, 1.5 * X[:,2])
Z = rand(m, 2)
Z = hcat(Z, 2 * X[:,1])

function colbycol(X::AbstractMatrix, Z::AbstractMatrix)
    m = size(X, 1)
    n = size(X, 2) + size(Z, 2)
    mm = zeros(m,0)
    R = zeros(0,0)
    R = qraddcol(mm,R,X[:,1])
    mm = hcat(mm, X[:,1])
    LinearIndependent = BitVector([true])
    for col ∈ 2:size(X, 2)
        tmp = qraddcol(mm,R,X[:,col])
        push!(LinearIndependent, tmp[end,end] > sqrt(eps()))
        if LinearIndependent[end]
            R = tmp
            mm = hcat(mm, X[:,col])
        end
    end
    for col ∈ 1:size(Z, 2)
        tmp = qraddcol(mm, R, Z[:,col])
        push!(LinearIndependent, tmp[end,end] > sqrt(eps()))
        if LinearIndependent[end]
            R = tmp
            mm = hcat(mm, Z[:,col])
        end
    end
    return mm, R.'R, LinearIndependent
end
function qrandupdate(X::AbstractMatrix, Z::AbstractMatrix)
    mm = hcat(X, Z)
    R = qrfact(mm)[:R]
    LinearIndependent = abs.(diag(R)) .> sqrt(eps())
    mm = mm[:, LinearIndependent]
    for col ∈ sort(find(.!LinearIndependent), rev = true)
        R = qrdelcol(R, col)
    end
    return mm, R.'R, LinearIndependent
end
function noupdates(X::AbstractMatrix, Z::AbstractMatrix)
    mm = hcat(X, Z)
    R = qrfact(mm)[:R]
    LinearIndependent = abs.(diag(R)) .> sqrt(eps())
    if any(LinearIndependent)
        mm = mm[:,LinearIndependent]
        R = qrfact(mm)[:R]
    end
    return mm, R.'R, LinearIndependent
end

colbycol(X, Z)[[1,3]] == qrandupdate(X, Z)[[1,3]] == noupdates(X, Z)[[1,3]]
colbycol(X, Z)[2] ≈ qrandupdate(X, Z)[2] ≈ noupdates(X, Z)[2]

@time colbycol(X, Z)
@time qrandupdate(X, Z)
@time noupdates(X, Z)