ラグランジュの未定乗数法~等式制約条件

準備

長方形の面積最大化

たとえば2変数の問題として、長方形の周囲長Lを一定として、その面積が最大となる長方形の形状と面積はどのようになるかを考える。この場合、長方形の辺の長さをx, yとすると、問題は以下のように表せる。

(1)    \begin{equation*} \mathrm{maximize} \quad S(x, y) = xy \quad \mathrm{subject~to} \quad x+y = \frac{L}{2} \end{equation*}

これは以下のように代数的に簡単に解けて、答えは正方形とわかる。

(2)    \begin{gather*} S = x \left( \frac{L}{2} - x \right) = -x^2 + \frac{L}{2} x = - \left( x - \frac{L}{4} \right)^2 + \frac{L^2}{16} \\ \max S = \frac{L^2}{16} \quad \mathrm{for} \: x = y = \frac{L}{4} \end{gather*}

ただし変数の数が増えたり、目的関数や制約条件が複雑になると、解析的に解くのが面倒になる。

Lgrangeの未定乗数法による解

解法から先に示す。Lagrangeの未定乗数法では、目的関数L(x, y)に対して以下の問題となる。

(3)    \begin{gather*} \mathrm{maximize} \quad L(x, y, \lambda) = S(x, y) - \lambda g(x, y) = xy - \lambda \left( x + y - \frac{L}{2} \right) \\ \mathrm{subject~to} \quad S(x, y) = xy, \quad g(x, y) = x + y - \frac{L}{2} \end{gather*}

L(x, y, \lambda)を最大化するために、x, y, \lambdaで偏微分した以下の方程式を設定する。

(4)    \begin{gather*} \frac{\partial L}{\partial x} = 0, \quad \frac{\partial L}{\partial y} = 0, \quad \frac{\partial L}{\partial \lambda} = 0 \end{gather*}

これを計算すると

(5)    \begin{gather*} y - \lambda = 0, \quad x - \lambda = 0, \quad x + y - \frac{L}{2}= 0 \\ \therefore \lambda = x = y = \frac{L}{4} \end{gather*}

Lagrangeの未定乗数法の一般形

一般には、変数\boldsymbol{x} = (x_1, \ldots, x_n)について、目的関数f(\boldsymbol{x})を制約条件g(\boldsymbol{x})=0の下で最大化/最小化する問題として与えられる。

(6)    \begin{align*} & \mathrm{maximize/minimize} \quad f(\boldsymbol{x}) \\ & \mathrm{subject~to} \quad g(\boldsymbol{x}) = 0 \end{align*}

この等式制約条件付き最大化/最小化問題は、以下のようにL(\boldsymbol{x}, \lambda)を導入して、連立方程式として表現される。

(7)    \begin{align*} & L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda g(\boldsymbol{x}) \\ & \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_1} = \cdots = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_n} = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial \lambda} = 0 \end{align*}

例題

例題1:平面と円

平面x + y - 1について、制約条件f(x, y) = x^2 + y^2の下での極値を求める。

(8)    \begin{align*} & \mathrm{minimize} \quad f(x, y) = x + y - 1 \\ & \mathrm{subject~to} \quad x^2 + y^2 = 1 \end{align*}

lagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(9)    \begin{equation*} L(x, y, \lambda) = x + y - 1 - \lambda (x^2 + y^2 - 1) \end{equation*}

(10)    \begin{align*} &\dfrac{\partial L}{\partial x} = 1 - 2 \lambda x = 0 \\ &\dfrac{\partial L}{\partial y} = 1 - 2 \lambda y = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x^2 - y^2 + 1 = 0 \end{align*}

この連立方程式を解くと以下のようになり、解として2つの極値を得るが、それらは最大値と最小値に相当する。

(11)    \begin{gather*} x = y = \frac{1}{2 \lambda} \quad \Rightarrow \quad \frac{1}{4 \lambda ^2} + \frac{1}{4 \lambda ^2} = 1 \quad \Rightarrow \quad \lambda = \pm \frac{1}{\sqrt{2}}\\ \therefore \; x = y = \pm \frac{1}{\sqrt{2}} \approx \pm 0.7071\\ \max f(x, y) = \sqrt{2} - 1 \approx 0.414\\ \min f(x, y) = -\sqrt{2} - 1 \approx -2.414 \end{gather*}

これを目的関数のコンターと制約条件の線で表すと以下の通り。

 

例題2:凸関数と直線

下に凸な関数f(x, y) = x^2 + y^2について、直線x + y = 1の制約条件下での最小値を求める。

(12)    \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 1 = 0 \end{align*}

ここでlagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(13)    \begin{equation*} L(x, y, \lambda) = x^2 + y^2 - \lambda (x + y - 1) \end{equation*}

(14)    \begin{align*} &\dfrac{\partial L}{\partial x} = 2x - \lambda = 0 \\ &\dfrac{\partial L}{\partial y} = 2y - \lambda = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x - y + 1 = 0 \end{align*}

この連立方程式を解くと以下のようになり、解は最小値1つとなる。

