Finding the closest integer fraction to a given random real between 0..1, given ranges of numerator and denominator

Question

Given two ranges of positive integers x: [1 ... n] and y: [1 ... m] and random real R from 0 to 1, I need to find the pair of elements (i,j) from x and y such that x_i / y_j is closest to R.

What is the most efficient way to find this pair?

Vouchsafe

What do you have so far? – Pane

I'm keeping Xi fixed and getting the closest Yi. I'm finding I'm not close enough. I know I can get closer by stepping Xi up and down and seeing what I get, but that seems gross. – Vouchsafe

It seems easy on first glance, but I think it might be hard. If there is not a perfect solution like 1/2 = .5, there may be multiple correct answers. Actually I guess in that case there are also multiple answers like 2/4. In the case where there are multiple answers, I want the biggest Xi and Yi in the range. – Vouchsafe

See this answer – Chukchi

Are x[] and y[] a list/array of numbers or a range of numbers? – Cline

Could this be posted on Mathematics version of StackOverflow...Especially regard some of the comments by Moron and YBungalow.... – Rms

possible duplicate of How can I turn a floating point number into the closest fraction represented by a byte numerator and denominator? – Caste

Answer 1

Using Farey sequence

This is a simple and mathematically beautiful algorithm to solve this: run a binary search, where on each iteration the next number is given by the mediant formula (below). By the properties of the Farey sequence that number is the one with the smallest denominator within that interval. Consequently this sequence will always converge and never 'miss' a valid solution.

In pseudocode:

input: m, n, R

a_num = 0, a_denom = 1
b_num = 1, b_denom = 1

repeat:
    -- interestingly c_num/c_denom is already in reduced form
    c_num = a_num + b_num
    c_denom = a_denom + b_denom

    -- if the numbers are too big, return the closest of a and b
    if c_num > n or c_denom > m then
        if R - a_num/a_denom < b_num/b_denom - R then
            return a_num, a_denom
        else
            return b_num, b_denom

    -- adjust the interval:
    if c_num/c_denom < R then
        a_num = c_num, a_denom = c_denom
    else
        b_num = c_num, b_denom = c_denom

goto repeat

Even though it's fast on average (my educated guess that it's O(log max(m,n))), it can still be slow if R is close to a fraction with a small denominator. For example finding an approximation to 1/1000000 with m = n = 1000000 will take a million iterations.

Saffron

-1: Why would you even expect this to work? Remember, the numerator and denominators are restricted. – Rosalinarosalind

@moron: it actually does work for restricted denominators and for numbers between 0 and 1 (which is the case here). in fact, i'm pretty sure it can be adapted for restricted numerators (which restrict the denominators, given the number). – Baalbek

I'm not sure about the O(log M) average case though, don't have time to analyze right now. – Saffron

@John: x = [5], y = [8], R = 3/5. This outputs 1 and stops (in step 3) which is not even a feasible solution. – Rosalinarosalind

I accepted because it looked the same as this - johndcook.com/blog/2010/10/20/best-rational-approximation/… – Vouchsafe

@John: I presume x and y are arrays of arbitrary positive numbers, while this answer assumes 1 <=x <= n and 1 <= y <= m. Which is it? Axn also had this question. – Rosalinarosalind

@John Shedletsky, I solve it exactly in better time and space order, think about algorithms instead of goggling similar ones. – Alikee

@ybung: Then I don't see the point of specifying x_i/y_j instead of just i/j in the question. Anyway, John will clarify I suppose. Probably my mistake. I agree, that with your interpretation, Farey sequences will give the right answer, not sure about your average case runtime though. – Rosalinarosalind

@DrXorile: obviously just using farey sequences won't give you the correct anwers. You also need to get the algorithm right. The code in the article is incorrect. Just run my pseudocode and got 17/28. You are welcome to find the differences. – Saffron

I implemented your algorithm and it does indeed return the correct result. Apologies. I deleted my comment. – Rickierickman

This algorithm works for a restricted range of denominators, because "the next Farey fraction between a and b" is the fraction with the smallest denominator which is strictly in between a and b. Hence this is the fraction which has the minimal denominator but improves the approximation. But why should it work for a given range of numerators? – Confer

@Echsecutor: because both increase monotonically, so when the first of them exceeds the bound then there's no point in looking further. – Saffron

Answer 2

The standard approach to approximating reals with rationals is computing the continued fraction series (see [1]). Put a limit on the nominator and denominator while computing parts of the series, and the last value before you break the limits is a fraction very close to your real number.

This will find a very good approximation very fast, but I'm not sure this will always find a closest approximation. It is known that

any convergent [partial value of the continued fraction expansion] is nearer to the continued fraction than any other fraction whose denominator is less than that of the convergent

but there may be approximations with larger denominator (still below your limit) that are better approximations, but are not convergents.

[1] http://en.wikipedia.org/wiki/Continued_fraction

Mesentery

I may be misunderstanding - I don't want a continued fraction as the answer, I want a single numerator and denominator. Are you saying that if I find the continued fraction then I have some sort of guarantee of optimality on a more simplified fraction? – Vouchsafe

What you probably want are the "Best rational approximations" (on the wikipedia page for the continued fraction), which is either a convergent to the continued fraction or has the final quotient of one of the convergents decreased by one. – Juster

