Extract ColVecs and RowVecs

[devmotion]

I couldn't find an existing issue so I open a new one here.

The ColVecs/RowVecs abstraction would be useful in other packages as well that otherwise don't want or need to depend on KernelFunctions. Thus I think we should extract it into a separate package. Maybe it could be named VecOfVecs or VectorsOfVectors? The latter might be too similar to https://github.com/JuliaArrays/ArraysOfArrays.jl though.

[davibarreira]

Indeed, it would be great to have this separately! If time is an issue, I might try contributing. By what I've seen from the code, it looks like the implementation is already quite compartmentalized.

[theogf]

Maybe be more explicit about what it does, like MatrixAsVecs or something like this?

[davibarreira]

Just my two cents about the name. I think VectorOfVectors or VecOfVecs is more aligned with how I think of it. I mean, if I have a matrix [1 2; 3 4] and someone say, "turn this into a vector of vectors" I would naturally write something like [[1,2],[3,4]]. But if you told me "turn this matrix into a vector", I'd do something like [1,2,3,4].

[devmotion]

I think MatrixAsVecs or MatrixAsVectors is a good name. It says what ColVecs and RowVecs do: take a matrix and view as a vector of vectors.

[theogf]

Should we go full explicit with MatrixAsVecOfVecs ? I think that this package is anyway always going to be used internally

[devmotion]

I'm not completely sure if this will be the case. If it is loaded only internally, then packages would be forced to reexport ColVecs and RowVecs it seems.

[theogf]

Yes, but that seems like a fair thing to do. We already export ColVecs/RowVecs after all and I don't think it should change.

[devmotion]

This would be different though for packages that only work with AbstractVector{<:AbstractVector} but not ColVecs or RowVecs explicitly, such as eg CalibrationErrors.

Nevertheless, I think this wouldn't be a strong argument for a shorter name - you can always use a short alias (in Julia >= 1.6 one can eg also use import MatrixAsVectors as MV).

[theogf]

@devmotion Since you mentioned ArraysOfArrays I had a look and they have nestedview which does the exact same thing as ColVecs

julia> using ArraysOfArrays
julia> A = rand(3, 5)
julia> x = nestedview(A)
5-element ArrayOfSimilarArrays{Float64, 1, 1, 2, Matrix{Float64}}:
 [0.2459596114320628, 0.1105183659893263, 0.11301322795465385]
 [0.3180768545993806, 0.04457502255219414, 0.9357462783278538]
 [0.42793803145344356, 0.3558943294674095, 0.8159156651542603]
 [0.12633273140535373, 0.9693234579440309, 0.13627099668977305]
 [0.7011049776409946, 0.16946481398138058, 0.9546857384793925]
julia> x isa AbstractVector{<:AbstractVector}
true