(15)    \begin{gather*} x = y = \frac{\lambda}{2} \quad \Rightarrow \quad \frac{\lambda}{2} + \frac{\lambda}{2} = 1 \quad \Rightarrow \quad \lambda = 1\\ \therefore \; x = y = \frac{1}{2} \quad , \quad \min f(x, y) = \frac{1}{2} \end{gather*}

これを目的関数のコンターと制約条件の線で表すと以下の通り。

なお、これを3次元で表示すると以下のようになる。青い曲面が目的関数で、赤い直線が制約条件となる。最適化問題は、制約条件を満たす曲面上の点(図中、赤い放物線)の最小値を求めることになる。

幾何学的説明

式(7)は以下のように書ける。

(16)    \begin{gather*} \left[ \begin{array}{c} \dfrac{\partial f}{\partial x_1} \\ \vdots \\ \dfrac{\partial f}{\partial x_n} \end{array} \right] = \lambda \left[ \begin{array}{c} \dfrac{\partial g}{\partial x_1} \\ \vdots \\ \dfrac{\partial g}{\partial x_n} \end{array} \right] \\ g(x_1, \ldots, x_n) = 0 \end{gather*}

さらにgradientで表すと

(17)    \begin{gather*} \nabla f = \lambda \nabla g \\ g(x_1, \ldots, x_n) = 0 \end{gather*}

すなわちこの式の解(x_1, \ldots, x_n)は、制約条件であるg(x_1, \ldots, x_n)=0を満足し、その曲線上にある。さらに解の点においてf(x_1, \ldots, x_n)の勾配ベクトルとゼロ平面上におけるg(x_1, \ldots, x_n)の勾配ベクトルが平行になる。これはゼロ平面上の解の点において制約条件の曲線とf(x_1, \ldots, x_n)のコンターのうち特定の曲線が接するのと同義であり、この点は停留点(stationary point)である。

つまりこのような停留点を発見する手順は、制約条件を満たす(制約条件の線上にある)点のうち、その点において目的関数のgradientと制約条件の関数のgradientが平行となる点を求めるということになる。再度これを式で表すと、

(18)    \begin{gather*} \nabla f(\boldsymbol{x}) = \lambda \nabla g(\boldsymbol{x}) \\ g(\boldsymbol{x}) = 0 \end{gather*}

となるが、これをLangrange関数L=f - \lambda gと定義したうえで各変数で偏微分したものをゼロと置いた方程式を解くと表現している。未定乗数λは停留点における目的関数のgradientと制約条件の関数のgradientの比を表している。

λの符号

λの符号に意味があるかどうか。

たとえば、以下の制約条件付き最適化問題を考える。

(19)    \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 2 = 0 \end{align*}

