平方根を使わずに距離を求める のバックアップ(No.3)

更新


公開メモ

平面状の2点間の距離を平方根を使わずに計算する

こちらの記事で、

http://d.hatena.ne.jp/nowokay/20120604

2点間の距離を平方根を使わず近似的に求める手法を知り、その精度について詳しく調べてみたのがこのページです。

最終的に、このようにすることで 約 1% の誤差で距離を近似できることが分かりました。

リンク先の内容を Mathematica で確かめる

正確な距離関数はこう:

d(x,y)=\sqrt{x^2+y^2}

LANG:Mathematica
distance[x_, y_] := Sqrt[x^2 + y^2]
Plot3D[distance[x, y], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large]

distance_wo_sqrt1.jpg

最大値、最小値を利用した近似はこう:

\tilde d(x,y;a,b)=a\,\mathrm{max}(|x|,|y|)+b\,\mathrm{min}(|x|,|y|)

LANG:Mathematica
distance2[x_, y_, a_, b_] := a Max[Abs[x], Abs[y]] + b Min[Abs[x], Abs[y]]
Plot3D[distance2[x, y, 1007/1024, 441/1024], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large]

distance_wo_sqrt2.jpg

プロットには、上記リンク先で紹介されていた係数 (a=1007/1024, b=441/1024) を使った。

近似値から真値を引いた誤差を表示すると:

e(x,y;a,b)=\tilde d(x,y;a,b)-d(x,y)

