平方根を使わずに距離を求める の履歴(No.5)

更新

公開メモ

平面上の２点間の距離を平方根を使わずに計算する†

こちらの記事で、

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

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

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

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

正確な距離関数はこう：

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

例えば点 (x_1,y_1) と点 (x_2,y_2) との間の距離を 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]

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

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

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

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

\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: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}}

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]

distance_wo_sqrt26.gif

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


区間を２つに分ける場合、次の式なら分母を 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 回

distance_wo_sqrt23.svg


区間を３分割すれば 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;

distance_wo_sqrt24.svg

区分法による最適化の方針†

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

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

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

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

例えば [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}]

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 になる。

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

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

distance_wo_sqrt28.svg

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

より広い範囲で見ると、

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

distance_wo_sqrt29.svg

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

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

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

distance_wo_sqrt30.svg

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

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

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

とはいえ†

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

３次元では†

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

でいいはず。

質問・コメント†




Counter: 16167 (from 2010/06/03), today: 6, yesterday: 6