(20)    \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(x + y - 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x - \lambda = 2y - \lambda = -x - y + 2 = 0\\ \lambda = 2, \; x=y=1 \end{gather*}

この問題の制約条件を変更して以下のようにした場合。

(21)    \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad - x - y + 2 = 0 \end{align*}

(22)    \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(- x - y + 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x + \lambda = 2y + \lambda = x + y - 2 = 0\\ \lambda = -2, \; x=y=1 \end{gather*}

このように、制約条件の正負を反転するとλの符号が逆になるが解は変わらない。これを表したのが以下の図。

等式制約条件の場合、制約条件の線上でgradientが平行になりさえすればよいので、λ符号(制約条件式の正負)には拘らなくてよい。

停留点が極致とならない例

以下の最適化問題の解をみてみる。

(23)    \begin{align*} & \mathrm{maxmize} \quad f(x, y) = x^3 + y^3 \\ & \mathrm{subject~to} \quad g(x, y) = x - y \end{align*}

(24)    \begin{align*} &L(x, y, \lambda) = x^3 + y^3 - \lambda(x - y) \\ &\left\{ \begin{array}{c} \dfrac{\partial L}{\partial x} = 3x^2 - \lambda = 0\\ \\ \dfrac{\partial L}{\partial y} = 3y^2 + \lambda= 0\\ \\ \dfrac{\partial L}{\partial \lambda} = x - y = 0 \end{array} \right. \end{align*}

第1式から第2式を引き、第3式を適用して、

(25)    \begin{gather*} \3x^2 - 3y^2 - 2\lambda = 0 \\ 3(x + y)(x - y) = 2\lambda = 0 \\ \therefore x = y = \lambda = 0 \end{gather*}

ここでz=f(x, y)の曲面と制約条件g(x, y)=0に対する曲面上の軌跡を描くと下図左のようになる。下図の右はt = x = yとして曲線を表したもので、x = y = 0で勾配は水平になっており、停留点ではあるが極大/極小となっていない。

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

ヘッセの標準形~点と平面の距離

2次元の場合

直線l、点Qの距離を考える。直線と各点の記号、座標を以下のように定義する。

(1)    \begin{equation*} l:ax + by + c = 0 \end{equation*}

Step-1:直線に直交するベクトル

まず、ベクトル(a, b)が直線lに直交することを示す。直線は以下のように媒介変数表示できて、ベクトル(u_x , u_y)は直線に平行なベクトル。

(2)    \begin{equation*} \left\{ \begin{array}{lll} x &=& u_x t + x_0 \\ y &=& u_y  t+ y_0 \end{array} \right. \end{equation*}

これを直線の式に代入して、

(3)    \begin{gather*} a (u_x t + x_0) + b (u_y t + y_0) + c = 0 \\ (a u_x + b u_y) t + (a x_0 + b y_0 + c) = 0 \end{gather*}

ここで任意のtに対して上式が成り立つことから、a u_x + b u_y = 0となり、ベクトル(a, b)は直線に垂直であることが示された。

別解:点と直線の距離をパラメーター(媒介変数)によって愚直に求める方法

Step-2:法線ベクトルとの平行条件による導出

与えられた点Pから直線lへの垂線の足をHとすると、\overrightarrow{PH} \parallel (a, b)なので、以下が成り立つ

(4)    \begin{equation*} \left| \overrightarrow{PH} \cdot (a, b) \right| = \left| \overrightarrow{PH} \right| \cdot | (a, b) | \end{equation*}

ここで(x_h, y_h)は直線l上にあることを考慮し、上式の左辺を以下のように変形できる。

(5)    \begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot (a, b) \right| &=& | (x_p - x_h ) a + (y_p - y_h) b | \\ &=& | a x_p + b y_p  - (a x_h + b y_h) | \\ &=& | a x_p + b y_p + c | \end{equation*}

これより

(6)    \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| a x_p + b y_p + c \right|}{\sqrt{a^2 + b^2}} \end{equation*}

3次元の場合

3次元平面の式

3次元空間内の平面は、たとえば以下のように表すことができる。

(7)    \begin{equation*} \pi : w_x x + w_y y + w_z z + w_0 = 0 \end{equation*}

一方、3次元平面上の点と法線ベクトル{\boldsymbol n} = (x_x, n_y, n_z)との直行条件から、以下のようにも表現できる。

(8)    \begin{gather*} {\boldsymbol n} ({\boldsymbol x} - {\boldsymbol x_0} ) = 0 \\ (n_x, n_y, n_z) \cdot (x - x_0, y - y_0, z - z_0) = 0 \\ n_x x + n_y y + n_z z + (-n_x x_0 -n_y y_0 -n_z z_0) = 0 \end{gather*}

上記2つの式より、ベクトル{\boldsymbol w} = (w_x, w_y, w_z)は平面に対する法線ベクトルであることがわかる。

法線ベクトルとの平行条件

この法線ベクトルがベクトル\left| \overrightarrow{PH} \right|と平行であることから、

(9)    \begin{equation*} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| = \left| \overrightarrow{PH} \right| \cdot \left| {\boldsymbol w} \right| \end{equation*}

上式の左辺は以下のように変形できる。

(10)    \begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| &=& \left| (x_p - x_h, y_p - y_h, z_p - z_h) \cdot (w_x, w_y, w_z) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p - (w_x x_h +w_y y_h + w_z z_h) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p + w_0 \right| \end{array} \end{equation*}

以上のことから、点{\boldsymbol x}_pから三次元平面\piへの距離については、以下で表される。

(11)    \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| w_x x_p + w_y y_p + w_z z_p + w_0 \right|}{| {\boldsymbol w} |} \end{equation*}

多次元の場合

n次元の超平面を以下の式で与える。

(12)   } \begin{equation*} {\boldsymbol w} \cdot {\boldsymbol x} + w_0 = 0 \; \Leftrightarrow \; w_0 + w_1 x_1 + \cdots + w_n x_n = 0 \end{equation*}

このとき、これまでと同様の考え方により、点{\boldsymbol x}_p (x_{p1}, \ldots , x_{pn})と上記の超平面との距離は以下で表される。

(13)    \begin{equation*} \left| \overrightarrow{PH} \right| = \frac{ {\boldsymbol w} \cdot {\boldsymbol x} + w_0 }{ \| {\boldsymbol w} \|} \end{equation*}

 

 

numpy – 行列(ndarray)

ベクトルと行列の定義

リテラル

ベクトルはnp.array()で引数にリストを指定して定義。

行列は同じくnp.array()で引数に二次元配列のリストを指定して定義。

単位行列

numpy.identity(n)でn×nの単位行列を生成。

転置

行列の転置にはtranspose()メソッドを使う。代替として.Tとしてもよい。

一次元配列で定義したベクトルにはtranspose()は効かない。列ベクトルに変換するにはreshape()メソッドを使う(reshape(行数, 列数))。

演算

定数倍

ベクトル・行列の定数倍は、各要素の定数倍。

加減

同じ要素数のベクトル、同じ次元の行列同士の下限は要素同士の加減

ベクトルの内積

同じ要素数のベクトルの内積(ドット積)はnp.dot()で計算。

     \begin{equation*} {\bf a} \cdot {\bf b} = \sum_{i=1}^n a_i b_i \end{equation*}

*演算子を使うと、要素ごとの積になる。

行列の積

行列同士の積もnp.dot()で計算。l×m行列とm×n行列の積はl×n行列になる。

     \begin{equation*} \left( \begin{array}{ccc} a_{11} & \cdots  & a_{1m} \\ \vdots & a_{ij} & \vdots \\ a_{l1} & \cdots & a_{lm} \\ \end{array} \right) \cdot \left( \begin{array}{ccc} b_{11} & \cdots & b_{1n} \\ \vdots & b_{jk} & \vdots \\ b_{m1} & \cdots & b_{mn} \\ \end{array}  \right) = \left( \sum_{j=1}^m a_{ij} b_{jk} \right) \end{equation*}

次元数が整合しないとエラーになる。

行ベクトルと行列の積は、ベクトルを前からかけてok。

     \begin{equation*} (1, 2, 3) \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array} \right) = (30, 36, 42) \end{equation*}

