Skip to content

Conversation

@DanielVandH
Copy link
Collaborator

@DanielVandH DanielVandH commented Nov 26, 2025

Now going to use this PR to get various other fixes for propagating CachedArrayStyle. Makes it easier for me when debugging but if it's preferred I can later break it up into smaller PRs. Once all the changes are in, I'll edit in all the new features here. Original PR text below.


This changes Vcat and Hcat to use the combined BroadcastStyle between a LazyArrayStyle and the combined style of its arguments. This fixes e.g. Vcat of cached arguments not being recognised as a CachedArrayStyle.

The implementation needs to use some recursion to get the BroadcastStyle of the arguments via tuple_type_broadcastlayout, but basically it's just doing the same as what you'd get from result_style(LazyArrayStyle{N}(), combine_style(A.args...)) where A is a Vcat/Hcat


@codecov
Copy link

codecov bot commented Nov 26, 2025

Codecov Report

❌ Patch coverage is 93.39623% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.76%. Comparing base (dc97ca8) to head (5c4ab23).

Files with missing lines Patch % Lines
src/cache.jl 91.93% 5 Missing ⚠️
src/lazyapplying.jl 88.88% 1 Missing ⚠️
src/lazybroadcasting.jl 83.33% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #391      +/-   ##
==========================================
- Coverage   95.94%   95.76%   -0.19%     
==========================================
  Files          17       17              
  Lines        3376     3467      +91     
==========================================
+ Hits         3239     3320      +81     
- Misses        137      147      +10     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@DanielVandH
Copy link
Collaborator Author

DanielVandH commented Nov 26, 2025

Similar cases missing a CachedArray style also show up in e.g. BlockArrays and LazyBandedMatrices:

julia> v = Accumulate(*, 1:10);

julia> (Base.BroadcastStyle  typeof).((v, unitblocks(v), BlockBroadcastArray(vcat, unitblocks(v), unitblocks(v))))
(LazyArrays.CachedArrayStyle{1}(), BlockArrays.BlockedStyle{1}(), LazyArrays.LazyArrayStyle{1}())

So I might have to eventually go and fix these, assuming it would be appropriate that these all return CachedArrayStyle{1}()? Similarly,

julia> Base.BroadcastStyle(typeof(vec(Vcat(v', v')))) # interlaced
Base.Broadcast.DefaultArrayStyle{1}()

will get a fix (this last one is now fixed in JuliaLinearAlgebra/ArrayLayouts.jl#268)

@DanielVandH DanielVandH changed the title Change BroadcastStyle of Vcat/Hcat to depend on BroadcastStyle of arguments Various fixes for propagating CachedArrayStyle correctly Nov 27, 2025
@DanielVandH DanielVandH marked this pull request as draft November 27, 2025 00:40
dlfivefifty pushed a commit to JuliaLinearAlgebra/ArrayLayouts.jl that referenced this pull request Nov 27, 2025
This allows e.g.

```julia-repl
julia> v = Accumulate(*, 1:10); Base.BroadcastStyle(typeof(vec(Vcat(v', v'))))
LazyArrays.LazyArrayStyle{2}()
```

(and, after JuliaArrays/LazyArrays.jl#391, this
would return `CachedArrayStyle{2}()`)
@dlfivefifty
Copy link
Member

I guess we don't need a CachedBandedArrayStyle because CachedStyle lowers to the data after resizing?

@DanielVandH
Copy link
Collaborator Author

I think that's right. I wasn't too sure exactly about cases like that (and e.g. the block matrices like I mentioned in my previous comment), but I think I reasoned similarly that CachedArrayStyle is fine

@DanielVandH
Copy link
Collaborator Author

DanielVandH commented Nov 27, 2025

This QRCompactWYQ type always seems to run into this issue. I'm not sure how to get around it properly this time, since there would need to be somewhere to call broadcastable earlier (which converts it into a Matrix), since by the time the calls get into using type information only it's too late. If LinearAlgebra just had

BroadcastStyle(::Type{<:LinearAlgebra.QRCompactWYQ}) = DefaultArrayStyle{2}()
BroadcastStyle(::Type{<:LinearAlgebra.AdjointQ}) = DefaultArrayStyle{2}()

it would be fine, but can't do that on our side without piracy. I don't really understand why they don't already do this.

Solutions like defining an internal _BroadcastStyle or overloading combine_styles seem like code smell. Have you had to deal with something similar before @dlfivefifty? An MWE for this is

julia> using LazyArrays, LinearAlgebra

julia> B = LazyArrays.CachedArray(randn(3, 3)); Q = qr(randn(3, 3)).Q; collect(B); Q * B  Q * B.data
ERROR: MethodError: no method matching ndims(::Type{LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}})
The function `ndims` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ndims(::Type{Union{}}, Any...)
   @ Base abstractarray.jl:276
  ndims(::Type{<:Ref})
   @ Base refpointer.jl:104
  ndims(::Type{<:Applied{<:Any, typeof(*), Args}}) where Args
   @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\linalg\mul.jl:41
  ...