However they do not have the option to use rows instead of columns (but that can be solved with nestedview(A')
One advantage is that they implement a lot of practical operations already (pop!, push! etc)
The additional dependencies would be:

Adapt
Requires
Statistics
UnsafeArrays

[davibarreira]

Are there performance issues when transposing the matrix or is this as efficient as RowVecs? I guess that if it's less efficient, submitting a PR to ArraysOfArrays.jl might be the faster solution.

[theogf]

using BenchmarkTools, KernelFunctions, ArraysOfArrays
A = rand(100, 10000)
B = rand(10000, 100)
x_col = ColVecs(A)
x_arr = nestedview(A)
f(x) = sum(x->sin.(x), x) # I added sin because ArofAr optimizes sum
@btime f($x_col)
7.353 ms (19999 allocations: 17.09 MiB)
@btime f($x_arr)
 7.599 ms (19999 allocations: 17.09 MiB)
x_row = RowVecs(B)
x_arr2 = nestedview(B')
@btime f($x_row)
  8.746 ms (19999 allocations: 17.09 MiB)
@btime f($x_arr2)
  8.479 ms (19999 allocations: 17.09 MiB)

[willtebbutt]

FWIW, I'd also want to ensure that AD works properly with ArraysOfArrays if we're considering adopting it. It looks to me like the inner constructors do quite a bit, so I'd want to be sure that they interoperate nicely with Zygote et al if we're giving up control of the types.

[theogf]

Kernel matrices benchmarks (which are unfair because we specialize on ColVecs and RowVecs)

@btime kernelmatrix(SqExponentialKernel(), $x_col)
  1.660 s (6 allocations: 1.49 GiB)
@btime kernelmatrix(SqExponentialKernel(), $x_arr)
  3.548 s (4 allocations: 1.49 GiB)
@btime kernelmatrix(SqExponentialKernel(), $x_row)
  1.679 s (8 allocations: 1.50 GiB)
@btime kernelmatrix(SqExponentialKernel(), $x_arr2)
  12.052 s (4 allocations: 1.49 GiB)

For AD

using Zygote, ForwardDiff, ReverseDiff
Zygote.gradient(X) do x
       sum(kernelmatrix(SqExponentialKernel(), nestedview(x)))
end
ERROR: 
Need an adjoint for constructor ArrayOfSimilarArrays{Float64, 1, 1, 2, Matrix{Float64}}. Gradient is of type Vector{Vector{Float64}}
ForwardDiff.gradient(X) do x
       sum(kernelmatrix(SqExponentialKernel(), nestedview(x)))
end
# Passes
ReverseDiff.gradient(X) do x
       sum(kernelmatrix(SqExponentialKernel(), nestedview(x)))
end
# Passes

[EDIT]
But! If I do :

using Distances
KernelFunctions.pairwise(d::SqEuclidean, x::VectorOfSimilarVectors) = Distances.pairwise(d, x.data, dims=2)
Zygote.gradient(X) do x
       sum(kernelmatrix(SqExponentialKernel(), nestedview(x)))
end
# This works!

And of course

@btime kernelmatrix(SqExponentialKernel(), $x_arr)
  1.700 s (6 allocations: 1.49 GiB)
@btime kernelmatrix(SqExponentialKernel(), $x_arr2)
  1.683 s (6 allocations: 1.49 GiB)

[willtebbutt]

Hmm probably the right way to fix this is to specialise some of our operations for the matrices in question so that a Vector{Vector{Float64}} cotangent isn't produced. (we ought to produce a Tangent instead).

[theogf]

Hmm probably the right way to fix this is to specialise some of our operations for the matrices in question so that a Vector{Vector{Float64}} cotangent isn't produced. (we ought to produce a Tangent instead).

@willtebbutt See edited reply above

[devmotion]

Since you mentioned ArraysOfArrays I had a look and they have nestedview which does the exact same thing as ColVecs

I know that it is similar, this is why I mentioned it 😄

As you already mentioned, it does not support RowVecs though.

Another major difference is that in contrast to ColVecs and RowVecs it is always a AbstractArray{<:Array}, i.e., the elements of the resulting array are Arrays and not general AbstractVectors (usually views) in our implementations.

[theogf]

Another major difference is that in contrast to ColVecs and RowVecs it is always a AbstractArray{<:Array}, i.e., the elements of the resulting array are Arrays and not general AbstractVectors (usually views) in our implementations.

Not totally exact (if I understood your sentence correctly) but weird nonetheless

julia> typeof(x)
ArrayOfSimilarArrays{Float64, 1, 1, 2, Matrix{Float64}}

julia> eltype(x)
Vector{Float64} (alias for Array{Float64, 1})

julia> typeof(x[1])
SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}

julia> x[1]
10-element view(::Matrix{Float64}, :, 1) with eltype Float64:
 0.3393166825154832
...

[devmotion]

One advantage is that they implement a lot of practical operations already (pop!, push! etc)

pop!, push!, append! are not supported for general arrays since you can't perform these operations with e.g. Matrix. You have to use eg. ElasticArrays as underlying data.

[theogf]

Oh you're right, it does not work. Since I saw this : https://github.com/JuliaArrays/ArraysOfArrays.jl/blob/a45514c4ab78a42b617972b7154319c87e738751/src/array_of_similar_arrays.jl#L306, I thought it would work out of the box

[devmotion]

In general, I assume at some point we can just use eachrow and eachcol and dispatch on something like Rows and Columns (JuliaLang/julia#32310). The PR tries to solve the same problem as RowVecs and ColVecs in base. It's an old PR though and hence I assumed we might not want to wait for it - but maybe it will be merged reasonably soon and we don't need a package for ColVecs and RowVecs anymore? 🤷‍♂️ Support of older Julia versions (not eachcol/eachrow but the Columns etc. types) would be added in Compat: JuliaLang/Compat.jl#663

[devmotion]

FYI JuliaLang/julia#32310 was merged 🙂 There's an open PR to Compat.jl, and when it is merged we could in principle start using RowSlices/ColumnSlices I think (with the eachcol/eachrow syntax on newer Julia versions).