行列と列ベクトルの積は、一次元配列のベクトルをreshape()で列ベクトルに変換してから。

     \begin{equation*} \left( \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{array} \right) \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right) = \left( \begin{array}{c} 14 \\ 32 \\ 50 \end{array} \right) \end{equation*}

なお、np.dot()の代わりに演算子@が使える。ベクトル同士なら内積、少なくともいずれか一つが行列なら行列積。

numpy.linalgパッケージ

行列式

行列式はnumpy.linalgパッケージのdet()関数で得られる。linalgは”linear algebra”の略で、慣例としてLAという名前で代替される。

逆行列

逆行列はnumpy.linalgパッケージのinv()関数で得られる。

固有値・固有ベクトル

正方行列の固有値と固有ベクトルを、eig()関数で得ることができる(行列の固有値・固有ベクトルの例題を用いた)。

結果は固有値が並んだベクトルと固有ベクトルが並んだ配列で、それぞれを取り出して利用する。なお、固有ベクトルはノルムが1となるように正規化されている。

注意が必要なのは固有ベクトルの方で、各固有ベクトルは配列の列ベクトルとして並んでいる。固有ベクトルを取り出す方法は2通り。

固有値に対応するサフィックスで列ベクトルを取り出す。この方法はnumpyの公式ドキュメントにも以下のように書かれている。

v(…, M, M) array
The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].

固有ベクトルの配列を転置して、行ベクトルの並びにする。

 

行列の固有値・固有ベクトル

概要

行列{\rm A}の固有値・固有ベクトルは以下で定義される。

(1)    \begin{equation*} \boldsymbol{Ax} = \lambda \boldsymbol{x} \end{equation*}

これを以下のように変形する。

(2)    \begin{equation*} (\boldsymbol{A} - \lambda \boldsymbol{I} ) \boldsymbol{x} = {\bf 0} \end{equation*}

この方程式が解をもつためには、以下の条件が必要。

(3)    \begin{equation*} | \boldsymbol{A} - \lambda \boldsymbol{I} | = 0 \end{equation*}

例題

以下の行列に対する固有値、固有ベクトルを求める。

(4)    \begin{equation*} \boldsymbol{A} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \end{equation*}

この行列に対する固有値方程式は以下の通り。

(5)    \begin{equation*} | \boldsymbol{A} - \lambda \boldsymbol{I} | = \left| \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) - \lambda \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) \right| = \left| \begin{array}{cc} 3 - \lambda & 1 \\ 2 & 4 - \lambda \end{array} \right| = 0 \end{equation*}

これを解くと、

(6)    \begin{align*} & (3 - \lambda) (4 - \lambda) - 2 = 0 \\ & \lambda ^2 - 7 \lambda + 10 = 0 \\ & (\lambda - 2)(\lambda - 5) = 0 \\ & \lambda = 2, \; 5 \end{align*}

次に、各固有値に対する固有ベクトルを求める。

まず\lambda = 2に対しては、

(7)    \begin{gather*} \left( \begin{array}{cc} 3 - 2 & 1 \\ 2 & 4 - 2 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{cc} 1 & 1 \\ 2 & 2 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \end{array} \right) \\ \Rightarrow \; y = -x \\ \therefore \; \boldsymbol{x} = (t, -t) \end{gather*}

確認してみると、

(8)    \begin{gather*} \boldsymbol{Ax} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} t \\ -t \end{array} \right) = \left( \begin{array}{c} 2t \\ -2t \end{array} \right) , \quad \lambda \boldsymbol{x} = 2 \left( \begin{array}{c} t \\ -t \end{array} \right) = \left( \begin{array}{c} 2t \\ -2t \end{array} \right) \end{gather*}

また\lambda = 5に対しては、

(9)    \begin{gather*} \left( \begin{array}{cc} 3 - 5 & 1 \\ 2 & 4 - 5 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{cc} -2 & 1 \\ 2 & -1 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \end{array} \right) \\ \Rightarrow \; y = 2x \\ \therefore \; \boldsymbol{x} = (t, 2t) \end{gather*}

こちらも確認してみると、

(10)    \begin{gather*} \boldsymbol{Ax} = \left( \begin{array}{cc} 3 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} t \\ 2t \end{array} \right) = \left( \begin{array}{c} 5t \\ 10t \end{array} \right) , \quad \lambda \boldsymbol{x} = 5 \left( \begin{array}{c} t \\ 2t \end{array} \right) = \left( \begin{array}{c} 5t \\ 10t \end{array} \right) \end{gather*}

なお、固有ベクトルを数値で表現する際、ノルムが1となるように正規化することが多い。

(11)    \begin{gather*} \boldsymbol{x} \rightarrow \frac{\boldsymbol{x}}{\| {\boldsymbol{x}} \|} \end{gather*}

