diff --git a/src/vector_of_array.jl b/src/vector_of_array.jl
index 95777e05..4000f3db 100644
--- a/src/vector_of_array.jl
+++ b/src/vector_of_array.jl
@@ -1045,6 +1045,14 @@ function Base.append!(
 end
 
 function Base.stack(VA::AbstractVectorOfArray; dims = :)
+    # Skipping the per-element `stack` is only valid where it returns an identical
+    # copy: it turns `BitArray`s into `Array{Bool}`, narrows abstract eltypes,
+    # and throws on 0-d and empty arrays.
+    ET = eltype(VA.u)
+    if ET <: Array{<:Number} && isconcretetype(ET) && isconcretetype(eltype(ET)) &&
+            ndims(ET) >= 1 && !any(isempty, VA.u)
+        return stack(VA.u; dims)
+    end
     return stack(stack.(VA.u); dims)
 end
 
diff --git a/test/Core/interface_tests.jl b/test/Core/interface_tests.jl
index 59f45857..ebdffe58 100644
--- a/test/Core/interface_tests.jl
+++ b/test/Core/interface_tests.jl
@@ -208,6 +208,22 @@ testva = VA[VA[ones(2, 2), 2ones(2, 2)], 3ones(2, 2, 2)]
 @test stack(testva) ==
     [1.0 1.0; 1.0 1.0;;; 2.0 2.0; 2.0 2.0;;;; 3.0 3.0; 3.0 3.0;;; 3.0 3.0; 3.0 3.0]
 
+numeric_stack = VectorOfArray([rand(2) for _ in 1:1000])
+stack(numeric_stack) # compile away type instability before measuring
+@test (@allocated stack(numeric_stack)) < 20_000 # final 2x1000 Float64 array alone is 16000 bytes
+@test stack(numeric_stack) == reduce(hcat, numeric_stack.u)
+
+stack_inputs = (
+    [rand(3) for _ in 1:2], [rand(2, 2) for _ in 1:3], [[1, 2], [3, 4]],
+    Vector{Real}[[1, 2.0], [3, 4.0]], [trues(3), falses(3)], [trues(2, 2), falses(2, 2)],
+)
+for u in stack_inputs, dims in (:, 1, 2)
+    @test typeof(stack(VectorOfArray(u); dims)) == typeof(stack(stack.(u); dims))
+    @test stack(VectorOfArray(u); dims) == stack(stack.(u); dims)
+end
+@test_throws BoundsError stack(VectorOfArray([fill(1.0), fill(2.0)]))
+@test_throws ArgumentError stack(VectorOfArray([Float64[], Float64[]]))
+
 # convert array from VectorOfArray/DiffEqArray
 t = 1:8
 recs = [rand(10, 7) for i in 1:8]