diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 000000000..700707ced
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 6771b3044..7bb83b0d9 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -5,11 +5,30 @@ on:
     - cron: 0 0 * * *
   workflow_dispatch:
 
+permissions:
+  contents: write
+  pull-requests: write
+
 jobs:
   CompatHelper:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - name: Check if Julia is already available in the PATH
+        id: julia_in_path
+        run: which julia
+        continue-on-error: true
+      - name: Install Julia, but only if it is not already available in the PATH
+        uses: julia-actions/setup-julia@v1
+        with:
+          version: '1'
+          arch: ${{ runner.arch }}
+        if: steps.julia_in_path.outcome != 'success'
+      - name: "Add the General registry via Git"
+        run: |
+          import Pkg
+          ENV["JULIA_PKG_SERVER"] = ""
+          Pkg.Registry.add("General")
+        shell: julia --color=yes {0}
       - name: "Install CompatHelper"
         run: |
           import Pkg
diff --git a/.github/workflows/benchmark-comment.yml b/.github/workflows/benchmark-comment.yml
index 770fad567..0e9d192d2 100644
--- a/.github/workflows/benchmark-comment.yml
+++ b/.github/workflows/benchmark-comment.yml
@@ -24,7 +24,7 @@ jobs:
       ${{ github.event.workflow_run.event == 'pull_request' &&
           github.event.workflow_run.conclusion == 'success' }}
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       # restore records from the artifacts
       - uses: dawidd6/action-download-artifact@v2
         with:
@@ -41,7 +41,7 @@ jobs:
           echo ::set-output name=body::$(cat ./pull-request-number.artifact)
       # check if the previous comment exists
       - name: find comment
-        uses: peter-evans/find-comment@v1
+        uses: peter-evans/find-comment@v3
         id: fc
         with:
           issue-number: ${{ steps.output-pull-request-number.outputs.body }}
diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index aef638939..41ad36a17 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -16,11 +16,11 @@ jobs:
     if: contains(github.event.pull_request.labels.*.name, 'performance critical')
     steps:
       # setup
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
         with:
           version: '1.7'
-      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-buildpkg@v1
       - name: install dependencies
         run: julia -e 'using Pkg; pkg"add PkgBenchmark BenchmarkCI@0.1"'
       # run the benchmark suite
@@ -47,7 +47,7 @@ jobs:
       - name: record pull request number
         run: echo ${{ github.event.pull_request.number }} > ./pull-request-number.artifact
       # save as artifacts (performance tracking (comment) workflow will use it)
