如何高效地对C++向量进行归一化

Question

我想知道如何在C++中高效地对向量进行归一化。 迄今为止，这就是我所拥有的内容。是否有一种方法可以使其更加高效，或者一次性完成？

std::array<float, MyClass::FEATURE_LENGTH> MyClass::normalize(const std::array<float, FEATURE_LENGTH>& arr) {
    std::array<float, MyClass::FEATURE_LENGTH> output{};
    double mod = 0.0;

    for (size_t i = 0; i < arr.size(); ++i) {
        mod += arr[i] * arr[i];
    }

    double mag = std::sqrt(mod);

    if (mag == 0) {
        throw std::logic_error("The input vector is a zero vector");
    }

    for (size_t i = 0; i < arr.size(); ++i) {
        output[i] = arr[i] / mag;
    }

    return output;
}

- cyrusbehr

---

你考虑过原地归一化向量吗？创建一个新的数组会导致内存分配，会比较耗费资源。 - OutOfBound

如果您的目标架构支持SIMD指令，这也可以提高性能。尝试使用“-march=native”编译和/或使用内部函数或矢量化库。 - OutOfBound

为什么这个问题被踩了？它很明确，不平凡，并且提供了一个自包含的示例。 - OutOfBound

@OutOfBound我唯一能想到的原因是这个问题更适合于Code Review。Stack Overflow非常适合处理不起作用的代码问题。Code Review旨在解决“我的代码可以运行，但如何让它变得更好？”的问题。请注意我链接了帮助页面。这是为了万一提问者想要删除这个问题并在Code Review上重新提问，他们应该先仔细阅读提问部分，以获得最好的回应。 - user4581301

@OutOfBound std::array 永远不会分配动态内存。 - smiling_nameless

Answer 1

根据问题的具体情况，有许多优化此算法实现的方法。

1. 对于所有循环，您可以使用SIMD向量化来增加吞吐量。
2. 如果您的向量非常宽，则可以使用多个线程计算幅度。每个线程将计算部分总和，然后一些串行代码将收集结果。
3. 如果您的值在范围内，则可以完全使用浮点数而不是双倍精度。
4. 您可以使用内置函数（例如x86上的RSQRTSS）或使用Quake's method计算幅度的倒数平方根。然后您将按比例缩放。

此外，通过与归一化融合操作，您可以获得更快的代码。假设您想要添加两个向量并规范化结果。您可以在单个传递中计算它们的总和和幅度，然后在第二个传递中进行缩放。

- Alex Reinking

Answer 2

在单次遍历中如何实现呢？显然，您需要使用所有项目来计算 mag ，并且在更新项目之前必须已经进行了计算。

由于执行除法可能比执行乘法要花费更多的时间，因此可能有一种优化方法是添加：

double mag_inv = 1.0 / mag;

然后你可以这样乘以项目：

output[i] = arr[i] * mag_inv;

如果向量已经被归一化的概率相对较高，您可能希望检查 mag 是否等于 1.0。

- Phil1970

---

（假设除法比乘法慢。对于我所知道的所有硬件来说都是如此。） - Rick James

因为 mag 是循环不变量，所以在使用 -Ofast 编译时，GCC会自动进行该更改。 - Alex Reinking

@AlexReinking 优化是编译器和选项相关的。一些编译器具有“快速数学”选项，这些选项不符合IEEE标准，但可以生成更快的代码。 - Phil1970

@Phil1970 - 是的，我知道。这就是-Ofast所启用的功能。我只是观察到你在答案中给出的转换（将除法转换为乘法）既不符合IEEE标准，也可以在设置了这些标志时自动执行。这很重要，因为如果OP正在使用快速数学编译，那么这将毫无作用。 - Alex Reinking

Answer 3

如果有人需要，这里是SIMD向量化代码的示例：

#include <immintrin.h> //header for SIMD functions

void Normalize(const float lpInput[4], float lpOutput[4]) {
    __m128 vInput = _mm_load_ps(lpInput); // load input vector (x, y, z, a)
    __m128 vSquared = _mm_mul_ps(vInput, vInput); // square the input values
    __m128 vHalfSum = _mm_hadd_ps(vSquared, vSquared); 
    __m128 vSum = _mm_hadd_ps(vHalfSum, vHalfSum); // compute the sum of values
    float fInvSqrt; _mm_store_ss(&fInvSqrt, _mm_rsqrt_ss(vSum)); // compute the inverse sqrt
    __m128 vNormalized = _mm_mul_ps(vInput, _mm_set1_ps(fInvSqrt)); // normalize the input vector
    _mm_store_ps(lpOutput, vNormalized); // store normalized vector (x, y, z, a)
}

为了正确编译它，您需要在编译器选项中启用SSE和AVX指令（对于gcc或clang，请使用-msse -mavx || 对于msvc，请使用/arch:sse /arch:avx）。

- Шериф Грей

---

这不是一个特别好的归一化例程，适用于任意长度的向量 - 这只处理4D情况。您假设输入和输出都是16字节对齐的（可能不是）。它使用rsqrt，所以会失去很多精度（对于表面法线可能还可以，但由于缺乏精度，在许多光线跟踪案例中将失败）。那段代码中没有任何需要AVX的内容。为什么在这里使用_mm_store_ss？_mm_cvtss_f32才是正确的方法。vSum中的每个元素都是相同的，因此_mm_rqsrt_ps(vSum)将避免存储和广播。_mm_dp_ps在这里是更好的选择。 - robthebloke

这是一个关于规范化固定长度向量的入门级问题，请放心。而且这仍然被认为是较高级别的代码，而不是直接的汇编代码。 O3优化将使用某种形式的广播单个值到整个向量而不是存储\设置。 - Шериф Грей