ガウスの消去法の前進消去を拡張すると正方行列の逆行列を求めることができる. これは,ガウス・ジョルダン法(Gauss-Jordan method),もしくは,掃出し法(sweeping-out method)と呼ばれる. まず,&ref(ls_gauss_jordan.eq1.gif,nolink,70%);の正方行列&ref(ls_gauss_jordan.eq2.gif,nolink,70%);を考える. &ref(ls_gauss_jordan.eq3.gif,nolink,70%);ならば&ref(ls_gauss_jordan.eq2.gif,nolink,70%);は正則行列(regular matrix)である (逆に&ref(ls_gauss_jordan.eq4.gif,nolink,70%);なら特異行列(singular matrix)). &ref(ls_gauss_jordan.eq1.gif,nolink,70%);の単位行列を&ref(ls_gauss_jordan.eq5.gif,nolink,70%);とすると,&ref(ls_gauss_jordan.eq2.gif,nolink,70%);が正則ならば,&ref(ls_gauss_jordan.eq6.gif,nolink,70%);を満足する行列&ref(ls_gauss_jordan.eq7.gif,nolink,70%);が1つだけ存在する. この&ref(ls_gauss_jordan.eq7.gif,nolink,70%);を逆行列と呼び&ref(ls_gauss_jordan.eq8.gif,nolink,70%);で表される. さて,ガウスの消去法で対象とした連立1次方程式&ref(ls_gauss_jordan.eq9.gif,nolink,70%);において,&ref(ls_gauss_jordan.eq10.gif,nolink,70%);を&ref(ls_gauss_jordan.eq7.gif,nolink,70%);に,&ref(ls_gauss_jordan.eq11.gif,nolink,70%);を&ref(ls_gauss_jordan.eq5.gif,nolink,70%);に置き換えることを考える. そうすると拡大行列は以下のようになる. #ref(ls_gauss_jordan.eq12.gif,nolink,70%) &ref(ls_gauss_jordan.eq6.gif,nolink,70%);であるので,両辺に&ref(ls_gauss_jordan.eq8.gif,nolink,70%);をかけると #ref(ls_gauss_jordan.eq13.gif,nolink,70%) となる. つまり,拡大行列の左側半分が単位行列となるように変換することで,右側半分に逆行列ができる. #ref(ls_gauss_jordan.eq14.gif,nolink,70%) このとき,逆行列&ref(ls_gauss_jordan.eq8.gif,nolink,70%);は, #ref(ls_gauss_jordan.eq15.gif,nolink,70%) となる. 前進消去を拡張し,これらの処理を行うことで逆行列を求めるのが, ガウス・ジョルダン法である. -ガウス・ジョルダン法の実装~ ***ガウス・ジョルダン法の実装 [#u219ca0b] ガウスの消去法における前進消去では上三角行列にするために, &ref(ls_gauss_jordan.eq16.gif,nolink,70%);の条件を満たす要素のみ対象として処理していたが, ガウス・ジョルダン法ではすべての要素に対象を広げる. ガウス・ジョルダンでの前身消去の式は以下である. #ref(ls_gauss_jordan.eq17.gif,nolink,70%) ガウス・ジョルダン法で逆行列を求めるコード例を以下に示す. #code(C){{ /*! * ガウス・ジョルダン法(ピボット選択なし) * @param[inout] A n×2nの拡張行列 * @param[in] n n元連立一次方程式 * @return 1:成功 */ int GaussJordan(vector< vector<double> > &A, int n) { // 拡張行列の右半分を単位行列にする for(int i = 0; i < n; ++i){ for(int j = 0; j < n; ++j){ A[i][j+n] = (i == j ? 1.0 : 0.0); } } // ガウス・ジョルダン法(Gauss-Jordan method)で逆行列計算 for(int k = 0; k < n; ++k){ double akk = A[k][k]; // 対角要素を1にするために,k行目のすべての要素をa_kkで割る for(int j = 0; j < 2*n; ++j){ A[k][j] /= akk; } // k列目の非対角要素を0にする for(int i = 0; i < n; ++i){ if(i == k) continue; double aik = A[i][k]; for(int j = 0; j < 2*n; ++j){ A[i][j] -= A[k][j]*aik; } } } return 1; } }}