LANG:Mathematica
Plot3D[distance2[x, y, 1007/1024, 441/1024] - distance[x, y], 
  {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]

distance_wo_sqrt3.jpg

誤差は -2% から +8% くらいまでの間に分布する。

どうやらこのパラメータは近似値が真の距離を大きく下回らないように選ばれているようです?

とはいえ、 \theta=n\pi/2 近辺では少々下回っているという謎仕様。。。

リンク先のコメント欄で紹介されている max + (3/8) min という非常に単純化された処理でも 誤差は -2.5% から +7% の間に収まるので、ちょっとこれはいただけない。

パラメータの最適化を考える

エラーの相対誤差の二乗を最小化してみる。

e_1(a,b)=\int_{10^{-5}}^1 dx\,\int_{10^{-5}}^1 dy\,|e(x,y;a,b)|^2/d(x,y)^2

LANG:Mathematica
error[a_, b_] := NIntegrate[
 (distance2[x, y, a, b] - distance[x, y])^2/
  distance[x, y]^2, {x, 10^(-5), 1}, {y, 10^(-5), 1}]
ContourPlot[error[a/1024, b/1024], {a, 500, 1500}, {b, 0, 1000}, ImageSize -> Large]

distance_wo_sqrt4.jpg

横軸が 1024 a、縦軸が 1024 b で、青い部分が誤差が小さいところ。

最小のエラーを与えるパラメータを求める。

LANG:Mathematica
FindMinimum[error[a/1024, b/1024], {a, 1000}, {b, 400}]
Out[]= {0.000552491, {a -> 963.831, b -> 419.924}}

プロットしてみる。

LANG:Mathematica
Plot3D[distance2[x, y, 964/1024, 420/1024] - distance[x, y], 
  {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]

distance_wo_sqrt5.jpg

誤差は -6% から +2.5% の範囲に分布する。

誤差がゼロの上下に分散したため、上記のパラメータを使うより絶対値が小さくなっている。

ただし、「誤差の二乗の平均値」を最小にしたため、 \theta=n\pi/4 の周辺の狭い範囲では大きなエラーが出るままだ。

最適値の周りでの誤差の変化を見てみる。

LANG:Mathematica
ContourPlot[error[a/1024, b/1024], 
  {a, 964 - 10, 964 + 10}, {b, 420 - 10, 420 + 10}, ImageSize -> Large]

distance_wo_sqrt6.jpg

左下から右上への成分と、左上から右下への成分とに分けられそう。

まず最小値付近を2次関数でフィッティングして、

LANG:Mathematica
Fit[Flatten[Table[{a, b, 
   error[(963.8306122312068` + a)/1024, (419.92395268555106` + b)/
     1024]}, {a, -1, 1}, {b, -1, 1}], 1], {a^2, a b, b^2, 1}, {a, b}]
Out[]= 0.000552491 + 7.48994*10^-7 a^2 + 6.61035*10^-7 a b + 2.04662*10^-7 b^2

係数を成分に分解します。

LANG:Mathematica
Eigensystem[{{7.48994, 6.61035/2}, {6.61035/2, 2.04662}}] // N
Out[]= {{9.04982, 0.486736}, {{-0.904343, -0.426806}, {0.426806, -0.904343}}}

それぞれの成分に相当する関数をプロットすると、

固有値 9.0 の主成分は、

distance_wo_sqrt7.jpg

固有値 0.5 の副次成分は、

distance_wo_sqrt8.jpg

のような形をしている。

主成分のみでほぼ完全に距離関数を復元しているのに対して、
副次成分は \theta=n\pi/2 \theta=n\pi/2+\pi/4 とで異なる符号を持ち、 一方を上げると一方が下がるために平均的な誤差への影響が小さいことが見て取れる。

誤差の最大値を最小化する

誤差の評価を、

e_2(a,b)=\mathrm{max}\Big(\ \big|\tilde d(x,y;a,b)-d(x,y)\big|\ \Big|\ \sqrt{x^2+y^2}=1\ \Big)

のようにして、誤差の最大値を最小化することを考える。

パラメータ a, b に対して誤差をプロットすると、

LANG:Mathematica
error2[a_, b_] := 
  Table[
    With[{x = Cos[t], y = Sin[t]}, 
      Abs[distance2[x, y, a, b] - distance[x, y]]
    ], {t, 0, Pi/4, Pi/4/1024}
  ] // Max
ContourPlot[error2[a/1024, b/1024], {a, 500, 1500}, {b, 0, 1000}, ImageSize -> Large]

distance_wo_sqrt9.jpg

拡大して、

LANG:Mathematica
ContourPlot[error2[a/1024, b/1024], 
  {a, 984 - 10, 984 + 10}, {b, 407 - 10, 407 + 10}, ImageSize -> Large, PlotPoints -> 100]
Plot3d[error2[a/1024, b/1024], 
  {a, 984 - 10, 984 + 10}, {b, 407 - 10, 407 + 10}, ImageSize -> Large, PlotPoints -> 100]

distance_wo_sqrt10.jpg distance_wo_sqrt10b.jpg

最小値を求める(始め a=1000, b=400 から探したところうまく見つからなかった)。

LANG:Mathematica
FindMinimum[error2[a/1024, b/1024], {a, 984}, {b, 407}]
Out[]= {0.0395663, {a -> 983.485, b -> 407.372}}
Plot3D[distance2[x, y, 983/1024, 407/1024] - distance[x, y], 
   {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 80]

distance_wo_sqrt11.jpg

(983/1024)\mathrm{max}+(407/1024)\mathrm{min} という近似で誤差は \pm 4 % の範囲に収まった。

※ 下で見るように、 (123/128)\mathrm{max}+(51/128)\mathrm{min} でもほぼ 4% に収まる

もうちょっとちゃんと考えてみる

リンク先でも指摘されているように、

&math( \begin{cases} x=r\cos\theta\\ y=r\sin\theta \end{cases} );

とすると、

d(x,y)=\sqrt{x^2+y^2}=|r|

&math( \tilde d(x,y;a,b) &=a\,\mathrm{max}(|x|,|y|)+b\,\mathrm{min}(|x|,|y|)\\ &=|r|\,\big[a\,\mathrm{max}(|\cos\theta|,|\sin\theta|)+b\,\mathrm{min}(|\cos\theta|,|\sin\theta|)\big]\\ );

となるから、相対的な近似精度を見るには r=1 と固定して \theta の関数として扱えば良い。

そこで、

d^*(\theta;a,b)=a\,\mathrm{max}(|\cos\theta|,|\sin\theta|)+b\,\mathrm{min}(|\cos\theta|,|\sin\theta|)

を、上記で求めた「誤差の最大値を最小化するパラメータ」に対してプロットすると、

LANG:mathematica
d2[t_, a_, b_] := a Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]] + b Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]
Plot[{1, d2[t, 983/1024, 407/1024]}, {t, 0, 2}, ImageSize -> Large]

