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

更新

公開メモ

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

こちらの記事で、

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

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

最終的に、このようにすることで 約 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 近辺では少々下回っているという謎仕様。。。

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

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

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]

distance_wo_sqrt10.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

誤差は \pm 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

このようになって、これらの線形結合で最も 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;
  int d = 2*max<5*min ?
    864*max+569*min : 1016*max+190*min;
  return (d+512)>>10;
}

これで「減点からの距離」が測れるとは驚きだ。

質問・コメント†




Counter: 18259 (from 2010/06/03), today: 15, yesterday: 2