# Making multivariate Kolmogorov-Smirnov benchmarks

**URL:** <https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114>\
**Category:** Performance\
**Created:** [January 24, 2022, 12:46pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114 "2022-01-24T12:46:11Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 24, 2022, 12:46pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/1 "2022-01-24T12:46:11Z")

</div>

Hey,

I’m trying to produce empirical multivariate Kolmogorov-Smirnov benchs to later benchmark an estimator. I have the following code for the moment :

```julia
# The goal of this file is to compute kolmogorov-smirnov p-values of stuff. 
using ProgressMeter
function get_ks_dist(newD,baseD)
    N,d = size(newD)
    dist = 0
    for i in 1:N # index of the sample. 
        x_1 = 0
        x_2 = 0
        for k in 1:N
            x1k = true
            x2k = true
            for l in 1:d
                x1k &= baseD[k,l] < newD[i,l]
                x2k &= newD[k,l] < newD[i,l]
            end
            x_1 += x1k ? 1 : 0
            x_2 += x2k ? 1 : 0
        end
        dist = max(dist, x_1 - x_2, x_2 - x_1)
    end
    return dist/N
end
function make_ks_benchmark(M,get_sample)
    # D is a multivariate distribution. 
    base_sample = get_sample()
    dists = zeros(M)
    @showprogress for i in 1:M
        dists[i] = get_ks_dist(get_sample(),base_sample)
    end
    return sort(dists)
end

function get_ln_indep_sample()
    d = 10
    N = 10000
    reshape(exp.(randn(d*N)),(N,d))
end
bench = make_ks_benchmark(100,get_ln_indep_sample)

```

And it takes a lot of time. I’ve checked `@code_warntype`, and it seems like `get_ks_dist` is type-stable. I was wandering if passing the function as argument was the problem, or if it is just my computations that are crazy long. If so, what can be done ?

Thnaks for any idea 🙂

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 24, 2022, 1:02pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/2 "2022-01-24T13:02:36Z")

</div>

Iterating on consecutive elts seems to speed up the computation  
Note that I have changed N to 2000 for convenience

```nohighlight
function get_ks_dist(newD,baseD)
    d,N = size(newD)
    dist = 0
    @inbounds for i in 1:N # index of the sample. 
        x_1 = 0
        x_2 = 0
        for k in 1:N
            x1k = true
            x2k = true
            for l in 1:d
                x1k &= baseD[l,k] < newD[l,i]
                x2k &= newD[l,k] < newD[l,i]
            end
            x_1 += Int(x1k)
            x_2 += Int(x2k)
            # x_1 += x1k ? 1 : 0
            # x_2 += x2k ? 1 : 0
        end
        dist = max(dist, x_1 - x_2, x_2 - x_1)
    end
    return dist/N
end
function make_ks_benchmark(M,get_sample)
    # D is a multivariate distribution. 
    base_sample = get_sample()
    dists = zeros(M)
    @showprogress for i in 1:M
        dists[i] = get_ks_dist(get_sample(),base_sample)
    end
    return sort(dists)
end

function get_ln_indep_sample()
    d = 10
    N = 2000

    reshape(exp.(randn(d*N)),(d,N))
end
@time bench = make_ks_benchmark(100,get_ln_indep_sample)

```

---

<div class="post-metadata">

**Author:** ![Andrea\_Pagnani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andrea_pagnani/32/3136_2.png) [@Andrea\_Pagnani](https://discourse.julialang.org/u/Andrea_Pagnani)\
**Post date:** [January 24, 2022, 1:04pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/3 "2022-01-24T13:04:16Z")

</div>

In julia the “fast” index for matrices is the leftmost one. You should reshape the matrix such that the innermost loop is. Also an inbound could help

```julia
@inbounds for l in 1:d
         x1k &= baseD[l,k] < newD[l,i]
         x2k &= newD[l,k] < newD[l,i]
end

```

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 24, 2022, 1:10pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/4 "2022-01-24T13:10:43Z")

