logo标签列表

采用成对求和法，需要多少项才能得到一个明显错误的结果？

pythonalgorithmnumpyfloating-pointsum
25

使用给定种类的浮点数（比如float16），构建完全错误的求和是很简单的。例如，使用Python / NumPy：

import numpy as np

one = np.float16(1)
ope = np.nextafter(one,one+one)

np.array((ope,one,-one,-one)).cumsum()
# array([1.001, 2.   , 1.   , 0.   ], dtype=float16)

在这里，我们使用cumsum来强制执行朴素求和。如果由于其自身原因未被限制，numpy将使用不同的求和顺序，从而得出更好的答案：

np.array((ope,one,-one,-one)).sum()
# 0.000977

以上是基于取消操作的。为了排除这类示例，让我们只允许非负项。对于朴素求和，仍然可以给出非常错误的总和的例子。以下是每个相等于10^-4的10^4个相同项的总和：

以上是基于取消操作的。为了排除这类示例，让我们只允许非负项。对于朴素求和，仍然可以给出非常错误的总和的例子。以下是每个相等于10^-4的10^4个相同项的总和：

np.full(10**4,10**-4,np.float16).cumsum()
# array([1.0e-04, 2.0e-04, 3.0e-04, ..., 2.5e-01, 2.5e-01, 2.5e-01],
  dtype=float16)

最后一项多除以了4。

同样地，允许numpy使用成对求和会得到更好的结果：

np.full(10**4,10**-4,np.float16).sum()
# 1.0

可以构造出比成对求和更好的求和方式。选择小于1的分辨率eps，我们可以使用1、eps、0、eps、3x0、eps、7x0、eps、15x0等，但这需要使用非常多的项。

我的问题是：使用float16和仅使用非负项，需要多少项才能获得一个比成对求和结果高至少一倍的结果。

附加题: 使用“正数”而不是“非负数”的情况下，相同的问题。这是否可能？

啊，不错。确实疯狂。 - Mark Dickinson
@MarkDickinson 是的，我认为这解决了问题。我已经设计好了这个问题。我很感激你从你的思考中提供一个部分答案。 - Paul Panzer
稍微与这个问题相关的是：在进行Knuth求和时可以避免“灾难性抵消”。这本质上是通过模拟具有两倍精度的总和来跟踪精度损失。（参见“复杂浮点算术中的准确求和、点积和多项式评估”） - kvantour
1
FYI，我已经开始实现DP了。 - David Eisenstat
深度为1432（即2^1432项）足以使真实总和超过计算总和两倍。 - David Eisenstat
显示剩余6条评论
3个回答
12
深度为1432（因此有2^1432个术语）足以使真实总和超过计算总和两倍。
我有一个想法，可以确定少于两倍所需的项数。
我们使用动态规划来回答以下问题：给定深度d和目标浮点总和s，具有配对和s的2^d个非负float16的最大真实总和是多少？
让这个数量为T(d，s)。我们得到一个递归式。
T(0, s) = s,    for all s.
T(d, s) =            max            (T(d-1, a) + T(d-1, b)),    for all d, s.
          a, b : float16(a + b) = s

每一步的递归都需要循环大约 2^29 个组合（因为我们可以假设 a ≤ b，并且负浮点数和特殊值是不允许的），而所需的深度不会超过 Hans 和你的答案所说的 10^4 左右。对我来说似乎是可行的。

DP 代码：

#include <algorithm>
#include <cstdio>
#include <vector>

using Float16 = int;
using Fixed = unsigned long long;

static constexpr int kExponentBits = 5;
static constexpr int kFractionBits = 10;
static constexpr Float16 kInfinity = ((1 << kExponentBits) - 1)
                                     << kFractionBits;

Fixed FixedFromFloat16(Float16 a) {
  int exponent = a >> kFractionBits;
  if (exponent == 0) {
    return a;
  }
  Float16 fraction = a - (exponent << kFractionBits);
  Float16 significand = (1 << kFractionBits) + fraction;
  return static_cast<Fixed>(significand) << (exponent - 1);
}

