Should UpperTriangular and LowerTriangular really be separate types?

[timholy]

Currently we have two types, UpperTriangular and LowerTriangular. (Internally, there are also UnitUpperTriangular and UnitLowerTriangular, but let's ignore these for now.) The existence of these two types means that cholfact(A, :U) cannot be type-stable. We have the Val mechanism, but that seems to be agreed to be a little ugly.

Is there really a justification for these being two separate types?

immutable Triangular{T,A<:AbstractMatrix} <: AbstractMatrix{T}
    isupper::Bool
    data::A
end
Triangular(isupper::Bool, A::AbstractMatrix) = Triangular{eltype(A),typeof(A)}(isupper, A)

function Base.A_ldiv_B!(A::Triangular, b::StridedVecOrMat)
    if A.isupper
        return Base.LAPACK.trtrs!('U', 'N', 'N', A.data, b)
    else
        return Base.LAPACK.trtrs!('L', 'N', 'N', A.data, b)
    end
end

function manymul(A, b, n)
    c = similar(b)
    s = zero(eltype(c))
    for i = 1:n
        A_mul_B!(c, A, b)
        for k = 1:length(c)
            s += c[k]
        end
    end
    s
end

function manydiv(A, b, n)
    c = similar(b)
    s = zero(eltype(c))
    for i = 1:n
        copy!(c, b)
        A_ldiv_B!(A, c)
        for k = 1:length(c)
            s += c[k]
        end
    end
    s
end


A = Float64[1 2; 0 3]
U = UpperTriangular(A)
b = rand(2)

manymul(A, b, 1)
manymul(U, b, 1)
manydiv(Triangular(true, A), b, 1)
manydiv(U, b, 1)
println("warmup @time:")
@time 1
println("regular multiplication, 2x2")
gc()
@time manymul(A, b, 10^5)
println("upper triangular multiplication, 2x2")
gc()
@time manymul(U, b, 10^5)
println("manually-dispatched division, 2x2")
gc()
@time manydiv(Triangular(true, A), b, 10^5)
println("upper triangular division, 2x2")
gc()
@time manydiv(U, b, 10^5)

A = Float64[1 2 3 4;
     0 5 6 7;
     0 0 8 9;
     0 0 0 4]
U = UpperTriangular(A)
b = rand(4)

println("regular multiplication, 4x4")
gc()
@time manymul(A, b, 10^5)
println("upper triangular multiplication, 4x4")
gc()
@time manymul(U, b, 10^5)
println("manually-dispatched division, 4x4")
gc()
@time manydiv(Triangular(true, A), b, 10^5)
println("upper triangular division, 4x4")
gc()
@time manydiv(U, b, 10^5)

Results:

warmup @time:
  0.000003 seconds (148 allocations: 10.151 KB)
regular multiplication, 2x2
  0.012172 seconds (7 allocations: 272 bytes)
upper triangular multiplication, 2x2
  0.015659 seconds (7 allocations: 272 bytes)
manually-dispatched division, 2x2
  0.032639 seconds (100.01 k allocations: 7.630 MB)
upper triangular division, 2x2
  0.031850 seconds (100.01 k allocations: 7.630 MB)
regular multiplication, 4x4
  0.013056 seconds (7 allocations: 288 bytes)
upper triangular multiplication, 4x4
  0.017556 seconds (7 allocations: 288 bytes)
manually-dispatched division, 4x4
  0.033030 seconds (100.01 k allocations: 7.630 MB)
upper triangular division, 4x4
  0.035649 seconds (100.01 k allocations: 7.630 MB)

With these performance numbers, I'm skeptical that we need these to be separate types. If there were just a boolean flag, we could have cholfact(A, :U) be type-stable.

[timholy]

Reference: https://groups.google.com/d/msg/julia-users/mGV-V-FYzrA/nHOF-_T2BAAJ

[andreasnoack]

1. I agree that chol(A, Val{:L}) is unbbearable ugly
2. Only chol, not cholfact needs the Val wrapper because cholfact is not parametric on :U/:L anymore.
3. The reason we have Upper/LowerTriangular is not fast dispatch but to make the types closed under basic arithmetic. Without the upper/lower in the type we cannot have LowerTringular+-*/\LowerTriagnular -> LowerTriangular because LowerTriangular+-*\/UpperTriangular -> Matrix.

Right now my preferred solution is to remove the upper/lower argument to chol. There is really no reason to expose this to the user.

[timholy]

We can still have your point 3 if instead of calling it Triangular we have

immutable SpecialMatrix{T, A<:AbstractMatrix} <: AbstractMatrix{T}
    shape::Symbol    # :U, :L, or :F (for full)
    data::A
end

So sometimes a SpecialMatrix will be no different from Matrix, but when it is one can (manually) dispatch to particularly efficient algorithms.

That said, I think I like your plan of removing the upper/lower argument.

[tkelman]

JuliaLang/julia#8240 would also be closed under arithmetic, and could maybe even still use the lapack triangular routines for the bandwidth = entire dimension cases by offsetting lda by 1. The SpecialMatrix proposed above is not generic enough to capture anything other than full or triangular, and for structures that don't fit in as special cases of banded, I'd rather use the type system and existing dispatch mechanisms than try to recreate them with flag fields.

[timholy]

SpecialMatrix was just intended for illustration: the point is that using flags is sometimes better than dispatch because it avoids a source of type-instability.

I'd rather use the type system and existing dispatch mechanisms than try to recreate them with flag fields

A year ago I might have agreed without hesitation, but now I'm less sure, particularly in situations like this one. Consider: checking a single flag is O(1), compared to the O(N^2) or O(N^3) of most linear algebra operations, and is much cheaper than runtime dispatch (which is what we rely on now in cases of type-instability). It also reduces the amount of time spent on JIT-compilation, which remains one of the main annoyances of writing and using big applications in julia. The downside is that this is basically C and not julia anymore.

Probably the best fix would be to have the compiler better at handling type-unstable code, and most likely it's best to see if we can hold out for this. A workaround might be a set of tools for making manual dispatch a little less awful. I haven't even begun to think what that might look like.

[toivoh]

I agree in general that we should consider using flags instead of creating more types, but in this case I think that would be problematic when it comes to type stability. The reason is that LL = L and UU = U, but LU and UL produce general matrices, so triangular matrices are not closed under matrix multiplication, you need to know at the type level whether the matrices are upper or lower triangular.

If we would subsume both types into a single BandedMatrix type I suppose it could work out, but the compiler couldn't know when one had actually become a full matrix.

[toivoh]

Another possibility of we get TransposedMatrix would be to define only one of upper and lower and use a transposed view for the other. (I realize this doesn't address your original concern.)

[andreasnoack]

@timholy With JuliaLang/julia#13680 merged, can we close this one or do you still think we need to consider this?