上の例で固有値ベクトルを正規化すると以下の通り。

(12)    \begin{gather*} \frac{(t, -t)}{\sqrt{t^2 + t^2}} = \frac{(1, -1)}{\sqrt{2}} \approx (0.7071, -0.7071) \\ \frac{(t, 2t)}{\sqrt{t^2 + 4t^2}} = \frac{(1, 2)}{\sqrt{5}} \approx (0.4472, 0.8944) \end{gather*}

 

ユークリッドの互除法

概要

ユークリッドの互除法(Euclidean Algorithm)は、2つの自然数の最大公約数を求める手順。

2つの自然数a, b \; (a \ge b)の最大公約数(GCD: greatest common divisor)は、bと剰余r = a \mod bの最大公約数に等しいという性質を利用。数を順次割り込んでいき、剰余がゼロとなったときの除数が最大公約数となる。

(1)   }\begin{align*}a \mod b &= r_1 \quad (a \ge b) \\b \mod r &= r_2 \\&\vdots \\r_{n-1} \mod r_n &= 0 \\&\Downarrow \\\gcd(a, b) &= r_n\end{align*}

証明

以下の余りあり除算を考える。

(2)   \begin{gather*}a \mod b = r \Rightarrow a = bd + r\end{gather*}

a, bの公約数をmとすると、上式は以下のように変形され、mは余りrの公約数でもあることがわかる。

(3)   \begin{gather*}a = mp , \quad b = mq \\r = mp - mqd = m(p - qd)\end{gather*}

一方、b, rの公約数をnとすると、naの公約数であることがわかる。

(4)   \begin{gather*}b = ns , \quad r = nt \\a = nsd + nt = n(sd + t)\end{gather*}

これより、a, bの公約数の集合とb, rの公約数の集合は等しく、最大公約数も等しくなる。

(5)   \begin{equation*}a = bd+r \Rightarrow \gcd(a, b) = \gcd(b, r)\end{equation*}

計算例

143と91の最大公約数を求める。

(6)   \begin{align*}&143 \mod 91 = 52 \\&91 \mod 52 = 39 \\&52 \mod 39 = 13 \\&39 \mod 13 = 0 \\&\therefore \gcd(143, 91) = 13\end{align*}

再帰関数による実装

PythonとCLispの再帰関数による実装例は以下の通り。ただし、第1引数>第2引数を前提としており、エラー処理はしていない。

Python

CLisp

 

最大公約数・最小公倍数

約数

約数の定義

整数Nの約数(divisor, factor)とは、Nを割り切る整数(余りが生じない除数)。

整数mNの約数であるとき、m|Nと表し、ある整数aに対してN=maが成り立つことでもある。一般には自然数あるいは0以上の整数で考える。

通常はm \ne 0の条件を課すが、0も含める場合は、N=0の時に限り0が約数になる。

例えば12の約数は、以下の6個。

(1)   \begin{align*}&12 \div 1 = 12 \\&12 \div 2 = 6 \\&12 \div 3 = 4 \\&12 \div 4 = 3 \\&12 \div 6 =2 \\&12 \div 12 = 1\end{align*}

効率的な約数の求め方

m|Nであるとき、\frac Nm|Nである。これよりm = \frac Nmとして、\sqrt N以下の約数m_iを求め、あとはN/m_iを計算すれば、手間が半分で済む。

Nが平方数の場合は\sqrt Nも約数となり、約数の総数は奇数個、平方数でない場合は偶数異なる。

0、1の約数

0の約数は0 = maとなる整数mであり、0以上の全ての整数である。

1の約数は1 = maとなる整数mであり、1のみ。

1は全ての整数の約数。

素数の約数

素数の定義が「1と自身以外に約数を持たない数」なので、約数は2個。

公約数

2つの数の公約数は、それらの最大公約数の約数。

0、1との公約数

aと1の公約数は1のみ。

aと0の公約数は、aの約数全て(0の約数は0以上の全ての整数)。

公約数と剰余

a \mod b = rのとき、a, bの公約数はb, rの公約数でもある。

最大公約数

2つの数の最大公約数(greatest common divisor)を、\gcd(a, b)のように表す。

最大公約数を求める手順として、にユークリッドの互除法がある。

0、1との最大公約数

(2)   \begin{gather*}\gcd(a, 0) = a \\\gcd(a, 1) = 1\end{gather*}

最大公約数と剰余

被除数と除数の最大公約数は、除数と剰余の最大公約数でもある。

(3)   \begin{equation*}a \mod b = r \Rightarrow \gcd(a, b) = \gcd(b, r)\end{equation*}

倍数

倍数(multiple)とは、ある数a(整数に限らない)を整数倍した数である。

(4)   \begin{gather*}\cdots \; -3a,\; -2a,\; -a,\; 0,\; a,\; 2a,\; 3a,\; \cdots\end{gather*}

  • 0の倍数は0のみ
  • 0はすべての数の倍数
  • すべての数は自分自身の倍数
  • すべての整数は1と-1の倍数

最小公倍数

最小公倍数(least common multiple)とは、2つの整数の公倍数のうち、正で最小のもの。