Stacktrace:
  [1] Base.Broadcast.BroadcastStyle(::Type{LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}})
    @ Base.Broadcast .\broadcast.jl:103
  [2] tuple_type_broadcastlayout
    @ c:\Users\djv23\.julia\dev\LazyArrays.jl\src\lazyconcat.jl:445 [inlined]
  [3] applybroadcaststyle(::Type{ApplyArray{Float64, 2, typeof(*), Tuple{…}}}, ::ApplyLayout{typeof(*)})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\lazyapplying.jl:334
  [4] Base.Broadcast.BroadcastStyle(M::Type{ApplyArray{Float64, 2, typeof(*), Tuple{…}}})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\lazyapplying.jl:336
  [5] combine_styles(c::ApplyArray{Float64, 2, typeof(*), Tuple{LinearAlgebra.QRCompactWYQ{…}, CachedArray{…}}})
    @ Base.Broadcast .\broadcast.jl:433
  [6] combine_styles(c1::ApplyArray{Float64, 2, typeof(*), Tuple{…}}, c2::Matrix{Float64})
    @ Base.Broadcast .\broadcast.jl:438
  [7] broadcasted
    @ .\broadcast.jl:1353 [inlined]
  [8] broadcast_preserving_zero_d
    @ .\broadcast.jl:882 [inlined]
  [9] -(A::ApplyArray{Float64, 2, typeof(*), Tuple{LinearAlgebra.QRCompactWYQ{…}, CachedArray{…}}}, B::Matrix{Float64})
    @ Base .\arraymath.jl:8
 [10] isapprox(x::ApplyArray{…}, y::Matrix{…}; atol::Int64, rtol::Float64, nans::Bool, norm::typeof(norm))
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.12.1+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\generic.jl:1982
 [11] isapprox(x::ApplyArray{Float64, 2, typeof(*), Tuple{LinearAlgebra.QRCompactWYQ{…}, CachedArray{…}}}, y::Matrix{Float64})
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.12.1+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\generic.jl:1978
 [12] top-level scope
    @ REPL[100]:1
Some type information was truncated. Use `show(err)` to see complete types.

…/BoundsError issue with cacheddata storing a Vcat+scalar
@dlfivefifty
Copy link
Member

I think

defining an internal _BroadcastStyle

is probably the best way around this?

I have played around in the past with making a type so that Q's are <: AbstractMatrix (I think AbstractQ is kinda poorly thought out btw):

JuliaLinearAlgebra/MatrixFactorizations.jl#69

@DanielVandH
Copy link
Collaborator Author

I'm running into a bit of an annoying problem. The problem stems from trying to support cacheddata(::ApplyArray) and cacheddata(::BroadcastArray), simply rewrapping the arguments in maybe_cacheddata calls back into another (Apply/Broadcast)Array. This is one of the last things we might need for our timing stuff since we're losing a lot of the CachedArrayStyle to ApplyArray/BroadcastArrays currently.

First I'll just explain a new function conforming_resize!, then some problems with it.

