楕円同士の接触判定と衝突判定

ググっても出てこなかったので。

2つの楕円が接している(内接 or 外接)かどうか判定する方法についてです。ついでに衝突判定もできます。

衝突判定だけしたい方

以下で説明する方法でも判定自体はできますが、非常に非効率です。悪いことは言いません。GJK法などを使いましょう。凸同士なので簡単にできます。

どうしても接触を判定したい方

心して読み進めてください。

事の発端

もちろん楕円が接するか判定する方法を知っているので、書くことにしました。

楕円の表現方法

楕円とはいわゆる2次曲線の一つで、円を線形に引き伸ばしたり潰したりすることで得られます。楕円の表現はいくつか考えられますが、ここでは正の対角行列 S と回転行列 R を用いて、原点中心の単位円を RS で線形変換(+必要に応じて平行移動)したものとして定義します。単位円は S で各軸方向に潰れて楕円になり、R で回転します。

S の1行1列成分が回転前の x 軸方向の半径、2行2列成分が回転前の y 軸方向の半径を表します。楕円の回転角の情報は回転行列 R が持っています。これで必要十分に楕円を表現することができます

以下、2つの楕円を考えます。楕円1を x_1 を中心として対角行列 S_1 と回転行列 R_1 で表される楕円とします。同様に楕円2を x_2, S_2, R_2 で表される楕円とします。

扱いやすいよう変形する

このままでは非常に扱いづらいので、2つの楕円の関係を保ったままできるだけ単純な形に持っていくことを考えます。

実は、あらゆる縮退しないアフィン変換[1]は楕円を保ちます。すなわち、好き勝手に空間全体に線形変換をかけても2つの楕円は依然として楕円の姿を保ちます。

楕円1を単位円に持っていく

さて、単位円を拡大縮小・回転して得られるのが楕円でしたから、適当な線形変換を空間全体にかけることで、楕円1だけでも単位円に持っていくことが可能であることが分かります。これは単純に、空間全体に R_1S_1 の逆行列 (R_1S_1)^{-1}=S_1^{-1}R_1^T を適用することで達成されます。

さて、このとき楕円2はどうなるかというと、単位円を R_2S_2 で変換したあとでさらに S_1^{-1}R_1^T という変換を受けたわけですから、全体として単位円を S_1^{-1}R_1^TR_2S_2 という行列で変換(+平行移動)した楕円になっています。

ここで困ったことが起きます。ここでの楕円の定義は「単位円を対角行列で拡大縮小したあと、回転行列で回転させたもの」でしたから、変換後の楕円2を楕円として認めるためには、S_1^{-1}R_1^TR_2S_2 もそのような行列の積(RS の形)として表現可能でなくてはなりません。

実は、そのような分解は不可能です。ある変換 A が回転行列 R と対角行列 S を用いて A=RS の形で表されるならば、自身の転置行列との積 A^TA

(RS)^TRS &= S^TR^TRS\\
&= SR^{-1}RS\\
&= S^2

と、対角行列にならなければなりませんが、S_1^{-1}R_1^TR_2S_2 を転置して掛けてみたところで

(S_1^{-1}R_1^TR_2S_2)^T(S_1^{-1}R_1^TR_2S_2) &= S_2^TR_2^TR_1S_1^{-T}S_1^{-1}R_1^TR_2S_2 \\
&= S_2R_2^TR_1S_1^{-2}R_1^TR_2S_2

で詰まってしまい、対角行列になる気配などありません。

好き勝手に空間全体に線形変換をかけても2つの楕円は依然として楕円の姿を保ちます。

これは嘘だったのでしょうか。

特異値分解

確かに任意の正則な線形変換 A は、単位円をある楕円に写します。しかし、単位円を同じ楕円に写すような線形変換は他にも多数存在します

例えば、ある楕円を表現する行列が RS であるとします。ここでもう一つ回転行列 R' を取り出して、RSR' という変換を考えてみます。実は、RSR' は単位円を RS と同じ楕円に写します。理由は単純で、単位円に最初にかかる変換である R' は単位円の形を変えないからです。

さらに、このような形で構成される変換は、全ての縮退しない線形変換を尽くします。すなわち、任意の正則な線形変換 A について、ある正の対角行列 S と回転行列 R, R' が存在し、A=RSR' を満たします。これが行列の特異値分解と呼ばれるものです[2]

