Math for borrow checking NumPy arrays

[adamreichold]

I am looking for someone who can help us with the math underlying dynamic borrow checking for NumPy arrays. We are trying to make the rust-numpy crate enforce aliasing discipline when creating Rust references into the interior of NumPy arrays. As NumPy supports multiple views into a single allocation and also wrapping foreign array types (we do this for example to wrap ndarray arrays), this leads to the following general problem which I would call a linear equation in a constrained integer domain:

Given an array A with data pointer P_A, dimensionality N_A, dimensions D_{1,A}, ..., D_{N_A,A} (positive integers) and strides S_{1,A}, ..., S_{N_A,A} (non-zero integers) and an array B with data pointer P_B, dimensionality N_B, dimensions D_{1,B}, ..., D_{N_B,B} and strides S_{1,B}, ..., S_{N_B,B}, does the equation

P_A + S_{1,A} * I_{1,A} + ... + S_{N_A,A} * I_{N_A,A} = P_B + S_{1,B} * I_{1,B} + ... S_{N_B,B} * I_{N_B,B}

have a solution in the domain defined by 0 <= I_{i,a} < D_{i,a} for each i in {1,...,N_a} and a in {A,B}.

So basically, are there any two sets indices for A and B which will point to the same element in the underlying allocation.

(The main obstacle leading to this is that Rust's type system is not around when these views are created, in contrast to e.g. ndarray, so that we have to infer this post-facto via runtime checks instead of checking e.g. the indices that define the views themselves.)

I would be grateful if anybody had advice on how best tackle this? Or could point me to some literature which would help me to better understand where similar problems are addressed?

(I think there is an angle to consider this as an optimization problem by aiming to minimize the difference between the left-hand side and the right-side of the equation and then checking if it is actually zero. But IIRC integer linear programs with bounds often end up being NP? I also vaguely fear that this might end up being equivalent to satisfiability via a binary representation of the indices.)

[adamreichold]

I also suspect that the bounds are not really fundamental to the issue: The arrays do not alias if they arrays do not overlap, i.e the intervals defined by the largest and smallest element addresses do not overlap, or they overlap but interleave without touching which should then be true for all indices, not just those within bounds. At least I cannot currently think of another situation where there is no aliasing of multiple views into the same allocation.

(Strides can be negative, so the smallest address is not necessarily at index 0, ..., 0.) (Interleaving would for example be something like a NxMx3 dimensional array representing RGB data where the views array[:,:,0], array[:,:,1] and array[:,:,2] would represent the separate R, G and B channels.)

[adamreichold]

@adamreichold

We are trying to make the rust-numpy crate enforce aliasing discipline when creating Rust references into the interior of NumPy arrays.

Maybe this should better explained: We already do enforce aliasing discipline, but by over-approximation, i.e. we consider any two views into the same allocation as possibly conflicting. Now we want to refine this as to not reject valid uses where the views do not actually alias because they do not overlap or interleave properly.

[alice]

Well, in 1D case, first you note that they can only overlap if d_A and d_B differ by a multiple of gcd(s_a, s_b). Here you should note that zero is a multiple of everything (though it's an easy "yes" if they are equal).

Then, if that is the case, you rewrite as follows:

d_a + s_a i_a = d_b + s_b i_b

s_a i_a + (-s_b) i_a = d_b - d_a

You can use the extended euclidean algorithm to find a solution to the above equation (which may be out of bounds). Then, by Bézout's identity, all solutions to the above equation are of the following form:

i_a_new = i_a + F_a * k
i_b_new = i_b + F_b * k

for integers k where F_a = (d_b - d_a) * s_b / (gcd(s_a, s_b)^2) and F_b = (d_b - d_a) * s_a / (gcd(s_a, s_b)^2), and where i_a and i_b is the solution you got from the extended euclidean algorithm. You can translate the above into four inequalities on k that are required to hold for the solution to be within bounds for both arrays by rewriting this:

0   <= i_a + F_a * k
L_a >= i_a + F_a * k
0   <= i_b + F_b * k
L_b >= i_b + F_b * k

into this:

-i_a / F_a <= k
(L_a - i_a) / F_a >= k
-i_b / F_b <= k
(L_b - i_b) / F_b >= k

Note that the fractions in the second form can result in non-integers. You should handle this by rounding up when the inequality is <= and rounding down when the inequality is >= since you know k is an integer. Here up/down means towards +/- infinity. Don't round towards zero when rounding down.

Anyway, the above equations will give a lower bound and upper bound on k. If the lower bound is greater than the upper bound, then there are no solutions, otherwise there are solutions.

Note that I use d_a and d_b for the data pointers and L_a and L_b for the lengths of the arrays to avoid using the same character for two things.

[adamreichold]

@alice

Well, in 1D case, first you note that they can only overlap if d_A and d_B differ by a multiple of gcd(s_a, s_b) .

I wonder whether this could be made into a general argument for rejecting overlap by considering the GCD of the smallest stride of each array. I think this is how one would argue in the RGB example.

@alice

Note that I use d_a and d_b for the data pointers and L_a and L_b for the lengths of the arrays to avoid using the same character for two things.

Good point, I will update the initial post to use P_A and P_B for the data pointers.

[alice]

I found this question which gives some ideas on how you might be able to extend the 1D case:

(Diophantine?) Equations With Multiple Variables

It certainly gives a simple partial solution if you also disallow overlap when it is out-of-bounds. That said, when taking gcd of many different things, you very quickly just end up with 1 (rejecting all) unless the strides are carefully chosen.

[adamreichold]

I feel a bit uncouth: I think my first reply should have been "Thank you for the detailed discussed of the one-dimensional case!". Sorry for forgetting that.

I will try to understand the Stack Exchange discussion and I think "Diophantine equation" will be useful starting point for a literature search. Thank you for that as well!

[adamreichold]

@alice

It certainly gives a simple partial solution if you also disallow overlap when it is out-of-bounds.

Since the aim is to refine an over-approximation, I think checking the GCD of all strides if the data ranges overlap would already be an improvement.

As written above, I also suspect that the strides coming from the same underlying array somehow constraint the problem so that we have only those two cases without aliasing: Either the data ranges are disjoint or they never alias, whether in bounds or not. At least, I cannot think of an example where the array views would alias only out of bounds.

[alice]

@adamreichold

At least, I cannot think of an example where the array views would alias only out of bounds.

What if you have two 1D arrays both with stride 1, but which are disjoint?

More generally, if any of your strides is 1, then the check always fails.

[adamreichold]

@alice

What if you have two 1D arrays both with stride 1, but which are disjoint?

The example was meant restricted to the case that the data ranges are disjoint, i.e.

@adamreichold

Either the data ranges are disjoint or they never alias, whether in bounds or not.

My current plan is to track the data range that the views cover and if that is disjoint, considering them non-conflicting. Otherwise, check if the GDC of the strides divides the difference data pointer and if not consider them non-conflicting. Otherwise, I consider them conflicting.

The last case might include situations where they meet only out of bounds and they do not really conflict, but it is a refinement over considering only the data ranges as I currently do in WIP: Refine dynamic borrow checking to track data ranges. by adamreichold · Pull Request #299 · PyO3/rust-numpy · GitHub.

[adamreichold]

@adamreichold

My current plan is to track the data range that the views cover and if that is disjoint, considering them non-conflicting. Otherwise, check if the GDC of the strides divides the difference data pointer and if not consider them non-conflicting. Otherwise, I consider them conflicting.

This seems to work nicely. Thank you again for your quick help!

[alice]

You're welcome!

[adamreichold]

@adamreichold

As written above, I also suspect that the strides coming from the same underlying array somehow constraint the problem so that we have only those two cases without aliasing: Either the data ranges are disjoint or they never alias, whether in bounds or not. At least, I cannot think of an example where the array views would alias only out of bounds.

So turns out this is not correct when a view is created by slicing with a step size that does not divide the dimension along that axis:

import numpy as np
a = np.zeros((10, 10))
a.strides # (80, 8)

v1 = a[:,::3]
v1.strides # (80, 24)

v2 = a[:,1::3]
v2.strides # (80, 24)

where v1 and v2 will have data pointers differing by 8 and the GCD of their strides is also 8. The smallest solution to 80*a+24*b-80*c-24*d=8 is a=0,b=7,c=2,d=0 which is however out of bounds for v1.

In contrast

v3 = a[:, ::2]
v3.strides # (80, 16)

v4 = a[:, 1::2]
v4.strides # (80, 16)

has a data pointer difference of 8 again but a GCD of the strides of 16 which means the GCD criterion alone already rules out conflicts.

So we either live we these false conflicts for now or we recursively solve the Diophantine equations using the extended Euclidean algorithm...

[alice]

Unfortunately, I can't think of an easy way out.