Directsum

[martincornejo]

Implements directsum (https://en.wikipedia.org/wiki/Matrix_addition#Direct_sum) and unicode operator ⊕. Currently it is limited to matrices but it could be extended to support multiple dimentsions (similar to the current behavior between Vector ⊕ Matrix)

[mcabbott]

For vectors, the direct sum is vcat(x, y) instead of cat(A, B; dims=(1,2)) for matrices. Is there a clear generalisation to all dimensions?

The present function seems to treat any mix of vectors & matrices as if they were matrices:

julia> A = [1 2 3; 4 5 6]; B = [7 8 9; 0 10 20];

julia> directsum([1,2,3], [77, 88, 99])
6×2 Matrix{Int64}:
 1   0
 2   0
 3   0
 0  77
 0  88
 0  99

julia> directsum(A, [77, 88, 99])
5×4 Matrix{Int64}:
 1  2  3   0
 4  5  6   0
 0  0  0  77
 0  0  0  88
 0  0  0  99

julia> @btime directsum($A, $B);  # this PR
  min 1.246 μs, mean 1.329 μs (11 allocations, 640 bytes)

julia> @btime cat($A, $B; dims=Val((1,2)));
  min 477.990 ns, mean 505.022 ns (15 allocations, 672 bytes)

[martincornejo]

For vectors, the direct sum is vcat(x, y) instead of cat(A, B; dims=(1,2))

Ah, you're right!

The present function seems to treat any mix of vectors & matrices as if they were matrices

Should mixing matrix and vector be undefined? And should in that case the user explicitly then provide a 1-column matrix instead?

julia> @btime cat($A, $B; dims=Val((1,2)));

Ah, that is a neat trick I just learned! cat over multiple dimensions 😊

Is there a clear generalisation to all dimensions?

To be honest, I'm not sure. But I remember reading somewhere (hopefully not taking it out of context), that it should follow $\dim A \oplus B = \dim A + \dim B$

In that sense, it can be generalized.

julia> AA = [A;;;A]; B = [B;;;B];

julia> size(cat(AA, BB, dims=(1,2,3))) == size(AA) .+ size(BB)
true

julia> size(cat(AA, BB, dims=(2,3))) == size(AA) .+ size(BB)
false

Maybe something like this

function directsum2(A::AbstractArray{T,N}, B::AbstractArray{S,M}) where {T,N,S,M}
    maxdim = max(N, M)
    return cat(A, B, dims=1:maxdim)
end

julia> directsum2(rand(2), rand(2))
4-element Vector{Float64}:
 0.7422427057509033
 0.8062369093673359
 0.484623525272175
 0.856149610025547

julia> directsum2(rand(2,2), rand(2,2))
4×4 Matrix{Float64}:
 0.247132  0.586941  0.0       0.0
 0.113415  0.250725  0.0       0.0
 0.0       0.0       0.369446  0.13091
 0.0       0.0       0.933767  0.595314

julia> @btime directsum2($A, $B);
  1.570 μs (28 allocations: 1.00 KiB)

But it seems there is a time penalty for generalizing)