Spektre - 7 months ago 54
C++ Question

# How approximation search works

[Prologue]

This Q&A is meant to explain more clearly the inner working of my approximations search class which I first published here

I was requested for more detailed info about this few times already (for various reasons) so I decided to write Q&A style topic about this which I can easily reference in the future and do not need to explain it over and over again.

[Question]

How to approximate values/parameters in Real domain (
`double`
) to achieve fitting of polynomials,parametric functions or solve (difficult) equations (like transcendental) ?

Restrictions

• real domain (
`double`
precision)

• C++ language

• configurable precision of approximation

• known interval for search

• fitted value/parameter is not strictly monotonic or not function at all

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(n.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

1. 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... }`

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

3. remember point `aa=x(i)` with minimal distance/error `ee`

4. stop when `x(i)>a1`

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

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

[notes]

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

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:

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

Source (Stackoverflow)