RFC:Add generic Cholesky decomposition and make Cholesky parametric on matrix type

[andreasnoack]

There are still some methods that need to be adjusted for this to work completely, but I think most of the functionality should be there.

[mlubin]

👍

[ViralBShah]

Nice!

[andreasnoack]

I think it is ready for review now so also time for some advertising. The Cholesky is generic both in the sense that you can use e.g. BigFloats

julia> A=big(randn(2,2));
julia> Apd=A'A;
julia> cholfact(Apd)
Cholesky{BigFloat,Array{BigFloat,2},:U} with factor:
2x2 Triangular{BigFloat,Array{BigFloat,2},:U,false}:
 1.857774050831274188486920823465157015411979051018802312773493386099520585257146e+00  …  1.76048668500619868309132341992563358432046517656847649557289132350834427503267e-01 
 0e+00                                                                                    5.829390443598957440055707197674804936840111279347445154902807401513834005139188e-01

and in the sense that you can factorise block matrices

julia> A=randn(4,4);
julia> Apd=A'A;
julia;Abl=Matrix{Float64}[Apd[(1:2)+2*(i-1), (1:2) + 2*(j-1)] for i = 1:2, j = 1:2]
2x2 Array{Array{Float64,2},2}:
 2x2 Array{Float64,2}:
  2.501     -0.861278
 -0.861278   3.16783   …  2x2 Array{Float64,2}:
 1.28158    0.655392
 0.625918  -3.63314 
 2x2 Array{Float64,2}:
 1.28158    0.625918
 0.655392  -3.63314        2x2 Array{Float64,2}:
 4.42572   1.74632
 1.74632  10.6957     

julia> cholfact!(copy(Abl))
Cholesky{Array{Float64,2},Array{Array{Float64,2},2},:U} with factor:
2x2 Triangular{Array{Float64,2},Array{Array{Float64,2},2},:U,false}:
 2x2 Array{Float64,2}:
 1.58145  -0.544611
 0.0       1.69447   …  2x2 Array{Float64,2}:
 0.81038    0.414424
 0.629849  -2.01091 
 2x2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0                         2x2 Array{Float64,2}:
 1.83638  1.45778
 0.0      2.08686      

because it isn't assumed that the elements commute under multiplication and therefore it also works for quaternions

julia> Q = Quaternions.Quaternion{Float64}[Quaternions.Quaternion(randn(4)...) for i = 1:2,j=1:2]
2x2 Array{Quaternion{Float64},2}:
 -0.0171977-0.886764i-0.427188j-1.4489k  …     0.192185-0.1924i+0.145574j+0.520387k
   0.156593-0.253024i+0.44728j+1.07369k     -0.0338218+1.12859i+0.244676j-0.362112k

julia> Qpd=Q'Q;

julia> cholfact(Qpd)
Cholesky{Quaternion{Float64},Array{Quaternion{Float64},2},:U} with factor:
2x2 Triangular{Quaternion{Float64},Array{Quaternion{Float64},2},:U,false}:
 4.50987+0.0i+0.0j+0.0k  -0.270315+0.172501i-0.38301j+0.227746k
     0.0+0.0i+0.0j+0.0k                  1.53041+0.0i+0.0j+0.0k

[ViralBShah]

The block matrix stuff is really cool. Perhaps it can be leveraged for a multi threaded implementation.

[andreasnoack]

I would like to try something like that.

[tknopp]

I heard multithreading, so I have to link #6741 ;-)
@andreasnoackjensen would be interesting if the parapply function in #6741 (see https://github.com/tknopp/julia/blob/jlthreading/base/threading.jl and https://github.com/JuliaLang/julia/pull/6741/files#diff-96e9e68057e887be8a21bb1275df1473R89 for an examplary parallel matrix vector multiplication) is flexible enough to address the needs for a parallel Cholesky implementation.

[simonbyrne]

It occurs to me that if UpLo == d, we don't need to triu!/tril!.

[ViralBShah]

@tknopp We will need some ways to manage dependencies within threads. In the distributed case, we have RemoteRefs that can be used. There should be a way to start out with something simple though. See the parallel implementation in examples/hpl.jl, which should be possible to convert into a Cholesky implementation.

[andreasnoack]

Rebased. I think this one is ready for merge, so I'll do so soon if no objections are made.

[jiahao]

👍

Also looking forward to the LDL' variant ;)

[IainNZ]

This just broke quite a few packages, not sure if that was intentional?

[IainNZ]

(also congratulations for having the first PR to break a nontrivial number of packages 🍰, I was wondering what would be the first one)

[timholy]

And private code, too!

But that's fine.

[ivarne]

Should probably add a NEWS.md entry if the breakage was intentional.

[andreasnoack]

Sorry for the silence and delay here. I have had some days in bed.

@IainNZ Thank you for pointing this out. It is still not possible to run the package tests on a pull request right?

@ivarne You are right. Done in f0a87e1

[IainNZ]

It is possible, but needs to be done manually