たとえば36と56の最小公倍数は504。

最大公約数と最小公倍数の積

2つの数a, bの積は、それらの最大公約数と最小公倍数の積に等しい。

(5)   \begin{equation*}ab = \gcd(a, b) \cdot \operatorname{lcm}(a, b)\end{equation*}

 

剰余と合同式

商と剰余

剰余(余り)について

a \div bの商がq、余りがrのとき、以下のように表される。d:divisor、q:quotient、r:reminderの意味。

(1)    \begin{eqnarray*} &a = dq + r \\ &a \mod d = r \end{eqnarray*}

余りの定義を「割る数未満の自然数あるいは0」とすると、

(2)    \begin{eqnarray*} 8 &=& 3 \times 2 + 2 \\ 7 &=& 3 \times 2 + 1 \\ 6 &=& 3 \times 2 + 0 \\ 5 &=& 2 \times 2 + 2 \\ 4 &=& 2 \times 2 + 1 \\ 3 &=& 2 \times 2 + 0 \end{eqnarray*}

割られる数が負の場合には、

(3)    \begin{eqnarray*} -9 &=& 3 \times (-3) + 0 \\ -8 &=& 3 \times (-3) + 1 \\ -7 &=& 3 \times (-3) + 2 \\ -6 &=& 3 \times (-2) + 0 \\ -5 &=& 3 \times (-2) + 1  \\ -4 &=& 3 \times (-2) + 2 \end{eqnarray*}

割る数が負の場合には、

(4)    \begin{eqnarray*} 9 &=& -3 \times (-3) + 0 \\ 8 &=& -3 \times (-3) + 1 \\ 7 &=& -3 \times (-3) + 2 \\ 6 &=& -3 \times (-2) + 0 \\ 5 &=& -3 \times (-2) + 1 \\ 4 &=& -3 \times (-2) + 2 \\ \end{eqnarray*}

負数の余り

余りとして負の数を認めることもできる。ただしその場合、商と余りの組み合わせが1つとは限らない。

(5)    \begin{eqnarray*} 10 &=& -3 \times (-3) + 1 \\ 10 &=& -3 \times (-4) -2 \end{eqnarray*}

余りの定義(要件)は以下の2通りがあり、いずれを採用するかは任意。

  • 余りを割る数の絶対値より小さい0以上とする(0 \le r < |d|)
  • 余りの絶対値が割る数の絶対値より小さい数とする(0 \le |r| < |d|)

特別な場合の余り

割る数が1あるいは-1のときは、余りは常に0。

(6)    \begin{eqnarray*} a &=& 1 \times a + 0 \\ a &=& -1 \times (-a) + 0 \end{eqnarray*}

割られる数が1のときの余りは1、割られる数が-1なら余りは割る数の絶対値から1を現じた値(余りを正と定義した場合)。

(7)    \begin{eqnarray*} 1 &=& b \times 0 + 1 \\ -1 &=& b \times 1 + b - 1 \quad (b > 0) \\ -1 &=& b \times (-1)1 + (-b - 1) \quad (b < 0) \end{eqnarray*}

合同式

合同式の定義

整数a, bを正の整数dで割った余りが等しいとき、以下のように表記し、「a, bdを法として合同である」という。

(8)    \begin{eqnarray*} a &\equiv& b \mod{d}\\ a &\equiv& b \pmod d \end{eqnarray*}

これは次のようにも表現できる。

(9)    \begin{equation*} a - pd = b - qd = r \end{equation*}

合同式の例

(10)    \begin{eqnarray*} 7 &\equiv& 4 \bmod 3\\ 6 &\equiv& 4 \bmod 2\\ 6 &\equiv& 0 \bmod 2\\ 5 &\equiv& 1 \bmod 2 \\ 5 &\equiv& (-1) \bmod 3 \end{eqnarray*}

合同式の性質

以下、合同式の\mod dを省略する。

合同式の和

(11)    \begin{equation*} a_1 \equiv b_1 , \;  a_2\equiv b_2 \; \Rightarrow \; a_1 + a_2 \equiv b_1 + b_2 \end{equation*}

【証明】

 \begin{eqnarray} a_1 - p_1 d &=& b_1 - q_1 d \\ a_2 - p_2 d &=& b_2 - q_2 d \\ &\Downarrow& \\ (a_1 + a_2) - (p_1 + p_2)d &=& (b_1 + b_2) - (q_1 + q_2) d

合同式の積

(12)    \begin{equation*} a_1 \equiv b_1 , \;  a_2\equiv b_2 \; \Rightarrow \; a_1 a_2 \equiv b_1 b_2 \end{equation*}

【証明】

(13)    \begin{eqnarray*} a_1 - p_1 d &=& b_1 - q_1 d \\ a_2 - p_2 d &=& b_2 - q_2 d \\ &\Downarrow& \\ (a_1 - p_1 d)(a_2 - p_2 d) &=& (b_1 - q_1 d)(b_2 - q_2 d) \\ &\Downarrow& \\ a_1 a_2 - (a_1 p_2 + a_2 p_1 + p_1 p_2 d) d &=& b_1 b_2 - (b_1 q_2 + b_2 q_1 + q_1 q_2 d) d \end{eqnarray*}

