Does Julia recognize block, banded, and block-banded matrices?

Specific Domains Numerics
,
1

For a simple example of why I need this, let’s say I want to compute the entropy of a multivariate Gaussian. One can do this as follows:

\text{Entropy} \propto -\frac{1}{2}*\log|-H|

where -H is the negative Hessian at the maximum of the Gaussian distribution. I.e., we need to calculate the log-determinant of the negative Hessian. One can calculate this in a stable way using the Cholesky decomposition, and in theory this can be done in O(T) time but only if you exploit the structure of the matrix (in my case this would be a banded matrix).

2

See BandedMatrices.jl

3

BandedMatrices.jl, BlockArrays.jl, BlockBandedMatrices.jl. Also LazyBandedMatrices.jl is available.

4

Do you know if these packages are compatible with other LinearAlgebra routines? It seems Banded.jl has some methods that override things like the \ operator, but for something like:

cholesky(my_banded_matrix)

or some other general method. I guess in the same way using Hermitian(some_matrix) from LinearAlgebra.jl provides speed-ups if the specific structure can be exploited?

5

Using the interfaces from ArrayLayouts.jl, yes they should e.g.

julia> using BandedMatrices, LinearAlgebra, ArrayLayouts

julia> A = Symmetric(brand(10,10,2,2) + 10I);

julia> cholesky(A) # using a banded chol
Cholesky{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
U factor:
10×10 UpperTriangular{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 3.17267  0.105034  0.144806    ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅       3.19379   0.0523953  0.251397   ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅        3.26506    0.151593  0.0750854   ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅         3.26117   0.0102327  0.0144174   ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅        3.16206    0.181666   0.107728    ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅         3.21687    0.0727688  0.117991    ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅         3.23698    0.0547971  0.191308    ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅         3.19253    0.0834566  0.278381
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅          ⋅         3.16049    0.229518
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅          ⋅          ⋅         3.21334

julia> ArrayLayouts.cholesky_layout(MemoryLayout(A), axes(A), A)
Cholesky{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
U factor:
10×10 UpperTriangular{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 3.17267  0.105034  0.144806    ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅       3.19379   0.0523953  0.251397   ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅        3.26506    0.151593  0.0750854   ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅         3.26117   0.0102327  0.0144174   ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅        3.16206    0.181666   0.107728    ⋅          ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅         3.21687    0.0727688  0.117991    ⋅          ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅         3.23698    0.0547971  0.191308    ⋅
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅         3.19253    0.0834566  0.278381
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅          ⋅         3.16049    0.229518
  ⋅        ⋅         ⋅          ⋅         ⋅          ⋅          ⋅          ⋅          ⋅         3.21334

julia> @which ArrayLayouts.cholesky_layout(MemoryLayout(A), axes(A), A)
cholesky_layout(::Union{HermitianLayout{<:AbstractBandedLayout}, SymmetricLayout{<:AbstractBandedLayout}}, ax, A; ...)
     @ BandedMatrices C:\Users\danjv\.julia\packages\BandedMatrices\WmDeb\src\symbanded\BandedCholesky.jl:96

Same goes for the other packages.

1 Like
6

Oh great thanks! Wasn’t aware of the existence of ArrayLayouts.jl. Will give it a go.

7

(Note that you don’t need to interact with ArrayLayouts.jl directly, just showing that the methods you don’t see overloaded are because of this interface, so cholesky(A) will just work with the most efficient method.)

1 Like