No method to construct factorization objects from pre-existing factors

[simonbyrne]

If I have an upper-triangular matrix U, there is no method by which I can wrap it in a CholeskyDense type easily. An example use case is sampling from a Wishart distribution using the Bartlett decomposition.

There was some discussion in #1374, where the best alternative was to construct the type with a 1 by 1 array, and then replace factor:

C = chold!(eye(1), true)
C.LR = U

However this seems suboptimal, as it still requires calling the lapack code, the result of which is discarded.

One solution might be to remove the computation from the constructor, and move it to a convert call. This may be relevant to the discussion of convert and the idiom of constructors in #1470.

[ViralBShah]

How about having a constructor along the lines of CholeskyDense(M::Matrix, ischol::Bool)?

[dmbates]

@ViralBShah I think that method will conflict with the current constructor where the Bool value is upper.

@simonbyrne Calling Lapack code for the Cholesky decomposition of eye(1) only requires evaluating sqrt(1) and I can't imagine that will make a noticeable impact on execution times. However, I do agree that the method described is inelegant.

In writing things the way I did, I was trying to ensure that the eventual call to *potrs was controlled through a single function and all the other constructors went through it. I guess it is a matter of whether one considers a factorization object as being the factorization of a matrix of a certain type or just something that behaves like a factorization.

If you are particularly interested in the Bartlett decomposition why not create a Bartlett decomposition type? Tearing apart the logic of the Cholesky decomposition types to accommodate this one case seems a bit drastic to me.

[simonbyrne]

In writing things the way I did, I was trying to ensure that the eventual call to *potrs was controlled through a single function and all the other constructors went through it.

I agree that it makes sense to keep the lapack call in one place, it is just that the constructor may not be the best place for it. It seems that there is still some debate about the distinction between constructors and convert methods, and I thought that this may be relevant.

I guess ultimately it comes down to whether one would ever want to construct a factorization (without creating a new type) by a method other than a call to the lapack functions, and this was one possible example.

[timholy]

Perhaps one solution would be a separate constructor with syntax C = CholeskyDense(Float64). This would conceptually be like calling it with an empty matrix (which currently gives an error).

However, I'm not entirely convinced that this issue is a problem: if you have the separate factors, can't you just use them directly? Is there a circumstance where bundling them into a factorization object makes your life significantly easier? That said, I'm not opposed to having this, since I don't think this would break anything or prevent future extensions.

[simonbyrne]

Is there a circumstance where bundling them into a factorization object makes your life significantly easier?

Well the main use is that it would allow the use of the methods of the factorization, without having to reimplement them. For example, to get a sample from an inverse Wishart, you can just call inv, rather than Lapack.potri!.

Admittedly, this is not a huge problem (one of the fantastic things about julia is that calling Lapack.potri! is actually quite easy), but I do think that consistency with other constructors (such as Tridiagonal) would be useful.

[toivoh]

I agree with Simon. If you bother to make the Cholesky factorization easy to use, you might as well make it easy to use in this case too.

[simonbyrne]

Resolved by f6bfc32