Heterogeneous blocked matrices

Specific Domains Numerics
1

I am struggling to come up with a sufficiently specific method signature to handle the Cholesky decomposition of heterogeneous blocked matrices, by which I mean that I have a square matrix of AbstractMatrix{Float64}, say, but the lower-level matrices can be of different storage types: dense, sparse, diagonal, etc.

I can define the union of all the possible storage types

typealias Block{T} Union{Diagonal{T}, Diagonal{StaticArrays.MMatrix},
    Hermitian{T, Matrix{T}}, LowerTriangular{T, Matrix{T}}, Matrix{T}}

but I don’t think AbstractMatrix{Block{T}} is meaningful in a method signature. Or at least I can’t seem to get it to work.

Currently I am patterning my code after the (unreleased) LinearAlgebra package by @andreasnoack. I want to create a method for his cholBlocked! function that will apply to my blocked matrices or square views of them and is more specific than AbstractMatrix. I need to somehow characterize AbstractMatrix{AbstractMatrix{Float64}}.

The gory details are in https://github.com/dmbates/MixedModels.jl/tree/LowerCholesky Look for the string cholBlocked! in the src tree.

2

If I’m understanding your post correctly, I think you’re just butting up against parametric invariance. You just need to do f{T<:Block}(B::AbstractMatrix{T}). A concrete example here is Diagonal{Matrix{Float64}}:

julia> D = Diagonal([[1 2; 3 4], [5 6; 7 8]])
2×2 Diagonal{Array{Int64,2}}:
 [1 2; 3 4]      ⋅
     ⋅       [5 6; 7 8]

julia> isblock{T<:Block}(x::AbstractArray{T}) = true
       isblock(x) = false
isblock (generic function with 2 methods)

julia> isblock(D)
true
1 Like
3

Thanks for the reply, Matt, but I don’t think that will work.

I should have been more specific in my explanation. It is the blocks that are heterogeneous, not the wrapper. The array of blocks is dense, square and lower triangular. Blocks on the diagonal are square and, by design, decreasing in size. The largest block, in the [1,1] position, is always diagonal. The trailing diagonal blocks are dense but often their size is negligible compared to the [1,1] block.

If I understand correctly, a signature like

f{T <: Block}(x::AbstactArray{T}) end

applies to a homogeneous array of blocks of a single type, T, which is one of the types in the Block union.

For the time being I will forgo trying to emulate the elegant recursive cholBlocked! function in @andreasnoack package and instead use explicit loops.

4

Ah, I think I understand. It’ll depend upon how you construct the blocked matrix. If you just do something like this:

julia> A = reshape([Diagonal([1.,2.,3.,4.]), zeros(2,4), zeros(4,2), Diagonal([5.,6.])], 2, 2)
2×2 Array{AbstractArray{Float64,2},2}:
 [1.0 0.0 0.0 0.0; 0.0 2.0 0.0 0.0; 0.0 0.0 3.0 0.0; 0.0 0.0 0.0 4.0]  [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
 [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]                                    [5.0 0.0; 0.0 6.0]

Then the eltype is widened beyond your Union to simply be a AbstractArray{Float64,2}. Depending upon how you do it, the eltype might even get widened the whole way to simply Any. So instead of using Block, you something like isblock{T<:AbstractMatrix}(A::AbstractMatrix{T}) = true should do the trick.