Isn't it very unintuitive that `cholesky()` stores data in the *upper* triangular part?

General Usage
1

A cholesky factor for a symmetric matrix A is usually described as a lower triangular matrix L such that A = LL'. But in Julia, the actual numerical value of L is stored in the upper triangular portion of cholesky(A).factors, that is, it is transposed.

Why?

This makes it awkward to actually access L. For example, how do I write an efficient half-vectorization routine? It feels like I need to throw around a bunch of transposes and UpperTriangular calls, and even then it’s still not clear how to access data in column major form.

2

You can ask for cholesky(A).L. More generally, Cholesky supports storing its factors either way — see the uplo field.

3 Likes
3

Thank you, but cholesky(A).L seems to be copying data internally, which is something I want to avoid.

I now see the uplo field, but have no idea how to actually change that. For example,

using LinearAlgebra
x = randn(5, 5);
sigma = x'*x;
L = cholesky(Symmetric(sigma));
L.uplo # returns 'U'

I guess L.uplo == 'U' means the data is actually stored in UpperTriangular(L.factors), but the documentation does not mention how to create a Cholesky type with L.uplo == 'L'.

The cholesky.jl file just seems a bit too complicated for me to just manually call the constructor of Cholesky <: Factorization directly.

4 
julia> C = cholesky(Hermitian([2 1; 1 2], :L))
Cholesky{Float64, Matrix{Float64}}
L factor:
2×2 LowerTriangular{Float64, Matrix{Float64}}:
 1.41421    ⋅ 
 0.707107  1.22474

julia> C.uplo
'L': ASCII/Unicode U+004C (category Lu: Letter, uppercase)
4 Likes