二分法の処理アルゴリズムと応用例

簡単には解けない複雑な計算式の解を求める手法の一つに二分法という手法があります。

これは、2点間の間に必ず解がある場合、2点間の半分の位置の値を求め、解が小さい方にあるのか?大きい方にあるのか?を確認し、解がある方のさらに半分の位置の値を求め...と解のある位置を追い込む手法です。

 

説明が難しいので具体例で示します。

下図は y = x2 – 10 のグラフですが、y=0 となるxの値を求める場合の例です。

この例では、普通に解を求める事ができますが、あえて二分法で解きます。

 

y=0 となるxの値は、-4 ~ -2 と 2 ~ 4 の間にありそうな事はわかります。

今回は 2 ~ 4 の間にある解を求めます。

 

2と4の中間の3の時のyの値を計算すると、yの値は負になるので、y=0となるxの値は
3 ~ 4 の間にあることがわかります。

 

さらに、3と4の中間の3.5の時のyの値を計算すると、yの値は正になるので、
y=0となるxの値は3 ~ 3.5 の間にあることがわかります。

 

さらに、3と3.5の中間の3.25の時のyの値を計算すると、yの値は正になるので、
y=0となるxの値は3 ~ 3.25 の間にあることがわかります。

 

・・・と、この処理を2点間の範囲が十分小さくなるまで繰り返すことで、y=0 となる値を求める手法が二分法となります。

 

二分法の考え方は非常に簡単ですが、個人的には数式の解を求めるために使うことは、ほとんど無く、二分法の考え方を応用することはよくあります。

 

その一例を示したいと思います。

 

プログラムのエラーの場所を特定する場合

プログラム中のどこかでエラーが発生しているときに、どこで発生しているのか?さっぱり見当が付かない場合に二分法の考え方を応用します。

 

まず、エラーが発生しているプログラムにおいて、プログラムシーケンスのだいたい真ん中あたりとなる位置にブレークポイントを置いてデバック実行を行います。

もし、ブレークポイントまでエラーが無ければ、ブレークポイントの後半の処理でエラーが発生している事になります。

 

次に後半の処理部分のさらに真ん中あたりにブレークポイントを置き、デバッグ実行をしたときに、後半の中の前半の部分でエラーが発生していれば、エラーは全体の後半の中の前半部分で発生することがわかります。

 

これを繰り返すと、エラーが発生している位置は、全体の1/2⇒1/4⇒1/8・・・と追い込む事ができるので、やみくもにエラーの場所を探そうとするよりは、バグの位置を特定しやすくなるかと思います。

 

 

ラインセンサカメラで撮影した画像の縦横比を合わせる場合

ラインセンサカメラ(スキャナのような物)で撮影した画像は、何も考えずに撮影すると、縦横比が合わなく、丸い物を撮影しても楕円状になってしまいます。

これをまじめに縦横比を合わせるには、画像の横方向の撮影分解能(mm/画素)を画像から算出し、カメラのスキャンレート、エンコーダの撮影分解能、被写体の移動速度などから縦横比の合う値を算出します。

 

しかし、計算するのも少し面倒なので、ただ、撮影するだけの場合、私は画像を見ながら被写体の搬送速度を調整して縦横比を調整しています。

(実際の装置の場合、搬送速度は装置の生産能力に係わる部分なので、調整することは仕様的に難しい場合が多く、カメラのスキャンレートやエンコーダの撮影分解能などを調整して縦横比を合わせます。)

 

以下、具体例です。

送り方向に関しては、パルスモータで搬送しているとして、画像が縦長、横長になるパルス速度を見つけます。

(画像が縦長になるのは搬送速度が遅く、画像が横長になるのは搬送速度が速い事を意味します。)

 

ここでは、パルス速度1000Hzで搬送した時に縦長、2000Hzで撮影した時に横長になったとすると、縦横比が合う搬送速度は1000~2000Hzの間にある事がわかります。

