为啥在梯度噪声发生器中从 Mersenne twister 切换到其他 PRNG 会产生不好的结果？

Posted 2023-02-14

【中文标题】为啥在梯度噪声发生器中从 Mersenne twister 切换到其他 PRNG 会产生不好的结果？

【英文标题】：Why does switching from Mersenne twister to other PRNGs in Gradient Noise Generator give bad results?为什么在梯度噪声发生器中从 Mersenne twister 切换到其他 PRNG 会产生不好的结果？

【发布时间】：2017-12-20 14:00:27

【问题描述】：

我一直在尝试创建一个通用的梯度噪声生成器（它不使用哈希方法来获取梯度）。代码如下：

class GradientNoise 
    std::uint64_t m_seed;
    std::uniform_int_distribution<std::uint8_t> distribution;
    const std::array<glm::vec2, 4> vector_choice = glm::vec2(1.0, 1.0), glm::vec2(-1.0, 1.0), glm::vec2(1.0, -1.0),
                                                    glm::vec2(-1.0, -1.0);

public:
    GradientNoise(uint64_t seed) 
        m_seed = seed;
        distribution = std::uniform_int_distribution<std::uint8_t>(0, 3);
    

    // 0 -> 1
    // just passes the value through, origionally was perlin noise activation
    double nonLinearActivationFunction(double value) 
        //return value * value * value * (value * (value * 6.0 - 15.0) + 10.0);
        return value;
    

    // 0 -> 1
    //cosine interpolation
    double interpolate(double a, double b, double t) 
        double mu2 = (1 - cos(t * M_PI)) / 2;
        return (a * (1 - mu2) + b * mu2);
    

    double noise(double x, double y) 
        std::mt19937_64 rng;
        //first get the bottom left corner associated
        // with these coordinates
        int corner_x = std::floor(x);
        int corner_y = std::floor(y);

        // then get the respective distance from that corner
        double dist_x = x - corner_x;
        double dist_y = y - corner_y;

        double corner_0_contrib; // bottom left
        double corner_1_contrib; // top left
        double corner_2_contrib; // top right
        double corner_3_contrib; // bottom right

        std::uint64_t s1 = ((std::uint64_t(corner_x) << 32) + std::uint64_t(corner_y) + m_seed);
        std::uint64_t s2 = ((std::uint64_t(corner_x) << 32) + std::uint64_t(corner_y + 1) + m_seed);
        std::uint64_t s3 = ((std::uint64_t(corner_x + 1) << 32) + std::uint64_t(corner_y + 1) + m_seed);
        std::uint64_t s4 = ((std::uint64_t(corner_x + 1) << 32) + std::uint64_t(corner_y) + m_seed);


        // each xy pair turns into distance vector from respective corner, corner zero is our starting corner (bottom
        // left)
        rng.seed(s1);
        corner_0_contrib = glm::dot(vector_choice[distribution(rng)], dist_x, dist_y);

        rng.seed(s2);
        corner_1_contrib = glm::dot(vector_choice[distribution(rng)], dist_x, dist_y - 1);


        rng.seed(s3);
        corner_2_contrib = glm::dot(vector_choice[distribution(rng)], dist_x - 1, dist_y - 1);


        rng.seed(s4);
        corner_3_contrib = glm::dot(vector_choice[distribution(rng)], dist_x - 1, dist_y);


        double u = nonLinearActivationFunction(dist_x);
        double v = nonLinearActivationFunction(dist_y);


        double x_bottom = interpolate(corner_0_contrib, corner_3_contrib, u);
        double x_top = interpolate(corner_1_contrib, corner_2_contrib, u);
        double total_xy = interpolate(x_bottom, x_top, v);
        return total_xy;
    
;

然后我生成一个 OpenGL 纹理以显示如下：

int width = 1024;
int height = 1024;
unsigned char *temp_texture = new unsigned char[width*height * 4];
double octaves[5] = 2,4,8,16,32;