「あらゆる(正則な)線形変換による円の像が楕円になる」ことをもって特異値分解の存在性に対する直観的な理解が可能になります。逆もまた然りです。この2つは密接に関連しています。

楕円の表現を見つける

特異値分解の存在によって、S_1^{-1}R_1^TR_2S_2=R_{12}S_{12}R' を満たすような正の対角行列 S_{12} と回転行列 R_{12}, R' が存在することが分かりました。このうち R_{12}S_{12} を取り出せば、変換後の楕円2の姿(=長軸と短軸の長さはいくらか、どの程度回転しているか)が明らかになります。

これは次のようにして求めることができます。まず自身の転置行列との積を求めます。

(S_1^{-1}R_1^TR_2S_2)(S_1^{-1}R_1^TR_2S_2)^T &= (R_{12}S_{12}R')(R_{12}S_{12}R')^T \\
&= R_{12}S_{12}R'R'^TS_{12}^TR_{12}^T \\
&= R_{12}S_{12}^2R_{12}^T

R' が消えて R_{12} によって S_{12}^2 に対角化される行列になりました。逆説的に、この行列を対角化すれば R_{12}S_{12}^2 が得られることになります。対角化の方法はさすがに省略しますが、2次正方行列の場合はこちらのページが参考になります。

マルペケさんによる楕円同士の衝突判定の近似の解説ページでも同様の変換が解説されていますが、④でやっている変換は本質的には特異値分解に他なりません

変換後の楕円2の軸を揃え原点へ

変換後の楕円2の正体が分かったところで、この楕円の軸が平行になるよう R_{12}^{-1}=R_{12}^T によって回転し、さらに中心が原点に来るよう全体を平行移動します。

これで扱いやすくするための変換は完了しました。状況を整理しましょう。

原点には軸に沿った変換後の楕円2があり、x 軸方向の半径は S_{12} の1行1列成分、y 軸方向の半径は S_{12} の2行2列成分です。

さらに、単位円に変換された元楕円1が、空間のどこかに浮いています。最終的に空間全体にかけた変換は、平行移動を除くと R_{12}^{-1}(R_1S_1)^{-1}=R_{12}^TS_1^{-1}R_1^T であるため、この単位円の中心座標は R_{12}^TS_1^{-1}R_1^T(x_1-x_2) で求められます。

変数の再定義

状況が簡単になったので、改めて変数を置き直します。以下、「変換後の楕円2」を単に「楕円2」、「単位円になった元楕円1」を「円1」と呼ぶことにします。また、対角行列 S_{12}単に S と表記します。

楕円2の x 軸, y 軸方向の半径をそれぞれ a, b とします。また、円1の中心座標を (x_c,y_c) とします。これを点c と呼びます。

接触判定を行う

あとは「2つの曲線 x^2/a^2+y^2/b^2=1(x-x_c)^2+(y-y_c)^2=1 は接するか?」という問題になります。

ここからのアプローチはいくつかあると思いますが、自分はベクトルと行列で解ける問題に傾きや角度を持ち込むのがあまり好きではない[3]なので、三角関数を使わない方針で解くことにします。

いきなり単位円と接するか調べるのではなく、点c を中心とする円が楕円2と接するような半径を求め、それが1であるかどうかを調べるという方針にします。円と接するような楕円上の点では、楕円のどちらの方向に動いても円の中心との距離が小さく(または大きく)なるためです[4]

次の最適化問題を考えます。

\text{minimize}\quad&(x-x_c)^2+(y-y_c)^2\\
\text{subject to}\quad&\frac{x^2}{a^2}+\frac{y^2}{b^2}=1,\ x,y\in\mathbb R

すなわち、楕円2上で点c との距離を最小化する問題を考えます。すると、目的関数 (x-x_c)^2+(y-y_c)^2 が楕円2上で極値を取るような場合に限り、点c を中心とする円と接することが分かります。

また、楕円の配置によっては以下のように最大4つの潜在的な接触点が発生します。

極値を求める

ではそのような半径を求めていきます。先程の最適化問題における制約を、単位ベクトルを用いて書き直します。楕円2上の任意の点はある単位ベクトル n=(n_x,n_y) を用いて Sn と表せるので、「楕円2上に点が存在する」という制約は「ベクトル n のノルムが1である」という制約に焼き直すことができます。

\text{minimize}\quad&\lVert Sn-c\rVert^2\\
\text{subject to}\quad&\lVert n\rVert^2=1,\ n\in\mathbb R^2