-      - uses: actions/upload-artifact@v2
+      - uses: actions/upload-artifact@v4
         with:
           name: Benchmarking
           path: ./*.artifact
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 9c9a4a0bd..628e5517d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -35,21 +35,12 @@ jobs:
           - 'MultiOutput'
           - 'Others'
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:
diff --git a/.github/workflows/doc_preview_cleanup.yml b/.github/workflows/doc_preview_cleanup.yml
index bc29462c0..94b9f70ad 100644
--- a/.github/workflows/doc_preview_cleanup.yml
+++ b/.github/workflows/doc_preview_cleanup.yml
@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v2
+        uses: actions/checkout@v4
         with:
           ref: gh-pages
       - name: Delete preview and history + push changes
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index a339654a2..90b31a14a 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -17,7 +17,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml
index f83b38f55..3f39a3fd3 100644
--- a/.github/workflows/format.yml
+++ b/.github/workflows/format.yml
@@ -13,8 +13,8 @@ jobs:
   format:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
         with:
           version: 1
       - run: |
diff --git a/Project.toml b/Project.toml
index 617f6e56c..24da5ddcc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.10.56"
+version = "0.10.63"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -25,13 +25,14 @@ ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 ChainRulesCore = "1"
 Compat = "3.7, 4"
 CompositionsBase = "0.1"
-Distances = "0.10"
+Distances = "0.10.9"
 FillArrays = "0.10, 0.11, 0.12, 0.13, 1"
 Functors = "0.1, 0.2, 0.3, 0.4"
 IrrationalConstants = "0.1, 0.2"
 LogExpFunctions = "0.2.1, 0.3"
 Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10, 1, 2"
+Statistics = "1"
 StatsBase = "0.32, 0.33, 0.34"
 TensorCore = "0.1"
 ZygoteRules = "0.2"
diff --git a/README.md b/README.md
index 837559bac..a70b77fcf 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # KernelFunctions.jl
 
-![CI](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/workflows/CI/badge.svg?branch=master)
+[![CI](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/actions/workflows/ci.yml/badge.svg?branch=master)](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/actions/workflows/ci.yml?query=branch%3Amaster)
 [![codecov](https://codecov.io/gh/JuliaGaussianProcesses/KernelFunctions.jl/branch/master/graph/badge.svg?token=rmDh3gb7hN)](https://codecov.io/gh/JuliaGaussianProcesses/KernelFunctions.jl)
 [![Documentation (stable)](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable)
 [![Documentation (latest)](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliagaussianprocesses.github.io/KernelFunctions.jl/dev)
diff --git a/docs/Project.toml b/docs/Project.toml
index a4792157b..264509c96 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,4 +12,5 @@ Documenter = "0.27"
 KernelFunctions = "0.10"
 Kronecker = "0.4, 0.5"
 PDMats = "0.11"
+Statistics = "1"
 julia = "1.3"
diff --git a/src/KernelFunctions.jl b/src/KernelFunctions.jl
index 63205b5bf..80711b4ec 100644
--- a/src/KernelFunctions.jl
+++ b/src/KernelFunctions.jl
@@ -48,7 +48,7 @@ export tensor, ⊗, compose
 
 using Compat
 using ChainRulesCore: ChainRulesCore, Tangent, ZeroTangent, NoTangent
-using ChainRulesCore: @thunk, InplaceableThunk
+using ChainRulesCore: @thunk, InplaceableThunk, ProjectTo, unthunk
 using CompositionsBase
 using Distances
 using FillArrays
diff --git a/src/basekernels/fbm.jl b/src/basekernels/fbm.jl
index 422400ac4..18dbdea14 100644
--- a/src/basekernels/fbm.jl
+++ b/src/basekernels/fbm.jl
@@ -57,7 +57,7 @@ end
 
 function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, x::AbstractVector)
     modx = _mod(x)
-    pairwise!(K, SqEuclidean(), x)
+    pairwise!(SqEuclidean(), K, x)
     K .= _fbm.(modx, modx', K, κ.h)
     return K
 end
@@ -72,7 +72,7 @@ end
 function kernelmatrix!(
     K::AbstractMatrix, κ::FBMKernel, x::AbstractVector, y::AbstractVector
 )
-    pairwise!(K, SqEuclidean(), x, y)
+    pairwise!(SqEuclidean(), K, x, y)
     K .= _fbm.(_mod(x), _mod(y)', K, κ.h)
     return K
 end
diff --git a/src/basekernels/matern.jl b/src/basekernels/matern.jl
index cb051ee0f..c009bf596 100644
--- a/src/basekernels/matern.jl
+++ b/src/basekernels/matern.jl
@@ -41,7 +41,8 @@ MaternKernel(; nu::Real=1.5, ν::Real=nu, metric=Euclidean()) = MaternKernel(ν,
 
 function _matern(ν::Real, d::Real)
     if iszero(d)
-        return one(d)
+        c = -ν / (ν - 1)
+        return one(d) + c * d^2 / 2
     else
         y = sqrt(2ν) * d
         b = log(besselk(ν, y))
diff --git a/src/chainrules.jl b/src/chainrules.jl
index 3b52860dd..549373876 100644
--- a/src/chainrules.jl
+++ b/src/chainrules.jl
@@ -111,36 +111,111 @@ end
 
 function ChainRulesCore.rrule(s::Sinus, x::AbstractVector, y::AbstractVector)
     d = x - y
-    sind = sinpi.(d)
-    abs2_sind_r = abs2.(sind) ./ s.r
+    abs2_sind_r = (sinpi.(d) ./ s.r) .^ 2
     val = sum(abs2_sind_r)
-    gradx = twoπ .* cospi.(d) .* sind ./ (s.r .^ 2)
+    gradx = π .* sinpi.(2 .* d) ./ s.r .^ 2
     function evaluate_pullback(Δ::Any)
-        return (r=-2Δ .* abs2_sind_r,), Δ * gradx, -Δ * gradx
+        r̄ = -2Δ .* abs2_sind_r ./ s.r
+        s̄ = ChainRulesCore.Tangent{typeof(s)}(; r=r̄)
+        return s̄, Δ * gradx, -Δ * gradx
     end
     return val, evaluate_pullback
 end
 
-## Reverse Rules SqMahalanobis
+function ChainRulesCore.rrule(
+    ::typeof(Distances.pairwise), d::Sinus, x::AbstractMatrix; dims=2
+)
+    project_x = ProjectTo(x)
+    function pairwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, dims)
+        x̄ = collect(zero(x))
+        r̄ = zero(d.r)
+        if dims == 1
+            for j in 1:n, i in 1:n
+                xi = view(x, i, :)
+                xj = view(x, j, :)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- xj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- xj) .^ 2 ./ d.r .^ 3
+                x̄[i, :] += ds
+                x̄[j, :] -= ds
+            end
+        elseif dims == 2
+            for j in 1:n, i in 1:n
+                xi = view(x, :, i)
+                xj = view(x, :, j)
+                ds = twoπ .* Δ[i, j] .* sinpi.(xi .- xj) .* cospi.(xi .- xj) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- xj) .^ 2 ./ d.r .^ 3
+                x̄[:, i] .+= ds
+                x̄[:, j] .-= ds
+            end
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄))
+    end
+    return Distances.pairwise(d, x; dims), pairwise_pullback
+end
 
 function ChainRulesCore.rrule(
-    dist::Distances.SqMahalanobis, a::AbstractVector, b::AbstractVector
+    ::typeof(Distances.pairwise), d::Sinus, x::AbstractMatrix, y::AbstractMatrix; dims=2
 )
-    d = dist(a, b)
-    function SqMahalanobis_pullback(Δ::Real)
-        a_b = a - b
-        ∂qmat = InplaceableThunk(
-            X̄ -> mul!(X̄, a_b, a_b', true, Δ), @thunk((a_b * a_b') * Δ)
-        )
-        ∂a = InplaceableThunk(
-            X̄ -> mul!(X̄, dist.qmat, a_b, true, 2 * Δ), @thunk((2 * Δ) * dist.qmat * a_b)
-        )
-        ∂b = InplaceableThunk(
-            X̄ -> mul!(X̄, dist.qmat, a_b, true, -2 * Δ), @thunk((-2 * Δ) * dist.qmat * a_b)
-        )
-        return Tangent{typeof(dist)}(; qmat=∂qmat), ∂a, ∂b
+    project_x = ProjectTo(x)
+    project_y = ProjectTo(y)
+    function pairwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, dims)
+        m = size(y, dims)
+        x̄ = collect(zero(x))
+        ȳ = collect(zero(y))
+        r̄ = zero(d.r)
+        if dims == 1
+            for j in 1:m, i in 1:n
+                xi = view(x, i, :)
+                yj = view(y, j, :)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- yj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- yj) .^ 2 ./ d.r .^ 3
+                x̄[i, :] .+= ds
+                ȳ[j, :] .-= ds
+            end
+        elseif dims == 2
+            for j in 1:m, i in 1:n
+                xi = view(x, :, i)
+                yj = view(y, :, j)
+                ds = π .* Δ[i, j] .* sinpi.(2 .* (xi .- yj)) ./ d.r .^ 2
+                r̄ .-= 2 .* Δ[i, j] .* sinpi.(xi .- yj) .^ 2 ./ d.r .^ 3
+                x̄[:, i] .+= ds
+                ȳ[:, j] .-= ds
+            end
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄)), @thunk(project_y(ȳ))
+    end
+    return Distances.pairwise(d, x, y; dims), pairwise_pullback
+end
+
+function ChainRulesCore.rrule(
+    ::typeof(Distances.colwise), d::Sinus, x::AbstractMatrix, y::AbstractMatrix
+)
+    project_x = ProjectTo(x)
+    project_y = ProjectTo(y)
+    function colwise_pullback(z̄)
+        Δ = unthunk(z̄)
+        n = size(x, 2)
+        x̄ = collect(zero(x))
+        ȳ = collect(zero(y))
+        r̄ = zero(d.r)
+        for i in 1:n
+            xi = view(x, :, i)
+            yi = view(y, :, i)
+            ds = π .* Δ[i] .* sinpi.(2 .* (xi .- yi)) ./ d.r .^ 2
+            r̄ .-= 2 .* Δ[i] .* sinpi.(xi .- yi) .^ 2 ./ d.r .^ 3
+            x̄[:, i] .+= ds
+            ȳ[:, i] .-= ds
+        end
+        d̄ = ChainRulesCore.Tangent{typeof(d)}(; r=r̄)
+        return NoTangent(), d̄, @thunk(project_x(x̄)), @thunk(project_y(ȳ))
     end
-    return d, SqMahalanobis_pullback
+    return Distances.colwise(d, x, y), colwise_pullback
 end
 
 ## Reverse Rules for matrix wrappers
@@ -150,8 +225,11 @@ function ChainRulesCore.rrule(::Type{<:ColVecs}, X::AbstractMatrix)
     function ColVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
         return error(
             "Pullback on AbstractVector{<:AbstractVector}.\n" *
-            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
-            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `ColVecs`",
+            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`,\n" *
+            "or because some external computation has acted on `ColVecs` to produce a vector of vectors." *
+            "In the former case, to solve this issue overload `kernelmatrix(_diag)` for your kernel for `ColVecs`." *
+            "In the latter case, one needs to track down the `rrule` whose pullback returns a `Vector{Vector{T}}`," *
+            " rather than a `Tangent`, as the cotangent / gradient for `ColVecs` input, and circumvent it."
         )
     end
     return ColVecs(X), ColVecs_pullback
@@ -162,8 +240,9 @@ function ChainRulesCore.rrule(::Type{<:RowVecs}, X::AbstractMatrix)
     function RowVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
         return error(
             "Pullback on AbstractVector{<:AbstractVector}.\n" *
-            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
-            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `RowVecs`",
+            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`,\n" *
+            "or because some external computation has acted on `RowVecs` to produce a vector of vectors." *
+            "If it is the former, to solve this issue overload `kernelmatrix(_diag)` for your kernel for `RowVecs`",
         )
     end
     return RowVecs(X), RowVecs_pullback
diff --git a/src/distances/pairwise.jl b/src/distances/pairwise.jl
index 8b5cb43e7..ab925fba6 100644
--- a/src/distances/pairwise.jl
+++ b/src/distances/pairwise.jl
@@ -6,11 +6,11 @@ end
 
 pairwise(d::PreMetric, X::AbstractVector) = pairwise(d, X, X)
 
-function pairwise!(out::AbstractMatrix, d::PreMetric, X::AbstractVector, Y::AbstractVector)
+function pairwise!(d::PreMetric, out::AbstractMatrix, X::AbstractVector, Y::AbstractVector)
     return broadcast!(d, out, X, permutedims(Y))
 end
 
-pairwise!(out::AbstractMatrix, d::PreMetric, X::AbstractVector) = pairwise!(out, d, X, X)
+pairwise!(d::PreMetric, out::AbstractMatrix, X::AbstractVector) = pairwise!(d, out, X, X)
 
 function pairwise(d::PreMetric, x::AbstractVector{<:Real})
     return Distances_pairwise(d, reshape(x, :, 1); dims=1)
@@ -20,14 +20,14 @@ function pairwise(d::PreMetric, x::AbstractVector{<:Real}, y::AbstractVector{<:R
     return Distances_pairwise(d, reshape(x, :, 1), reshape(y, :, 1); dims=1)
 end
 
-function pairwise!(out::AbstractMatrix, d::PreMetric, x::AbstractVector{<:Real})
-    return Distances.pairwise!(out, d, reshape(x, :, 1); dims=1)
+function pairwise!(d::PreMetric, out::AbstractMatrix, x::AbstractVector{<:Real})
+    return Distances.pairwise!(d, out, reshape(x, :, 1); dims=1)
 end
 
 function pairwise!(
-    out::AbstractMatrix, d::PreMetric, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
+    d::PreMetric, out::AbstractMatrix, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
 )
-    return Distances.pairwise!(out, d, reshape(x, :, 1), reshape(y, :, 1); dims=1)
+    return Distances.pairwise!(d, out, reshape(x, :, 1), reshape(y, :, 1); dims=1)
 end
 
 # Also defines the colwise method for abstractvectors
diff --git a/src/kernels/kernelsum.jl b/src/kernels/kernelsum.jl
index 77709502c..d76735962 100644
--- a/src/kernels/kernelsum.jl
+++ b/src/kernels/kernelsum.jl
@@ -45,6 +45,7 @@ Base.length(k::KernelSum) = length(k.kernels)
 
 _sum(f, ks::Tuple, args...) = f(first(ks), args...) + _sum(f, Base.tail(ks), args...)
 _sum(f, ks::Tuple{Tx}, args...) where {Tx} = f(only(ks), args...)
+_sum(f, ks::AbstractVector, args...) = sum(k -> f(k, args...), ks)
 
 (κ::KernelSum)(x, y) = _sum((k, x, y) -> k(x, y), κ.kernels, x, y)
 
diff --git a/src/matrix/kernelmatrix.jl b/src/matrix/kernelmatrix.jl
index e778f79eb..0b0f851d6 100644
--- a/src/matrix/kernelmatrix.jl
+++ b/src/matrix/kernelmatrix.jl
@@ -134,7 +134,7 @@ kernelmatrix_diag(κ::Kernel, x::AbstractVector, y::AbstractVector) = map(κ, x,
 
 function kernelmatrix!(K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector)
     validate_inplace_dims(K, x)
-    pairwise!(K, metric(κ), x)
+    pairwise!(metric(κ), K, x)
     return map!(x -> kappa(κ, x), K, K)
 end
 
@@ -142,7 +142,7 @@ function kernelmatrix!(
     K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector, y::AbstractVector
 )
     validate_inplace_dims(K, x, y)
-    pairwise!(K, metric(κ), x, y)
+    pairwise!(metric(κ), K, x, y)
     return map!(x -> kappa(κ, x), K, K)
 end
 
diff --git a/src/utils.jl b/src/utils.jl
index 0942fa5f7..29087259b 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -107,11 +107,11 @@ end
 function pairwise(d::PreMetric, x::ColVecs, y::AbstractVector{<:AbstractVector{<:Real}})
     return Distances_pairwise(d, x.X, reduce(hcat, y); dims=2)
 end
-function pairwise!(out::AbstractMatrix, d::PreMetric, x::ColVecs)
-    return Distances.pairwise!(out, d, x.X; dims=2)
+function pairwise!(d::PreMetric, out::AbstractMatrix, x::ColVecs)
+    return Distances.pairwise!(d, out, x.X; dims=2)
 end
-function pairwise!(out::AbstractMatrix, d::PreMetric, x::ColVecs, y::ColVecs)
-    return Distances.pairwise!(out, d, x.X, y.X; dims=2)
+function pairwise!(d::PreMetric, out::AbstractMatrix, x::ColVecs, y::ColVecs)
+    return Distances.pairwise!(d, out, x.X, y.X; dims=2)
 end
 
 """
@@ -178,11 +178,11 @@ end
 function pairwise(d::PreMetric, x::RowVecs, y::AbstractVector{<:AbstractVector{<:Real}})
     return Distances_pairwise(d, x.X, permutedims(reduce(hcat, y)); dims=1)
 end
-function pairwise!(out::AbstractMatrix, d::PreMetric, x::RowVecs)
-    return Distances.pairwise!(out, d, x.X; dims=1)
+function pairwise!(d::PreMetric, out::AbstractMatrix, x::RowVecs)
+    return Distances.pairwise!(d, out, x.X; dims=1)
 end
-function pairwise!(out::AbstractMatrix, d::PreMetric, x::RowVecs, y::RowVecs)
-    return Distances.pairwise!(out, d, x.X, y.X; dims=1)
+function pairwise!(d::PreMetric, out::AbstractMatrix, x::RowVecs, y::RowVecs)
+    return Distances.pairwise!(d, out, x.X, y.X; dims=1)
 end
 
 # Resolve ambiguity error for ColVecs vs RowVecs. #346
diff --git a/test/Project.toml b/test/Project.toml
index 5ec8f4529..a8a9e9928 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,7 @@
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -34,4 +36,5 @@ ReverseDiff = "1.2"
 SpecialFunctions = "0.10, 1, 2"
 StableRNGs = "1"
 StaticArrays = "1"
+Statistics = "1"
 Zygote = "0.6.38"
diff --git a/test/basekernels/periodic.jl b/test/basekernels/periodic.jl
index fb149dff5..540947b1b 100644
--- a/test/basekernels/periodic.jl
+++ b/test/basekernels/periodic.jl
@@ -15,7 +15,6 @@
     TestUtils.test_interface(PeriodicKernel(; r=[0.9, 0.9]), ColVecs{Float64})
     TestUtils.test_interface(PeriodicKernel(; r=[0.8, 0.7]), RowVecs{Float64})
 
-    # test_ADs(r->PeriodicKernel(r =exp.(r)), log.(r), ADs = [:ForwardDiff, :ReverseDiff])
-    @test_broken "Undefined adjoint for Sinus metric, and failing randomly for ForwardDiff and ReverseDiff"
+    test_ADs(r -> PeriodicKernel(; r=exp.(r)), log.(r))
     test_params(k, (r,))
 end
diff --git a/test/basekernels/sm.jl b/test/basekernels/sm.jl
index 5cbc37243..d9037aadd 100644
--- a/test/basekernels/sm.jl
+++ b/test/basekernels/sm.jl
@@ -57,5 +57,6 @@
     end
 
     # test_ADs(x->spectral_mixture_kernel(exp.(x[1:3]), reshape(x[4:18], 5, 3), reshape(x[19:end], 5, 3)), vcat(log.(αs₁), γs[:], ωs[:]), dims = [5,5])
-    @test_broken "No tests passing (BaseKernel)"
+    # No tests passing (BaseKernel)
+    @test_broken false
 end
diff --git a/test/basekernels/wiener.jl b/test/basekernels/wiener.jl
index 9dd60ba43..621b3dcbc 100644
--- a/test/basekernels/wiener.jl
+++ b/test/basekernels/wiener.jl
@@ -44,5 +44,6 @@
     TestUtils.test_interface(k2, x0, x1, x2)
     TestUtils.test_interface(k3, x0, x1, x2)
     # test_ADs(()->WienerKernel(i=1))
-    @test_broken "No tests passing"
+    # No tests passing
+    @test_broken false
 end
diff --git a/test/chainrules.jl b/test/chainrules.jl
index 03c2c3b1f..f2a4ae3a5 100644
--- a/test/chainrules.jl
+++ b/test/chainrules.jl
@@ -3,8 +3,6 @@
     x = rand(rng, 5)
     y = rand(rng, 5)
     r = rand(rng, 5)
-    Q = Matrix(Cholesky(rand(rng, 5, 5), 'U', 0))
-    @assert isposdef(Q)
 
     compare_gradient(:Zygote, [x, y]) do xy
         Euclidean()(xy[1], xy[2])
@@ -21,11 +19,23 @@
     compare_gradient(:Zygote, [x, y]) do xy
         KernelFunctions.Sinus(r)(xy[1], xy[2])
     end
-    if VERSION < v"1.6"
-        @test_broken "Chain rule of SqMahalanobis is broken in Julia pre-1.6"
-    else
-        compare_gradient(:Zygote, [Q, x, y]) do Qxy
-            SqMahalanobis(Qxy[1])(Qxy[2], Qxy[3])
+    @testset "rrules for Sinus(r=$r)" for r in (rand(3),)
+        dist = KernelFunctions.Sinus(r)
+        @testset "$type" for type in (Vector, SVector{3})
+            test_rrule(dist, type(rand(3)), type(rand(3)))
+        end
+        @testset "$type1, $type2" for type1 in (Matrix, SMatrix{3,2}),
+            type2 in (Matrix, SMatrix{3,4})
+
+            test_rrule(Distances.pairwise, dist, type1(rand(3, 2)); fkwargs=(dims=2,))
+            test_rrule(
+                Distances.pairwise,
+                dist,
+                type1(rand(3, 2)),
+                type2(rand(3, 4));
+                fkwargs=(dims=2,),
+            )
+            test_rrule(Distances.colwise, dist, type1(rand(3, 2)), type1(rand(3, 2)))
         end
     end
 end
diff --git a/test/distances/pairwise.jl b/test/distances/pairwise.jl
index 486097f52..6635579cb 100644
--- a/test/distances/pairwise.jl
+++ b/test/distances/pairwise.jl
@@ -11,10 +11,10 @@
 
     @test KernelFunctions.pairwise(d, x, y) ≈ pairwise(d, X, Y; dims=2)
     @test KernelFunctions.pairwise(d, x) ≈ pairwise(d, X; dims=2)
-    KernelFunctions.pairwise!(K, d, x, y)
+    KernelFunctions.pairwise!(d, K, x, y)
     @test K ≈ pairwise(d, X, Y; dims=2)
     K = zeros(Ns[1], Ns[1])
-    KernelFunctions.pairwise!(K, d, x)
+    KernelFunctions.pairwise!(d, K, x)
     @test K ≈ pairwise(d, X; dims=2)
 
     x = randn(rng, 10)
@@ -24,6 +24,6 @@
     K = zeros(10, 11)
     @test KernelFunctions.pairwise(d, x, y) ≈ pairwise(d, X, Y; dims=1)
     @test KernelFunctions.pairwise(d, x) ≈ pairwise(d, X; dims=1)
-    KernelFunctions.pairwise!(K, d, x, y)
+    KernelFunctions.pairwise!(d, K, x, y)
     @test K ≈ pairwise(d, X, Y; dims=1)
 end
diff --git a/test/kernels/kernelsum.jl b/test/kernels/kernelsum.jl
index 4b6f30f94..eb64aeee8 100644
--- a/test/kernels/kernelsum.jl
+++ b/test/kernels/kernelsum.jl
@@ -2,21 +2,24 @@
     k1 = LinearKernel()
     k2 = SqExponentialKernel()
     k = KernelSum(k1, k2)
-    @test k == KernelSum([k1, k2]) == KernelSum((k1, k2))
+    kvec = KernelSum([k1, k2])
+    @test k == kvec == KernelSum((k1, k2))
     for (_k1, _k2) in Iterators.product(
         (k1, KernelSum((k1,)), KernelSum([k1])), (k2, KernelSum((k2,)), KernelSum([k2]))
     )
         @test k == _k1 + _k2
+        @test kvec == _k1 + _k2
     end
-    @test length(k) == 2
-    @test repr(k) == (
+    @test length(k) == length(kvec) == 2
+    @test repr(k) ==
+        repr(kvec) ==
         "Sum of 2 kernels:\n" *
         "\tLinear Kernel (c = 0.0)\n" *
         "\tSquared Exponential Kernel (metric = Euclidean(0.0))"
-    )
 
     # Standardised tests.
     test_interface(k, Float64)
+    test_interface(kvec, Float64)
     test_interface(ConstantKernel(; c=1.5) + WhiteKernel(), Vector{String})
     test_ADs(x -> KernelSum(SqExponentialKernel(), LinearKernel(; c=exp(x[1]))), rand(1))
     test_interface_ad_perf(2.4, StableRNG(123456)) do c
diff --git a/test/runtests.jl b/test/runtests.jl
index caf43cb91..e054b992a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,7 @@
 using KernelFunctions
 using AxisArrays
+using ChainRulesCore
+using ChainRulesTestUtils
 using Distances
 using Documenter
 using Functors: functor
diff --git a/test/test_utils.jl b/test/test_utils.jl
index 8367fbd6d..aa8c81e1f 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -384,7 +384,7 @@ function test_zygote_perf_heuristic(
             @test_broken fwd[1] == fwd[2]
         end
         if passes[3]
-            @test pb[1] == pb[2]
+            @test abs(pb[1] - pb[2]) ≤ 1
         else
             @test_broken pb[1] == pb[2]
         end
diff --git a/test/transform/selecttransform.jl b/test/transform/selecttransform.jl
index b0e3bfa81..a6d382462 100644
--- a/test/transform/selecttransform.jl
+++ b/test/transform/selecttransform.jl
@@ -104,17 +104,41 @@
     end
 
     @testset "$(AD)" for AD in [:ReverseDiff]
-        @test_broken ga = gradient(AD, A) do a
-            testfunction(ta_row, a, 2)
+        @test_broken let
+            gx = gradient(AD, X) do x
+                testfunction(tx_row, x, 2)
+            end
+            ga = gradient(AD, A) do a
+                testfunction(ta_row, a, 2)
+            end
+            gx ≈ ga
         end
-        @test_broken ga = gradient(AD, A) do a
-            testfunction(ta_col, a, 1)
+        @test_broken let
+            gx = gradient(AD, X) do x
+                testfunction(tx_col, x, 1)
+            end
+            ga = gradient(AD, A) do a
+                testfunction(ta_col, a, 1)
+            end
+            gx ≈ ga
         end
-        @test_broken ga = gradient(AD, A) do a
-            testfunction(ta_row, a, B, 2)
+        @test_broken let
+            gx = gradient(AD, X) do x
+                testfunction(tx_row, x, Y, 2)
+            end
+            ga = gradient(AD, A) do a
+                testfunction(ta_row, a, B, 2)
+            end
+            gx ≈ ga
         end
-        @test_broken ga = gradient(AD, A) do a
-            testfunction(ta_col, a, C, 1)
+        @test_broken let
+            gx = gradient(AD, X) do x
+                testfunction(tx_col, x, Z, 1)
+            end
+            ga = gradient(AD, A) do a
+                testfunction(ta_col, a, C, 1)
+            end
+            gx ≈ ga
         end
     end
 
diff --git a/test/utils.jl b/test/utils.jl
index 42784548a..a7cabbd05 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -47,10 +47,10 @@
         @test vcat(DX, DY) isa ColVecs
         @test vcat(DX, DY).X == hcat(X, Y)
         K = zeros(N, N)
-        KernelFunctions.pairwise!(K, SqEuclidean(), DX)
+        KernelFunctions.pairwise!(SqEuclidean(), K, DX)
         @test K ≈ pairwise(SqEuclidean(), X; dims=2)
         K = zeros(N, N + 1)
-        KernelFunctions.pairwise!(K, SqEuclidean(), DX, DY)
+        KernelFunctions.pairwise!(SqEuclidean(), K, DX, DY)
         @test K ≈ pairwise(SqEuclidean(), X, Y; dims=2)
 
         let
@@ -105,10 +105,10 @@
         @test vcat(DX, DY) isa RowVecs
         @test vcat(DX, DY).X == vcat(X, Y)
         K = zeros(D, D)
-        KernelFunctions.pairwise!(K, SqEuclidean(), DX)
+        KernelFunctions.pairwise!(SqEuclidean(), K, DX)
         @test K ≈ pairwise(SqEuclidean(), X; dims=1)
         K = zeros(D, D + 1)
-        KernelFunctions.pairwise!(K, SqEuclidean(), DX, DY)
+        KernelFunctions.pairwise!(SqEuclidean(), K, DX, DY)
         @test K ≈ pairwise(SqEuclidean(), X, Y; dims=1)
 
         let