Comparing how to make diagonal block matrices

compare_block_diagonal.jl

using BlockDiagonals
using StatsBase
using BenchmarkTools
using Test
D = 2000
N = 100
K = 20
G = rand(D, N)
X = rand(D, N)
indices = vcat(1, sort(sample(2:D-1, K-1, replace = false)), D+1)
length(indices)
function dense_to_block(G, X, id)
    W = mean(eachcol(G) .* transpose.(eachcol(X)))
    Ws = [W[id[i]:id[i+1]-1, id[i]:id[i+1]-1] for i in 1:(length(id)-1)]
    return BlockDiagonal(Ws)
end



function direct_block(G, X, id)
    Ws = [mean(eachcol(G[id[i]:id[i+1]-1, :]) .* transpose.(eachcol(X[id[i]:id[i+1]-1, :]))) for i in 1:(length(id)-1)]
    return BlockDiagonal(Ws)
end

blockdiag = dense_to_block(G, X, indices)
blockdiag2 = direct_block(G, X, indices)
@test blockdiag ≈ blockdiag2

@btime dense_to_block($G, $X, $indices)
# 938.824 ms (443 allocations: 5.99 GiB)
@btime direct_block($G, $X, $indices)
# 73.298 ms (6765 allocations: 621.29 MiB)