解きやすそうな形になってきました。さらにラグランジュの未定乗数法を適用し、制約を消します。

\text{minimize}\quad&\sup_{\lambda\in\mathbb R}\left(\lVert Sn-c\rVert^2-\lambda(\lVert n\rVert^2-1)\right)\\
\text{subject to}\quad&n\in\mathbb R^2

追加の変数(a.k.a. 未定乗数) \lambda を用意し、\lambda を自由に動かして値を最大化させた結果を最小化する問題になりました。n が単位ベクトルである制約は消えましたが、\lVert n\rVert\neq1 となる場合は \lambda を自由にとれば発散してしまうので、最適化問題の解においては自動的に制約が満たされます

さて、元々求めたかったものは最適化問題の解ではなく、目的関数の極値(から得られる円の半径)でした。実はこれも、変換後の目的関数の偏微分から得ることができます。

上限を取る前の目的関数を \lambdan で偏微分します。

\frac\partial{\partial\lambda}\left(\lVert Sn-c\rVert^2-\lambda(\lVert n\rVert^2-1)\right)
&= -(\lVert n\rVert^2-1) \\
\nabla_n\left(\lVert Sn-c\rVert^2-\lambda(\lVert n\rVert^2-1)\right)
&= 2S^2n-2Sc-2\lambda n\\
&= 2(S^2-\lambda I)n-2Sc

ベクトルによるスカラー関数の微分は、以下の公式が役立ちます。

\nabla_x(x^TAx) &= (A+A^T)x\\
\nabla_x(x^Tb) &= b

\lambdan 双方による偏微分の結果が0になるとき、元の問題の極値[5]が求まります。書き直すと以下のようになります。

\lVert n\rVert^2 &= 1 \\
(S^2-\lambda I)n &= Sc

ここでよく考えると、行列 S^2-\lambda I が特異になる場合というのは、点c が座標軸上に存在し、かつ楕円2と2点で同時に接する場合であるということが分かります。

この場合、求まった \lambda に対して一意な単位ベクトル n を定めることができないので、行列が正則でないという内容と直観的にも合致しますね。

このようなケースについては簡単に判定と処理ができる[6]ため、以下では S^2-\lambda I は正則であると仮定します。すると次のように n を消去できます[7]

\left\lVert \frac{S}{S^2-\lambda I}c\right\rVert^2 = c^T\frac{S^2}{(S^2-\lambda I)^2}c = 1\\

展開すると以下のようになります。

\frac{a^2}{(a^2-\lambda)^2}x^2+\frac{b^2}{(b^2-\lambda)^2}y^2=1

分母を消すと、4次方程式が出てきます。極値の個数の最大値に一致していますね。

a^2(b^2-\lambda)^2x_c^2+b^2(a^2-\lambda)^2y_c^2=(a^2-\lambda)^2(b^2-\lambda)^2

ここから先は人力では辛いものがあるので、MATLAB に係数を計算させます。すると次のような方程式であることが分かります。

a_0+a_1\lambda+a_2\lambda^2+a_3\lambda^3+a_4\lambda^4=0
a_0 &= a^2b^2(a^2b^2 - a^2y_c^2 - b^2x_c^2)\\
a_1 &= -2a^2b^2(a^2 + b^2 - x_c^2 - y_c^2)\\
a_2 &= a^4 + b^4 - a^2x^2 - b^2y^2 + 4a^2b^2\\
a_3 &= -2(a^2 + b^2)\\
a_4 &= 1

4次方程式は解析的に解けるので、解きます[8]。今回は虚数解には興味がないので、虚数解は無視します。

得られた各実数解 \lambda に対して、点c を中心とする円が接するために必要な半径が出てきます。

\lVert Sn-c\rVert &= \left\lVert\frac{S^2}{S^2-\lambda I}c-c\right\rVert\\
&= \left\lVert\left(\frac{S^2}{S^2-\lambda I}-I\right)c\right\rVert\\
&= \left\lVert\frac{\lambda I}{S^2-\lambda I}c\right\rVert\\
&= \sqrt{\frac{\lambda^2}{(a^2-\lambda)^2}x_c^2 + \frac{\lambda^2}{(b^2-\lambda)^2}y_c^2}

この値が 1 になるような \lambda が見つかれば、2つの楕円は接しています。お疲れさまでした。

楕円の衝突判定を行う

