Return `UpperTriangular` from `chol`

[bshall]

To be consistent with the methods from Base I was wondering if the cholesky decomposition of a static array should return an UpperTriangular type? I know the AbstractTriangular types have setindex! and other mutating methods defined for them which may be an issue?

[andyferris]

Yeah - consistency with Base is definitely a primary consideration of this package.

I did have an idea once to implement our own triangular types, with the idea of keeping only (n^2 + n)/2 elements on the stack, as a micro-optimization. Not sure if that would ever by worthwhile. As it is, I'm not sure if the triangular type is simply slower than an SMatrix with extra zeros - probably solvable using ifelse inside getindex.

OTOH we can add more specializations to methods on UpperTriangular{<:Any, <:StaticMatrix} and LowerTriangular{<:Any, <:StaticMatrix} as necessary.

[c42f]

Personally I wouldn't bother worrying about the extra elements. There's a good chance that for really small matrices they're never even stored on the stack anyway, just passed in registers.

So wrapping in an UpperTriangular type makes a lot of sense to me.

Xref regarding cholfact - #192 (comment)

[andyferris]

So wrapping in an UpperTriangular type makes a lot of sense to me.

OK - let's do this.

One way to optimize around branching and sub2ind for immutable arrays is to fill in the zeros correctly upon construction, and mainly use the UpperTriangular tag for dispatch purposes (and allow fast, direct linear getindex for speed).

[c42f]

use the UpperTriangular tag for dispatch purposes and allow fast, direct linear getindex for speed

Sounds reasonable to me.

[KristofferC]

FWIW we store only the lower triangular part and use non linear indexing (with a branch) at Tensors.jl (https://github.com/KristofferC/Tensors.jl/blob/master/src/indexing.jl#L5). For the generated code, I have never seen a branch though, so the whole thing gets optimized away nicely.

[c42f]

@KristofferC Right it makes sense that the branches disappear if the whole thing is unrolled.

Ok, a good reason to have special storage would be if people want to store a lot of triangular static arrays inside an array. In that case it's a big win to only store the upper/lower part.

[c42f]

It seems we did this somewhere along the way:

julia> cholesky(@SMatrix[1 2; 2 4])
Cholesky{Float64,SArray{Tuple{2,2},Float64,2,4}}
U factor:
2×2 UpperTriangular{Float64,SArray{Tuple{2,2},Float64,2,4}}:
 1.0  2.0
  ⋅   0.0