合同式の商

a, \; dが互いに素(a \perp d)のとき、以下が成り立つ。

(14)    \begin{equation*} ab \equiv ac \bmod d \; \Rightarrow \; b \equiv c \bmod d \quad ( {\rm where} \; a \perp d ) \end{equation*}

【証明】

(15)    \begin{eqnarray*} ab \equiv ac \bmod d \; &\Rightarrow& \; ab = dp + r , \; ac = dq + r \\ &\Rightarrow& a(b - c) = d(p + q) \\ \end{eqnarray*}

ここでa, \; bは互いに素なので、b-cdの倍数となる。

(16)    \begin{eqnarray*} b - c = dk &\Rightarrow& \; b = dm + s , \; ac = dn + s \\ &\Rightarrow& b \equiv c \bmod d \\ \end{eqnarray*}

合同式の冪乗

(17)    \begin{equation*} a \equiv b \bmod d \; \Rightarrow \; a^n \equiv b^n \bmod d \end{equation*}

【証明】

(18)    \begin{equation*} a \equiv b \bmod d \quad \Leftrightarrow \quad \left\{ \begin{array}{l} a = pd + r \\ b = qd + r \end{array} \right. \end{equation*}

(19)    \begin{eqnarray*} a^n &=& (pd + r)^n = \sum_{k=0}^n \dbinom{n-k}{k} (pd)^{n-k} r^k \\ &=& pd \sum_{k=0}^{n-1} \dbinom{n-k}{k} (pd)^{n-k-1} r^k + r^n \\ b^n &=& (qd + r)^n = \sum_{k=0}^n \dbinom{n-k}{k} (qd)^{n-k} r^k \\ &=& qd \sum_{k=0}^{n-1} \dbinom{n-k}{k} (qd)^{n-k-1} r^k + r^n \end{array} \end{eqnarray*}

 

二項定理

二項定理の表現

(1)    \begin{align*} (x + y)^n &= \sum_{k=0}^n \binom{n}{k} x^{n-k} y^k \\ &= x^n + n x^{n-1}y + \cdots + \binom{n}{k} x^k y^{n-k} + \cdots + y^n \end{align*}

具体例

(2)    \begin{align*} (x + y)^2 &= x^2 + 2xy + y^2 \\ (x + y)^3 &= x^3 + 3x^2 y +3x y^2 + x y^3 \\ (x + y)^4 &= x^4 + 4x^3 y + 6x^2 y^2 + 4x y^3 + y^4 \end{align*}

証明

数学的帰納法で証明する。

n = 0のとき、

(3)    \begin{equation*} (x + y)^0 = \binom 00 x^0 y^0 = 1 \end{equation*}

n = 1のとき、

(4)    \begin{equation*} (x + y)^1 = \binom 10 x^1 y^0 + \binom 11 x^0 y^1 = x + y \end{equation*}

ここでn-1のときに以下が成り立つとする。

(5)    \begin{equation*} (x + y)^{n-1} = \sum_{k=0}^{n-1} \binom{n-1}{k} x^{n-1-k} y^k \end{equation*}

このときnに関しては、

(6)    \begin{eqnarray*} (x + y)^n &=& (x + y) \sum_{k=0}^{n-1} \binom {n-1}{k} x^{n-1-k} y^k \\ &=& \sum_{k=0}^{n-1} \binom {n-1}{k} (x^{n-k} y^k + x^{n-1-k} y^{k+1}) \end{eqnarray*}

これを右辺は以下のように展開できる。

 \begin{array}{rrrrrrrrrrrrrrr} & x^n & + & (n-1) x^{n-1} y & \cdots & \binom{n-1}{k} x^{n-k} y^k & \cdots & x y^{n-1} \\ + &&& x^{n-1} y & \cdots & \binom{n-1}{k-1} x^{n-k} y^k &\cdots & (n-1) x y^{n-1} & + & y^n \\ \hline & x^n & + & n x^{n-1} y & \cdots & \left( \binom{n-1}{k} - \binom{n-1}{k-1}\right) x^{n-k}{k} & \cdots & n x y^{n-1} & + & y^n \\ \end{array} \\

(7)    \begin{eqnarray*} \binom{n-1}{k} - \binom{n-1}{k-1} &=& \frac{(n-1)!}{(n-k-1)! k!} - \frac{(n-1)!}{(n-k)! (k-1)!} \\ &=& (n-1)! \frac{(n-k) - k}{(n-k)! k!} = \binom nk \end{eqnarray*}

これより、n \ge 0に対して(1)が証明された。

パスカルの三角形

数のような数による三角形をパスカルの三角形(Pascal’s triangle)と呼び、n行目の数の列が二項展開のn乗の係数となっている。

各項の計算の仕方は、一つ上の段の左右の数の和として求めていく(左端・右端の外側はゼロと考える)。

pascals-triangle

パスカルの三角形の各項が二項展開の計数となること、すなわちm段目のn項目の数をa_{m \cdot n}とし、これが\binom{m-1}{n-1}となることを、数学的帰納法で証明する。

(8)    \begin{eqnarray*} a_{2 \cdot 1} &=& \binom{1}{0} = 1 \\ a_{2 \cdot 2} &=& \binom{1}{1} = 1 \end{eqnarray*}

(9)    \begin{eqnarray*} a_{m \cdot n} &=& a_{m-1 \cdot n-1} + a_{m-1 \cdot n} = \binom{m-2}{n-2} + \binom{m-2}{n-1} \\ &=& \frac{(m-2)!}{(m-n)!(n-2)!} + \frac{(m-2)!}{(m-n-1)!(n-1)!} \\ &=& (m-2)! \frac{(n-1) + (m-n)}{(m-n)!(n-1)!} = \frac{(m-1)!}{(m-n)!(n-1)!} \\ &=& \binom{m-1}{n-1} \end{eqnarray*}

 

微分の公式

微分の定義

初等的な微分の定義は以下の通り。

(1)    \begin{equation*} f'(x) = \frac{d f(x)}{dx} = \lim_{\Delta x \rightarrow 0} \frac{f(x + \Delta x) - f(x)}{\Delta x} \end{equation*}

関数の積・商の微分

関数の積の微分

(2)    \begin{equation*} f(x) = u(x) v(x) \rightarrow f'(x) = u(x) v'(x) + u'(x) v(x) \end{equation*}

 

(3)    \begin{equation*} f'(x) = \lim_{\Delta x \rightarrow 0} \frac{u(x + \Delta x) v(x + \Delta x) - u(x) v(x)}{\Delta x} \end{equation*}

ここで

(4)    \begin{eqnarray*} u(x + \Delta x) &=& u(x) + u'(x) \Delta x + o({\Delta x}^2) \\ v(x + \Delta x) &=& v(x) + v'(x) \Delta x + o({\Delta x}^2) \end{eqnarray*}

したがって、

(5)    \begin{equation*} u(x + \Delta x) v(x + \Delta x) = u(x) v(x) + u(x) v'(x) \Delta x + u'(x) v(x) \Delta x + o({\Delta x}^2) \end{equation*}

これより以下を得る。

(6)    \begin{eqnarray*} f'(x) &=& \lim_{\Delta x \rightarrow 0} \frac{u(x) v'(x) \Delta x + u'(x) v(x) \Delta x + o({\Delta x}^2)}{\Delta x} \\ &=& u(x) v'(x) + u'(x) v(x) \end{eqnarray*}

関数の商の微分

(7)    \begin{equation*} f(x) = \frac{v(x)}{u(x)} \rightarrow f'(x) = \frac{u(x) v'(x) - u'(x) v(x)}{u(x) ^2} \end{equation*}

 

(8)    \begin{equation*} f(x) = \frac{v(x)}{u(x)} \Leftrightarrow v(x) = f(x) u(x) \end{equation*}

(9)    \begin{equation*} v'(x) = f'(x) u(x) + f(x) u'(x) = f'(x) u(x) + \frac{v(x) u'(x)}{u(x)} \end{equation*}

これより(7)を得る。

合成関数の微分

(10)    \begin{eqnarray*} y = f(u(x)) & \rightarrow & \frac{dy}{dx} = f'(u(x))u'(x) \\ &{\rm or}& \\ \frac{dy}{dx} &=& \frac{dy}{du} \frac{du}{dx} \end{eqnarray*}

 

u(x + \Delta x) - u(x) = \Delta u(x)とおくと、

(11)    \begin{eqnarray*} && lim_{x \rightarrow \infty} \frac{f(u(x + \Delta x)) - f(u(x))}{\Delta x} \\ &=& lim_{x \rightarrow \infty} \frac{f(\Delta u(x)) - f(u(x))}{\Delta u(x)} \frac{\Delta u(x)}{\Delta x} \\ &=& f'(u(x))u'(x) \end{eqnarray*}

逆関数の微分

(12)    \begin{eqnarray*} y = f^{-1}(x) & \rightarrow & \frac{dy}{dx} = \frac{1}{f'(x)} \\ &{\rm or}& \\ \frac{dx}{dy} &=& \frac{1}{\frac{dy}{dx}} \end{eqnarray*}

 

(13)    \begin{equation*} y = f^{-1}(x) \rightarrow f(y) = x \end{equation*}

とおいて、合成関数の微分より、

(14)    \begin{eqnarray*} f'(y) \frac{dy}{dx} = 1 \\ \therefore \frac{dy}{dx} = \frac{1}{f'(y)} \end{eqnarray*}

媒介変数表示

(15)    \begin{equation*} x = u(t), y = v(t) \rightarrow \frac{dy}{dx} = \frac{\frac{dy}{dt}}{\frac{dx}{dt}} = \frac{v'(t)}{u'(t)} \end{equation*}

 

\Delta u(t) = u(t + \Delta t) - u(t)\Delta v(t) = v(t + \Delta t) - v(t)とおいて、

(16)    \begin{equation*} \lim_{n \rightarrow \infty} \frac{\Delta v(t)}{\Delta u(t)} = \lim_{n \rightarrow \infty} \frac{v(t + \Delta t) - v(t)}{\Delta t} \frac{\Delta t}{u(t + \Delta t) - u(t)} = \frac{v'(t)}{u'(t)} \end{equation*}