Durand-Kerner算法多项式求根计算非确定性问题咨询

问题背景 #

使用Durand-Kerner算法计算多项式根时，相同输入下不同线程组的收敛速率存在差异，部分线程组甚至出现发散（推测产生NaN）。已排除同步bug（线程组仅包含一个warp，7<32），需找出代码中可能存在的非确定性操作。

代码实现 #

// Durand-Kerner算法实现
groupshared Complex roots[7];
groupshared float a[8];
groupshared float2 uv_sum[8];
groupshared float2 uv_max[8];
static const float PI = 3.14159265f;

[numthreads(7, 1, 1)]
void main( uint3 tid : SV_DispatchThreadID )
{
    uint index = tid.x;
    // z^7 + 3z^6 - 2z^5 + 10z^4 - 2z^3 + 8z^2 - z - 13
    float a_[8] = { 1.f, 3.f, -2.f, 10.f, -2.f, 8.f, -1.f, -13.f };
    a[index] = a_[index];
    if (index == 0)
    {
        a[7] = a_[7];
    }
    
    GroupMemoryBarrierWithGroupSync();
    
    float ui = 2 * pow(abs(a[index + 1]), 1.f / (index+1));
    float vi = 0.5f * pow(abs(a[7] / a[index]), 1.f / (7-index));
    uv_sum[index] = float2(ui, vi);
    uv_max[index] = float2(ui, vi);
    GroupMemoryBarrierWithGroupSync();
    
    for (unsigned int s = 8 / 2; s > 0; s >>= 1)
    {
        if (index < s)
        {
            uv_sum[index] += uv_sum[index + s];
            uv_max[index] = max(uv_max[index + s], uv_max[index]);
        }
        GroupMemoryBarrierWithGroupSync();
    }

    float2 cs;
    sincos(index * 2.f * PI / 7.f, cs.y, cs.x);
    float2 uv = uv_sum[0] / uv_max[0];
    cs *= (uv.x + uv.y) / 2;
    //Initial[index] = cs;
    roots[index] = Complex_(cs.x, cs.y);
    
    
    for (uint loop = 0; loop < 15; loop++)
    {
        GroupMemoryBarrierWithGroupSync();
        Complex zi = roots[index];
    
        Complex zipow[8];
        zipow[1] = zi;
        for (uint i = 2; i < 8; i++)
        {
            zipow[i] = zipow[i - 1] * zi;
        }
    
        //Complex Pzi = zip7 * a[0] + zip6 * a[1] + zip5 * a[2] + zip4 * a[3] + zip3 * a[4] + zip2 * a[5] + zi * a[6] + Complex_(a[7], 0);
        Complex Pzi = zipow[7] * a[0] + zipow[6] * a[1] + zipow[5] * a[2] + zipow[4] * a[3] + zipow[3] * a[4] + zipow[2] * a[5] + zi * a[6] + Complex_(a[7], 0);
        Complex d1 = (zi - roots[(index + 1) % 7]);
        Complex d2 = (zi - roots[(index + 2) % 7]);
        Complex d3 = (zi - roots[(index + 3) % 7]);
        Complex d4 = (zi - roots[(index + 4) % 7]);
        Complex d5 = (zi - roots[(index + 5) % 7]);
        Complex d6 = (zi - roots[(index + 6) % 7]);
        Complex product = d1 * d2 * d3 * d4 * d5 * d6;
        GroupMemoryBarrierWithGroupSync();
    
        Complex div = Pzi / product;
        Complex rooti = zi - div;
        roots[index] = rooti;
    }
    
    Complex root = roots[index];
    float theta = atan2(abs(root.imag), root.real) / PI;
    Output[index + 7 * tid.y] =
        int(theta * 1024);
}

可能的非确定性操作 #

• 浮点幂运算（pow）：GPU硬件对pow的实现存在细微精度差异，不同线程组的计算单元执行该指令时，结果会有极小偏差。Durand-Kerner是迭代算法，初始值的微小误差会被迭代放大，导致收敛速率不同甚至发散。
• 浮点除法与复数除法：当迭代过程中product（复数乘积）趋近于0时，Pzi / product会产生极大值或NaN。由于浮点精度的不确定性，部分线程组可能恰好触发这种极端情况，导致计算发散。此外，复数除法的硬件实现也可能存在精度波动。
• 归约阶段的浮点加法：初始化时的uv_sum归约计算中，浮点加法不满足结合律，不同线程组的加法执行顺序可能因硬件调度不同而变化，导致最终uv_sum[0]的结果存在细微差异，进而影响初始根的位置。
• 三角函数（sincos、atan2）：GPU对三角函数的硬件实现存在精度误差，不同计算单元的结果可能略有不同，同样会影响初始根的生成和最终结果的计算。

内容的提问来源于stack exchange，提问作者Tom Huntington