JuliaStats / PDMats.jl

Uniform Interface for positive definite matrices of various structures
Other
104 stars 43 forks source link

Avoid copying when accessing cholesky factors of PDMat #133

Closed sethaxen closed 3 years ago

sethaxen commented 3 years ago

The methods specialized on PDMat frequently access the U and L properties of the stored cholesky factor chol. Cholesky has an uplo parameter that determines whether it stores the upper or lower Cholesky factor. If chol.uplo=='U', then chol.U gives an UpperTriangular view of the stored factor, while chol.L allocates a transposed copy of this matrix. The opposite is true for chol.L and uplo == 'L'. To avoid unnecessary copies, the PDMat methods should check the uplo to decide whether to call, e.g. chol.U or chol.L'

sethaxen commented 3 years ago

As an example, consider the following two implementations of invquad:

# what PDMats currently uses
invquad(a::PDMat, x::AbstractVector) = sum(abs2, a.chol.L \ x)

# what I propose
function invquad(a::PDMat, x::AbstractVector)
    chol = a.chol
    return sum(abs2, (chol.uplo === 'L' ? chol.L : transpose(chol.U)) \ x)
end

Now let's benchmark for a 100x100 matrix:

julia> A = PDMat(exp(Symmetric(randn(100, 100))));

julia> x = randn(100);

julia> @btime invquad($A, $x)

The results are

  17.435 μs (4 allocations: 79.09 KiB)  # current version
  2.441 μs (3 allocations: 928 bytes)  # new version