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

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.

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

3 Likes

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.

``````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