Continued fractions do yield a rational approximation (by using the last converged with small enough numerator/denominator). But why should this be the best approximation to R in the given range of nominators/denominatirs? – Confer

Answer 3

Given that R is a real number such that 0 <= R <= 1, integers x: [1 ... n] and integers y: [1 ... m]. It is assumed that n <= m, since if n > m then x[n]/y[m] will be greater than 1, which cannot be the closest approximation to R.

Therefore, the best approximation of R with the denominator d will be either floor(R*d) / d or ceil(R*d) / d.

The problem can be solved in O(m) time and O(1) space (in Python):

from __future__ import division
from random import random
from math import floor

def fractionize(R, n, d):
    error = abs(n/d - R)
    return (n, d, error)  # (numerator, denominator, absolute difference to R)

def better(a, b):
    return a if a[2] < b[2] else b

def approximate(R, n, m):
    best = (0, 1, R)
    for d in xrange(1, m+1):
        n1 = min(n, int(floor(R * d)))
        n2 = min(n, n1 + 1) # ceil(R*d)
        best = better(best, fractionize(R, n1, d))
        best = better(best, fractionize(R, n2, d))
    return best

if __name__ == '__main__': 
    def main():
        R = random()
        n = 30
        m = 100
        print R, approximate(R, n, m)
    main()

Quijano

brute force is not always the best algorithm ;) – Confer

Answer 4

Prolly get flamed, but a lookup might be best where we compute all of the fractional values for each of the possible values.. So a simply indexing a 2d array indexed via the fractional parts with the array element containing the real equivalent. I guess we have discrete X and Y parts so this is finite, it wouldnt be the other way around.... Ahh yeah, the actual searching part....erm reet....

Rms

In my particular application, n and m are around 100,000. This makes pre-computation undesirable. I was hoping for some sort of hillclimb optimization. – Vouchsafe

Answer 5

The Solution: You can do this O(1) space and O(m log(n)) time:

there is no need to create any list to search,

The pseudo code may be is buggy but the idea is this:

r: input number to search.
n,m: the ranges.

for (int i=1;i<=m;i++)
{
    minVal = min(Search(i,1,n,r), minVal);
}

//x and y are start and end of array:
decimal Search(i,x,y,r)
{
   if (i/x > r)
      return i/x - r;

   decimal middle1 = i/Cill((x+y)/2); 
   decimal middle2 = i/Roof((x+y)/2);

   decimal dist = min(middle1,middle2)

   decimal searchResult = 100000;

   if( middle > r)
     searchResult = Search (i, x, cill((x+y)/2),r)
  else
     searchResult = Search(i, roof((x+y)/2), y,r)

  if  (searchResult < dist)
     dist = searchResult;

  return dist;
}

finding the index as home work to reader.

Description: I think you can understand what's the idea by code, but let trace one of a for loop: when i=1:

you should search within bellow numbers: 1,1/2,1/3,1/4,....,1/n you check the number with (1,1/cill(n/2)) and (1/floor(n/2), 1/n) and doing similar binary search on it to find the smallest one.

Should do this for loop for all items, so it will be done m time. and in each time it takes O(log(n)). this function can improve by some mathematical rules, but It will be complicated, I skip it.

Alikee

Any clever optimizations to do better than O(nm) space and O(nm lg (nm)) time? – Vouchsafe

No it is not. Especially not without proof. – Rosalinarosalind

@Moron, you want proof what? The algorithm as described above run in the specified order, and will get the best answer, for example for binary search you saying the proof, it finds exact match? no because the algorithm describes the trust, about the order, it's easy to proof it, if is there any ambiguity tell to describe it. – Alikee

I was responding to your comment to john. Not about your answer. – Rosalinarosalind

Answer 6

Rather than a completely brute force search, do a linear search over the shortest of your lists, using round to find the best match for each element. Maybe something like this:

best_x,best_y=(1,1)
for x in 1...n:
    y=max(1,min(m,round(x/R)))
    #optional optimization (if you have a fast gcd)
    if gcd(x,y)>1:
        continue

    if abs(R-x/y)<abs(R-bestx/besty):
        best_x,best_y=(x,y)
return (best_x,best_y)

Not at all sure whether the gcd "optimization" will ever be faster...

Directrix

How is this not "completely brute force"? – Confer

Answer 7

If the denominator of R is larger than m then use the Farey method (which the Fraction.limit_denominator method implements) with a limit of m to get a fraction a/b where b is smaller than m else let a/b = R. With b <= m, either a <= n and you are done or else let M = math.ceil(n/R) and re-run the Farey method.

def approx2(a, b, n, m):
    from math import ceil
    from fractions import Fraction
    R = Fraction(a, b)
    if R < Fraction(1, m):
        return 1, m
    r = R.limit_denominator(m)
    if r.numerator > n:
        M = ceil(n/R)
        r = R.limit_denominator(M)
    return r.numerator, r.denominator

>>> approx2(113, 205, 50, 200)
(43, 78)

It might be possible to just run the Farey method once using a limiting denominator of min(ceil(n/R), m) but I am not sure about that:

def approx(a, b, n, m):
    from math import ceil
    from fractions import Fraction
    R = Fraction(a, b)
    if R < Fraction(1, m):
        return 1, m
    r = R.limit_denominator(min(ceil(n/R), m))
    return r.numerator, r.denominator

Brindabrindell