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

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• バックアップ を表示
• プログラミング/平方根を使わずに距離を求める へ行く。
  • 1 (2018-03-26 (月) 17:15:34)
  • 2 (2018-03-26 (月) 23:15:29)
  • 3 (2018-03-27 (火) 02:27:36)
  • 4 (2018-03-27 (火) 12:25:59)
  • 5 (2018-03-28 (水) 17:09:05)
  • 6 (2018-03-28 (水) 17:09:05)
  • 7 (2018-03-29 (木) 14:03:17)
  • 8 (2018-03-29 (木) 22:11:14)
  • 9 (2018-03-29 (木) 22:11:14)
  • 10 (2018-04-03 (火) 10:27:36)
  • 11 (2018-04-03 (火) 10:27:36)
  • 12 (2019-04-10 (水) 13:06:58)
  • 13 (2019-04-10 (水) 13:06:58)

[[公開メモ]]

#contents

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

こちらの記事で、

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

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

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

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

正確な距離関数はこう：

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

&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%);

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

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

 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 の２倍の周期の三角波を

&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;
 }

* もっと粗い近似で構わない場合 [#u6b9ba84]

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

 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);


* それにしても [#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]

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

* ３次元では [#f64e777f]

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

でいいはず。

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

#article_kcaptcha

PukiWiki 1.5.0 Copyright © 2001-2006 PukiWiki Developers Team. License is GPL.
Based on "PukiWiki" 1.3 by yu-ji. Powered by PHP 7.4.33. HTML convert time: 0.003 sec.
Also based on Bootstrap, Start Bootstrap - SB Admin 2, highlight.js, Google fonts, Google recaptcha, Google chart, Bitly.