How Approximation Search Works

How approximation search works

Approximation search

This is analogy to binary search but without its restrictions that searched function/value/parameter must be strictly monotonic function while sharing the O(log(n)) complexity.

For example Let assume following problem

We have known function y=f(x) and want to find x0 such that y0=f(x0). This can be basically done by inverse function to f but there are many functions that we do not know how to compute inverse to it. So how to compute this in such case?

knowns

y=f(x) - input function
y0 - wanted point y value
a0,a1 - solution x interval range

Unknowns

x0 - wanted point x value must be in range x0=<a0,a1>

Algorithm

probe some points x(i)=<a0,a1> evenly dispersed along the range with some step da
So for example x(i)=a0+i*da where i={ 0,1,2,3... }
for each x(i) compute the distance/error ee of the y=f(x(i))
This can be computed for example like this: ee=fabs(f(x(i))-y0) but any other metrics can be used too.
remember point aa=x(i) with minimal distance/error ee
stop when x(i)>a1
recursively increase accuracy
so first restrict the range to search only around found solution for example:
```
a0'=aa-da;
a1'=aa+da;
```
then increase precision of search by lowering search step:
```
da'=0.1*da;
```
if da' is not too small or if max recursions count is not reached then go to #1
found solution is in aa

This is what I have in mind:

On the left side is the initial search illustrated (bullets #1,#2,#3,#4). On the right side next recursive search (bullet #5). This will recursively loop until desired accuracy is reached (number of recursions). Each recursion increase the accuracy 10 times (0.1*da). The gray vertical lines represent probed x(i) points.

Here the C++ source code for this:

//---------------------------------------------------------------------------
//--- approx ver: 1.01 ------------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _approx_h
#define _approx_h
#include <math.h>
//---------------------------------------------------------------------------
class approx
    {
public:
    double a,aa,a0,a1,da,*e,e0;
    int i,n;
    bool done,stop;

    approx()            { a=0.0; aa=0.0; a0=0.0; a1=1.0; da=0.1; e=NULL; e0=NULL; i=0; n=5; done=true; }
    approx(approx& a)   { *this=a; }
    ~approx()           {}
    approx* operator = (const approx *a) { *this=*a; return this; }
    //approx* operator = (const approx &a) { ...copy... return this; }

    void init(double _a0,double _a1,double _da,int _n,double *_e)
        {
        if (_a0<=_a1) { a0=_a0; a1=_a1; }
        else          { a0=_a1; a1=_a0; }
        da=fabs(_da);
        n =_n ;
        e =_e ;
        e0=-1.0;
        i=0; a=a0; aa=a0;
        done=false; stop=false;
        }
    void step()
        {
        if ((e0<0.0)||(e0>*e)) { e0=*e; aa=a; }         // better solution
        if (stop)                                       // increase accuracy
            {
            i++; if (i>=n) { done=true; a=aa; return; } // final solution
            a0=aa-fabs(da);
            a1=aa+fabs(da);
            a=a0; da*=0.1;
            a0+=da; a1-=da;
            stop=false;
            }
        else{
            a+=da; if (a>a1) { a=a1; stop=true; }       // next point
            }
        }
    };
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

This is how to use it:

approx aa;
double ee,x,y,x0,y0=here_your_known_value;
//            a0,  a1, da,n, ee  
for (aa.init(0.0,10.0,0.1,6,&ee); !aa.done; aa.step())
    {
    x = aa.a;        // this is x(i)
    y = f(x)         // here compute the y value for whatever you want to fit
    ee = fabs(y-y0); // compute error of solution for the approximation search
    }

in the rem above for (aa.init(... are the operand named. The a0,a1 is the interval on which the x(i) is probed, da is initial step between x(i) and n is the number of recursions. so if n=6 and da=0.1 the final max error of x fit will be ~0.1/10^6=0.0000001. The &ee is pointer to variable where the actual error will be computed. I choose pointer so there are not collisions when nesting this and also for speed as passing parameter to heavily used function creates heap trashing.

[notes]

This approximation search can be nested to any dimensionality (but of coarse you need to be careful about the speed) see some examples

Approximation of n points to the curve with the best fit
Curve fitting with y points on repeated x positions (Galaxy Spiral arms)
Increasing accuracy of solution of transcendental equation
Find Minimum area ellipse enclosing a set of points in c++
2D TDoA Time Difference of Arrival
3D TDoA Time Difference of Arrival

In case of non-function fit and the need of getting "all" the solutions you can use recursive subdivision of search interval after solution found to check for another solution. See example:

Given an X co-ordinate, how do I calculate the Y co-ordinate for a point so that it rests on a Bezier Curve

What you should be aware of?

you have to carefully choose the search interval <a0,a1> so it contains the solution but is not too wide (or it would be slow). Also initial step da is very important if it is too big you can miss local min/max solutions or if too small the thing will got too slow (especially for nested multidimensional fits).

Problem nesting approximation search algorithm

Well I ported my original 2D solution into 3D and then tried your values... Looks like I found the problem.

In 3D and 3 receivers there are many positions with the same TDoA values hence the first found is returned (which might not be the one you seek). To verify just print the simulated TDoA values for found solution and compare to your input TDoA values. Another option is to print the final optimization error variable for outer most loop (in your case ax.e0) after computation and see if it is near zero. If it is it means approximation search found solution...

To remedy your case you need at least 4 receivers... so just add one ... Here my updated 3D C++/VCL code:

//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Unit1.h"
#include "backbuffer.h"
#include "approx.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
backbuffer scr;
//---------------------------------------------------------------------------
// TDoA Time Difference of Arrival
// https://stackoverflow.com/a/64318443/2521214
//---------------------------------------------------------------------------
const int n=4;          // number of receivers
double recv[n][4];      // (x,y,z) [m] receiver position,[s] time of arrival of signal
double pos0[3];         // (x,y,z) [m] object's real position
double pos [3];         // (x,y,z) [m] object's estimated position
double v=299800000.0;   // [m/s] speed of signal (light)
double err=0.0;         // [m] error
double aps_e=0.0;       // [m] aprox search error
bool _recompute=true;
//---------------------------------------------------------------------------
void compute()
    {
    int i;
    double x,y,z,a,da,t0;
    //---------------------------------------------------------
    // init positions
    da=2.0*M_PI/(n);
    for (a=0.0,i=0;i<n;i++,a+=da)
        {
        recv[i][0]=2500.0+(2200.0*cos(a));
        recv[i][1]=2500.0+(2200.0*sin(a));
        recv[i][2]=500.0*sin(a);
        }
    // simulate measurement
    t0=123.5;                           // some start time
    for (i=0;i<n;i++)
        {
        x=recv[i][0]-pos0[0];
        y=recv[i][1]-pos0[1];
        z=recv[i][2]-pos0[2];
        a=sqrt((x*x)+(y*y)+(z*z));      // distance to receiver
        recv[i][3]=t0+(a/v);            // start time + time of travel
        }
    //---------------------------------------------------------
    // normalize times into deltas from zero
    a=recv[0][3]; for (i=1;i<n;i++) if (a>recv[i][3]) a=recv[i][3];
    for (i=0;i<n;i++) recv[i][3]-=a;
    // fit position
    int N=8;
    approx ax,ay,az;
    double e,dt[n];
              // min,   max,  step,recursions,&error
     for (ax.init( 0.0,5000.0,200.0,         N,   &e);!ax.done;ax.step())
      for (ay.init( 0.0,5000.0,200.0,         N,   &e);!ay.done;ay.step())
       for (az.init( 0.0,5000.0, 50.0,         N,   &e);!az.done;az.step())
        {
        // simulate measurement -> dt[]
        for (i=0;i<n;i++)
            {
            x=recv[i][0]-ax.a;
            y=recv[i][1]-ay.a;
            z=recv[i][2]-az.a;
            a=sqrt((x*x)+(y*y)+(z*z));  // distance to receiver
            dt[i]=a/v;                  // time of travel
            }
        // normalize times dt[] into deltas from zero
        a=dt[0]; for (i=1;i<n;i++) if (a>dt[i]) a=dt[i];
        for (i=0;i<n;i++) dt[i]-=a;
        // error
        e=0.0; for (i=0;i<n;i++) e+=fabs(recv[i][3]-dt[i]);
        }
    pos[0]=ax.aa;
    pos[1]=ay.aa;
    pos[2]=az.aa;
    aps_e=ax.e0;    // approximation error variable e for best solution
    //---------------------------------------------------------
    // compute error
    x=pos[0]-pos0[0];
    y=pos[1]-pos0[1];
    z=pos[2]-pos0[2];
    err=sqrt((x*x)+(y*y)+(z*z));        // [m]
    }
//---------------------------------------------------------------------------
void draw()
    {
    scr.cls(clBlack);
    int i;
    const double pan_x=0.0;
    const double pan_y=0.0;
    const double zoom=512.0/5000.0;
    double x,y,r=8.0;

    #define tabs(x,y){ x-=pan_x; x*=zoom; y-=pan_y; y*=zoom; }
    #define trel(x,y){ x*=zoom; y*=zoom; }

    scr.bmp->Canvas->Font->Color=clYellow;
    scr.bmp->Canvas->TextOutA(5, 5,AnsiString().sprintf("Error: %.3lf m",err));
    scr.bmp->Canvas->TextOutA(5,25,AnsiString().sprintf("Aprox error: %.20lf s",aps_e));
    scr.bmp->Canvas->TextOutA(5,45,AnsiString().sprintf("pos0 %6.1lf %6.1lf %6.1lf m",pos0[0],pos0[1],pos0[2]));
    scr.bmp->Canvas->TextOutA(5,65,AnsiString().sprintf("pos  %6.1lf %6.1lf %6.1lf m",pos [0],pos [1],pos [2]));

    // recv
    scr.bmp->Canvas->Pen->Color=clAqua;
    scr.bmp->Canvas->Brush->Color=clBlue;
    for (i=0;i<n;i++)
        {
        x=recv[i][0];
        y=recv[i][1];
        tabs(x,y);
        scr.bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);
        }
    // pos0
    scr.bmp->Canvas->Pen->Color=TColor(0x00202060);
    scr.bmp->Canvas->Brush->Color=TColor(0x00101040);
    x=pos0[0];
    y=pos0[1];
    tabs(x,y);
    scr.bmp->Canvas->Ellipse(x-r,y-r,x+r,y+r);

    // pos
    scr.bmp->Canvas->Pen->Color=clYellow;
    x=pos[0];
    y=pos[1];
    tabs(x,y);
    scr.bmp->Canvas->MoveTo(x-r,y-r);
    scr.bmp->Canvas->LineTo(x+r,y+r);
    scr.bmp->Canvas->MoveTo(x+r,y-r);
    scr.bmp->Canvas->LineTo(x-r,y+r);

    scr.rfs();

    #undef tabs(x,y){ x-=pan_x; x*=zoom; y-=pan_y; y*=zoom; }
    #undef trel(x,y){ x*=zoom; y*=zoom; }
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    Randomize();
    pos0[0]=2000.0;
    pos0[1]=2000.0;
    pos0[2]= 900.0;
    scr.set(this);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::tim_updateTimer(TObject *Sender)
    {
    if (_recompute)
        {
        compute();
        _recompute=false;
        }
    draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormClick(TObject *Sender)
    {
    pos0[0]=scr.mx1*5000.0/512.0;
    pos0[1]=scr.my1*5000.0/512.0;
    pos0[2]=Random()*1000.0;
    _recompute=true;
    }
//---------------------------------------------------------------------------

So as before just ignore the VCL stuff and port to your environment/language... Its configured so it does not compute for too long (less than 1 sec) and the error is <0.01 m

Here preview:

preview

also with printed the approximation error variable final value (abs sum of difference between solution's and input TDoA times in seconds) ... and pos0,pos values...

Beware even this is not fully safe as it has blind spots ... The receivers should not be at the same altitude and also maybe even not distributed evenly on circle because that can cause the TDoA duplicates...

In case you got additional info (like know that pos0 is at ground and have the terrain map) then it might be possible to do this with 3 receivers where the z coordinate will not be approximated anymore but computed from x,y instead...

PS. your nesting has a slight cosmetics "bug" as you have:

          ay.e = error; ax.e = error; az.e = error
          az.step()

        ay.step()

    ax.step()

I would feel much safer with (not sure though as I do not code in Python):

          az.e = error
          az.step()

        ay.e = error
        ay.step()

    ax.e = error
    ax.step()

Is Approximate String Matching / Fuzzy String Searching possible with BigQuery?

Unfortunately, approximate string matching is not supported. The closest you can get is by using regular expressions. Your best bet may be to normalize the data before it gets to BigQuery -- i.e transform "Rhodes USA" and "Rhodes, USA. " into the same string. I'll add a feature request bug for this support, however.

Big O, how do you calculate/approximate it?

I'll do my best to explain it here on simple terms, but be warned that this topic takes my students a couple of months to finally grasp. You can find more information on the Chapter 2 of the Data Structures and Algorithms in Java book.

There is no mechanical procedure that can be used to get the BigOh.

As a "cookbook", to obtain the BigOh from a piece of code you first need to realize that you are creating a math formula to count how many steps of computations get executed given an input of some size.

The purpose is simple: to compare algorithms from a theoretical point of view, without the need to execute the code. The lesser the number of steps, the faster the algorithm.

For example, let's say you have this piece of code:

int sum(int* data, int N) {
    int result = 0;               // 1

    for (int i = 0; i < N; i++) { // 2
        result += data[i];        // 3
    }

    return result;                // 4
}

This function returns the sum of all the elements of the array, and we want to create a formula to count the computational complexity of that function:

Number_Of_Steps = f(N)

So we have f(N), a function to count the number of computational steps. The input of the function is the size of the structure to process. It means that this function is called such as:

Number_Of_Steps = f(data.length)

The parameter N takes the data.length value. Now we need the actual definition of the function f(). This is done from the source code, in which each interesting line is numbered from 1 to 4.

There are many ways to calculate the BigOh. From this point forward we are going to assume that every sentence that doesn't depend on the size of the input data takes a constant C number computational steps.

We are going to add the individual number of steps of the function, and neither the local variable declaration nor the return statement depends on the size of the data array.

That means that lines 1 and 4 takes C amount of steps each, and the function is somewhat like this:

f(N) = C + ??? + C

The next part is to define the value of the for statement. Remember that we are counting the number of computational steps, meaning that the body of the for statement gets executed N times. That's the same as adding C, N times:

f(N) = C + (C + C + ... + C) + C = C + N * C + C

There is no mechanical rule to count how many times the body of the for gets executed, you need to count it by looking at what does the code do. To simplify the calculations, we are ignoring the variable initialization, condition and increment parts of the for statement.

To get the actual BigOh we need the Asymptotic analysis of the function. This is roughly done like this:

Take away all the constants C.
From f() get the polynomium in its standard form.
Divide the terms of the polynomium and sort them by the rate of growth.
Keep the one that grows bigger when N approaches infinity.

Our f() has two terms:

f(N) = 2 * C * N ^ 0 + 1 * C * N ^ 1

Taking away all the C constants and redundant parts:

f(N) = 1 + N ^ 1

Since the last term is the one which grows bigger when f() approaches infinity (think on limits) this is the BigOh argument, and the sum() function has a BigOh of:

O(N)

There are a few tricks to solve some tricky ones: use summations whenever you can.

As an example, this code can be easily solved using summations:

for (i = 0; i < 2*n; i += 2) {  // 1
    for (j=n; j > i; j--) {     // 2
        foo();                  // 3
    }
}

The first thing you needed to be asked is the order of execution of foo(). While the usual is to be O(1), you need to ask your professors about it. O(1) means (almost, mostly) constant C, independent of the size N.

The for statement on the sentence number one is tricky. While the index ends at 2 * N, the increment is done by two. That means that the first for gets executed only N steps, and we need to divide the count by two.

f(N) = Summation(i from 1 to 2 * N / 2)( ... ) = 
     = Summation(i from 1 to N)( ... )

The sentence number two is even trickier since it depends on the value of i. Take a look: the index i takes the values: 0, 2, 4, 6, 8, ..., 2 * N, and the second for get executed: N times the first one, N - 2 the second, N - 4 the third... up to the N / 2 stage, on which the second for never gets executed.

On formula, that means:

f(N) = Summation(i from 1 to N)( Summation(j = ???)(  ) )

Again, we are counting the number of steps. And by definition, every summation should always start at one, and end at a number bigger-or-equal than one.

f(N) = Summation(i from 1 to N)( Summation(j = 1 to (N - (i - 1) * 2)( C ) )

(We are assuming that foo() is O(1) and takes C steps.)

We have a problem here: when i takes the value N / 2 + 1 upwards, the inner Summation ends at a negative number! That's impossible and wrong. We need to split the summation in two, being the pivotal point the moment i takes N / 2 + 1.

f(N) = Summation(i from 1 to N / 2)( Summation(j = 1 to (N - (i - 1) * 2)) * ( C ) ) + Summation(i from 1 to N / 2) * ( C )

Since the pivotal moment i > N / 2, the inner for won't get executed, and we are assuming a constant C execution complexity on its body.

Now the summations can be simplified using some identity rules:

Summation(w from 1 to N)( C ) = N * C
Summation(w from 1 to N)( A (+/-) B ) = Summation(w from 1 to N)( A ) (+/-) Summation(w from 1 to N)( B )
Summation(w from 1 to N)( w * C ) = C * Summation(w from 1 to N)( w ) (C is a constant, independent of w)
Summation(w from 1 to N)( w ) = (N * (N + 1)) / 2

Applying some algebra:

f(N) = Summation(i from 1 to N / 2)( (N - (i - 1) * 2) * ( C ) ) + (N / 2)( C )

f(N) = C * Summation(i from 1 to N / 2)( (N - (i - 1) * 2)) + (N / 2)( C )

f(N) = C * (Summation(i from 1 to N / 2)( N ) - Summation(i from 1 to N / 2)( (i - 1) * 2)) + (N / 2)( C )

f(N) = C * (( N ^ 2 / 2 ) - 2 * Summation(i from 1 to N / 2)( i - 1 )) + (N / 2)( C )

=> Summation(i from 1 to N / 2)( i - 1 ) = Summation(i from 1 to N / 2 - 1)( i )

f(N) = C * (( N ^ 2 / 2 ) - 2 * Summation(i from 1 to N / 2 - 1)( i )) + (N / 2)( C )

f(N) = C * (( N ^ 2 / 2 ) - 2 * ( (N / 2 - 1) * (N / 2 - 1 + 1) / 2) ) + (N / 2)( C )

=> (N / 2 - 1) * (N / 2 - 1 + 1) / 2 = 

   (N / 2 - 1) * (N / 2) / 2 = 

   ((N ^ 2 / 4) - (N / 2)) / 2 = 

   (N ^ 2 / 8) - (N / 4)

f(N) = C * (( N ^ 2 / 2 ) - 2 * ( (N ^ 2 / 8) - (N / 4) )) + (N / 2)( C )

f(N) = C * (( N ^ 2 / 2 ) - ( (N ^ 2 / 4) - (N / 2) )) + (N / 2)( C )

f(N) = C * (( N ^ 2 / 2 ) - (N ^ 2 / 4) + (N / 2)) + (N / 2)( C )

f(N) = C * ( N ^ 2 / 4 ) + C * (N / 2) + C * (N / 2)

f(N) = C * ( N ^ 2 / 4 ) + 2 * C * (N / 2)

f(N) = C * ( N ^ 2 / 4 ) + C * N

f(N) = C * 1/4 * N ^ 2 + C * N

And the BigOh is:

O(N²)

Unit Testing Approximation Algorithms

You need to separate two concerns here. The quality of your approximation algorithms and the correctness of implementation of those algorithms.

Testing the quality of an approximation algorithm usually will not lend itself to unit testing methods used in software development. For example you would need to generate random problems that is representative of the real sizes of problems. You might need to do mathematical work to get some upper/lower bound to judge the quality of your algorithms for unsolvable large instances. Or use problem test sets that have known or best known solutions and compare your results. But in any case unit testing would not help you much in improving the quality of the approximation algorithms. This is where your domain knowledge in optimization and math will help.

The correctness of your implementation is where unit tests will be really useful. You can use toy sized problems here and compare known results (solving by hand, or verified through careful step by step debugging in code) with what your code generates. Having small problems is not only enough but also desirable here so that tests run fast and can be run many times during development cycle. These types of tests makes sure that overall algorithm is arriving at the correct result. It is somewhere between a unit test and an integration tests since you are testing a large portion of the code as a black box. But I have found these types of tests to be extremely useful in optimization domain. One thing I recommend doing for this type of testing is removing all randomness in your algorithms through fixed seeds for random number generators. These tests should always run in a deterministic way and give exactly the same result 100% of the time.
I also recommend unit testing at the lower level modules of your algorithms. Isolate that method that assigns weights to arcs on the graph and check if the correct weights are assigned. Isolate your objective function value calculation function and unit test that. You get my point.

One other concern that cuts both of these slices is performance. You cannot reliably test performance with small toy problems. Also realizing a change that degrades performance significantly for a working algorithm quickly is very desirable. Once you have a running version of your algorithms you can create larger test problems where you measure the performance and automate it to be your performance/integration tests. You can run these less frequently as they will take more time but at least will notify you early of newly introduced performance bottlenecks during refactoring or new feature additions to algorithms

Approximate Matching and Replacement in Text R

Some agrep examples:

agrep("lasy", "1 lazy 2")
agrep("lasy", "1 lazy 2", max = list(sub = 0))
agrep("laysy", c("1 lazy", "1", "1 LAZY"), max = 2)
agrep("laysy", c("1 lazy", "1", "1 LAZY"), max = 2, value = TRUE)
agrep("laysy", c("1 lazy", "1", "1 LAZY"), max = 2, ignore.case = TRUE)

agrep is in base. If you load stringdist you can calculate string-distance using Jarro-Winkler with (you guessed) stringdist or if you're lazy you can just use ain or amatch. For my purposes, I tend to use Damerau–Levenshtein (method="dl") more, but your mileage might vary.

Just make sure to read up on exactly how the algorithm's parameters work before you use it (i.e. set your p, q and maxDist values to levels that make sense given what you're doing)

Approximate a curve with a limited number of line segments and arcs of circles

the C1 requirement demands the you must have alternating straights and arcs. Also realize if you permit a sufficient number of segments you can trivially fit every pair of points with a straight and use a tiny arc to satisfy slope continuity.

I'd suggest this algorithm,

1 best fit with a set of (specified N) straight segments. (surely there are well developed algorithms for that.)

2 consider the straight segments fixed and at each joint place an arc. Treating each joint individually i think you have a tractable problem to find the optimum arc center/radius to satisfy continuity and improve the fit.

3 now that you are pretty close attempt to consider all arc centers and radii (segments being defined by tangency) as a global optimization problem. This of course blows up if N is large.

How Approximation Search Works