# Package/algorithm for barycentric grid interpolation

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:

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)
ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
@test ẑ ≈ 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.

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.

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).

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.

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. https://github.com/PetrKryslUCSD/FinEtools.jl can generate these kinds of meshes, with different templates for the cells. In particular, refer to the lines below.

5 Likes

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

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)
ẑ = mapreduce((v, α) -> v .* α, (v1, v2) -> v1 .+ v2, vs, αs)
@test ẑ ≈ 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