Python与C++矩阵乘法的浮点精度差异及符号异常问题

问题背景 #

我执行以下三步操作：

1. 定义一个4×4二维方阵
2. 计算该矩阵的转置
3. 原矩阵与转置矩阵相乘

分别用C++和Python实现后，输入矩阵完全相同，但乘法结果存在明显精度差异：非对角线元素不仅精度差距大，部分元素符号甚至反转。需要解决：

• 如何让两种语言得到相同精度的结果
• 符号反转的原因

原始代码 #

C++ 实现 #

#include <iostream>
#include <iomanip>  // For setw

int main() {
    float R[16] = {
        0.5, 0.63245553, -0.5, 0.31622777,
        0.5, 0.31622777, 0.5, -0.63245553,
        0.5, -0.31622777, 0.5, 0.63245553,
        0.5, -0.63245553, -0.5, -0.31622777
    };

    const int nRows = 4;
    const int nCols = 4;
    float result[nRows][nRows] = { 0 };

    // Perform matrix multiplication
    for (int i = 0; i < nRows; i++) {
        for (int j = 0; j < nRows; j++) {
            for (int k = 0; k < nCols; k++) {
                result[i][j] += R[i * nCols + k] * R[j * nCols + k];
            }
        }
    }

    // Print the result with left-aligned columns and a padding of 15 characters
    for (int i = 0; i < nRows; i++) {
        for (int j = 0; j < nRows; j++) {
            std::cout << std::left << std::setw(15) << result[i][j] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

Python 实现 #

import numpy as np


R = np.array([
    0.5, 0.63245553, -0.5, 0.31622777,
    0.5, 0.31622777, 0.5, -0.63245553,
    0.5, -0.31622777, 0.5, 0.63245553,
    0.5, -0.63245553, -0.5, -0.31622777
]).reshape(4, 4)

result = np.dot(R, R.T)

# Print the result
print(result)

运行结果 #

C++ 输出 #

1               -1.49012e-08    0               -7.45058e-09
-1.49012e-08    1               0               0
0               0               1               -1.49012e-08
-7.45058e-09    0               -1.49012e-08    1   

Python 输出 #

[[1.00000000e+00 2.77555756e-17 0.00000000e+00 5.32461714e-11]
 [2.77555756e-17 1.00000000e+00 5.32461852e-11 0.00000000e+00]
 [0.00000000e+00 5.32461852e-11 1.00000000e+00 2.77555756e-17]
 [5.32461714e-11 0.00000000e+00 2.77555756e-17 1.00000000e+00]]

原因分析 #

1. 数据类型精度不匹配 #

C++代码使用单精度浮点数（float），有效数字约6-7位；Python的NumPy默认使用双精度浮点数（float64），有效数字约15-17位。两者存储输入矩阵时的误差基础不同，累积到乘积结果后，精度差异被放大。

2. 符号反转的本质 #

符号反转是浮点数舍入误差的典型表现：输入矩阵中的非对角线元素理论值应为0，但原始小数（如0.63245553是√(2/5)的近似值）无法用二进制浮点数精确表示。累加多个接近抵消的数时，舍入误差的累积方向不同，导致C++和Python的结果符号相反，但绝对值都属于可接受的浮点误差范围（远小于1e-6）。

3. 计算实现的底层差异 #

NumPy的np.dot调用优化的BLAS/LAPACK库实现矩阵乘法，计算顺序和精度优化逻辑与手动编写的C++三层循环不同，进一步导致误差累积路径的差异。

解决方案 #

1. 统一数据类型 #

方案A：Python改用单精度

修改Python代码的数组定义，指定dtype=np.float32，与C++的float对齐：

R = np.array([
    0.5, 0.63245553, -0.5, 0.31622777,
    0.5, 0.31622777, 0.5, -0.63245553,
    0.5, -0.31622777, 0.5, 0.63245553,
    0.5, -0.63245553, -0.5, -0.31622777
], dtype=np.float32).reshape(4, 4)

方案B：C++改用双精度

将C++中的float替换为double，与NumPy默认的float64对齐：

double R[16] = {
    0.5, 0.63245553, -0.5, 0.31622777,
    0.5, 0.31622777, 0.5, -0.63245553,
    0.5, -0.31622777, 0.5, 0.63245553,
    0.5, -0.63245553, -0.5, -0.31622777
};
double result[nRows][nRows] = { 0 };

2. 统一计算逻辑（可选） #

如果需要完全一致的计算顺序，在Python中手动实现与C++相同的三层循环，避免NumPy的优化逻辑影响：

import numpy as np

R = np.array([
    0.5, 0.63245553, -0.5, 0.31622777,
    0.5, 0.31622777, 0.5, -0.63245553,
    0.5, -0.31622777, 0.5, 0.63245553,
    0.5, -0.63245553, -0.5, -0.31622777
], dtype=np.float32).reshape(4, 4)

nRows, nCols = 4, 4
result = np.zeros((nRows, nRows), dtype=np.float32)

for i in range(nRows):
    for j in range(nRows):
        for k in range(nCols):
            result[i][j] += R[i, k] * R[j, k]

print(result)

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