# Blockdiag for dense matrices

**URL:** <https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941>\
**Category:** Specific Domains\
**Tags:** linearalgebra\
**Created:** [January 29, 2020, 4:54pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941 "2020-01-29T16:54:32Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [January 29, 2020, 4:54pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/1 "2020-01-29T16:54:32Z")

</div>

Quite often I want to form dense block diagonal matrices from a few dense matrices. Is there an easy way to do this? There is `blockdiag` for sparse matrices:

```julia
help?> blockdiag
search: blockdiag

  blockdiag(A...)

  Concatenate matrices block-diagonally. Currently only implemented for sparse matrices.

```

Is there some reason that this hasn’t been implemented for dense matrices or is there some other way to achieve this? (preferably without using additional packages)

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [January 29, 2020, 5:43pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/2 "2020-01-29T17:43:08Z")

</div>

Apply `Matrix` to the result of `blockdiag`?

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [January 29, 2020, 6:38pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/3 "2020-01-29T18:38:20Z")

</div>

`blockdiag` only accepts sparse matrices as inputs. Converting the dense matrices to sparse matrices and then converting the result to a dense matrix seems inefficient and inconvenient.

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [January 29, 2020, 6:43pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/4 "2020-01-29T18:43:09Z")

</div>

Perhaps methods in the [`BlockDiagonals` package](https://github.com/invenia/BlockDiagonals.jl) would be more suitable?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [January 29, 2020, 11:16pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/5 "2020-01-29T23:16:07Z")

</div>

A simple test:

```julia
julia> v = [rand(2,2), rand(3,3)];

julia> @benchmark BlockDiagonal(v)
BenchmarkTools.Trial:
  memory estimate: 16 bytes
  allocs estimate: 1
  --------------
  minimum time: 84.702 ns (0.00% GC)
  median time: 91.884 ns (0.00% GC)
  mean time: 110.563 ns (7.64% GC)
  maximum time: 56.262 μs (99.83% GC)
  --------------
  samples: 10000
  evals/sample: 961

julia> @benchmark begin
       Matrix(blockdiag(sparse.(v)...))
       end
BenchmarkTools.Trial:
  memory estimate: 2.66 KiB
  allocs estimate: 26
  --------------
  minimum time: 5.117 μs (0.00% GC)
  median time: 6.183 μs (0.00% GC)
  mean time: 7.641 μs (17.00% GC)
  maximum time: 7.730 ms (99.83% GC)
  --------------
  samples: 10000
  evals/sample: 6

```

**Conclusion** : creating the block diagonal matrix using package BlockDiagonals with method BlockDiagonal is roughly 60 times faster than converting between dense and sparse matrices – for this simple case.

If I convert the BlockDiagonal matrix to a regular Matrix, most of the advantage is lost:

```julia
julia> @benchmark collect(BlockDiagonal(v))
BenchmarkTools.Trial:
  memory estimate: 7.17 KiB
  allocs estimate: 77
  --------------
  minimum time: 3.275 μs (0.00% GC)
  median time: 4.275 μs (0.00% GC)
  mean time: 6.096 μs (25.25% GC)
  maximum time: 6.393 ms (99.90% GC)
  --------------
  samples: 10000
  evals/sample: 8

```

Finally, it is of interest to see how much slower multiplication with a BlockDiagonal matrix is:

```julia
julia> M_bl = BlockDiagonal(v);

julia> M = collect(M_bl);

julia> b = rand(5);

julia> @benchmark M_bl*b
BenchmarkTools.Trial:
  memory estimate: 624 bytes
  allocs estimate: 7
  --------------
  minimum time: 402.000 ns (0.00% GC)
  median time: 418.995 ns (0.00% GC)
  mean time: 642.634 ns (19.39% GC)
  maximum time: 261.817 μs (99.77% GC)
  --------------
  samples: 10000
  evals/sample: 200

julia> @benchmark M*b
BenchmarkTools.Trial:
  memory estimate: 128 bytes
  allocs estimate: 1
  --------------
  minimum time: 132.351 ns (0.00% GC)
  median time: 135.640 ns (0.00% GC)
  mean time: 179.146 ns (7.58% GC)
  maximum time: 69.557 μs (99.69% GC)
  --------------
  samples: 10000
  evals/sample: 881

```

**Conclusion** : matrix-vector multiplication when the matrix is BlockDiagonal is 3 times _slower_ than when the matrix is a regular Matrix – for this simple example.

---

<div class="post-metadata">

**Author:** ![oxinabox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oxinabox/32/206603_2.png) [@oxinabox](https://discourse.julialang.org/u/oxinabox)\
**Post date:** [January 30, 2020, 12:12am UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/6 "2020-01-30T00:12:47Z")

</div>

You benchmarked with nonconst globals.  
That will not give realistic results.  
You need to either make them `const` or interpolate them into the benchmark expression.

See docs: [https://docs.julialang.org/en/latest/manual/performance-tips/#Avoid-global-variables-1](https://docs.julialang.org/en/latest/manual/performance-tips/#Avoid-global-variables-1)  
BenchmarkTools Readme: [GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language](https://github.com/JuliaCI/BenchmarkTools.jl#quick-start)

* * *

also thats a pretty small matrix, its not surprising to see dense matrix multiplcation work better.  
StaticArrays would likely be even faster still at that size.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [January 30, 2020, 8:05pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/7 "2020-01-30T20:05:24Z")

</div>

Making functions of the various operations didn’t change much. Using `BlockDiagonal` to create the matrix with `v = [rand(2,2), rand(3,3)]` was ca. 50-60 times faster than if I created a sparse matrix, created the block diagonal, and converted back to dense matrix.

With `v = [rand(200,200), rand(300,300)]` instead, using `BlockDiagonal` was some 1e5 faster. Probably because the wrapping of `sparse`, `blockdiagonal`, and `Matrix` function creates temporary arrays. So I’d need to be more careful for a fair comparison. On the other hand, using `BlockDiagonal` is extremely simple, and seems to be well optimized.

Wrt. multiplying a matrix of type `BlockDiagonal` with a dense vector, this was some 3 times _slower_ for the small matrix (5,5) compared to multiplying a dense matrix with a dense vector. However, for the large matrix (500,500), the computational times were relatively similar.

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [January 30, 2020, 8:30pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/8 "2020-01-30T20:30:05Z")

</div>

Hmm, `BlockDiagonals` seems great for exploiting the structure of block diagonal matrices, but I am simply looking for a convenient way of forming a (dense, relatively small) block diagonal matrix. I was hoping this could be done without using yet another package. It seems like a quite common thing to do (at least for me) and that it would be nice if there was a function for this in `Base`? In particular since there is the corresponding function for sparse matrices.

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [January 30, 2020, 9:02pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/9 "2020-01-30T21:02:54Z")

</div>

How’s this?

```julia
julia> blocks = [rand(i,i) for i in 1:3];
julia> cat(blocks...,dims=(1,2))
6×6 Array{Float64,2}:
 0.448769 0.0 0.0 0.0 0.0 0.0      
 0.0 0.744482 0.225303 0.0 0.0 0.0      
 0.0 0.463063 0.332314 0.0 0.0 0.0      
 0.0 0.0 0.0 0.688196 0.263436 0.896587 
 0.0 0.0 0.0 0.694443 0.403921 0.0913202
 0.0 0.0 0.0 0.486398 0.299153 0.526531 

```

---

<div class="post-metadata">

**Author:** ![oxinabox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oxinabox/32/206603_2.png) [@oxinabox](https://discourse.julialang.org/u/oxinabox)\
**Post date:** [January 30, 2020, 9:07pm UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/10 "2020-01-30T21:07:25Z")

</div>

> [@olof3](#):
>
> It seems like a quite common thing to do (at least for me) and that it would be nice if there was a function for this in `Base` ?

Julia, like TeX, is a language in which packages are the standard way to get functionality.  
It’s largely a historical artifact that sparse matrixes have a standard library.

---

<div class="post-metadata">

**Author:** ![olof3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olof3/32/23875_2.png) [@olof3](https://discourse.julialang.org/u/olof3)\
**Post date:** [February 3, 2020, 9:51am UTC](https://discourse.julialang.org/t/blockdiag-for-dense-matrices/33941/11 "2020-02-03T09:51:42Z")

</div>

@tobydriscoll, nice, yes it was something like that I was looking for, that looks as good as it gets without a dedicated method.

Actually stumbled upon a dedicated method. It seems to be a bit faster

```julia-repl
julia> blocks = [rand(i,i) for i in 1:3];
julia> @btime ControlSystems.blockdiag(blocks...)
  275.169 ns (3 allocations: 592 bytes)
julia> @btime cat(blocks...,dims=(1,2))
  9.080 μs (52 allocations: 2.31 KiB)

```

`blockdiag` is quite analogous to `LinearAlgebra.diagm` which is also a convenience function for creating dense matrices, so `LinearAlgebra` would be the appropriate place. But I take it that there are not so many others that have found a need for this.