bool Plus(Float16 a, Float16 b, Float16* c) {
  Fixed exact_sum = FixedFromFloat16(a) + FixedFromFloat16(b);
  int exponent = 64 - kFractionBits - __builtin_clzll(exact_sum);
  if (exponent <= 0) {
    *c = static_cast<Float16>(exact_sum);
    return true;
  }
  Fixed ulp = Fixed{1} << (exponent - 1);
  Fixed remainder = exact_sum & (ulp - 1);
  Fixed rounded_sum = exact_sum - remainder;
  if (2 * remainder > ulp ||
      (2 * remainder == ulp && (rounded_sum & ulp) != 0)) {
    rounded_sum += ulp;
  }
  exponent = 64 - kFractionBits - __builtin_clzll(rounded_sum);
  if (exponent >= (1 << kExponentBits) - 1) {
    return false;
  }
  Float16 significand = rounded_sum >> (exponent - 1);
  Float16 fraction = significand - (Float16{1} << kFractionBits);
  *c = (exponent << kFractionBits) + fraction;
  return true;
}

int main() {
  std::vector<Fixed> greatest0(kInfinity);
  for (Float16 a = 0; a < kInfinity; a++) {
    greatest0[a] = FixedFromFloat16(a);
  }
  for (int depth = 1; true; depth++) {
    auto greatest1 = greatest0;
    for (Float16 a = 1; a < kInfinity; a++) {
      Fixed greatest0_a = greatest0[a];
      for (Float16 b = a; b < kInfinity; b++) {
        Float16 c;
        if (!Plus(a, b, &c)) {
          continue;
        }
        Fixed& value = greatest1[c];
        value = std::max(value, greatest0_a + greatest0[b]);
      }
    }

    std::vector<double> ratios;
    ratios.reserve(kInfinity - 1);
    for (Float16 a = 1; a < kInfinity; a++) {
      ratios.push_back(greatest1[a] / static_cast<double>(FixedFromFloat16(a)));
    }
    std::printf("depth %d, ratio = %.17g\n", depth,
                *std::max_element(ratios.begin(), ratios.end()));
    greatest0.swap(greatest1);
  }
}

我会运行这个程序，并在完成后发布更新。

9

如果允许零，则需要大量的术语，这几乎是不可能的；如果不允许零，由于溢出，实际上是不可能的。维基百科总结了一些由Nicolas Higham引起的error bounds。由于所有项都是非负数，条件数为1，因此n个项的相对误差被限制为|En|/|Sn| ≤ ε log2 n / (1 - ε log2 n)，其中ε是机器精度。要偏差两倍，我们需要|En| ≥ |Sn|，这只有在ε log2 n ≥ 1/2时才可能，这等价于n ≥ 21/(2 ε) = 21024（对于float16）。

0

剩下的问题是，如果您允许在求和中添加零（*），那么总和是否如此尖锐，以至于您可以通过配对求和获得2的相对误差。

简单的答案是肯定的，通过使用指数数量的零填充cum-sum序列的错误部分，方法如下（其中a1，a2，a3，... an对于普通求和具有问题）：

a1,
a2,
a3, 0,
a4, 0, 0, 0,
a5, 0, 0, 0, 0, 0, 0, 0,
a6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
...

对于成对求和，它将生成相同的总和，并产生相同的舍入误差，而您仅需要2 **（n-1）个术语，而不是n个。因此，由于10 ** 4 个术语可以为正常求和生成4的因子，则2 **（10 ** 4-1）个术语可以为成对求和提供4的因子。

*：David Eistenstat的答案表明，在禁止零之前，总和将溢出并变得有问题。（我假设成对求和递归到最后。）

谢谢你的回答。但恐怕它并没有提供任何新的东西。你所提出的构造方式与问题中已经描述的方式相同（倒数第二段：“可以构造和…”），只是后者明显节省了初始斜坡的时间。我希望得到更好的界限，即比David Eisenstat提供的下界更大和/或比零填充更低的下界（或具有更好非零元素的零填充）。 - Paul Panzer
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接