for( int i = 0; i < height; i++)
    for(int j = 0; j < width; j++)
        double d_noise = 0;
        d_noise += temp_1.noise(j/octaves[0], i/octaves[0]);
        d_noise += temp_1.noise(j/octaves[1], i/octaves[1]);
        d_noise += temp_1.noise(j/octaves[2], i/octaves[2]);
        d_noise += temp_1.noise(j/octaves[3], i/octaves[3]);
        d_noise += temp_1.noise(j/octaves[4], i/octaves[4]);
        d_noise/=5;
        uint8_t noise = static_cast<uint8_t>(((d_noise * 128.0) + 128.0));
        temp_texture[j*4 + (i * width * 4) + 0] = (noise);
        temp_texture[j*4 + (i * width * 4) + 1] = (noise);
        temp_texture[j*4 + (i * width * 4) + 2] = (noise);
        temp_texture[j*4 + (i * width * 4) + 3] = (255);
    

效果不错：

但是 gprof 告诉我 Mersenne twister 占用了我 62.4% 的时间，并且随着纹理的增大而增长。没有其他人会花费如此多的时间。虽然 Mersenne twister 在初始化后很快，但我每次使用它时都对其进行初始化这一事实似乎使它变得相当慢。

这个初始化是 100% 需要的，以确保相同的 x 和 y 在每个整数点生成相同的梯度（因此您需要一个哈希函数或每次都为 RNG 播种）。

我尝试将 PRNG 更改为线性同余生成器和 Xorshiftplus，虽然两者的运行速度都快了几个数量级，但结果却很奇怪：

LCG（1次，使用前运行5次）

Xorshiftplus

经过一次迭代

经过 10,000 次迭代。

我试过了：

在利用输出之前运行生成器数次，这会导致执行缓慢或只是不同的工件。

使用初始种子后两次连续运行的输出再次播种 PRNG 并使用病房后的值。结果没有区别。

发生了什么？我该怎么做才能获得与 mersenne twister 相同质量的更快结果？

好的大更新：

我不知道为什么会这样，我知道它与使用的素数有关，但是在搞砸了一点之后，似乎以下工作：

第 1 步，将 x 和 y 值分别合并为种子（并与它们合并一些其他偏移值或附加种子值，该数字应该是素数/非平凡因子）

第 2 步，使用这两个种子结果将生成器再次播种回函数中（就像 geza 所说，制作的种子不好）

第 3 步，当得到结果时，而不是使用模数项 (4) 试图得到，或者 & 3，将结果模数以素数 first 然后应用 & 3。我不确定素数是梅森素数是否重要。

这是 prime = 257 和使用 xorshiftplus 的结果！ （注意我用的是 2048 乘 2048，其他的是 256 乘 256）

【问题讨论】：

顺便说一句，为什么 rng 是类成员，而不是自动变量？

您正在使用 PRNG 作为一个非常昂贵的哈希函数。尝试使用实际（加密？）散列函数。

@yurikilochek 我该怎么做？

@snb 什么不清楚？只需通过散列函数传递你的种子（或直接坐标），从结果中选择两个位来选择你的向量

@snb：你的最后一张（10000 次迭代）图片显示了什么？你的意思是，你播种了 xorshift，然后你生成并忽略了 10000 个数字，然后你使用了第 10001 个？然后，你甚至得到了这张带有这些图案的照片？

【参考方案1】：

众所周知，LCG 不适合您的目的。

Xorshift128+ 的结果很糟糕，因为它需要良好的播种。提供良好的播种会破坏使用它的全部目的。我不推荐这个。

但是，我建议使用整数哈希。例如，来自Bob's page 的一个。

这是该页面第一个哈希的结果，在我看来还不错，而且速度很快（我认为它比 Mersenne Twister 快得多）：

这是我为生成此代码而编写的代码：

#include <cmath>
#include <stdio.h>

unsigned int hash(unsigned int a) 
    a = (a ^ 61) ^ (a >> 16);
    a = a + (a << 3);
    a = a ^ (a >> 4);
    a = a * 0x27d4eb2d;
    a = a ^ (a >> 15);
    return a;


unsigned int ivalue(int x, int y) 
    return hash(y<<16|x)&0xff;


float smooth(float x) 
    return 6*x*x*x*x*x - 15*x*x*x*x + 10*x*x*x;


float value(float x, float y) 
    int ix = floor(x);
    int iy = floor(y);
    float fx = smooth(x-ix);
    float fy = smooth(y-iy);

    int v00 = ivalue(iy+0, ix+0);
    int v01 = ivalue(iy+0, ix+1);
    int v10 = ivalue(iy+1, ix+0);
    int v11 = ivalue(iy+1, ix+1);
    float v0 = v00*(1-fx) + v01*fx;
    float v1 = v10*(1-fx) + v11*fx;
    return v0*(1-fy) + v1*fy;


