近似搜索的工作原理
How approximation search works
[序言]
此Q&A 旨在更清楚地解释我首次在此处发布的近似搜索class 的内部工作原理
- Increasing accuracy of solution of transcendental equation
已经有几次(出于各种原因)要求我提供有关这方面的更多详细信息,所以我决定就此编写 Q&A 风格的主题,以便将来可以轻松参考不需要再解释一遍
[问题]
如何在实域 (double
) 中逼近 values/parameters 以实现多项式、参数函数的拟合或求解(困难的)方程(如超越方程)?
限制
- 真实域(
double
精度)
- C++ 语言
- 可配置的近似精度
- 搜索的已知间隔
- 拟合 value/parameter 不是严格单调的或根本不起作用
近似搜索
这类似于二分查找,但没有其限制,搜索 function/value/parameter 必须严格单调函数,同时共享 O(log(n))
复杂性。
例如假设以下问题
我们已经知道函数 y=f(x)
并且想要找到 x0
使得 y0=f(x0)
。这基本上可以通过 f
的反函数来完成,但是有很多函数我们不知道如何计算它的反函数。那么在这种情况下如何计算呢?
知名度
y=f(x)
- 输入函数
y0
- 通缉点数 y
值
a0,a1
-求解x
区间范围
未知数
x0
- 通缉点 x
值必须在 x0=<a0,a1>
范围内
算法
探测一些点 x(i)=<a0,a1>
沿着范围均匀分布一些点 da
例如 x(i)=a0+i*da
其中 i={ 0,1,2,3... }
对每个x(i)
计算y=f(x(i))
[=64=的distance/erroree
]
这可以这样计算,例如:ee=fabs(f(x(i))-y0)
但也可以使用任何其他指标。
记住点aa=x(i)
带最小distance/erroree
在x(i)>a1
时停止
递归提高准确性
所以首先将范围限制为仅在找到的解决方案周围搜索,例如:
a0'=aa-da;
a1'=aa+da;
然后通过降低搜索步长来提高搜索精度:
da'=0.1*da;
如果 da'
不是太小或者如果未达到最大递归计数则转到 #1
找到的解决方案在aa
这是我的想法:
左侧是图示的初始搜索(项目符号 #1、#2、#3、#4)。在右边下一个递归搜索(bullet #5)。这将递归循环,直到达到所需的精度(递归次数)。每次递归都会将准确度提高 10
倍 (0.1*da
)。灰色垂直线代表探测的 x(i)
个点。
这里是 C++ 源代码:
//---------------------------------------------------------------------------
//--- 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
//---------------------------------------------------------------------------
使用方法如下:
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
}
上面的rem中for (aa.init(...
是命名的操作数。 a0,a1
是探测 x(i)
的区间,da
是 x(i)
之间的初始步骤,n
是递归的次数。因此,如果 n=6
和 da=0.1
,则 x
拟合的最终最大误差将为 ~0.1/10^6=0.0000001
。 &ee
是指向将计算实际错误的变量的指针。我选择指针,这样嵌套时就不会发生冲突,而且速度也是如此,因为将参数传递给频繁使用的函数会造成堆垃圾。
[备注]
这个近似搜索可以嵌套到任何维度(但粗略的你需要注意速度)看一些例子
- Approximation of n points to the curve with the best fit
- Increasing accuracy of solution of transcendental equation
- Find Minimum area ellipse enclosing a set of points in c++
在非函数拟合和需要获得“所有”解决方案的情况下,您可以在找到解决方案后使用搜索间隔的递归细分来检查另一个解决方案。参见示例:
- Given an X co-ordinate, how do I calculate the Y co-ordinate for a point so that it rests on a Bezier Curve
你应该注意什么?
你必须仔细选择搜索区间 <a0,a1>
以便它包含解决方案但不要太宽(否则会很慢)。同样,初始步骤 da
非常重要,如果它太大,您可能会错过局部 min/max 解决方案,或者如果太小,事情就会变得太慢(特别是对于嵌套的多维拟合)。
结合secant (with bracketing, but see correction at the bottom) and bisection方法更好:
我们通过割线找到根近似值,并像二分法一样将根括在括号中。
始终保持区间的两条边,使得一条边的delta为负,另一条边为正,保证根在里面;而不是减半,使用正割法。
伪代码:
given a function f
given two points a, b, such that a < b and sign(f(a)) /= sign(f(b))
given tolerance tol
find root z of f such that abs(f(z)) < tol -- stop_condition
DO:
x = root of f by linear interpolation of f between a and b
m = midpoint between a and b
if stop_condition holds at x or m, set z and STOP
[a,b] := [a,x,m,b].sort.choose_shortest_interval_with_
_opposite_signs_at_its_ends
这显然在每次迭代时将间隔 [a,b]
减半,甚至更好;因此,除非该函数的行为非常糟糕(例如,sin(1/x)
接近 x=0
),否则它将很快收敛,对于每个迭代步骤,最多只对 f
进行两次评估.
我们可以通过检查 b-a
不会变得太小来检测不良行为情况(特别是如果我们使用有限精度工作,如双精度)。
update: apparently this is actually double false position 方法,如上面的伪代码所述,它是带括号的割线。即使在最病态的情况下,像二分法一样通过中间点增加它也能确保收敛。
[序言]
此Q&A 旨在更清楚地解释我首次在此处发布的近似搜索class 的内部工作原理
- Increasing accuracy of solution of transcendental equation
已经有几次(出于各种原因)要求我提供有关这方面的更多详细信息,所以我决定就此编写 Q&A 风格的主题,以便将来可以轻松参考不需要再解释一遍
[问题]
如何在实域 (double
) 中逼近 values/parameters 以实现多项式、参数函数的拟合或求解(困难的)方程(如超越方程)?
限制
- 真实域(
double
精度) - C++ 语言
- 可配置的近似精度
- 搜索的已知间隔
- 拟合 value/parameter 不是严格单调的或根本不起作用
近似搜索
这类似于二分查找,但没有其限制,搜索 function/value/parameter 必须严格单调函数,同时共享 O(log(n))
复杂性。
例如假设以下问题
我们已经知道函数 y=f(x)
并且想要找到 x0
使得 y0=f(x0)
。这基本上可以通过 f
的反函数来完成,但是有很多函数我们不知道如何计算它的反函数。那么在这种情况下如何计算呢?
知名度
y=f(x)
- 输入函数y0
- 通缉点数y
值a0,a1
-求解x
区间范围
未知数
x0
- 通缉点x
值必须在x0=<a0,a1>
范围内
算法
探测一些点
x(i)=<a0,a1>
沿着范围均匀分布一些点da
例如
x(i)=a0+i*da
其中i={ 0,1,2,3... }
对每个
x(i)
计算y=f(x(i))
[=64=的distance/erroree
]这可以这样计算,例如:
ee=fabs(f(x(i))-y0)
但也可以使用任何其他指标。记住点
aa=x(i)
带最小distance/erroree
在
时停止x(i)>a1
递归提高准确性
所以首先将范围限制为仅在找到的解决方案周围搜索,例如:
a0'=aa-da; a1'=aa+da;
然后通过降低搜索步长来提高搜索精度:
da'=0.1*da;
如果
da'
不是太小或者如果未达到最大递归计数则转到 #1找到的解决方案在
aa
这是我的想法:
左侧是图示的初始搜索(项目符号 #1、#2、#3、#4)。在右边下一个递归搜索(bullet #5)。这将递归循环,直到达到所需的精度(递归次数)。每次递归都会将准确度提高 10
倍 (0.1*da
)。灰色垂直线代表探测的 x(i)
个点。
这里是 C++ 源代码:
//---------------------------------------------------------------------------
//--- 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
//---------------------------------------------------------------------------
使用方法如下:
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
}
上面的rem中for (aa.init(...
是命名的操作数。 a0,a1
是探测 x(i)
的区间,da
是 x(i)
之间的初始步骤,n
是递归的次数。因此,如果 n=6
和 da=0.1
,则 x
拟合的最终最大误差将为 ~0.1/10^6=0.0000001
。 &ee
是指向将计算实际错误的变量的指针。我选择指针,这样嵌套时就不会发生冲突,而且速度也是如此,因为将参数传递给频繁使用的函数会造成堆垃圾。
[备注]
这个近似搜索可以嵌套到任何维度(但粗略的你需要注意速度)看一些例子
- Approximation of n points to the curve with the best fit
- Increasing accuracy of solution of transcendental equation
- Find Minimum area ellipse enclosing a set of points in c++
在非函数拟合和需要获得“所有”解决方案的情况下,您可以在找到解决方案后使用搜索间隔的递归细分来检查另一个解决方案。参见示例:
- Given an X co-ordinate, how do I calculate the Y co-ordinate for a point so that it rests on a Bezier Curve
你应该注意什么?
你必须仔细选择搜索区间 <a0,a1>
以便它包含解决方案但不要太宽(否则会很慢)。同样,初始步骤 da
非常重要,如果它太大,您可能会错过局部 min/max 解决方案,或者如果太小,事情就会变得太慢(特别是对于嵌套的多维拟合)。
结合secant (with bracketing, but see correction at the bottom) and bisection方法更好:
我们通过割线找到根近似值,并像二分法一样将根括在括号中。
始终保持区间的两条边,使得一条边的delta为负,另一条边为正,保证根在里面;而不是减半,使用正割法。
伪代码:
given a function f
given two points a, b, such that a < b and sign(f(a)) /= sign(f(b))
given tolerance tol
find root z of f such that abs(f(z)) < tol -- stop_condition
DO:
x = root of f by linear interpolation of f between a and b
m = midpoint between a and b
if stop_condition holds at x or m, set z and STOP
[a,b] := [a,x,m,b].sort.choose_shortest_interval_with_
_opposite_signs_at_its_ends
这显然在每次迭代时将间隔 [a,b]
减半,甚至更好;因此,除非该函数的行为非常糟糕(例如,sin(1/x)
接近 x=0
),否则它将很快收敛,对于每个迭代步骤,最多只对 f
进行两次评估.
我们可以通过检查 b-a
不会变得太小来检测不良行为情况(特别是如果我们使用有限精度工作,如双精度)。
update: apparently this is actually double false position 方法,如上面的伪代码所述,它是带括号的割线。即使在最病态的情况下,像二分法一样通过中间点增加它也能确保收敛。