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

更新

[[公開メモ]]

#contents

* 平面上の2点間の距離を平方根を使わずに計算する [#v08053c6]

こちらの記事で、

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

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

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

* リンク先の内容を Mathematica で確かめる [#y6b0a391]

正確な距離関数はこう:

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

例えば点 &math((x_1,y_1)); と点 &math((x_2,y_2)); との間の距離を
&math(d(x_2-x_1,y2-y_1)); として求められる。

この関数をグラフにすると、

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

&attachref(distance_wo_sqrt1.jpg,,50%);

最大値、最小値を利用した近似関数のグラフはこう。

&math(\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]

&attachref(distance_wo_sqrt2.jpg,,50%);

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

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

&math(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]

&attachref(distance_wo_sqrt3.jpg,,50%);

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

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

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

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

* パラメータの最適化を考える [#p8c95460]

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

&math(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]

&attachref(distance_wo_sqrt4.jpg,,50%);

横軸が 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]

&attachref(distance_wo_sqrt5.jpg,,50%);

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

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

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

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

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

&attachref(distance_wo_sqrt6.jpg,,50%);

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

まず最小値付近を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 の主成分は、

&attachref(distance_wo_sqrt7.jpg,,50%);

固有値 0.5 の副次成分は、

&attachref(distance_wo_sqrt8.jpg,,50%);

のような形をしている。

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

* 誤差の最大値を最小化する [#jf7830a7]

誤差の評価を、

&math(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]

&attachref(distance_wo_sqrt9.jpg,,50%);

拡大して、

 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]

&attachref(distance_wo_sqrt10.jpg,,50%);
&attachref(distance_wo_sqrt10b.jpg,,50%);

最小値を求める(始め &math(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]

&attachref(distance_wo_sqrt11.jpg,,50%);

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

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

* もうちょっとちゃんと考えてみる [#jf9e0004]

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

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

とすると、

&math(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]\\
);

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

そこで、

&math(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]

&attachref(distance_wo_sqrt12.svg);

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

どうしてこんな関数になるかを理解するため、
&math(\mathrm{max}(|\cos\theta|,|\sin\theta|)); と
&math(\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}]

&attachref(distance_wo_sqrt13.svg);

青が max、黄色が min。

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

max も min も &math(\theta=0\sim \pi/4); の範囲が反転しながら繰り返し現れるから、
誤差を最小化するには &math([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}]

&attachref(distance_wo_sqrt14.svg);

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

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

&math(\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}]

&attachref(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]
 ]

&attachref(distance_wo_sqrt16.svg);

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

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

* 計算量を減らす [#qbba8f2b]

計算量を増加させずに精度を上げるには、
リンク先と同様に 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
 ]

&attachref(distance_wo_sqrt19.jpg,,50%);
&attachref(distance_wo_sqrt18.jpg,,50%);

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

相対誤差を評価すると、

 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}
 ]

&attachref(distance_wo_sqrt17.svg);

