Package/algorithm for barycentric grid interpolation

Specific Domains Numerics
1

I am looking for a package or an algorithm (which I am happy to code into a package and make it available) that expresses a coordinate inside a simplex as a convex combination “near” of points with integer coordinates.

Formally, for a given integer K \in \mathbb{N}, consider the simplex

S = \{ x \in \mathbb{R}^M \mid x_i \ge 0\ \forall i, \sum_{i=0}^M x_i \le K \}

defined by the sum of coordinates being \le K, and the points with integer coordinates

I = S \cap \mathbb{N}^M

For an z \in S, I am looking for vertices v_i \in I, i = 0, \dots, M and weights \sum_{i = 0}^M \alpha_i = 1 such that

z = \sum_{i = 0}^M \alpha_i v_i

Ideally the v_i should be “near” z, but not necessarily the nearest if that’s difficult.

The problem is very easy for M = 2, but the naive algorithms I came up with fail for M = 3 — 2D is very easy to partition with self-similar simplexes like this:

2d
2d1093×1006 58.5 KB

But, importantly, for this algorithm the simplexes formed by the v_i for each z do not have to form a partition of S.

Here is some test code I would like to pass, if anyone wants to provide concrete code:

using LinearAlgebra, Test

f(z, K) = IMPLEMENTATION # returns a vector of vs and αs

ϵ = 1e-8 # just picked some arbitrary tolerance

for _ in 1:100
    M = rand(3:5)
    K = rand(2:10)
    z = normalize(rand(M), 1) .* K
    vs, αs = f(z, K)
    ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
    @test ẑ ≈ z atol = ϵ
    @test sum(αs) ≈ 1 atol = ϵ
    @test all(αs .≥ 0)
    for v in vs
        @test eltype(v) == Int # or all v are integers, can be stored as floats
        @test all(v .≥ 0)
        @test sum(v) ≤ K
    end
end

But, again, I appreciate any kind of suggestions.

2

Are you looking for something like
“Wasserstein Barycentric Coordinates:
Histogram Regression Using Optimal Transport”
Or Are you looking for the barycenter itself?

Barycenters under and optimal transport cost are relatively straightforward to calculate, especially using the sinkhorn algorithm.

3

I am looking for the solution of the problem above. I read the article, but I am not sure how it is related (sorry for being dense, this is not my area of expertise).

5

I found a nice thesis about these topics, Moore (1992). I may be able to employ the symmetric subdivision algorithm (Section 2.2.3) using the Kuhn triangulation, I just need to understand it first.

6

If I understand your problem, the solution may be a simplicial triangulation of “points” with integer coordinates. The problem is tricky with the Delaunay triangulation, since the triangulation is not unique. However, a regular mesh like this may be readily generated in 2 and 3 dimensions using a cell template. You mentioned Kuhn triangulation: in 3d that is one possibility. GitHub - PetrKryslUCSD/FinEtools.jl: Finite Element tools in Julia can generate these kinds of meshes, with different templates for the cells. In particular, refer to the lines below.

5 Likes
7

Thanks — I am may be able to generalize this to \mathbb{R}^N.

Is there a book/introductory article you would recommend about this topic? (symmetric subdivision of simplexes).

8

I think I figured it out. The key is to transform to what Coxeter calls R-simplexes, where the Kuhn triangulation will always be in the large simplex by construction.

I include a pedagogical implementation below if anyone is interested — a much nicer, optimized, non-allocating, and efficient version will be made available in a package soon (which accompanies a paper that motivated this whole thing) — the vs of course do not need to be constructed, the permutation is sufficient. I will of course link in both here. In the meantime,

###
### horribly inefficient pedagogical implementation
###

StoR(z) = reverse(cumsum(reverse(z))) # did I say horribily inefficient?

RtoS(x) = vcat(.- diff(x), x[end:end])

function f(z, K)
    # go to canonical R-simplex
    c = StoR(z)
    # fractional and integer parts
    fv = map(x -> ((f, i) = modf(x); (f, Int(i))), c)
    # “lower” corner
    v0 = last.(fv)
    # fractional part as a convex combination, Kuhn triangulation
    fi = sort!(collect(enumerate(first.(fv))); by = last, rev = true)
    # gridpoints mutated buffer for walking the vertices
    vv = copy(v0)
    vs = [(vv[fi[1]] += 1; copy(vv)) for fi in fi]
    push!(vs, v0)
    # weights
    αs = RtoS((last.(fi)))
    push!(αs, 1 - sum(αs))
    # back to inflated standard simplex
    map(RtoS, vs), αs
end

ϵ = 1e-8 # just picked some arbitrary tolerance

for _ in 1:1000
    M = rand(3:5)
    K = rand(2:10)
    z = normalize(rand(M), 1) .* (K * rand())
    vs, αs = f(z, K)
    ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
    @test ẑ ≈ z atol = ϵ
    @test sum(αs) ≈ 1 atol = ϵ
    @test all(αs .≥ 0)
    @test length(αs) == length(vs) == M + 1
    for v in vs
        @test eltype(v) == Int
        @test length(v) == M
        @test all(v .≥ 0)
        @test sum(v) ≤ K
    end
end

BTW, the book I have been looking for is

4 Likes