www.cnblogs.com

回应CSDN肖舸《做程序,要“专注”和“客观”》,实验比较各离散采样算法

2010-05-27 23:57 by Milo Yip, ... 阅读, ... 评论, 收藏, 编辑

自从肖舸在其CSDN博客上说“拒绝回答博客园等网站网友的问题”,实质上不单是拒绝回答,而且还删去包括一些网友及本人对于纯粹技术探讨的评论。当然每位博主都有自由这么做,但个人认为这对于社区的交流发展有负面影响。为了探讨这个技术问题,本人唯有把回应发表于博客园内。本文会阐述之前的论点,评论肖舸的实现,并进行了兩个实验比较不同算法、实现的优劣之处。

之前的“交流”

约一个月前,肖舸于《实际中常用的一个随机数产生器(分类别概率随机)》(下简称《实》)一文中,介绍了一个宣称“0bug”的实现。而实质上这实现至少有两个bug,其中一个已被网友发现,而肖舸也随后更新,不过没有公开向该网友致谢,连网友的id或名字也没写。本人还发现另一个会做成崩溃bug,稍后再说。

因为有感《实》欠缺背后理论,而且该实现的需求不清楚、代码难读,所以本人于5日后撰写博文《用JavaScript玩转游戏编程(一)掉宝类型概率》(下简称《用》 ),从统计学解释这问题,并以JavaScript实作互动的示范程序。文中提及:

这里用了线性搜寻(linear search),……另外,也可以用二分搜寻(binary search),那么复杂度会降低为O(lg N),……那么,还有没有更快的方法呢?答案是肯定的,例如别名方法(alias method)、近似方法等,有兴趣的读者可参考……当然,在N很小的情况下,线性搜寻和二分搜寻也足够。
笔者撰写本文,灵感来自这篇博文(指《实》)。其算法实际上是储存CDF的逆函数采样,利用空间和有限的CDF精确度,换取O(1)的时间复杂度。衡量N的大小、精确度、空间需求、缓存延迟后,或许该方法也能适合某些个别需求。但对于该文作者说N最大为100,二分搜寻只需最多7次迭代,因缓存问题可能二分搜寻更快。

之后,肖舸在另一篇博文《做程序,要“专注”和“客观”》中,当中有一段似乎是回应上文:

就好比我前面写的一篇博文《实》,我在文中的代码里,明明实现了O(1)的复杂度,但是就有人,为了攻击我本人和我的书《……》,专门撰文,说用其他办法,O(7)的复杂度也可以实现,我这个办法不值得提倡。
我晕,我们做算法优化,有时候,O(n)这个值,能减少1都是巨大的成功,因为程序是有循环的,循环次数是被乘数哦,这是乘法关系,这个核心算法复杂度减少1,放出去就是几千万甚至几亿的时钟开销,效率提升就是巨大的。很多时候,我做优化,都在为了减少这个1在努力。
不过,这毕竟是少数人,准确的讲,说这话的人不能算技术人员,因为针对到科学的,算法的,优化的问题上,一是一,二是二,不能带着个人感情讨论。这是技术人员,特别是程序员基本的职业道德。
这里,我希望广大程序员朋友一定要养成一个习惯,“客观”和“严谨”是程序员的基本职业修养,也是我们能在这个行业里面立足的根本,千万不要丢掉了。


因为文中说“O(7)”,我想应该是回应我写的7次迭代,如果不是, 就当本人对号入座吧。对于算法是否属于“科学”,大家可以想想。我自问只谈到技术上的意见,不知道为什么会说到“个人感情”,甚至推理至是否一个技术人员、有没有职业道德问题了。
也许是本人写得不够清楚,也许是读者有误解,无论如何,本人也在此解释一下。

复杂度与性能

首先,肖舸说“O(n)这个值”、“ 核心算法复杂度减少1” 、“O(7)”,都是不正确的说法。

研究算法的运行时间,最常见是采用Big-O表示法,例如O(1)、O(n)等。这表示法是指,算法在n接近无限大的时候,其运行时间的渐近复杂度(asymptotic complexity)的上界(upper bound)。举个例子,例如解决同一问题的两个算法﹐,它们的运行时间(例如单位是秒)为:

\begin{align*} t_A(n) &= 2n + 1 \\ t_B(n) &= 9 \end{align*}

用Big-O表示法,A算法是线性速度O(n),B算法是常数速度O(1)。也就是说n大到某个程度,算法A比算法B慢。那么n要有多大?只要联合两个函数,即可求出,当n > 4时,算法A比B慢。相反,当n<4时,算法A会比算法B快。如果有另一组算法,其方程求出来的n是少于或等于0,说明当中的其中一个算法永远比另一个优。

使用Big-O表示法的目的,是方便比较不同的算法。 Big-O会隐藏的算法实现时的系数(例如上面例子中的2、1、9)。而在实际应用中,设计或实现算法,还要考虑多个因素,例如

  • 算法的应用场合: 算法所需的精确度、时间要求(通常二者须此消彼长)。如算法分为几个步骤(如本问题可分割为初始化和采样),每个步骤的实际使用次数。另外亦要考虑输入数据的分布。
  • 内存的使用: 同时间内所需的内存最大值(可用Big-O表示法)、内存存取模式(access pattern)。
  • 算法实现的难度: 除了开发及维护时间,还影响程序是否健壮(robust),例如复杂的浮点运算会做成误差。
  • 硬件: 例如主频、核心数量、内存速度、特殊指令(如SIMD),或者不同架构的可编程系统(如GPU)、设备等
  • 软件: 会执行于什么操作系统、编译器、库等等,这些软件对算法是否有影响。

很多时候,解决同一问题有许多算法,并没有绝对的好坏之分(当然有时候也有些能完全取胜,有些完全无用),每个算法都其优点缺点。程序员须要分析实际应用,去选择及实现算法。理想地,更可以实验多个算法,用实际数据去作出比较。

评论肖舸的实现