distance_wo_sqrt12.svg

となる。(横軸は \pi を単位とした方位角)

どうしてこんな関数になるかを理解するため、 \mathrm{max}(|\cos\theta|,|\sin\theta|) \mathrm{min}(|\cos\theta|,|\sin\theta|) とを別々にプロットしてみると、

LANG:mathematica
Plot[{Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]}, {t, 0, 2}]

distance_wo_sqrt13.svg

青が max、黄色が min。

これらの線形結合で最も 1 に近い関数を作ったのが上記の結果になる。

max も min も \theta=0\sim \pi/4 の範囲が反転しながら繰り返し現れるから、 誤差を最小化するには [0,\pi/4] の範囲だけ考えれば良いことも、リンク先で指摘されている通りだ。

LANG:mathematica
Plot[{Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]}, {t, 0, 1/4}]

distance_wo_sqrt14.svg

このプロットから、このパラメータに対する誤差の最大値は約 4% であることが分かる。

さらに誤差を小さくするには min の2倍の周期の三角波を

\big|\mathrm{min}(x,y)-(49/128) \tilde d(x,y;a,b)\big|

のようにして作って、

LANG:mathematica
Plot[{
  Max[Abs[Cos[t Pi]], Abs[Sin[t Pi]]], 
  Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]],
  Abs[Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]]-49/128 d2[t,a/1024,b/1024]]
}, {t, 0, 2}]

distance_wo_sqrt15.svg

LANG:mathematica
With[{a = 938, b = 384}, 
  Plot[{1, 
    d2[t, a/1024, b/1024] + 
    55/256 Abs[
      Min[Abs[Cos[t Pi]], Abs[Sin[t Pi]]] - 
      49/128 d2[t, a/1024, b/1024]])
  }, {t, 0, 1/4}, 
  PlotRange -> {Full, {0.98, 1.02}}, GridLines -> Automatic]
]

distance_wo_sqrt16.svg

とすると、最大誤差を 1% 程度に抑えられる。

ただ、計算量が増えてしまうのでもう少し何とかならないか?

計算量を減らす

計算量を増加させずに精度を上げるには、 リンク先と同様に min と max の大小関係で分岐するのが良さそう。

LANG:mathematica
Plot3D[Sqrt[x^2 + y^2], {x, -1, 1}, {y, -1, 1}, ImageSize -> Large, PlotPoints -> 100]
Plot3D[Module[{max, min},
  max = Max[Abs[x], Abs[y]];
  min = Min[Abs[x], Abs[y]];
  If[max < (5/2) min,
     864 max + 569 min,
    1016 max + 190 min]/1024
  ], {x, -1, 1}, {y, -1, 1}, 
  ImageSize -> Large, PlotPoints -> 100
]

distance_wo_sqrt19.jpg distance_wo_sqrt18.jpg

左が正確な距離関数、右が近似関数。

相対誤差を評価すると、

LANG:mathematica
Plot[{1, Module[{x, y, max, min},
  x = Cos[t Pi]; y = Sin[t Pi];
  max = Max[x, y]; min = Min[x, y];
  If[max < (5/2) min,
     864 max + 569 min,
    1016 max + 190 min]/1024
  ]}, {t, 0, 1/4}, 
  GridLines -> Automatic, PlotRange -> {0.98, 1.02}
]

distance_wo_sqrt17.svg

横軸は \pi を単位とする方位角、縦軸は相対誤差です。

max と min とを混ぜ合わせる簡単な処理で 1% 程度の精度で距離関数を再現できることが分かります。

コードにするとこんな感じ。