</div>

That’s indeed a lot better, twice as fast on my machine. Thanks.

Edit: I also threaded the main loop. Surprisingly, even for `d=10`, the distribution does look a lot like the univariate one from Wikipedia :

```julia
function get_ks_dist(newD,baseD)
    d,N = size(newD)
    dist = 0
    @inbounds for i in 1:N # index of the sample. 
        x_1 = 0
        x_2 = 0
        for k in 1:N
            x1k = true
            x2k = true
            for l in 1:d
                x1k &= (baseD[l,k] < newD[l,i])
                x2k &= (newD[l,k] < newD[l,i])
            end
            x_1 += x1k ? 1 : 0
            x_2 += x2k ? 1 : 0
        end
        dist = max(dist, x_1 - x_2, x_2 - x_1)
    end
    return dist/N
end
function make_ks_benchmark(M,get_sample)
    # D is a multivariate distribution. 
    base_sample = get_sample()
    dists = zeros(M)
    p = Progress(M)
    Threads.@threads for i in 1:M
        dists[i] = get_ks_dist(get_sample(),base_sample)
        next!(p)
    end
    return sort(dists)
end

function get_ln_indep_sample()
    d = 10
    N = 1000

    reshape(exp.(randn(d*N)),(d,N))
end
@time bench = make_ks_benchmark(10000,get_ln_indep_sample);
histogram(bench)

```

![plot_2](https://global.discourse-cdn.com/julialang/original/3X/e/6/e611423805ae14e0a4da6da469228c4fcc119400.png)  
To compare to the plot there : [Kolmogorov–Smirnov test - Wikipedia](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test#Kolmogorov_distribution)  
My guess is that it is because the variables are independent in this simple example, and something else would append with a dependence structure, although i’m not an expert.

But it is still a lot of computations ! I’d like to put `N=10_000` instead of `N=1000` so I still lack an order of magnitude. In fact, switching from `N=1000` to `N=10000` makes my runtime jump from half a minut to 50mins, which is to be expected since the loop is quadratic. So two orders of magnitudes…

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 24, 2022, 2:52pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/5 "2022-01-24T14:52:31Z")

</div>

Depending on the value of `d` it may be faster to find a branchless (and simd friendly) way to determine the values of `x1k` and `x2k`.  
(kind of _sorting network_)

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 24, 2022, 2:54pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/6 "2022-01-24T14:54:36Z")

</div>

Hum… `d` goes from 1 to ~ 200 in my applications. I did not tried yet with `d` bigger than `20`, but I’ll have to soon. Thanks for the tip, i’ll take a look.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 24, 2022, 2:56pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/7 "2022-01-24T14:56:00Z")

</div>

Well if d is equal to 200 it is pretty sure that a shortcut solution will win: exit the loop before reaching the end

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [January 24, 2022, 2:57pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/8 "2022-01-24T14:57:29Z")

</div>

Which is why i used `all()` at the begining instead of this loop, as the `all` function does take these kind of shortcuts. For `d=10` it did not improve things, suprisingly, but i might try again for bigger dimensions.

---

<div class="post-metadata">

**Author:** ![jstrube](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jstrube/32/525_2.png) [@jstrube](https://discourse.julialang.org/u/jstrube)\
**Post date:** [June 26, 2022, 10:01pm UTC](https://discourse.julialang.org/t/making-multivariate-kolmogorov-smirnov-benchmarks/75114/9 "2022-06-26T22:01:47Z")

</div>

No Julia implementation, yet, but maybe this code might be useful to speed things up: [GitHub - pnnl/DDKS: A high-dimensional Kolmogorov-Smirnov distance for comparing high dimensional distributions](https://github.com/pnnl/DDKS). There’s also a paper at [[2106.13706] Accelerated Computation of a High Dimensional Kolmogorov-Smirnov Distance](https://arxiv.org/abs/2106.13706). Hope it helps.