Hi everyone!

I’ve been a Julia user for several years, but this is my first time posting here, so apologies if this has been asked somewhere before. I was curious about Julia’s implementation of Cholesky factorization — specifically accessing and using the factors of the Cholesky data structure.

Say I have some (positive definite) matrix A = L L^T, I perform the Cholesky factorization as

`C = cholesky(A)`

, and I want to use L for some calculations (e.g., a matrix-vector multiplication). I realized recently that, given some vector `v`

,

```
C.L*v
```

is quite a bit slower than

```
C.U'*v
```

because `C.L`

actually forms a copy of the underlying `factors`

matrix every time it’s called (and because storage of the factor in the upper part is the default behavior). The relevant code in `cholesky.jl`

is:

```
if d === :U
return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
elseif d === :L
return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors'))
elseif d === :UL
return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
else
```

I could imagine that many a user might naively call `C.L`

without realizing the significant performance hit, so I’m curious: why not return something like `LowerTriangular(Cfactors')`

instead of `LowerTriangular(copy(Cfactors'))`

(and similarly for `C.U`

if the factor is stored in the lower part)? Of course, there’s a difference in meaning by returning an adjoint, and maybe something would break (I can’t think of anything at the moment) for an `AbstractMatrix`

data structure that’s not just a simple array of floats. But for most use cases (e.g., multiplication, solves, changing coefficients, accessing coefficients, etc.), I would imagine this is advantageous. Are there some deep reasons why this choice of copying was made?

Thanks !