Does 3-arg `kernelmatrix_diag` require same length?

[simsurace]

The docs state that https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/blob/c97a246d9406765facceff9af1f65c38bc06a0ca/src/matrix/kernelmatrix.jl#L80-L83

However, testing this function for different-length output works for many kernels (did not try all though). So is the docstring wrong or is this an accidental feature? (incidentally, diag on a non-square matrix works, so the feature could in principle be useful)

[willtebbutt]

Hmm this definitely isn't an intended feature. Do you know what it computes when they're not the same length?

[simsurace]

Yes the same as what diag on the full matrix would give.

[simsurace]

Basically,

julia> x = rand(5); y = rand(4);

julia> kernelmatrix(SEKernel(), x, y)
5×4 Matrix{Float64}:
 0.913674  0.981759  0.925923  0.984204
 0.99906   0.848392  0.999942  0.979586
 0.998322  0.796401  0.995909  0.954717
 0.996092  0.869745  0.998438  0.987599
 0.997917  0.85858   0.999488  0.983593

julia> kernelmatrix_diag(SEKernel(), x, y)
4-element Vector{Float64}:
 0.9136737830751345
 0.8483918154068688
 0.9959090633450111
 0.9875987951145452

julia> kernelmatrix(SEKernel(), x, y) |> diag
4-element Vector{Float64}:
 0.9136737830751345
 0.8483918154068689
 0.9959090633450111
 0.9875987951145452

[simsurace]

553 showed that the unequal-length case errors for many kernel and input type combinations.