To properly get the CachedArrayStyle to propagate correctly, I needed to add support for cacheddata on ApplyArray and BroadcastArray, which creates an ApplyArray/BroadcastArray over the same arguments except with their cached data. This will also have to make sure all the cacheddata calls on the arguments of those arrays match all the other arguments (since e.g., BroadcastVector(+, [1, 2], [1, 2, 3]) as an output from a cacheddata when the arguments have differerent datasizes wouldn't make sense). The logic for this is implemented (which, admittedly, looks quite rough for the moment) in conforming_resize! here

LazyArrays.jl/src/cache.jl

Lines 622 to 701 in 78b3c33

##
# Tuples of potential cached arrays
# Assumes that the Tuples contain arrays of the same shape
# Some of the complication here is in making sure that we can handle things like mixing vectors and matrices that have one column
##
@inline _tuple_wrap(x::Int) = (x,) # Some CachedArrays mistakenly use Int instead of Tuple for vector sizes (e.g. https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl/blob/a29007f2815134180b8433fdab46a23acfcdcd01/src/neg1c.jl#L19 and https://github.com/DanielVandH/InfiniteRandomArrays.jl/blob/a859d565f5bd8278d3f3499cd26b62916b4867e0/src/vector.jl#L18, to name a couple)
@inline _tuple_wrap(x::Tuple) = x
@inline _datasize(x) = size(x)
@inline _datasize(::Number) = (1,) # Numbers have size (), but we treat them as size (1,) for resizing purposes
@inline _datasize(x::AbstractCachedArray) = x.datasize
@inline _datasize(x::SubArray) = _datasize(parent(x))
@inline _datasize(x::AdjOrTrans) = reverse(_datasize(parent(x)))
@inline _datasize(x::AdjOrTrans{<:Any, <:AbstractVector}) = (sz = only(_datasize(parent(x))); (min(1, sz), sz))
@inline _datasizes(::Tuple{}) = ()
@inline _datasizes(t::Tuple) = (_datasize(first(t)), _datasizes(Base.tail(t))...) # equivalent to _datasize.(t)
@inline _has_vector(::Tuple{}) = false # Numbers have size ()
@inline _has_vector(::Tuple{Any}) = true
@inline _has_vector(::Tuple{Any, Vararg}) = false
@inline has_vector(::Tuple{}) = false
@inline has_vector(t::Tuple) = _has_vector(first(t)) || has_vector(Base.tail(t))
@inline to_vector_size(sz::Tuple{Any}) = sz # Already 1-tuple
@inline to_vector_size(sz::Tuple{Any, Vararg}) = (prod(sz),) # Multi-tuple -> 1-tuple
function max_datasize(sizes::Tuple)
if has_vector(sizes)
normalized = map(to_vector_size, sizes)
return reduce(@inline((a, b) -> max.(a, b)), normalized)
else
return reduce(@inline((a, b) -> max.(a, b)), sizes)
end
end
@inline _expand_size(sz::Tuple{Any}, ::Tuple{Any}) = sz
@inline _expand_size(sz::Tuple{Any}, ::Tuple{Any, Vararg}) = (only(sz), 1)
@inline _expand_size(sz::Tuple{Any, Vararg}, ::Tuple{Any, Vararg}) = sz
@inline _effective_ndims(sz::Int) = 1
@inline function _effective_ndims(sz::Tuple)
length(sz) == 1 && return 1
length(sz) == 2 && minimum(sz) <= 1 && return 1
return length(sz)
end
@inline _all_same_ndims(::Tuple{}) = true
@inline _all_same_ndims(t::Tuple{Any}) = true
@inline function _all_same_ndims(t::Tuple)
first_sz = _arraysize(first(t))
first_ndim = _effective_ndims(first_sz)
_check_same_ndims_args(first_ndim, Base.tail(t))
end
@inline _check_same_ndims_args(::Int, ::Tuple{}) = true
@inline function _check_same_ndims_args(expected_ndim::Int, t::Tuple)
sz = _arraysize(first(t))
_effective_ndims(sz) == expected_ndim && _check_same_ndims_args(expected_ndim, Base.tail(t))
end
@inline _arraysize(x::Number) = (1,)
@inline _arraysize(x) = size(x)
function conforming_resize!(args::Tuple)
isempty(args) && return args
if !_all_same_ndims(args)
throw(ArgumentError("Cannot conform arrays with incompatible dimensions: $(map(_arraysize, args))"))
end
sizes = _datasizes(args)
sz = max_datasize(sizes)
_resize_each!(args, sizes, sz)
return args
end
@inline _resize_each!(::Tuple{}, ::Tuple{}, sz) = nothing
@inline function _resize_each!(args::Tuple, sizes::Tuple, sz)
arg = first(args)
orig_sz = first(sizes)
target_sz = _expand_size(sz, orig_sz)
resizedata!(arg, target_sz...)
_resize_each!(Base.tail(args), Base.tail(sizes), sz)
end

and called e.g. for ApplyArray here

function cacheddata(A::ApplyArray{T, N, F}) where {T, N, F}
args = arguments(A)
conforming_resize!(args)
ApplyArray{T, N}(A.f, map(maybe_cacheddata, args)...)
end

The code for conforming_resize! is a bit complex because of issues I was having with e.g. recognising that cache(rand(1, 10))' is a vector but its size looks more like a matrix. You can look at the tests here

@testset "conforming_resize!" begin
args = (cache(1:10), cache(1:10));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((0,), (0,));
LazyArrays.resizedata!(args[2], 4);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((4,), (4,));
args = (1:10, cache(1:10));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((10,), (10,));
args = (cache(1:10)', cache(1:10)');
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((0, 0), (0, 0));
LazyArrays.resizedata!(args[1], 1, 3);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((1, 3), (1, 3));
args = (cache(1:10)', LazyArrays.CachedArray(rand(1, 10)));
@test LazyArrays._datasizes(args) == ((0, 0), (0, 0));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((0, 0), (0, 0));
LazyArrays.resizedata!(args[1], 1, 4);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((1, 4), (1, 4));
LazyArrays.resizedata!(args[2], 1, 6)
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((1, 6), (1, 6));
args = (cache(1:10), LazyArrays.CachedArray(rand(1, 10))');
LazyArrays.resizedata!(args[1], 4);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((4,), (4, 1));
args = (cache(1:10), LazyArrays.CachedArray(rand(1, 10))', 1:10);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((10,), (10, 1), (10,));
args = (1, cache(1:2));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((1,), (1,)) # because scalars are treated as size (1,)
args = (rand(10, 10), LazyArrays.CachedArray(rand(10, 10)));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((10, 10), (10, 10));
args = (LazyArrays.CachedArray(rand(10, 10)), LazyArrays.CachedArray(rand(10, 10)));
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((0, 0), (0, 0));
LazyArrays.resizedata!(args[1], 3, 4);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((3, 4), (3, 4));
LazyArrays.resizedata!(args[2], 5, 6);
LazyArrays.conforming_resize!(args);
@test LazyArrays._datasizes(args) == ((5, 6), (5, 6));
@testset "conforming_resize! dimension mismatch" begin
args = (cache(1:10), cache(1:10))
@test LazyArrays.conforming_resize!(args) === args
args = (cache(1:10), LazyArrays.CachedArray(rand(10, 5)))
@test_throws ArgumentError LazyArrays.conforming_resize!(args)
args = (LazyArrays.CachedArray(rand(3, 4)), LazyArrays.CachedArray(rand(5, 2)))
@test LazyArrays.conforming_resize!(args) === args
args = (1, cache(1:10))
@test LazyArrays.conforming_resize!(args) === args
args = (cache(1:10), LazyArrays.CachedArray(rand(2, 3)), reshape(cache(1:8), 2, 2, 2))
@test_throws ArgumentError LazyArrays.conforming_resize!(args)
args = (reshape(cache(1:8), 2, 2, 2), reshape(cache(1:27), 3, 3, 3))
@test LazyArrays.conforming_resize!(args) === args
args = ()
@test LazyArrays.conforming_resize!(args) === args
args = (cache(1:10),)
@test LazyArrays.conforming_resize!(args) === args
end
end
end

to try and see all the cases I was supporting so far if you did want to understand why it has to be so complex (not that there isn't probably a better way, which hopefully we can discuss).

Some problems:

  1. This approach doesn't really work in the case of e.g.
julia> A, b = rand(2, 2), rand(2);

julia> M = BroadcastMatrix(*, cache(A), cache(b));

julia> LazyArrays.cacheddata(M)
ERROR: ArgumentError: Cannot conform arrays with incompatible dimensions: ((2, 2), (2,))
Stacktrace:
 [1] conforming_resize!
   @ c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:698 [inlined]
 [2] cacheddata(A::BroadcastMatrix{Float64, typeof(*), Tuple{Matrix{Float64}, Vector{Float64}}})
   @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\lazybroadcasting.jl:127
 [3] top-level scope
   @ REPL[167]:1

and I'm not seeing the way to fix that currently.

  1. The current implementation tries to resize to an infinite arrays in the case that one of the arguments is infinite and not a cached array:
julia> A = Vcat((1:∞)', cache(1:∞)');

julia> B = view(A, :, 1:10);

julia> LazyArrays.cacheddata(B) 
ERROR: ArgumentError: Cannot create infinite Array
Stacktrace:
  [1] Vector{Int64}(::UndefInitializer, ::Tuple{Infinities.InfiniteCardinal{0}})
    @ InfiniteArrays C:\Users\djv23\.julia\packages\InfiniteArrays\v6xnl\src\infarrays.jl:19
  [2] Vector{Int64}(x::LazyArrays.CachedArray{Int64, 1, Vector{Int64}, Zeros{Int64, 1, Tuple{InfiniteArrays.OneToInf{…}}}})
    @ Base .\array.jl:621
  [3] convert(::Type{Vector{Int64}}, a::LazyArrays.CachedArray{Int64, 1, Vector{Int64}, Zeros{Int64, 1, Tuple{…}}})
    @ Base .\array.jl:614
  [4] setproperty!(x::LazyArrays.CachedArray{…}, f::Symbol, v::LazyArrays.CachedArray{…})
    @ Base .\Base_compiler.jl:57
  [5] _vec_resizedata!(B::LazyArrays.CachedArray{…}, n::Infinities.InfiniteCardinal{…})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:254
  [6] resizedata!(::DenseColumnMajor, ::LazyArrays.LazyLayout, B::LazyArrays.CachedArray{…}, n::Infinities.InfiniteCardinal{…})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:266
  [7] resizedata!(B::LazyArrays.CachedArray{…}, mn::Infinities.InfiniteCardinal{…})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:229
  [8] resizedata!(A::Adjoint{Int64, LazyArrays.CachedArray{…}}, m::Int64, n::Infinities.InfiniteCardinal{0})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:237
  [9] _resize_each!
    @ c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:711 [inlined]
 [10] _resize_each!
    @ c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:712 [inlined]
 [11] conforming_resize!(args::Tuple{Adjoint{…}, Adjoint{…}})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:702
 [12] cacheddata(A::ApplyArray{Int64, 2, typeof(vcat), Tuple{Adjoint{…}, Adjoint{…}}})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\lazyapplying.jl:233
 [13] cacheddata(V::SubArray{Int64, 2, ApplyArray{…}, Tuple{…}, false})
    @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\cache.jl:617
 [14] top-level scope
    @ REPL[172]:1
Some type information was truncated. Use `show(err)` to see complete types.

this example fails because it wants to create cacheddata(A) which in turn wants to expand the cacheddata for the second argument to infinity since the first argument is a non-cached argument with infinite length. This is the most difficult part for me.

There are a few ideas I have, like trying to reimplement cacheddata specifically for views of BroadcastArray/ApplyArray. This is all starting to feel very over-engineered however, so I wanted to see if you had any ideas @dlfivefifty before I add even more code that probably wouldn't hold up to even more intense debugging anyway.

…dispatch for _bc_resizecacheddata! rather than a MemoryLayout dispatch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants