现在的位置: 首页 > 综合 > 正文

实验比较各离散采样算法

2017年11月05日 ⁄ 综合 ⁄ 共 21179字 ⁄ 字号 小 中 大 ⁄ 评论关闭

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

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

自从肖舸在其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)。举个例子，例如解决同一问题的两个算法﹐，它们的运行时间(例如单位是秒)为:

用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()):

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#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++呢。

乱数产生器不均分布

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

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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):

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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):

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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]附录的代码实作:

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
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)的匀均分布。

?
1
2
3
4
5
6
7
8
9
10
11
template <typename T>
class RandomCRT
{
public:
    RandomCRT(unsigned
seed = 0) {
        srand(seed);
    }
 
    T
operator()() {
        return rand()
* ((T)1 / (RAND_MAX + 1));
    }
};
?
1
2
3
4
5
6
7
8
9
10
11
12
13
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的总误差:

?
1
2
3
4
5
6
7
8
9
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;
}

当中是采样等于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的单位为字节):

sampler  N  init  sample  error  memory 
Tony 4 0.000664 0.027172 0.081564 11008
Linear<CRT, double> 4 0 0.036899 0.001122 44
Binary<CRT, double> 4 0 0.037687 0.000393 44
Alias<CRT, double> 4 0.000901 0.039043 0.001285 96
Linear<LCG, double> 4 0 0.012864 0.001358 44
Binary<LCG, double> 4 0 0.015508 0.001358 44
Alias<LCG, double> 4 0.000901 0.020245 0.001269 96
Linear<LCG, float> 4 0 0.014556 0.001358 28
Binary<LCG, float> 4 0 0.013205 0.001358 28
Alias<LCG, float> 4 0.000388 0.019169 0.001269 56
Tony 8 0.000866 0.02719 0.114413 11008
Linear<CRT, double> 8 0 0.041352 0.001867 76
Binary<CRT, double> 8 0 0.044668 0.001192 76
Alias<CRT, double> 8 0.000613 0.040483 0.002028 160
Linear<LCG, double> 8 0 0.021303 0.001793 76
Binary<LCG, double> 8 0 0.024576 0.001793 76
Alias<LCG, double> 8 0.000605 0.01771 0.001696 160
Linear<LCG, float> 8 0 0.02296 0.001793 44
Binary<LCG, float> 8 0 0.020713 0.001793 44
Alias<LCG, float> 8 0.000626 0.016415 0.001696 88
Tony 16 0.001118 0.02716 0.100027 11008
Linear<CRT, double> 16 0 0.046342 0.002474 140
Binary<CRT, double> 16 0 0.049653 0.001773 140
Alias<CRT, double> 16 0.000965 0.040003 0.003115 288
Linear<LCG, double> 16 0 0.026674 0.002771 140
Binary<LCG, double> 16 0 0.033027 0.002771 140
Alias<LCG, double> 16 0.000772 0.019638 0.002637 288
Linear<LCG, float> 16 0 0.029094 0.002771 76
Binary<LCG, float> 16 0 0.02711 0.002771 76
Alias<LCG, float> 16 0.000757 0.018408 0.002637 152
Tony 32 0.00168 0.027125 0.117478 11008
Linear<CRT, double> 32 0 0.054437 0.003796 268
Binary<CRT, double> 32 0 0.060345 0.003682 268
Alias<CRT, double> 32 0.001328 0.04101 0.00407 544
Linear<LCG, double> 32 0 0.0378 0.004307 268
Binary<LCG, double> 32 0 0.038143 0.004307 268
Alias<LCG, double> 32 0.001088 0.021245 0.004277 544
Linear<LCG, float> 32 0 0.03846 0.004307 140
Binary<LCG, float> 32 0 0.039474 0.004307 140
Alias<LCG, float> 32 0.001063 0.020093 0.004277 280
Tony 64 0.002716 0.027327 0.125566 11008
Linear<CRT, double> 64 0 0.068234 0.006185 524
Binary<CRT, double> 64 0 0.071522 0.005066 524
Alias<CRT, double> 64 0.002016 0.041144 0.005638 1056
Linear<LCG, double> 64 0 0.051763 0.005387 524
Binary<LCG, double> 64 0 0.048 0.005387 524
Alias<LCG, double> 64 0.001813 0.021425 0.004949 1056
Linear<LCG, float> 64 0 0.052001 0.005387 268
Binary<LCG, float> 64 0 0.049192 0.005387 268
Alias<LCG, float> 64 0.001906 0.020285 0.004949 536
Linear<CRT, double> 128 0 0.103351 0.008263 1036
Binary<CRT, double> 128 0 0.079406 0.007868 1036
Alias<CRT, double> 128 0.003739 0.04089 0.008274 2080
Linear<LCG, double> 128 0 0.087546 0.006949 1036
Binary<LCG, double> 128 0 0.05506 0.006949 1036
Alias<LCG, double> 128 0.003442 0.021405 0.006866 2080
Linear<LCG, float> 128 0 0.087284 0.006943 524
Binary<LCG, float> 128 0 0.056546 0.006943 524
Alias<LCG, float> 128 0.003467 0.020354 0.006866 1048
Linear<CRT, double> 256 0 0.172883 0.011585 2060
Binary<CRT, double> 256 0 0.087077 0.011734 2060
Alias<CRT, double> 256 0.006625 0.041404 0.011714 4128
Linear<LCG, double> 256 0 0.156308 0.010972 2060
Binary<LCG, double> 256 0 0.063352 0.010972 2060
Alias<LCG, double> 256 0.006188 0.021585 0.011057 4128
Linear<LCG, float> 256 0 0.157396 0.010976 1036
Binary<LCG, float> 256 0 0.064066 0.010968 1036
Alias<LCG, float> 256 0.006339 0.020347 0.011057 2072