ここで、二分法の考え方を応用して、1500Hzで撮影した時に画像が横長になったとすると、縦横比が合う搬送速度は1000~1500Hzの間となります。

 

1250Hzで撮影した時に画像が横長になったとき、1000~1250Hzの間、

1125Hzで撮影した時に画像が縦長になったとき、1125~1250Hzの間、

1187.5Hzで撮影した時に・・・と繰り返していくと、縦横比の合う搬送速度を見つけ出す事ができます。

 

 

 

このように、二分法を知っていると、数式を解くだけでなく、効率的に目的の場所や値を求めることが出来るようになります。

 

使える数学へ戻る

標準偏差計算の注意点

標準偏差の説明としては、

   標準偏差はデータのバラツキを表す

というのが多いでしょうか?

その性質からデータの誤差などの指標としても用いられる事が多くあります。

しかしながら、 標準偏差=バラつき とだけ覚えていると、実際に標準偏差を用いる時に、思わぬ落とし穴にハマる場合があります。

 

そこで、標準偏差の式を改めて眺めてみると

$$\sigma =\sqrt { \frac { 1 }{ n } \sum _{ i=1 }^{ n }{ { \left( { x }_{ i }-\overline { x } \right) }^{ 2 } } } $$

となっています。

この式が何の計算をしているのか?を言葉で表すと

全データの平均(\(\overline { x } \))と各データ(\({ x }_{ i }\))の差の2乗の平均の平方根

を計算しています。

この計算を図示すると下図のようになります。

 

平均とデータの差の二乗の平均をルートした値が標準偏差です。

 

それでは、具体的に標準偏差を計算するとどうなるのか?を見ていきたいと思います。

 

問題1)AとBのデータの標準偏差は、どちらが大きいでしょうか?

答え)B

平均とデータの差が全体的に大きいBのデータの方が標準偏差が大きくなります。

 

問題2)AとBのデータの標準偏差は、どちらが大きいでしょうか?

AとBのデータのばらつきは、ほぼ同じでデータの個数が異なる場合

答え)AとBもほぼ同じ

バラつきの具合が同じ場合、データの個数が変化しても、標準偏差は平均とデータの差の二乗の平均のルートを計算しているので、標準偏差の値はほぼ同じになります。
ただし、データ個数の少ないBの方がデータの個数が少ないため、1つ1つのデータに平均の値が影響されやすくなり、標準偏差の値が各データの値に依存しやすく(標準偏差の値がバラつきやすく)なります。

 

ここまでは、標準偏差=バラつき と思っていても問題にならない場合ですが、次からは注意が必要な場合の例を示します。

 

問題3)AとBのデータの標準偏差は、どちらが大きいでしょうか?

AとBのデータのばらつきは、ほぼ同じでデータの個数が異なります。
ただし、1点だけ大きく平均からズレたデータが含まれます。

答え)B

AとBのデータのバラつき具合は、一見、同じように見えますが、平均とデータの差の二乗の値の合計をデータ個数で割る際に、データ個数の少ないBの方が値としては大きくなるため、標準偏差はBの方が大きくなります。

 

問題4)AとBのデータの標準偏差は、どちらが大きいでしょうか?

Aのデータは前半25個のデータと、後半25個のデータは同じ。
BのデータはAのデータの前半25個のみ

答え)AとBも同じ

一見、Aのデータの方が大きな誤差を含む点が2点(Bのデータは1点)あるため、Aの方が標準偏差の値が大きくなるように見えなくもありませんが、全データの平均と各データの差の2乗の平均の値(分散という)が同じになるため、標準偏差も同じになります。

 

問題5)AとBのデータの標準偏差は、どちらが大きいでしょうか?

Bのデータは、Aのデータに比べ、細かいデータのバラつきは小さいですが、全体的にデータが傾いています。

答え)B

全データの平均と各データの差の2乗の平均の値がBの方が大きくなるため、標準偏差もBの方が大きくなります。

 

問題6)AとBのデータの標準偏差は、どちらが大きいでしょうか?