副産物として、楕円の衝突判定を厳密に行う手法が得られます。先程得られた半径のうち1より小さいものが存在すれば、円1よりも小さい半径で楕円2と接することになるので、円1は楕円2と重なっていることになります。

ただし、点c がそもそも楕円2の内部にある場合はこの判定法は成立しないので、c_x^2/a^2+c_y^2/b^2\lt1 を満たす場合は早期に重なっていると判定して終了します。もちろんコーナーケースの判定も忘れないようにしてください。

しかし、冒頭で述べたように、この方法は楕円の衝突判定を行うにはあまりに非効率です。衝突だけ調べたいのであれば、楕円の凸性を活かした反復手法による判定を行うことを強く推奨します。

追記: 衝突判定をしたい場合は円と楕円に変換してから Signed Distance Function を使う という方法もありますが、こちらも解析的に求まる方は4次方程式を直接解いているため、数値的に非常に安定するわけではないようです。

プログラム

実装にあたっては、Python など線形代数の各種操作が簡単にできる言語とライブラリを使うことをお勧めします。でないと特に4次方程式の求解などの非本質な部分で苦しむことになります。

苦しんだ場合の実装例全文を載せておきます(Haxe によるもの)。コーナーケースの処理は書かれていないので注意してください。

using Lambda;

// 楕円の定義
typedef Ellipse = {
    var x:Float;
    var y:Float;
    var a:Float; // x軸の半径
    var b:Float; // y軸の半径
    var angle:Float; // 角度
}

function isZero(x:Float):Bool {
    return x < 1e-6 && x > -1e-6;
}

// ベクトル
class Vec2 {
    public var x:Float;
    public var y:Float;

    public function new(x:Float, y:Float) {
        this.x = x;
        this.y = y;
    }

    public function dot(v:Vec2):Float {
        return x * v.x + y * v.y;
    }

    public function length():Float {
        return Math.sqrt(dot(this));
    }

    public function normalize():Vec2 {
        final l = 1 / length();
        return new Vec2(x * l, y * l);
    }
}

// 複素数を扱う用
@:forward(r, i)
abstract Complex(ComplexData) {
    public function new(real:Float, imag:Float) {
        this = new ComplexData(real, imag);
    }

    @:from
    static function fromFloat(f:Float):Complex {
        return new Complex(f, 0);
    }

    @:op(A + B)
    static function opAdd(a:Complex, b:Complex):Complex {
        return new Complex(a.r + b.r, a.i + b.i);
    }

    @:op(A - B)
    static function opSub(a:Complex, b:Complex):Complex {
        return new Complex(a.r - b.r, a.i - b.i);
    }

    @:op(A * B)
    static function opMul(a:Complex, b:Complex):Complex {
        return new Complex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
    }

    @:op(A / B)
    static function opDiv(a:Complex, b:Complex):Complex {
        final d = b.r * b.r + b.i * b.i;
        return new Complex((a.r * b.r + a.i * b.i) / d, (a.i * b.r - a.r * b.i) / d);
    }

    @:op(-A)
    static function opNeg(a:Complex):Complex {
        return new Complex(-a.r, -a.i);
    }

    @:op(A * B)
    @:commutative
    static function opMulFloat(a:Complex, b:Float):Complex {
        return a * b;
    }

    public function arg():Float {
        return this.r < 0 && isZero(this.i) ? Math.PI : Math.atan2(this.i, this.r);
    }

    public function abs():Float {
        return Math.sqrt(this.i * this.i + this.r * this.r);
    }

    public function pow(to:Float):Complex {
        return cast ComplexData.polar(Math.pow(abs(), to), arg() * to);
    }

    public function isReal():Bool {
        return isZero(this.i);
    }

    public function toString():String {
        return this.r + " + " + this.i + "i";
    }
}

// 補助データ
private class ComplexData {
    public var r:Float;
    public var i:Float;

    public function new(r:Float, i:Float) {
        this.r = r;
        this.i = i;
    }

    public static function polar(r:Float, arg:Float):ComplexData {
        return new ComplexData(Math.cos(arg) * r, Math.sin(arg) * r);
    }

    public function toString():String {
        return r + (i < 0 ? " - " + -i : " + " + i) + "i";
    }
}

// 行列
class Mat2 {
    public var e00:Float;
    public var e01:Float;
    public var e10:Float;
    public var e11:Float;

    public function new(e00:Float, e01:Float, e10:Float, e11:Float) {
        this.e00 = e00;
        this.e01 = e01;
        this.e10 = e10;
        this.e11 = e11;
    }

