R -- 对矩阵列进行排序后，某些列的列平均值现在有所不同......为什么？ [重复]

Question

因此，我想绘制根据大小排序时各列的累积和。这最终产生了一个看起来不像我期望看到的情节。我的期望可能是错误的，但可以肯定的是，在对几列进行排序后，列的平均值已经改变。我一生都无法弄清楚为什么。

事实上，这是一个虚假的数据练习，所以你可以自己看看：

df.exercise = .2
set.seed(78455823)

toad = matrix(rt(1000 * 4000, df.exercise), ncol = 1000)

for(i in 1:1000){
  toad.sort[, i] = toad[order(abs(toad[, i])), i]
}

which(colMeans(toad) != colMeans(toad.sort))
which(colSums(toad) != colSums(toad.sort))

loop.sums = loop.sums.sort = numeric(1000)

for(i in 1:1000){
 loop.sums.sort[i] = sum(toad.sort[, i])
 loop.sums[i] = sum(toad[, i])
}

which(loop.sums.sort != loop.sums)
which(loop.sums.sort != colSums(toad.sort))
which(loop.sums != colSums(toad))

还有一种情况是

colSums(toad) / 4000

不会导致

colMeans(toad)

。

那么，有谁知道为什么会发生分歧？

Answer 1

发现得好！浮点算术的怪癖之一是，与传统数学不同，运算符

+

关心数据呈现的顺序（及其符号和大小）。

有一些简单的方法可以解决这个问题 - 始终将总和转化为更高精度的变量，该变量在尾数中有足够的空间来处理结果总和。将实数或整数放入双精度累加器，将双精度数转换为适当的长双精度数（真正的 80 位必需的 - 有时

long double == double

，例如 MSC 编译器或 MS 兼容模式下的 Intel 是 64 位）。

对数组进行预排序，首先处理最小值将给出最准确的常规结果。相反，对数组进行预排序以便首先对最大的数组求和将得到最不准确的结果。

可以采用各种方案，例如递归地对成对的有序值求和（根据需要用零填充）。这可以很好地降低错误，并且比简单求和更快。

最终的方法将获得更好的结果（假设您的优化编译器不会破坏它）是Kahan Summation，它会向前推进一个补偿值

c

，代表每次加法中产生的错误。它利用给定

a, b

的精确方法来计算一对

c, d

，使得

a + b == c + d

准确无误。这是它的简单 C 代码版本（从伪代码参考转换而来 - 可能存在拼写错误）。

double KahanSum(double *x, int n)
{
  float sum = 0.0;
  float c = 0.0;   // running compensation for lost low-order bits.
  // sum is potentially big, y small, so low-order digits of y are lost.

  for (int i = 0; i < n; i++)
  {
     float t, y;       
     var y = x[i] - c;    // c is zero the first time around.
     t = sum + y;         // (t - sum) cancels the high-order part of y;
    
     c = (t - sum) - y;   // subtracting y recovers negative (low part of y)
     sum = t;
     // Algebraically, c should always be zero (but in IEEE754 FP math it isn't). 
     // Beware of overly-aggressive optimizing compilers! --ffast-math breaks this code
     // The lost low part will be added to y in a fresh attempt.
  }
  return sum;
}

我过去曾使用过这个，我认为完全按照上面的表述进行操作。但现在在这里输入它作为答案我想知道为什么返回值不是

sum - c

。在我看来，循环终止是在未应用最终待定补偿

c

的情况下发生的。其他人怎么看？

在同一个 Wiki 页面上，Nuemaier 和 Klein 提出了一些较新的、更复杂的算法，具有更好的校正能力，但速度较慢。 Kahan 的方法值得更多人了解（对于不提供硬件 80 位浮点支持的编译器特别有用）。