Bのデータは、Aのデータに比べ、細かいデータのバラつきは小さいですが、全体的にデータが歪んでいます。

答え)B

問題5と同様に、全データの平均と各データの差の2乗の平均の値がBの方が大きくなるため、標準偏差もBの方が大きくなります。

 

まとめ

問題1~6に示したように、人の目で見たデータのバラつきの大きさと、標準偏差の値の大小とは、必ずしも一致しない場合があります。

そのため、 標準偏差=バラつき とだけ認識するのではなく、標準偏差の計算である

標準偏差 = 全データの平均と各データの差の2乗の平均の平方根

として、計算方法を言葉で覚えておく事をおススメします。

 

それでは、人の感じるバラつきと一致するような値を求めるには、どうすれば良いのか?を考えると、決まりきった計算方法は知らないのですが、問題3のようにデータの個数が異なるデータを比較する場合は、標準偏差の計算をデータ全体で計算するのではなく、部分的に標準偏差を計算すると、ある領域ではAとBのバラつきは同じだが、他の領域ではAの方がバラつきが小さい事が分かります。

問題5、問題6のようにデータ全体に傾きや歪みがある場合は、ガウシアンフィルタなどのノイズ除去処理を強めにかけて(下図の紫の線)、そのノイズ除去された値と各データの値の差(下図の赤の点)に関して標準偏差を計算すると、データ全体の傾きや歪みに影響される事なく、バラつきを求める事ができます。

 

 

 

使える数学へ戻る

最小二乗法をExcelで解く

エクセルではグラフの機能で近似曲線の追加というのはありますが、任意の式で近似したい場合や、近似式の係数を取得したい場合にはエクセルの近似曲線の機能だけでは物足りない場合があります。
そんなときにエクセルでまじめに最小二乗法で解く方法を紹介します。

 

まずは、最小二乗法を解くとき、逆行列を用いるところまではご理解願います。
詳しくはこちらで紹介しています。

 

最小二乗法
一般式による最小二乗法

 

また、Excelで行列を計算する方法はこちらを参照下さい。

 

行列の積、逆行列、転置行列の計算

 

以下では円の最小二乗法を例にとって説明をしたいと思います。

 

円の最小二乗法を解くときの円の一般式を

 

X2 + Y2 + AX + BY + C = 0

 

としたときの行列の式は

 

 

です。
まず、行列の式に使われているの値を求めます。
はエクセル関数の「SUM」関数を用いればよいので、簡単に求まります。

 

 

まず、行列の式の各値を各セルに書いておきます。
次に逆行列を求めるのですが、エクセル関数に「MINVERSE」というのがあります。

 

逆行列を表示したい複数のセル(下記の例では3×3のセル)を選択します。
次に

 

=MINVERSE(逆行列を求める複数のセル)

 

と入力します。ここでは文字を入力するだけでEnterキーは押さないこと

 

 

 

で、次が最大のポイント!
通常ではEnterキーを押したいところですが、行列の計算では

 

Ctrlキー  + Shiftキー + Enterキー

 

を同時に押します。
これで、逆行列の値が求まります。

 

 

あとは行列の積が計算できれば、求めたいA,B,Cの値を求めることができるのですが、
行列の積の計算にはエクセル関数の「MMULT」という関数を用います。

 

逆行列を計算したときと同様に行列の積の結果を表示したいセル複数のセルを選択し、

 

=MMULT(行列1,行列2)

 

を入力します。ここでも最後はEnterキーだけではなく

 

Ctrlキー  + Shiftキー + Enterキー

 

を押すので注意して下さい。

 

 

これで求めたかった未知数A,B,Cが求まります。

 

 

今回は円の近似なので、A,B,Cの値を用いで円の中心座標(a,b)、円の半径rを求めるには

 

a = -A / 2
b = -B / 2
r = SQRT(a2 + b2 – C)

 

なので、

 

 

となり、円の近似式を求めることができます。

 