    public static function scale(x:Float, y:Float):Mat2 {
        return new Mat2(x, 0, 0, y);
    }

    public static function rotation(angle:Float):Mat2 {
        final s = Math.sin(angle);
        final c = Math.cos(angle);
        return new Mat2(c, -s, s, c);
    }

    public function mul(m:Mat2):Mat2 {
        return new Mat2(e00 * m.e00 + e01 * m.e10, e00 * m.e01 + e01 * m.e11, e10 * m.e00 + e11 * m.e10, e10 * m.e01 + e11 * m.e11);
    }

    public function mulVec(v:Vec2):Vec2 {
        return new Vec2(e00 * v.x + e01 * v.y, e10 * v.x + e11 * v.y);
    }

    // 転置
    public function transpose():Mat2 {
        return new Mat2(e00, e10, e01, e11);
    }

    // 対角化のための行列を計算する
    public function diagonalize():Mat2 {
        if (isZero(e01) && isZero(e10))
            return new Mat2(1, 0, 0, 1);
        final b = e00 + e11;
        final c = e00 * e11 - e01 * e10;
        final d = Math.sqrt(b * b - 4 * c);
        final l1 = 0.5 * (b - d);
        final l2 = 0.5 * (b + d);
        if (isZero(e01)) {
            final c0 = new Vec2(l1 - e11, e10).normalize();
            final c1 = new Vec2(l2 - e11, e10).normalize();
            return new Mat2(c0.x, c1.x, c0.y, c1.y);
        } else {
            final c0 = new Vec2(e01, l1 - e00).normalize();
            final c1 = new Vec2(e01, l2 - e00).normalize();
            return new Mat2(c0.x, c1.x, c0.y, c1.y);
        }
    }

    // 逆行列
    public function invert():Mat2 {
        final d = 1 / (e00 * e11 - e01 * e10);
        return new Mat2(e11 * d, -e01 * d, -e10 * d, e00 * d);
    }

    public function toString():String {
        return "[" + e00 + ", " + e01 + "; " + e10 + ", " + e11 + "]";
    }
}

function touching(e1:Ellipse, e2:Ellipse):Bool {
    // 各楕円の回転行列
    final r1 = Mat2.rotation(e1.angle);
    final r2 = Mat2.rotation(e2.angle);
    // 各楕円の拡縮回転行列
    final rs1 = r1.mul(Mat2.scale(e1.a, e1.b));
    final rs2 = r2.mul(Mat2.scale(e2.a, e2.b));
    final invRS1 = rs1.invert();
    // 楕円2の楕円1に対する相対変換行列
    final relT = invRS1.mul(rs2);
    // 特異値分解で楕円2の相対変換後の軸の長さを求める
    final relTT = relT.mul(relT.transpose());
    final relR2 = relTT.diagonalize();
    final sing2 = relR2.transpose().mul(relTT).mul(relR2);
    final singX = Math.sqrt(sing2.e00);
    final singY = Math.sqrt(sing2.e11);
    // 相対変換後の楕円2の拡縮行列
    final relS2 = Mat2.scale(singX, singY);
    // 相対変換後の楕円2の長軸と短軸に沿った座標系から見た楕円1の相対座標
    final relX1 = relR2.transpose().mul(invRS1).mulVec(new Vec2(e1.x - e2.x, e1.y - e2.y));

    final x = relX1.x;
    final y = relX1.y;
    final a = relS2.e00;
    final b = relS2.e11;
    final x2 = x * x;
    final y2 = y * y;
    final a2 = a * a;
    final b2 = b * b;

    // 4次方程式の係数
    final c0 = a2 * b2 * (a2 * b2 - a2 * y2 - b2 * x2);
    final c1 = -2 * a2 * b2 * (a2 + b2 - x2 - y2);
    final c2 = a2 * a2 + b2 * b2 - a2 * x2 - b2 * y2 + 4 * a2 * b2;
    final c3 = -2 * (a2 + b2);
    final c4 = 1;

    // 未定乗数
    final lambdas = realRootsQuart(c0, c1, c2, c3, c4);
    final rs = lambdas.map(l -> {
        final p = l / (a2 - l);
        final q = l / (b2 - l);
        Math.sqrt(p * p * x2 + q * q * y2);
    });
    // 0.95~1.05 の範囲にいれば接触とする
    return rs.exists(r -> r > 0.95 && r < 1.05);
}