unsigned char pic[1024*1024];

int main() 
    for (int y=0; y<1024; y++) 
        for (int x=0; x<1024; x++) 
            float v = 0;

            for (int o=0; o<=9; o++) 
                v += value(x/64.0f*(1<<o), y/64.0f*(1<<o))/(1<<o);
            

            int r = rint(v*0.5f);

            pic[y*1024+x] = r;
        
    

    FILE *f = fopen("x.pnm", "wb");
    fprintf(f, "P5\n1024 1024\n255\n");
    fwrite(pic, 1, 1024*1024, f);
    fclose(f);

如果您想了解散列函数的工作原理（或者更好的是，了解一个好的散列具有哪些属性），请查看 Bob 的页面，例如 this。

【讨论】：

你的渐变在哪里？

我不想让这段代码更复杂。它只有八度和。检查哈希质量就足够了。

我很困惑，这是您用于创建上述图像的确切代码吗？我很困惑这是如何在没有渐变的情况下工作的：/

@snb：是的，这是（你可以试试，它是可编译的。例如，可以用 gimp 查看生成的 x.pnm）。你不需要渐变来得到这样的图片，只需要八度相加。

我以为你需要渐变；那么渐变实际上有什么帮助呢？看起来这个结果很好。该代码使其看起来像随机值噪声。

【参考方案2】：

您（在不知不觉中？）实现了 PRNG 非随机模式的可视化。看起来很酷！

除了 Mersenne Twister，您测试的所有 PRNG 似乎都不适合您的目的。由于我自己没有做进一步的测试，我只能建议尝试和测量进一步的 PRNG。

【讨论】：

查看我更新的帖子，显然他们确实满足，但对素数有些小题大做，我就是不明白为什么！

【参考方案3】：

众所周知，LCG 的随机性对其参数的选择很敏感。特别是，period of a LCG 是相对于 m 参数的 - 最多它将是 m（您的主要因素）并且对于许多值它可以更小。

同样，需要仔细选择参数才能获得long period from Xorshift PRNG。

您已经注意到，一些 PRNG 提供了良好的程序生成结果，而另一些则没有。为了查明原因，我将排除 proc gen 的内容并直接检查 PRNG 输出。可视化数据的一种简单方法是构建灰度图像，其中每个像素值都是（可能缩放的）随机值。对于基于图像的东西，我发现这是一种查找可能导致视觉伪影的东西的简单方法。您看到的任何工件都可能导致您的 proc gen 输出出现问题。

另一种选择是尝试Diehard tests 之类的方法。如果上述图像测试未能发现任何问题，我可能会使用它来确保我的 PRNG 技术值得信赖。

【参考方案4】：

请注意，您的代码播种 PRNG，然后从 PRNG 生成一个伪随机数。您发现 xorshift128+ 中的非随机性的原因是 xorshift128+ 在更改其状态之前简单地添加了种子的两半（并使用结果 mod 2⁶⁴ 作为生成的数字）（查看其source code)。这使得 PRNG 与散列函数有很大不同。

【参考方案5】：

您所看到的是 PRNG 质量的实际演示。 Mersenne Twister 是性能最好的 PRNG 之一，它通过了 DIEHARD 测试。要知道，生成一个随机数并不是一项简单的计算任务，因此寻求更好的性能必然会导致质量不佳。 LCG 被认为是有史以来设计的最简单和最差的 PRNG，它清楚地显示了您图片中的二维相关性。 Xorshift 生成器的质量很大程度上取决于位数和参数。它们肯定比 Mersenne Twister 差，但有些 (xorshift128+) 可能工作得很好，可以通过 BigCrush 电池组 TestU01 测试。

换句话说，如果您要进行重要的物理建模数值实验，您最好继续使用 Mersenne Twister，众所周知，它可以在速度和质量之间取得良好的平衡，并且它包含在许多标准库中。在不太重要的情况下，您可以尝试使用 xorshift128+ 生成器。要获得最终结果，您需要使用加密质量的 PRNG（此处未提及的任何一个都可能用于加密目的）。

【讨论】：

不要夸大你的情况，LCG 不是有史以来设计的最简单和最差的 PRNG。冯·诺依曼真正介入其中的是middle-square method，它具有定点行为和一些种子值的非常短的周期。