hvcat?

[chriscoey]

I have a big block matrix with multiple rows and columns of heterogeneous blocks (sparse, dense, structured, LinearMaps, UniformScaling...). It is too large to fit in memory and I don't want to convert any of the blocks. I want to be able to use it with IterativeSolvers - will that one day be possible with LazyArrays?

[tkf]

You can use BlockArrays.mortar to do something similar

julia> B = mortar(reshape([
           Ones(2, 2), Zeros(2, 2),
           spzeros(2, 2), ones(2, 2),
       ], 2, 2))
4×4 BlockArray{Float64,2,AbstractArray{Float64,2}}:
 1.0  1.0  │  0.0  0.0
 1.0  1.0  │  0.0  0.0
 ──────────┼──────────
 0.0  0.0  │  1.0  1.0
 0.0  0.0  │  1.0  1.0

julia> B.blocks[1, 1]
2×2 Ones{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}:
 1.0  1.0
 1.0  1.0

julia> B.blocks[2, 1]
2×2 Zeros{Float64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}:
 0.0  0.0
 0.0  0.0

julia> B.blocks[1, 2]
2×2 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> B.blocks[2, 2]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

There will be an indirection due to AbstractArray but I guess for such large size problem this does not matter. I'm not sure if it works with LinearMaps, though.

[dlfivefifty]

Hvcat would be great to have, any PR adding it would be greatly appreciated.

[tkf]

@dlfivefifty Can we actually get rid of Vcat and Hcat by allowing StaticArray in BlockArray's blocks? This way, I guess we get the full lazy cat for "free"?

[dlfivefifty]

I think a better approach to avoid contaminating packages is to replace them with ApplyArray(hcat, ...). I wonder if this already works for hvcat?

[tkf]

Ah, yes, a good point. Nonhomogeneous case needs different treatment usually.

ApplyArray(hcat, ...)

I think you'd want a special container for matrix case to efficiently use, e.g., Cij += Aij * Bij for each block ij to do C = A * B for the whole matrix. Using @~ hcat(vcat(...), vcat(...), ...) or even @~ [A B; C D] as a surface syntax sounds great but I think dispatching on a complex tree structure like this would be a nightmare. So maybe add instantiate phase like Broadcasted does? The idea is that lazy hcat-of-vcats would be internally converted to an object that supports efficient mul-add etc.

[chriscoey]

♥️ everything @tkf just said. Ideally no allocations for mul! with hvcat

[tkf]

Here is an idea about supporting @~ [A B; C D]: https://discourse.julialang.org/t/surface-ast-to-surface-ast-lowering-desugaring/28067

[dlfivefifty]

Just to clarify my comment above: I wasn't proposing combining hcat and vcat, but rather support lazy hvcat via:

julia> A = randn(2,2);

julia> hvcat((2,2),A,A,A,A) # equivalent to [A A; A A]
4×4 Array{Float64,2}:
 0.232976  -0.0681057  0.232976  -0.0681057
 0.808114  -0.674706   0.808114  -0.674706 
 0.232976  -0.0681057  0.232976  -0.0681057
 0.808114  -0.674706   0.808114  -0.674706 

julia> ApplyArray(hvcat, (2,2), A,A,A,A) # should match above
ERROR: MethodError: no method matching ndims(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(hvcat),Tuple{Tuple{Int64,Int64},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}})
Closest candidates are:
  ndims(::Base.Generator) at generator.jl:53
  ndims(::Number) at number.jl:67
  ndims(::Type{#s75} where #s75<:Number) at number.jl:68
  ...
Stacktrace:
 [1] ApplyArray{Any,N,F,Args} where Args<:Tuple where F where N(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(hvcat),Tuple{Tuple{Int64,Int64},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}}) at /Users/solver/Projects/LazyArrays.jl/src/lazyapplying.jl:127
 [2] ApplyArray(::LazyArrays.Applied{LazyArrays.DefaultApplyStyle,typeof(hvcat),Tuple{Tuple{Int64,Int64},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}}) at /Users/solver/Projects/LazyArrays.jl/src/lazyapplying.jl:128
 [3] ApplyArray(::Function, ::Tuple{Int64,Int64}, ::Vararg{Any,N} where N) at /Users/solver/Projects/LazyArrays.jl/src/lazyapplying.jl:132
 [4] top-level scope at REPL[9]:1

All we have to do is add size and ndims and the basic functionality works:

julia> Base.ndims(::LazyArrays.Applied{<:Any,typeof(hvcat)}) = 2

julia> Base.size(::LazyArrays.Applied{<:Any,typeof(hvcat)}) = (4,4)

julia> ApplyArray(hvcat, (2,2), A,A,A,A)
4×4 ApplyArray{Any,2,typeof(hvcat),Tuple{Tuple{Int64,Int64},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}}:
 0.232976  -0.0681057  0.232976  -0.0681057
 0.808114  -0.674706   0.808114  -0.674706 
 0.232976  -0.0681057  0.232976  -0.0681057
 0.808114  -0.674706   0.808114  -0.674706 

Of course, this is using generic fall back so will be slow until various operations like materialize(::MulAdd{HvcatLayout}) are overloaded.

[dlfivefifty]

I could in theory add generic fallbacks ndims(A::Applied) = ndims(materialize(A)) and size(A::Applied) = size(materialize(A))...

[tkf]

I wasn't proposing combining hcat and vcat

Ah, sorry hcat-of-vcat thing was due to my misunderstanding of how [A B; C D] work.

[dlfivefifty]

I think this has sort of been taken care of.