このエクセルによる円の最小二乗法のファイルはこちら↓

LSM_Circle.xls

 

また、一般的なn次式による最小二乗法を解いたのファイルはこちら↓

n_LeastSquare.xls

 

使える数学へ戻る

 

一般式による最小二乗法(円の最小二乗法)

n次曲線による最小二乗法につていは説明しましたが、円や楕円、その他の一般式についても最小二乗法による近似は可能です。
今回は円の最小二乗法を例にとって説明しますが、他の一般式についても要領は同じです。

 

まず初めに近似する一般式を作成します。

円の場合、円の中心座標(a,b)、円の半径 r とすると

 

(X – a)2 + (Y – b)2 = r2 ・・・①

 

となります。

n次曲線による最小二乗法ではX座標がXiの時のY座標Yiと近似曲線上のY座標f(Xi)との差の二乗が最小になるように計算していましたが、円の場合、X座標がXiの時にY座標が存在しなかったり、Y座標が2点存在してしまう場合もあります。

 

そこで、今回はY座標の差の2乗の総和を求めるのではなく、①式を=0の式に変形し、
点の座標(Xi,Yi)を代入し、その式の2乗の総和を求めます。
(n次式近似でやった計算も結局は近似式=0になるように変形していた訳ですが...)

 

①式を=0になるように変形し、2乗の総和は

 

