[WIP] cholesky

[willtebbutt]

A first attempt at implementing the cholesky factorisation for Kronecker matrices, will close #10 if completed. Attempting to implement this has highlighted a few aspects of the codebase that we really should sort out:

1. GeneralizedKroneckerProduct <: AbstractMatrix{Number} - this was bad causes type stability issues. If you don't parametrise GeneralisedKroneckerProduct{T}, you wind up constructing things like Cholesky{Number, ...}. This required quite a few changes throughout the codebase. Seems to work fine (tests pass), but of course I might have missed something. This also seems to resolve an issue whereby attempting to display UpperTriangular(d) for some KroneckerProduct d threw an error.
2. We're not doing a good job of handling UpperTriangular and LowerTriangular. For Cholesky stuff to really work, we need this to work. The approach I've taken here is to cheat and make UpperTriangular of an AbstractKroneckerProduct return an AbstractKroneckerProduct of UpperTriangular matrices. While this makes sense mathematically I'm not sure that it's a great move semantically. It's also been causing some stack overflow errors that I've not managed to track down.

I also took the (highly opinionated) liberty of renaming the files in test so that they don't begin with the word test - this word being redundant as they're in the test folder. Please feel free to suggest that I revert this :)

Let me know what you think :)

[MichielStock]

1. Thanks for pointing this out. @Beramos and I tried to improve it at the hackathon at Juliacon, but I think we did not 100% got the hang of it. General parametric type is a good idea, because ideally Kronecker.jl should work with other stuff than numbers, such as strings.
2. This should work I think and is the approach I would have taken. I do worry that the LowerTriangular of a Kronecker product is not the Kronecker product of the LT of the individual matrices. The latter itself is lower triangular though (even sparser).

I will try out your branch and see I can resolve the errors.

[MichielStock]

I fixed an error here by changing it into

function getproperty(C::KroneckerCholesky, d::Symbol)
    Cfactors = getfield(C, :factors)
    Cuplo    = getfield(C, :uplo)
    info     = getfield(C, :info)

    Cuplo != 'U' && throw(NotImplementedError(""))

    if d == :U
        return UpperTriangular(Cfactors)
    elseif d == :L
        return LowerTriangular(copy(Cfactors'))
        #return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
    elseif d == :UL
        return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
    else
        return getfield(C, d)
    end
end

I also needed my own version of copy which provides shallow copies of Kronecker products. This gives the right answer.

[MichielStock]

I am not sure this is correct. If I understand it correctly, the main performance of Cholesky decomposition originates from the fact that it is easy to solve a triangular system. The vec trick kind of destroys this structure...

[willtebbutt]

Hmm. I admit that I hadn't really thought very hard about this when writing it, but I'm pretty sure that you can use the vec trick to compute backsolves. Note that

vec(inv(L2) X inv(L1)) = (inv(L1) ⊗ inv(L2)) vec(X) = inv(L1 ⊗ L2) vec(X)

so we should be able to implement backsolving with kronecker matrices efficiently provided that the individual backsolves are efficient.

[MichielStock]

So you would leave it as is for the moment?

[willtebbutt]

This should work I think and is the approach I would have taken. I do worry that the LowerTriangular of a Kronecker product is not the Kronecker product of the LT of the individual matrices. The latter itself is lower triangular though (even sparser).

You're right. I shouldn't have been doing the above. The cholesky factorisation of a kronecker product is the kronecker product of the cholesky factors though, which is what I would ultimately like to exploit. Probably requires some care to make this work correctly.

[dkarrasch]

Why do you require that both Kronecker factors have the same eltype T1? I don't think that's necessary.

[dkarrasch]

Should this be

chol_P, chol_Q = cholesky(P; check=check), cholesky(Q; check=check)

? Otherwise, the kwarg doesn't enter anywhere.

[dkarrasch]

UpperTriangular is the constructor of the UpperTriangular type. So I guess this should return an object of that type, and you may wish to consider returning

UpperTriangular(UpperTriangular(A) ⊗ UpperTriangular(B))

? Otherwise, subsequent code will not benefit from upper-triangularizing an AbstractKroneckerProduct, it will be "invisible" from the type point of view.

[willtebbutt]

Closing in favour of #50