// 4次方程式の実数解を列挙する
function realRootsQuart(a0:Float, a1:Float, a2:Float, a3:Float, a4:Float):Array<Float> {
    // 元の係数
    final o0 = a0;
    final o1 = a1;
    final o2 = a2;
    final o3 = a3;
    final o4 = a4;

    // 誤差を緩和するために解のスケーリングを行う
    final scale = (Math.pow(Math.abs(a0 / a4), 1 / 4) + Math.pow(Math.abs(a1 / a4), 1 / 3) + Math.pow(Math.abs(a2 / a4), 1 / 2) +
        Math.abs(a3 / a4)) * 0.25;
    a1 *= scale;
    a2 *= scale * scale;
    a3 *= scale * scale * scale;
    a4 *= scale * scale * scale * scale;

    // モニック化
    a0 /= a4;
    a1 /= a4;
    a2 /= a4;
    a3 /= a4;
    a4 = 1;

    final alpha = -3 * a3 * a3 / 8 + a2;
    final beta = a3 * a3 * a3 / 8 - a3 * a2 / 2 + a1;
    final gamma = -3 * a3 * a3 * a3 * a3 / 256 + a2 * a3 * a3 / 16 - a3 * a1 / 4 + a0;

    var cands:Array<Complex> = [];
    if (isZero(beta)) {
        final p:Complex = -a3 / 4;
        final r:Complex = alpha * alpha - 4 * gamma;
        final q = r.pow(0.5);
        cands.push(p + ((-alpha + q) * 0.5).pow(0.5));
        cands.push(p + ((-alpha - q) * 0.5).pow(0.5));
        cands.push(p - ((-alpha + q) * 0.5).pow(0.5));
        cands.push(p - ((-alpha - q) * 0.5).pow(0.5));
    } else {
        final p:Complex = -alpha * alpha / 12 - gamma;
        final q:Complex = -alpha * alpha * alpha / 108 + alpha * gamma / 3 - beta * beta / 8;
        final r = -q / 2 + (q * q / 4 + p * p * p / 27).pow(0.5);
        final u = r.pow(1 / 3);
        final y = -5 / 6 * alpha + (q.abs() >= u.abs() ? -q.pow(1 / 3) : u - p / (3 * u));
        final w = (alpha + 2 * y).pow(0.5);
        final s = -a3 / 4;
        final t = -(3 * alpha + 2 * y);
        final v = 2 * beta / w;
        cands.push(s + (w + (t - v).pow(0.5)) * 0.5);
        cands.push(s + (w - (t - v).pow(0.5)) * 0.5);
        cands.push(s + (-w + (t + v).pow(0.5)) * 0.5);
        cands.push(s + (-w - (t + v).pow(0.5)) * 0.5);
    }
    cands = cands.map(c -> c * scale);

    // 実数解のみを取り出す
    var res = cands.filter(a -> a.isReal()).map(a -> a.r);

    // 精度が不十分なことがあるので、ニュートン法を回して精度を高める
    res = res.map(r -> {
        for (iter in 0...8) {
            final t = o0 + o1 * r + o2 * r * r + o3 * r * r * r + o4 * r * r * r * r;
            final d = o1 + 2 * o2 * r + 3 * o3 * r * r + 4 * o4 * r * r * r;
            if (!isZero(d))
                r -= t / d;
        }
        r;
    });
    return res;
}

変な場所があったりしたら Twitter などで教えてもらえると助かります。


  1. 線形変換と平行移動を組み合わせたもの ↩︎
  2. ここでは \mathbb R^2 上正則な変換のみを考えましたが、特異値分解自体はあらゆる行列に対して考えることが可能です。 ↩︎
  3. 往々にしてシンプルさを欠いたり特異点が生じてロバストでない方法になったりするからです ↩︎
  4. 例外は楕円2が点c を中心とする円である場合です。 ↩︎
  5. 正確には停留点ですが、今回の問題では曲率の符号が変化しないため常に極値になります ↩︎
  6. 邪魔だった定数項が消えるため、全体として2次方程式を解くだけで済みます ↩︎
  7. S\lambda I も対角行列で積が可換なので実数のような扱いをしています ↩︎
  8. 複素数の四則演算と平方根・立方根をまともに考える必要があるので、実は残念ながらここで三角関数が出てきます…… ↩︎

このエントリーをはてなブックマークに追加