誤差のない浮動小数点数の sum

仕事で float32 で sum を計算したら結構な誤差が出るということがあって、それ自体は float64 にすれば満足な精度にはなったのだけど、そもそも計算の方法を工夫すれば誤差は抑えられる気がして、試してみた。

まず、浮動小数点数の誤差は大きい数と小さい数を足し合わせたときに発生するので、要するに小さい数は小さい数同士、大きい数は大きい数同士で足せれればよいのではないかと考え、絶対値でソートしてみた。

func mysum(values []float64) float64 {
  sorted := make([]float64, len(values))
  copy(sorted, values)
  sort.Slice(sorted, func(i, j int) bool {
    return math.Abs(sorted[i]) < math.Abs(sorted[j])
  })
  var result float64
  for _, x := range sorted {
    result += x
  }
  return result
}

これでも何も考えずに足すよりはよいのかもしれないが、以下のような凶悪な入力を与えるとあっさり破綻し、1.7484020030031562e+86 のようなめちゃくちゃな結果になる。

func main() {
  repeat := 10000
  part := []float64{1, 1e100, 1, -1e100}
  input := make([]float64, 0, len(part)*repeat)
  for i := 0; i < repeat; i++ {
    input = append(input, part...)
  }
  {
    // shuffle
    rand.Seed(time.Now().UnixNano())
    n := len(input)
    for i := n - 1; i >= 0; i-- {
      j := rand.Intn(i + 1)
      input[i], input[j] = input[j], input[i]
    }
  }

  fmt.Println("mysum", mysum(input))
}

ソートした上で分割統治で足し合わせながら絞り込んでいくという方法も試してみた。こちらも変わらず全然ダメ。

func mysum2(values []float64) float64 {
  sorted := make([]float64, len(values))
  copy(sorted, values)
  sort.Slice(sorted, func(i, j int) bool {
    return math.Abs(sorted[i]) < math.Abs(sorted[j])
  })
  for step := 1; step <= len(sorted); step *= 2 {
    for i := 0; i < len(sorted)-step; i += step * 2 {
      sorted[i] += sorted[i+step]
      sorted[i+step] = 0
    }
  }
  return sorted[0]
}

調べてみると、Stack Overflow のやりとりを発見。論文があるらしい (.ps ファイルで読めない)。

コードを見るに、ソートは特にしていない。ただ、足し算をしたときに誤差が発生したらそれも足していっているのは分かる。よさそうなのは感覚的に分かるが、これで完璧になる理屈はちゃんと論文を読まないと分からなさそう。

Go に移植もしてみた。実際誤差が出ないことも確認。すばらしい。

https://github.com/mono0x/go-msum

func Sum(values []float64) float64 {
  var partials []float64
  for _, x := range values {
    i := 0
    for _, y := range partials {
      if math.Abs(x) < math.Abs(y) {
        tmp := x
        x = y
        y = tmp
      }
      hi := x + y
      lo := y - (hi - x)
      if lo != 0.0 {
        partials[i] = lo
        i++
      }
      x = hi
    }
    if i < len(partials) {
      partials[i] = x
      partials = partials[0 : i+1]
    } else {
      partials = append(partials, x)
    }
  }
  var result float64
  for _, x := range partials {
    result += x
  }
  return result
}

カハンの加算アルゴリズムというものも有名らしい。

func kahanSum(values []float64) float64 {
  sum := 0.0
  c := 0.0
  for _, x := range values {
    y := x - c
    t := sum + y
    c = (t - sum) - y
    sum = t
  }
  return sum
}

こちらは今回試した凶悪な入力に対しては正しい結果を出せなかったが、こちらのほうが速そうに見えるので、適用しようとしている問題で十分な精度が出るのであれば、使ってみてもいいのかもしれない。

ただ、最初に触れたようにまずは float の精度を上げられるなら上げてから考えたほうがよいだろう。