这个是肖舸在《实》公布的实现代码(本人加入了#ifdef MILO_HOTFIX及GetMemory()):

#define CTonyRandomArea_TOKEN_MAX 100           //最大类型数   
#define CTonyRandomArea_TOKEN_AREA_MAX 10000    //类型数组单元数,精确到小数点后两位   
//输入最大100个元素的数组,每个数组表示每类占有的百分比,内部自带百分比调整。   
//即如果外部输入的数字之和不是整数100,内部会根据百分比,自动调整其比例,使总和=100.0   
//然后内部建立10000个单元的类型数组,根据传入的每种类型的比例,在类型数组中批量填充对应的类型值   
//总之,类型数组中每种类型的数量,占据的比例正好是输入的百分比   
//最后,在0~10000中取随机,然后在对应的类型数组单元中取类型值,即为返回的类型   
class CTonyRandomArea   
{   
public:   
    CTonyRandomArea(double* pTokenPercentArray,char cTokenCount)   
    {   
        m_nTokenCount=cTokenCount;   
        if(CTonyRandomArea_TOKEN_MAX<m_nTokenCount)   
            m_nTokenCount=CTonyRandomArea_TOKEN_MAX;   
        int i=0;   
        for(i=0;i<m_nTokenCount;i++)   
        {   
            m_dTokenPercentArray[i]=*(pTokenPercentArray+i);   
        }   
        //动态调整内部的值   
        //有时候试验人员,测得几个状态出现的数字,可能懒得再计算成百分比   
        //程序帮忙自动计算   
        double dNumberCount=0;   
        for(i=0;i<m_nTokenCount;i++)   
        {   
            dNumberCount+=m_dTokenPercentArray[i];   
        }   
        if(100.0!=dNumberCount)   
        {   
            for(i=0;i<m_nTokenCount;i++)   
            {   
                m_dTokenPercentArray[i]/=dNumberCount; 
                m_dTokenPercentArray[i]*=100; 
            }   
        }   
        //以小数点后两位精度,开始计算在10000个总单元中,每种类型对应的数量。   
        for(i=0;i<m_nTokenCount;i++)   
        {   
            m_sTokenPercentArray[i]=(short)(m_dTokenPercentArray[i]*100);   
        }   
  
        //按比例填充类型数组   
        int j=0;   
        int nFillMin=0;   
        int nFillMax=0;   

#ifdef MILO_HOTFIX
		for(i=0;i<CTonyRandomArea_TOKEN_AREA_MAX;i++)
        {   
            m_cTokenPercentArrayAreaUp[i]=0;
        }
#else
		for(i=0;i<m_nTokenCount;i++)
		{
			m_cTokenPercentArrayAreaUp[i]=-1;   
		}
#endif
  
        for(i=0;i<m_nTokenCount;i++)   
        {   
            nFillMax=nFillMin+m_sTokenPercentArray[i];   
            for(j=nFillMin;j<nFillMax;j++)   
            {   
                m_cTokenPercentArrayAreaUp[j]=i;   
            }   
            nFillMin=nFillMax;   
        }   
//      m_cTokenPercentArrayAreaUp[CTonyRandomArea_TOKEN_AREA_MAX-1]=i-1;   
    }   
    ~CTonyRandomArea(){}   
    void PrintfInfo(void)   
    {   
        int i=0;   
        double dNumberCount=0;   
        for(i=0;i<m_nTokenCount;i++)   
        {   
            dNumberCount+=m_dTokenPercentArray[i];   
            printf("%d = %f\n",i,m_dTokenPercentArray[i]);   
        }   
        printf("All = %f\n",dNumberCount);   
  
        //打印10000个单元的类型分布,看看排得对不对   
        //这段打印起来太长,一般隐掉,需要了再打印   
//      int nOutMax=400;   
//      int nInMax=25;      //二者相乘为10000   
//      int j=0;   
//      for(i=0;i<nOutMax;i++)   
//      {   
//          printf("%05d - ",i*nInMax);   
//          for(j=0;j<nInMax;j++)   
//          {   
//              printf("%d ",m_cTokenPercentArrayAreaUp[i*nInMax+j]);   
//          }   
//          printf("\n");   
//      }   
    }   
  
    //取类型数组对应单元的值   
    char GetType(int nTypeIndex)    //输入参数0~10000   
    {   
        if(10000<=nTypeIndex) nTypeIndex=9999;   
        if(0>nTypeIndex) nTypeIndex=0;   
        return m_cTokenPercentArrayAreaUp[nTypeIndex];   
    }   
    //真实的工作函数,利用输入的概率来产生随机type   
    char GetRandomType(void)   
    {   
        return GetType(GetRandomBetween(0,10000));   
    }   

	// MILO ADD {
	size_t GetMemory() const {
		return sizeof(*this);
	}
	// } MILO

private:   
    double m_dTokenPercentArray[CTonyRandomArea_TOKEN_MAX];             //比例数组   
    short m_sTokenPercentArray[CTonyRandomArea_TOKEN_MAX];              //比例数组,短整型,1~10000,相当于精确到小数点后两位   
    char m_nTokenCount;   
    char m_cTokenPercentArrayAreaUp[CTonyRandomArea_TOKEN_AREA_MAX];    //类型数组   
};
////这是测试代码   
//void TestCTonyRandomArea(void)   
//{   
//    double dTokenPercentArray[100];   
//  
//    dTokenPercentArray[0]=12.0;   
//    dTokenPercentArray[1]=40.0;   
//    dTokenPercentArray[2]=40.0;   
//    dTokenPercentArray[3]=7.0;   
//    dTokenPercentArray[4]=1.0;   
//  
//    CTonyRandomArea Area1(dTokenPercentArray,5);   
//    Area1.PrintfInfo();   
//  
//    int i=0;   
//    for(i=0;i<20;i++)   
//    {   
//        printf("RandType = %d\n",Area1.GetRandomType());   
//    }   
//}

崩溃Bug

在测试这个实现时,发现GetRandomType()返回的值会>=N,而N为概率数组的大小。这样有机会做成程序崩溃。例如在代码中,把返回值作为索引,填进数组就会写到不合法的地址。

调试之后,发现程序的bug位于三个地方。第一,21-36行检查概率总和是否等于100个百分点,否则会自动进行正规化(normalization)。正规化后,浮点小数只能表示近似值。第二,在37-41行会把浮点小数乘上10000,再转做整数。因为这里是下取整,计算出来的整数总和有机会少于10000。于是m_cTokenPercentArrayAreaUp数组最后的几个元素未被填入。第三,第54-57行对m_cTokenPercentArrayAreaUp数组进行初始化,但用错循环次数,应为CTonyRandomArea_TOKEN_AREA_MAX而不是m_nTokenCount。

例如使用{ 3, 5, 7, 11, 13 } 作为输入,CTonyRandomArea_TOKEN_AREA_MAX数组就会填入前9998个元素,最后两个未被初始化。

这三个错处要完全更正会改动太多,不如重写。所以本人只是作最少改动,原整地初始化CTonyRandomArea_TOKEN_AREA_MAX数组为0。

用肉眼检验正确性?

测试函数TestCTonyRandomArea()的作用就是列印20个采样。该程序员或测试人员是用肉眼观测那20个采样,看看是否大概和输入的概率分布接近么?或是记下每个返回值的频率,笔算其频率除以20后的百分比?

这个问题本质上就是能编写测试代码去自动测试的。以下的实验就会测试各个实现的特性,包括其误差。

浪费内存

m_dTokenPercentArray和m_sTokenPercentArray都没必要作为成员变量。前者在PrintfInfo()函数还有一点用途,后者只在构造函数使用。

const-ness

构造函数的指针没加const。使用者用的时候要const_cast。 C语言的函数都要适当地使用const指针(可参考strlen原型),何况是C++呢。

乱数产生器不均分布

以下是其的乱数产生器代码:

inline int _GetNot0(void)   
{   
    int nRet=rand();   
    if(!nRet) nRet++;   
    return nRet;   
}   
inline int GetRandomBetween(int nBegin,int nEnd)   
{   
    int n=_GetNot0();   
    int nBetween=0;   
    if(0>nBegin) nBegin=-nBegin;   
    if(0>nEnd) nEnd=-nEnd;   
    if(nBegin>nEnd)   
    {   
        nBetween=nEnd;   
        nEnd=nBegin;   
        nBegin=nBetween;   
    }   
    else if(nBegin==nEnd)   
        nEnd=nBegin+10;   
    nBetween=nEnd-nBegin;   
    n=n%nBetween;   
    n+=nBegin;   
    return n;   
}

这个在书评里已写了很多,原来只需三句代码都可以写成这样,更把乱数的分布变得更不平均,不再多说了。

算法实现

本来以为《用》一文中的说法,已足够解释。以上的情况,加上本人未试过实现别名方法(alias method),便做了二个实验去比较及分析各种算法。实验有其局限性,只是在某个客观环境(计算机、编译器等),对某个算法实现的测量。但本实验的结果可能可以支持我在《用》所提及的观点。

逆变换方法(线性搜寻)

这个算法于《用》有详细说明,以下是其C++代码,复杂度是O(N):

template <typename T, typename RandomGenerator>
class InverseLinear {
public:
	InverseLinear(const T* cdf, size_t N, RandomGenerator& random) : mCDF(cdf), mN(N), mRandom(random) {
		assert(N > 0);
		assert(cdf[N - 1] == (T)1.0);
	}

	size_t operator()() {
		T y = mRandom();
		for (size_t x = 0; x < mN - 1; x++)
			if (y < mCDF[x])
				return x;

		return mN - 1;
	}

	size_t GetMemory() const {
		return sizeof(*this) + mN * sizeof(T);
	}

private:
	const T* mCDF;
	const size_t mN;
	RandomGenerator mRandom;
};

逆变换方法(二元搜寻)

这个算法用二元搜寻代替线性搜寻,复杂度降低为O(lg N):

template <typename T, typename RandomGenerator>
class InverseBinary {
public:
	InverseBinary(const T* cdf, size_t N, RandomGenerator& random) : mCDF(cdf), mN(N), mRandom(random) {
		assert(N > 0);
		assert(cdf[N - 1] == (T)1.0);
	}

	size_t operator()() {
		T y = mRandom();
		size_t left = 0, right = mN;

		while (left < right - 1) {
			size_t mid = (left + right) / 2;
			assert(mid - 1 < mN);
			if (mCDF[mid - 1] < y)
				left = mid;
			else
				right = mid;
		}

		return left;
	}

	size_t GetMemory() const {
		return sizeof(*this) + mN * sizeof(T);
	}

private:
	const T* mCDF;
	const size_t mN;
	RandomGenerator mRandom;
};

必须注意这个算法和一般的二元搜寻条件有些差异。

别名方法

别名方法能提供O(1)的复杂度,但初始化过程比较复杂。参照[1]附录的代码实作:

template <typename T, typename RandomGenerator>
class Alias {
public:
	Alias(const T* pdf, size_t N, RandomGenerator& random) : mN((T)N), mRandom(random), mAliasTable(N) {
		assert(N > 0);
		std::vector<size_t> smallBag, largeBag;
		smallBag.reserve(N);
		largeBag.reserve(N);

		const T stdWeight = (T)1 / N;
		for (size_t i = 0; i < N; i++) {
			mAliasTable[i].prob = pdf[i];

			if (pdf[i] > stdWeight)
				largeBag.push_back(i);
			else
				smallBag.push_back(i);
		}

		while (!largeBag.empty() && !smallBag.empty()) {
			const size_t sm = smallBag.back();
			const size_t lg = largeBag.back();
			smallBag.pop_back();
			largeBag.pop_back();

			mAliasTable[lg].prob -= stdWeight - mAliasTable[sm].prob;
			mAliasTable[sm].prob *= (T)N;
			mAliasTable[sm].index = lg;

			if (mAliasTable[lg].prob > stdWeight)
				largeBag.push_back(lg);
			else
				smallBag.push_back(lg);
		}

		for (std::vector<size_t>::iterator itr = smallBag.begin(); itr != smallBag.end(); ++itr)
			mAliasTable[*itr].prob = (T)1;

		for (std::vector<size_t>::iterator itr = largeBag.begin(); itr != largeBag.end(); ++itr)
			mAliasTable[*itr].prob = (T)1;
	}

	__forceinline size_t operator()() {
		T u = mRandom() * mN;
		size_t i = (size_t)u;
		return u - i > mAliasTable[i].prob ? mAliasTable[i].index : i;
	}

	size_t GetMemory() const {
		return sizeof(*this) + sizeof(AliasItem) * mAliasTable.capacity();
	}

private:
	T mN;
	RandomGenerator mRandom;

	struct AliasItem {
		T prob;
		size_t index;
	};
	std::vector<AliasItem> mAliasTable;
};

伪随机数产生器

Visual C++的C运行时库(C runtime library, CRT)中,提供的rand()是线程安全的,因而有额外的开销(需要存取TLS)。而且,其精确度有限(只有15-bit)。因此,在实验中使用了两个伪乱数产生器,一个采用rand(),一个采用最普通的线性同余产生器(Linear congruential generator, LCG)。两个产生器都是产生半合区间[0,1)的匀均分布。

template <typename T>
class RandomCRT {
public:
	RandomCRT(unsigned seed = 0) {
		srand(seed);
	}

	T operator()() {
		return rand() * ((T)1 / (RAND_MAX + 1));
	}
};
template<typename T>
class RandomLCG {
public:
	RandomLCG(unsigned seed = 0) : mSeed(seed) {}

	T operator()() {
		mSeed = 214013 * mSeed + 2531011;
		return mSeed * ((T)1 / 4294967296);
	}
	
private:
	unsigned mSeed;
};

实验一: 方法

主要的测试方式如下。设N为4, 8, 16, ..., 256,用乱数产生一组N个元素的pdf,之后对每个算法,执行其初始化1,000次,采样1,000,000次,分别量度时间。此外,用采样来计算频表,之后用以下算式比较频表和pdf的总误差:

Error(x, h)=\sum_{0}^{N-1}{ \| h_i - P(X=x_i) \| }
template <typename T>
T ComputeError(const unsigned* histogram, const T* pdf, size_t N, size_t sampleCount) {
	T error = 0.0;
	for (size_t i = 0; i < N; i++) {
		T delta = (T)histogram[i] / sampleCount - pdf[i];
		error += abs(delta);
	}
	return error;
}

当中h_i是采样等于i的出现的次数除以总采样次数。

有些算法分别采用float和double两种类型,亦会使用CRT和LCG两种伪乱数产生器,所有测试组合如下:

  • Tony
  • Linear<CRT, double>
  • Binary<CRT, double>
  • Alias<CRT, double>
  • Linear<LCG, double>
  • Binary<LCG, double>
  • Alias<LCG, double>
  • Linear<LCG, float>
  • Binary<LCG, float>
  • Alias<LCG, float>

实验环境:

  • Intel i7 920 @2.67Hz
  • Microsoft Windows 7 64-bit

编译器及主要编译参数:

  • Microsoft Visual Studio 2008 9.0.30729.1 SP
  • Platform: Win32
  • Optimization: Maximize Speed (/O2)
  • Inline Function Expansion: Any Suitable (/Ob2)
  • Enable Intrinsic Function: Yes (/Oi)
  • Favor Size of Speed: Favor Fast Code (/Ot)
  • Enable C++ Exception (No)
  • Runtime Library: Multi-threaded DLL (/MD)
  • Buffer Security Check: No (/GS-)
  • Enable Enhanced Instruction Set: Streaming SIMD Extensions 2 (/arch:SSE2)
  • Floating Point Model: Fast (/fp:fast)
  • Default Char Unsigned: Yes (/J)
  • Enable Run-Time Type Info: No (/GR-)

宏:

  • #define _SECURE_SCL 0

实验一: 结果及分析

以下结果皆使用Google Visualization API显示,因此在RSS或转载中可能无法显示。我初次试用这API,其HTML和JavaScript代码是由C++测试程序中自动生成的。

原始数据表如下(init和sample的单位为秒,memory的单位为字节):

初始化时间

纵轴是时间,单位为秒。注意横轴为对数比例。

除了Tony和Alias,其余算法的初始化都接近零。 Tony和Alias都是以线性上升。大部份情况下,Tony的初始化性能较低。

采样时间

最好成绩的是Alias的LCG版本,接近常数性能,其float和double版本差异不大。 Tony亦为常数性能,比Alias慢30%左右。

一如所料,Linear和Binary的性能为O(N)和O(lg N)(因横轴为对数比例,对数在图里会接近直线,例如橙线的Binary)。值得注意的是,Linear和Binary的LCG版本在N=4时,超越Alias和Tony。这是由于这两个算法简单,系数低,所以有机会超越O(1)的对手。

需要较高性能的话,别用CRT的rand()。

准确度

如果要说Tony算法的主要问题,就在于其准确度远远差于其他算法的实现(在数百倍以上)。就算使用整数的百份比,用那N=5的例子,Tony的误差约为0.08,其他算法则介乎0.001至0.002。

内存用量

Tony算法所使用的内存是常量11008个字节。而其他的算法为O(N),Alias因为要多存一个别名索引,其内存是其他除Tony以外的两倍(事实上,Alias的索引可以因N的大小而采用unsigned short或unsigned char,以减低开销)。

这里亦要留意一点,Linear和Binary其实只储存cdf的常数指针、N和Random,实际上只有12字节。 cdf数组可供多个采样器共用,所以Linear和Binary的内存开销是非常少的。

N的范围

Tony算法的另一缺点是只支持至N=100。相反,其他的算法可支持较大的N,例如用Tony所需的內存,Alias的float版本可支持N=1375。

实验二: 方法

上一个实验,量度各个算法的采样器在不同的N时的性能。但是在应用上,很多时候需要同时使用多个采样器。肖舸在《实》中提及:

实话跟你讲吧,这段代码是有前提的,我们要做5000万条记录,中间有20万个设备的记录,每个设备的采样频率不一样,我要并发模拟,你再想想我写这么复杂有道理没?

那么就实际测试一下。建立M个采样器(M设为1、10、...、100000),对每个采样器轮流采样,M个采样器合共采样1百万次。而N则取Tony算法可接受的最大值100。测试的循环代码是这样的:

const size_t interleavedSample = sampleCount / M;

size_t unused = 0; // prevent optimized out
LARGE_INTEGER start, end;
QueryPerformanceCounter(&start);
for (size_t i = 0; i < interleavedSample; i++)
	for (size_t j = 0; j < M; j++)
		unused += (*samplers[j])();
QueryPerformanceCounter(&end);
cout << unused << endl;

从理论上来说,用一个采样器去做一百万次采样,和用M个采样器合共采样一百万次,在效能上应该没有分别。但在实际的测试中又是否这样呢?

因为上一个实验已测量准确度,以及那些算法组合在这个情况下比较好,就只测试以下的组合:

  • Tony
  • Binary<LCG, float>
  • Alias<LCG, float>

实验二:结果及分析

原始数据如下:

用图表显示采样性能(横轴为M個采样器,纵轴是采样1百万次时间,单位为秒):

可以发现,Alias在整个测试范围里,表现也是最佳的。若按理论内说,三条线应该都是平的,或接近平的。为什么O(1)的Tony算法却在M=1000以后,输给O(lg N)的Binary?而且,如果用多个采样器是有额外开销的话,为何Alias和Binary不是一直随M上升而线性递增,反而会是个U形,有个甜蜜点(sweet spot)?

Tony算法除了准确度问题,最大问题是其内存开销。当同时使用多个采样器,其存取的内存超过CPU缓存空间,就会严重影响性能。这可以作为例子证实在某些情况下,我在《用》里说的假设:

但对于该文作者说N最大为100,二分搜寻只需最多7次迭代,因缓存问题可能二分搜寻更快。

至于甜蜜点,在实验之前也没预计过。我认为这也是由缓存做成的。因为缓存本身也有平行运作的特性,分散存取整个缓存空间,会比集中存取一小块缓存更有效率。观看上表,可发现甜蜜点的内存使用大约在40K-80K,这似乎刚能对上i7的64K L1缓存。但本人对硬件了解不足,还望各位网友指导。

那么在M=100,000之后,最终Tony是否又会再超越Binary呢?从理论上来说,只要两个算法都使用超过缓存大小的内存(以i7 920来说是8MB的L3),缓存就会慢慢失去效果。可是我没有再尝试,因为在M=100,000,Tony算法已使用了超过1.1GB内存,而Binary才约41MB。回想起来,肖舸说要做20万个设备,每个有不同的采样率。如果要同时使用,Tony算法要2.2GB内存啊,而Binary才82MB,后者比前者简单、准确、快。

总结

本文使用实验比较了4个离散随机变量的产生算法。实验量度的主要四种数据,可以帮助读者选择合适的算法。在于N比较大的情况下,建议使用别名方法,因为其采样性能最高、内存比Tony少、准确度比Tony高,只是初始化性能稍慢。而本人参照[1]的别名方法实现,也可以尝试进一步优化,例如用一个alloca()代替两个vector,减少配置内存的开销。

本人认为,使用类似Tony的算法,在某些情况是有用的。例如,可以设计一个采样器,其输入为整数,代表每个类别的相对比重。例如输入整数数组{ 1, 2, 3, 1, 2 },表示概率{ 1/9, 2/9, 1/3, 1/9, 2/9 },那么只需要9个元素的表就能实现快速的采样。这样的输入,大概十多行代码就能实现,而且准确度非常高,也不需浮点运算。缺点(或特征)是其复杂度与数组的元素总和成线性关系。这个实现留给读者尝试吧。

希望本文也能带出几个要点。在实际应用中,算法不能单靠Big-O时间复杂度来选择。以采样器为例,需同时考虑初始化和采样的时间。例如每次建立一个采样器之后,会执行50次采样,那么就可以比较各种算法的实际执行时间,也许Binary会比Alias好。另外,程序的准确性通常比性能更重要(当然也有例外,例如游戏中的粒子系统效果就不需要很准确),但无论要求如何,也应谨慎对待准确性。除了嵌入式设备有内存限制,因现代计算机的内存速度和缓存缺失(cache miss),编程时应尽量减少内存的开销。最后,也请注意不同算法,代码的复杂程度不尽相同。以上所说,皆是选择算化的重要考虑因素。

本人是否为一名技术人员、是否有职业道德,就不在此回应了。

有与趣的读者可下载本实验的源代码包(不保証0bug)。此外,C++ TR1有部份关于随机数的实现。

参考

更新

  • 2010-05-28: 肖舸在其CSDN博客上发表《请博客园立即停止侵权行为的公开信》。本人认为本文只是适当引用肖舸的文章,为介绍、评论其文章及说明问题,并在发表时已加上原文名称、链结及原作者名字。因此,按法律并不需得到原作者同意。
  • 2010-05-29: 推友@DOTIAN做了一个关于别名方法(alias method)的详细分析。本文主要在性能比较,没有介绍别名方法原理,读者可参考这文章。
  • 2010-06-02: 學生大本營中"編程之美"老師轉載本文。
  • 2010-06-06: 肖舸刪去其CSDN博客所有文章。此非吾所願。
  • 2010-06-10: 肖舸的CSDN博客和學生大本營帪戶(被?)關閉。CSDN要求其他幾位學生大本營老師刪除有關
  • 2010-06-11: 學生大本營中"編程之美"老師轉載的本文被刪。