初始化时间

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

Initialization TimeTonyLinear<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>481632641282560.0000.0020.0040.0060.008

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

采样时间

Sampling TimeTonyLinear<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>481632641282560.000.050.100.150.20

最好成绩的是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()。

准确度

ErrorTonyLinear<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>481632641282560.000.040.080.120.16

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

内存用量

MemoryTonyLinear<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>4816326412825603,0006,0009,00012,000

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。测试的循环代码是这样的:

?
1
2
3
4
5
6
7
8
9
10
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>

实验二：结果及分析

原始数据如下:

sampler  M  init  sample  -  memory 
Tony 1 0.000016 0.027413 0 11008
Binary<LCG, float> 1 0 0.057389 0 412
Alias<LCG, float> 1 0.000005 0.022055 0 824
Tony 10 0.000085 0.02733 0 110080
Binary<LCG, float> 10 0.000003 0.046024 0 4120
Alias<LCG, float> 10 0.000041 0.017833 0 8240
Tony 100 0.000654 0.029593 0 1100800
Binary<LCG, float> 100 0.000025 0.042469 0 41200
Alias<LCG, float> 100 0.000283 0.015489 0 82400
Tony 1000 0.00923 0.055332 0 11008000
Binary<LCG, float> 1000 0.000251 0.049636 0 412000
Alias<LCG, float> 1000 0.002857 0.018702 0 824000
Tony 10000 0.101176 0.093712 0 110080000
Binary<LCG, float> 10000 0.002474 0.060276 0 4120000
Alias<LCG, float> 10000 0.02918 0.030357 0 8240000
Tony 100000 0.850493 0.118865 0 1100800000
Binary<LCG, float> 100000 0.025231 0.105215 0 41200000
Alias<LCG, float> 100000 0.292381 0.047954 0 82400000

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

Sampling time with N=100, M samplersTonyBinary<LCG,float>Alias<LCG,float>1101001000100001000000.000.030.060.090.12

可以发现，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: 學生大本營中"編程之美"老師轉載的本文被刪。
返回
【上篇】
【下篇】

作者:

抱歉!评论已关闭.

返回首页

Copyright © 2013-2018 学步园  保留所有权利.
软文销售 QQ客服：2265327166 点击这里给我发消息（其他合作也可洽谈）
必威体育
必威电竞