Finding a common step value between an array of non-uniformly distributed points

I have an array of values that do not have a uniform step value between them. I'm looking to find the largest step value such that a linspace array from 0 to the maximum value of the array with that step value will include all the points of the original array,

I feel like it must have something to do with common factors, but I can't quite figure it out.

I've listed some examples below. Any help would be appreciated!

A = [0 0.1 0.5];
Astep = 0.1; % step value required
Arange = 0:Astep:max(A);
Aidx = ismember(A,Arange);

B = [0 0.1 0.5 0.75];
Bstep = 0.05; % step value required
Brange = 0:Bstep:max(B);
Bidx = ismember(B,Brange);

C = [0 0.1 0.5 0.75 5.12];
Cstep = 0.01; % step value required
Crange = 0:Cstep:max(C);
Cidx = ismember(C,Crange);

correctstep = all([Aidx Bidx Cidx])

correctstep = logical
   1

Answer 1

[John D'Errico]

This is generally impossible to do exactly, IF your numbers are floating point. Let me repeat, effectively impossible. If your numbers are integers, then it is trivial. Just compute the greatest common divisor of the diff between elements. FOR EXAMPLE...

V = [0 5 10 20 35 40];

Clearly, the stride you would need is 5. We can see that here:

dV = diff(V)

dV = 1×5

5 5 10 15 5

But can we do it programmatically? The problem is, gcd does not work to find the gcd of an entire set of numbers, at least not for doubles. But gcd does work for a sym vector. So we can do this:

stride = double(gcd(sym(dV)))

stride = 5

0:stride:max(V)

ans = 1×9

0 5 10 15 20 25 30 35 40

And that works nicely. HOWEVER, you cannot do the same thing if your numbers are floating point. For example, suppose we had:

format long g
V = [0,sort(rand(1,5))]

V = 1×6

0 0.488391892303935 0.627274353393812 0.719333733143841 0.97495042449865 0.993659562773364

dV = diff(V)

dV = 1×5

0.488391892303935 0.138882461089878 0.0920593797500285 0.255616691354809 0.0187091382747141

Clearly gcd will not work on floating point numbers anyway. (If you look carefully at the standard GCD algorithm, you would understand why not.) And if it did, you would get garbage there, where the gcd of those differences is probably as small as eps.

But even suppose your numbers are nice, and "well-behaved"?

V = linspace(0,1,7)

V = 1×7

0 0.166666666666667 0.333333333333333 0.5 0.666666666666667 0.833333333333333 1

dV = diff(V)

dV = 1×6

0.166666666666667 0.166666666666667 0.166666666666667 0.166666666666667 0.166666666666667 0.166666666666667

And while those numbers may look identical, they are not. For example, consider this:

numel(unique(dV))

ans = 
     4

So there were 4 distinct differences found there, all of them varying by an amount on the order of eps. Again, any GCD operation will fail there. Sorry. And this is why I said your question is essentially at least technically very difficult to solve. Yes, a good coder who completely understands floating point arithmetic could probably make some headway. So in theory, I could probably cook up some code. It would not be anything trivial to write, and it would probably not be perfect, and would definitely require a tolerance. Thinking about the idea, I might do it using some sort of quadratic programming code, designed to find an optimal stride, with a tolerance. Hmm. That won't work, because this is a discrete thing. Ugh.

[Hassaan]

@John D'Errico Good observation.

[Jack]

Would this be possible to do with a minimum step size to circurmvent the issues with floating points? Say 0.01, and you would round everything beforehand to this precision before calculating LCM?

And might the problems with Muhammad's answer above be fixed by dividing the array by the largest values such that the whole array is less than or equal to 1 and so there is no integer issue any more...

Just some thoughts.

Answer 2

[Hassaan]

% Test the function with your arrays
A = [0, 0.1, 0.5];
B = [0, 0.1, 0.5, 0.75];
C = [0, 0.1, 0.5, 0.75, 5.12];

A_step = findSmallestStep(A);
B_step = findSmallestStep(B);
C_step = findSmallestStep(C);

disp(['A_step: ', num2str(A_step)]);  % Expected to be around 0.1

A_step: 0.1

disp(['B_step: ', num2str(B_step)]);  % Expected to be around 0.05

B_step: 0.05

disp(['C_step: ', num2str(C_step)]);  % Expected to be around 0.01

C_step: 0.01

function stepValue = findSmallestStep(arr)
    % Calculate differences between consecutive elements
    diffs = diff(arr);

    % Convert differences to fractions and find LCM of denominators
    lcm_den = 1;
    for i = 1:length(diffs)
        [~, den] = rat(diffs(i));
        lcm_den = lcm(lcm_den, den);
    end

    % The step is the reciprocal of the LCM of denominators
    stepValue = 1 / lcm_den;
end

This function first finds the differences between consecutive elements. Each difference is then converted to its fractional representation using rat, and the LCM of all denominators is computed. The required step is the reciprocal of this LCM, ensuring it's the smallest step size that can accommodate all intervals in the array.

---------------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[Jack]

This is great and just what I was looking for, thank you! If I wanted to fix the maximum number of decimal places of the step, I assume that would be done in the tolerance input for rat.

[Dyuman Joshi]

@Jack, The code posted above does not work for all cases -

% Test the function with your arrays
A = [5 7.5 12.5 50]

A = 1×4

5.0000 7.5000 12.5000 50.0000

B = [23 69 230]

B = 1×3

23 69 230

C = [0.7 2.8 49]

C = 1×3

0.7000 2.8000 49.0000

A_step = findSmallestStep(A);
B_step = findSmallestStep(B);
C_step = findSmallestStep(C);

disp(['A_step: ', num2str(A_step)]);  %Should be 2.5

A_step: 0.5

disp(['B_step: ', num2str(B_step)]);  %Should be 23

B_step: 1

disp(['C_step: ', num2str(C_step)]);  %Should be 0.7

C_step: 0.1

function stepValue = findSmallestStep(arr)
    % Calculate differences between consecutive elements
    diffs = diff(arr);

    % Convert differences to fractions and find LCM of denominators
    lcm_den = 1;
    for i = 1:length(diffs)
        [~, den] = rat(diffs(i));
        lcm_den = lcm(lcm_den, den);
    end

    % The step is the reciprocal of the LCM of denominators
    stepValue = 1 / lcm_den;
end

[John D'Errico]

Consider the trivial case, where A is entirely integer. Surely it should work properly there?

stepValue = findSmallestStep([0 5 10 20 35 40])

stepValue = 1

And of course, it sees a step of 1. That might not be that terrible, but for larger integers, it still fails.

stepValue = findSmallestStep([0 500000 1000000 2000000 3500000 4000000])

stepValue = 1

And here that very small step seems utterly wrong.

As I pointed out in my answer, this is not as trivial a problem as you might think. You could probably patch the integer case by a pre-test to see if all of the elements are integer. Then you could use the direct gcd code to find the best stride.

Another idea might be to find a good approximation to the stride, then converting the problem to an integer one. For example...

A = [0 2 3 5 7 11]*2/7;

Here by observation, I would argue the stride should be 2/7. How might we find this?

Adiff = diff(A);
Admin = min(Adiff);
Ahat = A/Admin

Ahat = 1×6

0 2 3 5 7 11

stride = double(gcd(sym(Ahat)))*Admin

stride = 0.2857

A problem is, sometimes, that might not result in exact integers. Or maybe there is some case where some of the elements of Ahat are not themselves integer. So we might apply round to the result, before passing it to GCD.

stride = double(gcd(sym(round(Ahat))))*Admin

stride = 0.2857

However, that scheme also fails, against this next case. See that Ahat is not integer.

A = [0 3 5 7 11]*2/7;
Adiff = diff(A);
Admin = min(Adiff);
Ahat = A/Admin

Ahat = 1×5

0 1.5000 2.5000 3.5000 5.5000

And now GCD would fail, and we surely do not want to apply round there. Yes, I could probably still find a patch for that case too. Maybe you could apply findSmallestStep to that vector, in the case where my trick does not result in tegers.

function stepValue = findSmallestStep(arr)
% Calculate differences between consecutive elements
diffs = diff(arr);
% Convert differences to fractions and find LCM of denominators
lcm_den = 1;
for i = 1:length(diffs)
    [~, den] = rat(diffs(i));
    lcm_den = lcm(lcm_den, den);
end
% The step is the reciprocal of the LCM of denominators
stepValue = 1 / lcm_den;
end

You need to look for all of these special cases to write robust code.