∑{(Xi – a)2 + (Yi – b)2 – r22 = 0 ・・・②

 

となります。

ただし、このまま未知数a,b,rに関して②式を偏微分してもa,b,rは4次関数となるので、各未知数の偏微分=0の時が最小とは限りません。

 

そこで、②式を展開し、下の式のように置きます。

 

∑{Xi2 + Yi2 + AXi + BYi + C}2 = 0 ・・・③

 

ただし

A = -2a ・・・④
B = -2b ・・・⑤
C =  a2 + b2 – r2 ・・・⑥

 

③式に関してA,B,Cについて偏微分すると

 

∂/∂A = A∑Xi2 + B∑XiYi + C∑Xi + ∑Xi3 + ∑XiYi2 = 0 ・・・⑦
∂/∂B = A∑XiYi + B∑Yi2 + C∑Yi + ∑Xi2Yi + ∑Yi3 = 0 ・・・⑧
∂/∂C = A∑Xi + B∑Yi + C∑1  + ∑Xi2 + ∑Yi2 = 0 ・・・⑨

 

⑦、⑧、⑨式を行列を用いて解くと

 

 

となり、逆行列を求めるとA,B,Cが求まります。

さらに④、⑤、⑥式より円の中心座標(a,b)、円の半径 rが求まります。

  \(a=-\frac{A}{2}\)

  \(b=-\frac{B}{2}\)

  \(r=\sqrt{a^{2}+b^{2}-C}\)

 

使える数学へ戻る

 

関連記事

最小二乗法の計算方法

最小二乗法の計算方法

最小自乗法という人もいますが、私は最小二乗法派です。

最小二乗法とは?

点の集まりから、近似直線(曲線)などの近似式を解くための手法です。
単に直線や二次曲線、三次曲線などに限らず、最小二乗法の応用範囲は広くあります。

最小二乗法は、言葉の通りに2乗を最小にする方法なのですが、この2乗と言うのが、方程式上の理論値と近似するデータとの差の2乗の合計を最小にする方法となります。と言っても分かりづらいので、具体的な最小二乗法の例を示します。

 

最小二乗法の計算方法

今回は説明を簡単にするため、f(X)= aX + b の直線でデータを近似する場合で説明します。

 

 

X = Xi のとき、直線上のY座標は

$$Y_{i}=f(X_{i})=aX_{i}+b$$

より、点と近似直線との差は

$$Y_{i}-f(X_{i})=Y_{i}-(aX_{i}+b)$$

となります。

 

この差の合計が一番小さい時にもっとも点の集まりを直線で近似できていることとなりますが、
そのままこの差の総和を計算すると差の値がプラスの時とマイナスの時があるため、例えば下図

 

 

のように、正の差と負の差が同じようにある場合、差の合計が小さくなってしまうため、差を2乗してから合計した値が最小になるようにします。

 

よって、差の2乗の合計は

$$\sum(Y_{i}-(aX_{i}+b))^{2}\\
=\sum(a^{2}X_i^2+2aX_{i}b+b^{2}-2aX_{i}Y_{i}-2bY_{i}+Y_i^2 )$$

・・・式①

となり、この総和の値が最小となるaとbを求めれば近似直線を求められることとなります。

 

この式を良くみると、aとbの二次関数でできているため、この二次関数が最小となる値を求めればよい。
二次関数が最小となるのは微分した値 = 0 の場合であるから、式①をaとbで偏微分した値が
0(ゼロ)となる値を求めればよい。

 

式①をaとbでそれぞれ偏微分すると

$$\frac{\partial}{\partial a}=2a\sum X_i^2+2b\sum X_{i}-2\sum X_{i}Y_{i}=0\\
\frac{\partial}{\partial b}=2a\sum X_i+2b\sum 1-2\sum Y_{i}=0$$

より

$$a\sum X_i^2+b\sum X_{i}-\sum X_{i}Y_{i}=0\\
a\sum X_i+b\sum 1-\sum Y_{i}=0$$

となり、未知数がaとbで式が2本あることから、aとbの連立法的式を解くことで未知数a,bを求めることができます。

 

連立方程式を行列で表すと

$$\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\frac{1}{\sum X_i^2\sum 1-\sum X_{i}\sum X_{i}}\left(\begin{array}{c}\sum 1 & -\sum x_i\\ -\sum x_i & \sum X_i^2\end{array}\right)\left(\begin{array}{c}\sum X_iY_i\\ \sum Y_i\end{array}\right)$$

となり、行列を解くと未知数a,bが求まり、近似直線を求めることができます。

 

ちなみに、二次式(aX2 + bX + c)で近似すると

$$\left(\begin{array}{c}a\\ b \\ c\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

三次式(aX3 + bX2 + cX + d)で近似すると

$$\left(\begin{array}{c}
a\\
b\\
c\\
d
\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^6 & \sum X_i^5 & \sum X_i^4 & \sum X_i^3\\
\sum X_i^5 & \sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^4 & \sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^3 & \sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^3Y_i\\
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

と規則的に変化するため、n次式近似のプログラムも簡単に作ることができます。
ただし、行列の部分を二次元配列と考えると、近似する次数(n)が上がって行くと、配列の先頭に値が追加されていくのは、扱いにくいので、n次式近似のプログラムを作る場合は、近似式を

$$Y=a+bX+cX^2+dX^3 \cdot\cdot\cdot\cdot$$

として、行列の表現は

$$\left(\begin{array}{c}a\\ b\\ c\\d\\:
\end{array}\right)=
\left(\begin{array}{c}
\sum 1 & \sum X_i & \sum X_i^2 & \sum X_i^3 & .. \\
\sum X_i & \sum X_i^2 & \sum X_i^3 & \sum X_i^4 & .. \\
\sum X_i^2 & \sum X_i^3 & \sum X_i^4 & \sum X_i^5 & ..\\
\sum X_i^3 & \sum X_i^4 & \sum X_i^5 & \sum X_i^6 & ..\\
: & : & : & : & :
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum Y_i\\
\sum X_iY_i \\
\sum X_i^2Y_i\\
\sum X_i^3Y_i\\
:
\end{array}\right)$$

として、n次式で近似した最小二乗法を求めます。

 

この計算をエクセルで行ったものはExcelによる最小二乗近似からダウンロードできます。
逆行列を求めるプログラムにはガウスの消去法などが用いられます。

 

最小二乗法での注意点

三次式(a+bX+cX2+dX3)の近似を例にとってみると、近似処理をする過程でXに関する計算式はX~X6まで計算されます。

 

ここで、例えばXの値の範囲が

 

0.01~1000

 

だったとすると、Xに関する計算では

 

0.000000000001~1000000000000000000

 

までの計算がされます。
逆行列の計算ではこれらの値に関して足し算や掛け算がなされる訳ですからコンピュータの
計算誤差(桁落ち、情報落ちなど?)が発生する場合があります。
これを防ぐために、最小二乗法のプログラムでは近似式をそのまま使うのではなく、Xの値から
Xの平均値(X)を引いた計算式

$$Y=a+b(X-X_b)+c(X-X_b)^2+d(X-X_b)^3$$

で処理を行うのが一般的です。

 

最小二乗法の関連記事

一般式による最小二乗法(円の最小二乗法)

最小二乗法をExcelで解く

最小二乗法の最適化(高速化)

ロバスト推定法(Tukey’s biweight)

 

使える数学へ戻る

 

任意点周りの回転移動(アフィン変換)

回転行列では原点周りに点を回転させますが、任意の点(C、C周りに回転させたい場合にはどうするのか?

 

 

これまでの知識を少し応用することで、意外と簡単に求めることができます。

 

まず、回転する点を回転中心座標が原点と一致するように点を移動させます。

 

次に移動した点を原点周りに回転移動させます。

 

 

回転移動後、点を原点から元の回転中心位置へ移動させます。

 

 

これで、任意点(C、C周りに点を回転移動させることができました。

 

この 原点へ移動→原点周りの回転→元へ戻す の一連の処理を行列であらわすと

 

 

行列部分を整理すると

 

 

となり、任意点(C、C周りに点を回転移動させる行列を求めることができます。

 

今回は回転移動について説明しましたが、拡大縮小についても、任意点を基点に拡大縮小する場合についても同様です。

この考え方は、二次元の座標だけでなく、三次元の場合も使えます。。

参考

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

 

使える数学へ戻る

 

回転行列、拡大縮小行列、平行移動行列(三次元座標の場合)

二次元座標(X,Y座標)の場合のアフィン変換行列についてはこちらで説明しましたが、今回は三次元座標(X,Y,Z座標)のアフィン変換となります。

三次元座標の場合、まず座標軸の定義、回転方向の定義を明確に覚えます。

この座標は右手座標系と呼ばれます。
フレミングの法則のときのように右手親指人差し指中指をそれぞれ
直交するようにします。
このとき親指から順に親指がX軸人差し指がY軸中指がZ軸の方向と
なります。
回転方向は電流と磁界の向きと同じように電流軸の向き磁界回転方向
に相当します。(右ねじの法則と同じです。)

 

回転行列

三次元の回転行列の前に二次元の回転行列のおさらいです。
二次元の回転行列は以下の通りとなります。

 

 

これをベースに三次元座標の場合では、回転する軸の正の方向から原点の方向を見たときに、X軸、Y軸はそれぞれ何軸に相当するのか?を考えれば、二次元座標のXやYの変数の置き換えで導き出すことができます。
行列変換しない軸に関しては単位行列でそのまま残します。

 

【X軸周りの回転】

【Y軸周りの回転】

【Z軸周りの回転】

拡大縮小行列

点(x, y, z)を原点に関してX軸方向に、Y軸方向に、Z軸方向にZする行列は

 

平行移動行列

点(x, y, z)をX軸方向に、Y軸方向に、Z軸方向にZだけ移動する行列は

 

補足

三次元の座標変換に関して検索すると座標変換は下記のように

行ベクトルで表記される場合もあるのですが、変換行列の値が変わるので、
混同しないようご注意下さい。
この表現はマイクロソフトがお得意で、DirectX(Direct3D)や.NETのアフィン変換でしか使われないので、特に必要の無い場合は覚えない方が無難です。

 

関連記事

二次元座標のアフィン変換についてはこちら↓にまとめています。

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

 

使える数学へ戻る

 

回転行列、拡大縮小行列、平行移動行列

本記事は、私がアフィン変換を勉強し始めた当初の記事になります。
今では、3×3行列の同次座標行列と呼ばれる行列しか用いておらず、こちらの方が断然おススメなので、下記ページを参照ください。

アフィン変換(平行移動、拡大縮小、回転、スキュー行列)

以下は、2×2行列を使ったアフィン変換の説明です。

回転行列

点(x, y)を原点まわりに反時計方向にθ度回転する行列は

 

拡大縮小行列

点(x, y)を原点に関してX軸方向に、Y軸方向にする行列は

 

平行移動行列

点(x, y)をX軸方向に、Y軸方向にだけ移動する行列は

 

 

ただし、平行移動だけ行列の足し算になると、扱いにくい場合があるので3×3行列を用いて以下のように表す場合もあります。
というより、こちらを使う方が便利です。(私はこちらしか使いません。)

【回転行列】

【拡大縮小行列】

【平行移動行列】

 

とすることで、すべての座標変換を行列の積で扱うことができます。

この行列を同次座標行列と言います。
詳細はこちらを参照ください。

 

参考まで...

個人的には回転行列を覚えるのは苦手で、SinとCosが逆になっりマイナスのつける位置を間違ったりしていたのですが、次のように考えることで少しは覚えやすくなりました。

下図のように

点(1,0)をθ度回転すると(Cosθ、Sinθ)
点(0,1)をθ度回転すると(-Sinθ、Cosθ)

に移動することはすぐにわかります。

このことを行列で表現すると
点(1,0)が(Cosθ、Sinθ)になることから

点(0,1)が(-Sinθ、Cosθ)になることから

という事がわかります。
これを合わせて表現すると

となり、回転行列が求まります。
この計算を何回か繰り返すと、そのうち覚えると思います。

 

使える数学へ戻る

 

行列で連立方程式を解く

行列を使って連立方程式を解く方法を紹介します。

 

【計算例】

3つの方程式

 

 

を行列であらわすと

 

 

となります。この行列を逆行列を使ってX、Y、Zに関して解くと

 

 

となり連立方程式を行列で解くことができます。

 

 

…というのは普通すぎて面白くないので、3点からなる円の方程式を行列で解く方法を紹介します。

 

中心(a、b)、半径 r の円の方程式

 

(X – a)2 + (Y – b)2 = r2

 

を展開して

 

X2 + Y2 – 2aX – 2bY + a2 + b2 – r2 =  0

 

となり、A = 2a、B = 2b、C = r2 – a2 – b2 とおくと、上式は

 

AX + BY + C  =   X2 + Y2

 

となります。

 

ここで、円上の3点(X0、Y0)、(X1、Y1)、(X2、Y2を代入すると

 

AX0 + BY0 + C  =  X02 + Y02
AX1 + BY1 + C  =  X12 + Y12
AX2 + BY2 + C  =  X22 + Y22

 

の3本の式が成り立ち、これを行列で表現すると、

 

 

となり、A、B、Cに関して行列を解くと

 

 

となり、A、B、Cが求まることから、A = 2a、B = 2b、C = r2 – a2 – b2 より
中心半径を求めることができます。

 

行列を解く部分は

ガウスの消去法

のページを参照下さい。

 

使える数学へ戻る

 

ガウスの消去法

ガウス(Gauss)の消去法は連立一次方程式を解くのに用いられます。

 

基本的な方針は、下記のような連立一次方程式

 

 

を行列であらわすと、

 

 

となりますが、対角成分を全て、左下の成分をになるように、行の入替えや足し算、引き算などを行い、下記の行列になるように調整します。

 


(?の部分は任意の値で可)

 

これまでの処理を前進消去といいます。
ここで、最後の行の部分(3の部分)に着目すると、答えが確定しています。
この値を使って最後から2行目の値を計算すると答えが算出できます。
このように最後の行から順々に計算すると答えが全て計算することができます。
この処理は後退代入といいます。
この前進消去と後退代入の処理を合わせてガウスの消去法といいます。

 

さらに、前進消去のときに対角成分を1にするときの割り算の計算のときにおいて、
ピボット選択を行ったGauss-Jordan法でも説明しましたが、

 

  • 0(ゼロ)で割ることはできない。
  • 値を絶対値の小さい値で割ると、値に誤差が含まれる場合、計算結果に大きく影響が出てしまう。

     

    という性質を考慮します。

     

    ということで、以下、具体的な計算例を示します。

     

    連立一次方程式


     

    のうち、まず最初に1の項のどれか1つを1にするのですが、割り算の計算で誤差が少なくなるように、1の係数が一番大きい2行目の式を一番上に持ってきます。

     

     

    次に1行目の式の1の係数がになるように1行目の式全体を1の係数()で割ります。

     

     

    次に1行目以下の式の1の項が消えるように、
    2行目-1行目×23行目-1行目を計算します。

     

     

    ここで、3行目の式の3の項(行列で表現すると対角成分)が消えてしまったので、
    このままだと対角成分をにできないので、2行目3行目を入れ替えます。

     

     

    次に2行目の式の2の係数をにするように、2行目の式全体に3/2を掛けます。

     

     

    次に2行目以下の2の項が消えるように3行目+2行目 × 2/3を計算します。

     

     

    3行目の式のの係数はすでになので、前進消去はこれで終了です。
    この3本の式を行列であらわすと

     

     

    というように、対角成分が全てがで左下の成分がになったことが分かります。

     

    次に後退代入です。

     

    一番下の行の式はすでに答えが確定(X3 = 3)しているのがわかります。
    一番下の答えを用いて、下から2行目の式の答えを計算、
    一番下、下から2行目の答えを用いて下から3行目の式の答えを計算していきます。
    (今回の例だけ特別に、下から2行目の式の答えが確定してしまっています。)

     

     

    こうして全ての答えが求まります。

     

     

    この一連の処理がガウスの消去法です。

     

    このガウスの消去法を用いると、連立方程式を求めたり、最小二乗法の未知数を
    求めることができます。

     

    使える数学へ戻る

     

    逆行列(Gauss-Jordan法)

    2×2行列の逆行列

    行列

    の逆行列は

    となります。
    ただし、ad-bc = 0 のとき、逆行列は存在しません。

     

    3×3以上の行列の逆行列

    逆行列を解く手法はいくつかありますが、ここでは比較的分かり易いGauss-Jordan法を紹介します。

     

    Gauss-Jordan法では行列の右側に単位行列を付けたして、行ごとに掛け算、足し算、引き算を行い、行列の左側が単位行列になるように計算を繰り返し、最後に右側に残った行列が逆行列となります。

     

    といっても分かりづらいと思うので、具体的な計算例は以下の通りです。

     

    行列

    の右側に単位行列を追加します。

    1行1列目の要素が1となるように1行目2で割ります

    1列目の要素が(1 0 0)となるように

    [2行目] = [2行目]ー[1行目]
    [3行目] = [3行目]ー[1行目]×4

    を計算します。

    2行2列目の要素が1となるように2行目2倍します。

    2列目の要素が(0 1 0)となるように

    [1行目] = [1行目]ー[2行目]×3/2
    [3行目] =
    [3行目]+[2行目]

    を計算します。

    ここで、3行3列目の要素はすでに1なので、3列目の要素が(0 0 1)となるように

    [1行目] = [1行目]+[3行目]×2
    [2行目] = [2行目]ー[3行目]×2

    を計算します。

    これで、左側が単位行列となり、右側に出来た行列が求める逆行列となります。

     

    ただ、このままの方法では、求める行列の対角要素(行数と列数の同じ場所)に0(ゼロ)がある場合は対角要素を1に出来ない(0で割れない)ので、ここにピボット選択という手法を導入します。
    このピボット選択についてはピボット選択を行ったGauss-Jordan法にて紹介しています。

     

    使える数学へ戻る