LANG:cpp
int distance2(int x, int y)
{
    int absx = abs(x);
    int absy = abs(y);
    int max = absx > absy ? absx : absy;
    int min = absx > absy ? absy : absx;

    int result;
    if( 2 * max < 5 * min ) {
        result = 864 * max + 569 * min;
    } else {
        result = 1016 * max + 190 * min;
    }

    return (result + 512) / 1024; // 四捨五入
}

掛け算を使わないで書くとこんな感じ。

LANG:cpp
int distance3(int x, int y)
{
    int absx = abs(x);
    int absy = abs(y);
    int max = absx > absy ? absx : absy;
    int min = absx > absy ? absy : absx;

    int result;
    if( (max << 1) < ((min << 2) + min) ) { // 2 * max < 5 * min
        // 864 * max + 569 * min
        // 864 = 0x360 = 0b001101100000
        // 569 = 0x239 = 0b001000111001
        result = ((max << 10) - (max << 7) - (max << 5)) + 
                 ((min << 9) + (min << 6) - (min << 3) + min);
    } else {
        // 1016 * max + 190 * min
        // 1016 = 0x3f8 = 0b001111111000
        // 190  = 0xbe  = 0b000010111110
        result = ((max << 10) - (max << 3)) +
                 ((min << 8) - (min << 6) - (min << 1));
    }

    return (result + 512) >> 10; // 四捨五入
}

こちらで

https://code.sololearn.com/

以下のコードを試したところ、 1000\le x,y\le 10000 の範囲で 相対誤差の最大値は 1.07576% でした。

LANG:cpp
#include <math.h>
#include <stdio.h>

double distance(int x, int y)
{
    return sqrt((double)x*x+(double)y*y);
}

int distance2(int x, int y)
{
    int absx = abs(x);
    int absy = abs(y);
    int max = absx > absy ? absx : absy;
    int min = absx > absy ? absy : absx;

    int result;
    if( 2 * max < 5 * min ) {
        result = 864 * max + 569 * min;
    } else {
        result = 1016 * max + 190 * min;
    }

    return (result + 512) / 1024;
}

int main() {
    int x, y;
    double max = 0;
    for(x=1000; x<10000; x+=10) {
        for(y=1000; y<10000; y+=10) {
            double d = distance(x, y);
            double d2 = distance2(x, y);
            if (max < abs(d-d2) / d)
                max = abs(d-d2) / d;
        }
    }
    printf("%lg", max);
    return 0;
}

もっと粗い近似で構わない場合

リンク先のコメント欄にある、

LANG:c
int distance(int x, int y) {
  x = abs(x); y = abs(y);
  int max = x>y ? x:y;
  int min = x>y ? y:x;
  return max + (min + min + min) >> 3; // max + (3/8)*min
}

これだけでも、誤差は -2.5% から 7% 程度に収まる。

\mathrm{max}+(3/8)\mathrm{min}
distance_wo_sqrt20.svg


分母を 32 にすれば -5% ~ 6% の範囲にできる。

\mathrm{max}+(11/32)\mathrm{min}
distance_wo_sqrt21.svg


max にも係数を付ければ分母が 128 の段階で誤差をほぼ 4% 以内に収められる。

(123/128)\mathrm{max}+(51/128)\mathrm{min}
distance_wo_sqrt22.svg

それにしても

LANG:c
int distance(int x, int y) {
  x = abs(x); y = abs(y);
  int max = x>y ? x:y;
  int min = x>y ? y:x;
  int d = 2*max<5*min ?
    864*max+569*min : 1016*max+190*min;
  return (d+512)>>10;
}

これで「原点からの距離」が 1% の精度で測れるとは驚きだ。

とはいえ

本当に必要なのは「距離の逆数」の場合が多いんだよなあ。。。

3次元では

LANG:c
int distance3d(int x, int y, int z)
{
  return distance2d(distance2d(x, y), z);
}

でいいはず。

質問・コメント





Counter: 14709 (from 2010/06/03), today: 4, yesterday: 3