横軸は &math(\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/

以下のコードを試したところ、&math(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;
 }

* ビット幅などいろいろ工夫の余地は残る [#r129c2f4]

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

 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% 程度に収まる。

&math(\mathrm{max}+(3/8)\mathrm{min}); ~
&attachref(distance_wo_sqrt20.svg);

~

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

&math(\mathrm{max}+(11/32)\mathrm{min}); ~
&attachref(distance_wo_sqrt21.svg);

~

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

&math((123/128)\mathrm{max}+(51/128)\mathrm{min});~
&attachref(distance_wo_sqrt22.svg);

区間を分けない限り、これがほぼ最良なのは、この周辺での誤差分布をプロットすれば分かる。

 LANG:mathematica
 ContourPlot[
   Max[Table[
     Module[{x, y, max, min},
       x = Cos[t Pi]; y = Sin[t Pi];
       max = Max[x, y]; min = Min[x, y];
       Abs[(a max + b min)/128 - 1]
     ], 
     {t, 0, 1/4, 1/400}
   ]],
   {a, 123 - 2, 123 + 2}, {b, 51 - 2, 51 + 2},
   PlotPoints -> 400, PlotLegends -> Automatic,
   GridLines -> Automatic, Method -> {"GridLinesInFront" -> True}
 ]
 LANG:mathematica
 FindMinimum[
   Max[Table[
     Module[{x, y, max, min},
       x = Cos[t Pi]; y = Sin[t Pi];
       max = Max[x, y]; min = Min[x, y];
       Abs[(a max + b min)/128 - 1]
     ], 
     {t, 0, 1/4, 1/400}
   ]],
   {{a, 123}, {b, 51}}
 ]
 Out[]= {0.0396404, {a -> 122.926, b -> 50.9694}}

&attachref(distance_wo_sqrt25.png);

最小値は (a,b)=(122.926, 50.9694) の時の 0.03964

どうして三角形状の谷になっているかというと・・・

 Table[
   Module[{a, b},
     a = 122.92602775001741` + 3 Cos[tt]; 
     b = 50.96935967667149` + 3 Sin[tt];
     GraphicsRow[{
       Show[%110, ListPlot[{{a, b}}, PlotStyle -> White]],
       Plot[{1, Module[{x, y, max, min}, 
         x = Cos[t Pi]; y = Sin[t Pi];
         max = Max[x, y]; min = Min[x, y];
         (a max + b min)/128]
         }, {t, 0, 1/4}, PlotRange -> {0.9, 1.1}
       ]}, ImageSize -> Large
     ]
   ], {tt, 0, 2 Pi, Pi/16}
 ];
 Export["animation.gif", %%, "DisplayDurations" -> 0.2]

&attachref(distance_wo_sqrt26.gif);

パラメータを右上に動かすと中央部分、左や左下へ動かすと両端部分のエラーが増加する。

~

区間を2つに分ける場合、次の式なら分母を 512 にできて、精度は 1024 を分母とした上記のものより
むしろほんの少しだけ良くなる。ただし、加減算の回数が少し増える。

 LANG:cpp
 int d = 5 * max < 12 * min ? 430 * max + 287 * min : 507 * max + 101 * min;
 return (d+ 256) >> 9;
 //   5 = 0b0101 加減算 2 回
 //  12 = 0b1100 加減算 2 回
 // 430 = 0x1ae = 0b000110101110  加減算 3 回
 // 287 = 0x11f = 0b000100011111  加減算 2 回
 // 507 = 0x1f8 = 0b000111111000  加減算 1 回
 // 101 = 0x065 = 0b000001100101  加減算 3 回

&attachref(distance_wo_sqrt23.svg);

~

区間を3分割すれば 0.5% の精度が出る。

 LANG:cpp
 int d = 4 * max > 15 * min ? 510 * max +  67 * min :
         4 * max >  7 * min ? 476 * max + 195 * min :
                              409 * max + 312 * min ;
 return (d+ 256) >> 9;

&attachref(distance_wo_sqrt24.svg);

* 区分法による最適化の方針 [#p4f2d6ff]

区間を分けた場合、ベストな結果を与えるのは区間の両端の値が等しい時なので、
ここから a と b との関係が決まる。

つまり、区間 &math([t_1, t_2]); での最良値を与える a, b の間には、

&math(d^*(t_1;a,b)=d^*(t_2;a,b));

が成り立つので、ここから &math(b = A a); のような線形関係を導ける。

例えば &math([0,0.1]); に対してなら、

 LANG:mathematica
 Solve[d2[0, a, b] == d2[0.1, a, b], b][[1]] // N
 Out[]= {b -> 0.158384 a}

この制約の上で a を動かした際の誤差は

 LANG:mathematica
 maxerrorinsegment[t1_, t2_, a_] :=
   Module[{solution},
     solution = Solve[d2[t1, a, b] == d2[t2, a, b], b][[1]] // N;
     FindMaximum[d2[t, a, b /. solution], {t, (t1 + t2)/2}][[1]] //
     Max[
       Abs[d2[t1, a, b /. solution] - 1],
       Abs[#1 - 1]
     ] &
   ]
 Plot[maxerrorinsegment[0, 0.1, a], {a, 0.8, 1.2}]

&attachref(distance_wo_sqrt27.svg);

のようにほぼ一次関数になるので、

 LANG:mathematica
 optimizeinsegment[t1_, t2_] :=
   FindMinimum[
     maxerrorinsegment[t1, t2, a] // Hold,
     {a, 1}
   ] // Quiet
 optimizeinsegment[0, 0.1]
 Out[]= {0.00619396, {a -> 0.993806}}

のように最適化するのは簡単。

例えば [0, 0.1] の区間での最適値は a = 0.993806 で、そのときの最大誤差は 0.00619396 になる。

これを使って、区間を2つに区切る位置に対する左右区間の誤差の大きさを見てみると、

 LANG:mathematica
 Plot[{
     optimizeinsegment[0, t][[1]], 
     optimizeinsegment[t, 1/4][[1]]
   }, {t, 0.12, 0.13}, GridLines -> Automatic
 ]

&attachref(distance_wo_sqrt28.svg);

目の子で探した上記のパラメータがほぼ最適値になっていることを確認できる。

より広い範囲で見ると、

 LANG:mathematica
 Plot[{
     optimizeinsegment[0, t][[1]], 
     optimizeinsegment[t, 1/4][[1]]
   }, {t, 0, 1/4}
 ]

&attachref(distance_wo_sqrt29.svg);

のようになっていて、このグラフは区間をどれほど小さくすればどれだけの精度で近似できるかの指標を与える。

見やすくするために対数プロットすると、

 LANG:mathematica
 LogPlot[optimizeinsegment[0, t][[1]], {t, 0.005, 1/4}, 
   GridLines -> Automatic]

&attachref(distance_wo_sqrt30.svg);

横軸を N 分割の分割数 N にすると、

 LANG:mathematica
 ListLogPlot[Table[{n, optimizeinsegment[0, 1/4/n][[1]]}, {n, 1, 20}], GridLines -> Automatic]

&attachref(distance_wo_sqrt31.svg);

精度を1桁よくするのに分割数を約3倍にする必要がある感じだ。

* それにしても [#u4233dc8]

 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% の精度で測れるとは驚きだ。

* とはいえ [#tc1de7de]

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

* 3次元では [#f64e777f]

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

でいいはず。

* 質問・コメント [#uabcd372]

#article_kcaptcha

Counter: 14711 (from 2010/06/03), today: 3, yesterday: 0