2017-07-10

再開します

最後の更新から一年以上が経過していましたが、再開したいと思います

ここ一年ぐらいあまり勉強できてなかったというのもあるので、インプットしつつ、このブログでアウトプットできればいいなあと

ゆるやかに続けていきたいです

2016-03-04

Gauss-Seidel法

前回はヤコビ法について書きましたが、今回はそこから一歩推し進めたGauss-Seidel(ガウス・ザイデル)法について書きます。ヤコビ法では収束までの反復数が多いというのが問題でした。これを解決するために考えだされたのがガウス・ザイデル法です。簡単に違いを述べると、ヤコビ法は一回前の反復で得られた近似解を使って新しいより良い近似解を作ろうというものでしたが、ガウス・ザイデル法は更新が終わった近似解の値、つまり新しい近似解と古い近似解の両方を使って更新を行っていくというものです。
Ax=bを解く上で、k反復目の近似解を $x^k$ とします。例えばk+1反復目の近似解のi行目を求めることを考えます。その際、1～i-1行目の近似解に関しては既に更新が終わっているのでk+1反復目のものを使い、i+1行目からN行目までの近似解に関してはまだ更新が終わってないのでk反復目までのものを用います。最も新しい近似解を随時使っていくことによって収束を早めることが可能となります。
プログラムは下のような形になります。

void solve_gauss_seidel(real *dmat, real *x, real *y, int M, int N) {
  
  int i, j, k, it;
  int M, N;
  real old_x, r; //residual

  /* Set initial solution */
  for (i = 0; i < N; i++) {
    x[i] = 0;
  }
  
  /* Iteration */
  for (it = 0; it < MAXIT; it++) {
    r = 0;
    for (i = 0; i < M; i++) {
      old_x = x[i];
      x[i] = y[i];
      for (j = 0; j < N; j++) {
	if (i != j) {
	  x[i] -= dmat[i * M + j] * x[j];
	}
      }
      x[i] /= dmat[i * M + i];
      r += fabs(x[i] - old_x);
    }

    if (r < EPS) {
      break;
    }

    cout << it << ": ";
    print_vector(x, M);

  }

}

前回のヤコビ法の実装では古い近似解の値を残しておきましたが、ガウス・ザイデル法では随時更新していくことから残して置く必要は無いため、消えています。手元にある適当なテストデータで実験してみたところ、収束までの反復数も削減されており、ヤコビ法と比較して収束のスピードが上がっていることを確認しました。ただ、近似解の行間に依存関係があるため、並列化が行いづらいという問題があります。

SOR法

がっつり書くつもりは無いですが、更に収束のスピードを上げる手法としてSOR(Successive Over-Relaxation)法(逐次緩和法)というものがあります。これはガウス・ザイデル法で得られた近似解をそのままその反復での解にするのではなく、その値とひとつ前の近似解を組み合わせたものを新たな近似解とする、というものです。k+1反復目でガウス・ザイデル法によって得られた解を $x^*_{k+1}$ とし、k反復目での近似解を $x_k$ とすると、新たな近似解をパラメータ $\omega$ を用いて
$x_{k+1}=(1-\omega)x_{k} + \omega x^*_{k+1}$
としたものが、SOR法です。
パラメータ $\omega$ は $0<\omega<2$ の範囲で設定しなければ収束はしないそうです。ただうまく設定してあげればガウス・ザイデル法よりも収束が高速化されます。ちなみに $\omega=1$ だと単なるガウス・ザイデル法になります。ここらへんのパラメータ調整的な話はきりが無い気がするので、SOR法についてこのぐらいで終わりにしたいと思います。

2016-03-04

Jacobi法

ここまで連立一次方程式を解く手法としてガウスの消去法をベースとしたものを見てきましたが、ここからいくつかは反復法と呼ばれるものを見ていきたいと思います。LU分解などでは行列の値を直接操作することによって、解を求めていました。このような手法での計算量は $O(N^3)$ であり、行列のサイズが大きくなると解くのが困難となってきます。また、係数行列が疎行列(行列の内部の要素のほとんどがゼロであるような行列)であるような場合には、直接法ではフィルイン(ゼロだった要素がゼロでなくなる、つまり疎行列であることの利点を活かしづらくなる)が発生してしまいます。反復法というのはある適当な初期解を設定した後、繰り返し計算を行うことによって、求めるべき解に近づけていくというものです。正しい解を得られるか、十分に解との差が小さくなるまで処理を行います。反復法には大きく分けて、定常反復法と非定常反復法の二つがあります。今回書くJacobi法は定常反復法の方になります。反復法として広く知られている共役勾配法は非定常反復法にカテゴライズされます。では早速Jacobi法について、見ていきたいと思います。

Jacobi法の流れ

Jacobi(ヤコビ)法の大まかなイメージとしては、繰り返し処理を行う中で、ひとつ前の繰り返しで得られた近似解をもとに新たな近似解を生成するというものです。
数式的に表すと、Ax=bを解く上で、まず係数行列Aを
$A=L + D + U$
という形で分解をします。なお、Lは下三角行列(対角成分もゼロ)、DはAの対角成分を並べた対角行列、Uは上三角行列(対角成分はゼロ)となります。そして、k回目の反復で得られた近似解 $x_k$ からk+1回目の近似解を次のように定義するのがヤコビ法です。
$x^{k+1} = D^{-1}(-(L + U)x^k + b)$
ここでDの逆行列についてですが、対角行列ですので、単純にDの各要素の逆数を取ったものがDの逆行列となります。数式のイメージとしてはxのi行目の値の更新を行う際には、ひとつ前の近似解のi行目以外のデータをもとにi行目の方程式を満たすように解を求めてあげるというようなものになります。繰り返し処理を行っていくと最終的に近似解が固まるようになります。つまり、]
$x^* = D^{-1}(-(L + U)x^* + b)$
となります。
この式の両辺にDを掛けて項を移動すれば $Ax^*=b$ となり、 $x^*$ がこの方程式の解であることがわかります。
このように流れとしては非常に単純でわかりやすいものとなっています。

Jacobi法のプログラム

Jacobi法のプログラムが下のものになります。中身の処理自体に関しては上に書いた流れのとおりです。追加して書くべきこととしては、最大の反復回数を設定してあげることと、どこで反復を打ち切るのかについてです。最大の反復回数に関しては対象とする行列のサイズに併せて設定してあげればいいと思われます。どこで反復を打ち切るのか、つまり、近似解と解との誤差が十分に小さいと判断するのはどうやって行うかについてですが、ここでは新たに得られた近似解とその前の反復で得られた解について、差の絶対値をとったものを採用し、それが十分に小さいか(例えば、EPS=1.0e-5、EPSは精度を表すイプシロンの略)を判定し、小さければ反復終了、大きければ反復を継続というプログラムになっています。数式で表すと
$|x^{k}-x^{k-1}| < \mathrm{EPS}$
というような感じになります。

void solve_jacobi(real *dmat, real *x, real *y, int M, int N) {
  
  int i, j, k, it;
  real r; //residual
  real *x_old;
  x_old = (real *)malloc(sizeof(real) * N);

  /* Set initial solution */
  for (i = 0; i < N; i++) {
    x[i] = 0;
    x_old[i] = 0;
  }
  
  /* Iteration */
  for (it = 0; it < MAXIT; it++) {
    
    for (i = 0; i < M; i++) {
      x[i] = y[i];
      for (j = 0; j < N; j++) {
	if (i != j) {
	  x[i] -= dmat[i * M + j] * x_old[j];
	}
      }
      x[i] /= dmat[i * M + i];
    }

    r = 0;
    for (i = 0; i < N; i++) {
      r += fabs(x[i] - x_old[i]);
    }
    
    if (r < EPS) {
      break;
    }

    cout << it << ": ";
    print_vector(x, M);
    for (i = 0; i < N; i++) {
      x_old[i] = x[i];
    }

  }

}

ヤコビ法の利点としては各行の処理を独立して行えるので並列化を行いやすいというところですね。ただ収束までの反復数が多いというのが難点です。
この次は順番的に反復数の削減を可能にしたガウス・ザイデル法についてかな。

2016-03-04

Cholesky分解

今回は特定の条件を満たした時に使えるLU分解として、Cholesky(コレスキー)分解について書きます。特定の条件というのは、方程式Ax=bについて、係数行列Aが正定値対称行列である、というものです。対称行列というのは転置を行っても同じ行列になる、つまり $A^T = A$ となるものです。正定値行列というのは行列と同じサイズの任意のベクトルxを持ってきて常に、 $x^TAx > 0$ となる行列Aを言います。このような行列はすべての固有値が正であるなど特殊な性質を持っています。このような正定値対称行列に対して例えばLU分解を行ってあげると、 $L^T = U$ となります。つまり、Lが求まれば自動的にUも求められるということになります。LU分解ではLとUのそれぞれを計算して求めていましたが、このような特殊な性質を用いてあげることによって、計算量を半分近くまで減らすことが可能となります。これがコレスキー分解の目標となります。

コレスキー分解の流れとプログラム

実質的にやっていることはLU分解とほとんど変わらず、行列の左上からひたすら値の更新を行っていくというものになります。
①対象となる行を選択
②その行の対角成分に着目
③分解後の対角成分を求める
④対角成分と同じ列の要素を更新する(下三角行列Lにおいて、着目している対角成分の下側を更新する<=>上三角行列Uにおいて、着目している対角成分の右側を更新する)
⑤①から④を繰り返す

③や④の操作に関しては、対称行列であるという特徴を活かすと、うまく計算量を削減することが可能になります。

プログラムは下のようになります。

void cholesky_decompose(real *dmat, int M, int N) {

  int i, j, k;
  real vii, vji;

  for (i = 0; i < M; i++) {
    vii = dmat[i * N + i];
    for (j = 0; j < i; j++) {
      vii -= dmat[i * N + j] * dmat[i * N + j];
    }
    dmat[i * N + i] = sqrt(vii);

    for (j = i + 1; j < M; j++) {
      vji = dmat[j * N + i];
      for (k = 0; k < i - 1; k++) {
        vji -= dmat[j * N + k] * dmat[i * N + k];
      }
      dmat[j * N + i] = vji / dmat[i * N + i];
    }
  }

}

計算式にルート計算が入るので、計算のコストが大きいというのはあります。また、正定値性を持つ行列なら問題は無いのですが、ルートの中が正でなければ実数空間での処理は出来ないのでその点は問題です。
通常のコレスキー分解を改良することで、ルートを使わず、かつ対象範囲を対称行列(正定値性は問わない)まで拡張した修正コレスキー分解というのも存在して、こちらのほうがよく使われているらしいです。ざっくり調べた感じでは、修正コレスキー分解では対角成分だけ抜き出して考えるというもののようです。つまり、対角行列Dを持ってきて、
$A = L^TDL$
と分解するというものです。従来のコレスキー分解同様Lに関してはLU分解の半分の計算量ですみ、対角成分の計算に関してはルート計算が発生しないという強みがあります。

ここまではLU分解とCholesky分解という二つの代表的な連立一次方程式の解法を見てきました。この二つはどちらも直接法と呼ばれる方法に分類されていて、行列を直接操作することによって、解を厳密に求めていくというものになります。これに対して、反復法と呼ばれる異なるアプローチで解を求めにいく方法があります。今後はそちらの方を見ていく予定です。

2016-03-03

LU分解

前回はガウスの消去法について書きましたが、今回はそこから一歩進めた方法であるLU分解について記します。
簡単に説明するとある方程式Ax=bの係数行列Aについて、下三角行列Lと上三角行列Uを用いて
${\displaystyle A=LU }$
という形で分解してあげるのがLU分解となります。あとは方程式のほうに戻ってあげるとLUx=bとなりますので、まずUx=cとでも置いて、Lc=bを解きます。cが求まった後にUx=cを解くと、無事にxを求めることが出来ます。Lは下三角行列なので、Lc=bを解く操作を前進代入、Uは上三角行列なので、Ux=cを解く操作を後進代入といいます。このように分解をすると何が嬉しいのかについてですが、ガウスの消去法では方程式の右辺ベクトルを含めて求解を進めていきます。そのため、行列が同じでも右辺ベクトルが異なるような方程式が現れた場合にはもう一度最初から解く必要が出てきます。このような行列が同じで右辺ベクトルだけが異なるというようなシチュエーションは多いらしく、ガウスの消去法では非常にコストが大きいという問題があります。LU分解を用いて予め行列Aを分解しておくことによって、右辺ベクトルが変わった際に前進代入と後進代入をさくっと行うだけで解を求めることができるようになります。そのため、方程式を解く上で非常に有用な手段であるといえます。また、この分解操作をどのように行うかはこの後書いていく予定ですが、この分解操作をうまくやれば行列行列積を使って行うことが出来るようです。行列行列積はメモリアクセス量に対して演算が多く、キャッシュで再利用させるなどうまく処理してあげることによって高い性能を出すことが可能な演算です。なので、単純に入れ替えだけを行う操作に比べて高い性能を出すことがLU分解では出来るようです。

LU分解のアルゴリズム

LU分解のアルゴリズムとして、right-looking、left-looking、Crout法の3つがあります。ここではそれぞれどのようなアルゴリズム化を見ていきたいと思います。併せてプログラムもここで記します。

Right-looking

大まかな流れとしては
①対象となる行を決める
②対象となる行の対角成分を見る
③対角成分をもとに、その成分と同じ列の要素を更新(なお対象となるのはその成分より下の行の要素)
④対象となる行の要素を用いて、対角成分より右下の要素を更新していく
⑤この操作をすべての行に対して行う

対象となる対角成分の右部分を参照しながら右下部分の更新を行うことからright-lookingと呼ばれています。なお、Lの対角成分が1となることを仮定しています。

プログラムで表すと下のような形になります。ちなみに行列はMxNとします。

void lu_decompose_right(real *dmat, int M, int N) {
  
  int i, j, k;
  real temp, vij;

  for (i = 0; i < 1; i++) {
    temp = 1.0 / dmat[i * N + i];
    for (j = i + 1; j < M; j++) {
      dmat[j * N + i] *= temp;
    }
    for (j = i + 1; j < N; j++) {
      vij = dmat[i * N + j];
      for (k = i + 1; k < M; k++) {
	dmat[k * N + j] -=dmat[k * N + i] * vij;
      }
    }
  }
}

Left-looking

Right-lookingでは右側を見て右側を更新していたのに対して、left-lookingでは左側を見て更新を行っていきます。これもLの対角成分が1になることを仮定しています。
こちらもおおまかな流れを書くと
①対象となる行を決める
②その行の対角成分を見る
③その対角成分と同じ列の上の部分と左側の部分を元にして上の部分を更新する
④最後に対角成分と同じ列の下の部分を対角成分を元にして更新
⑤これを繰り返す

④はright-lookingと同じですが、①～③の流れは完全に別物になっています。

プログラムは以下のようになります。

void lu_decompose_left(real *dmat, int M, int N) {

  int i, j, k;
  real vji, vii;

  for (i = 0; i < M; i++) {
    for (j = 0; j < i; j++) {
      vji = dmat[j * N + i];
      for (k = j + 1; k < N; k++) {
	dmat[k * N + i] -= dmat[k * N + j] * vji;
      }
    }
    vii = 1.0 / dmat[i * N + i];
    for (k = i + 1; k < N; k++) {
      dmat[k * N + i] *= vii;
    }
  }
}

今回は行方向に連続になるように行列のメモリ配置を取っていましたが、列方向に連続になるように取ってあげれば、更新の部分がメモリ連続になって良さそうだなあとか思ったり。

Crout法

Crout法では右だけとか左だけとではなく、様々なところを参照しながら更新を行っていきます。また、Crout法ではUの対角成分を1とします。
流れとしては
①対象となる行を決める
②その行の対角成分を見つける
③その対角成分の左下と右上を元に、対角成分と同じ行と列の更新を行う
④これをすべての行について繰り返す

プログラムは下のようになります。

void lu_decompose_crout(real *dmat, int M, int N) {

  int i, j, k;
  real vii;  

  for (i = 0; i < M; i++) {
    for (k = i; k < M; k++) {
      for (j = 0; j < i; j++) {
	dmat[k * N + i] -= dmat[k * N + j] * dmat[j * N + i];
      }
    }
    vii = dmat[i * N + i];
    for (j = i + 1; j < N; j++) {
      for (k = 0; k < i; k++) {
	dmat[i * N + j] -= dmat[k * N + j] * dmat[i * N + k];
      }
      dmat[i * N + j] /= vii;
    }
  }
}

Crout法では更新する際に2箇所の参照領域が必要となります。そのため、行列全体を分割して並列に処理を行うような場合には、参照領域を求めてプロセス間で通信が発生してしまいます。そのため、並列計算にはあまり向いてはいません。ただ、ノード内でメモリが共有されているような状態でしたら問題はありません。

ここまでで3つのアルゴリズムの紹介が終わりました。どうやら並列計算とかを念頭に置くと、right-lookingが最も良いらしいですね。あと、最初に行列行列積がうんぬんと言いましたが、どうやらブロック形式のアルゴリズムで用いられるらしいです。今回は力尽きたのでここまでとしますが、やる気が出た時にでもブロック形式は勉強したいなあというところです。

2016-03-01

ガウスの消去法

連立一次方程式を解くというような話になった時に一番最初に出てくるであろうガウスの消去法から取り上げていきたいと思います。おそらくこれは中高生のころに習った方程式の解き方みたいな感じなので理解はしやすいです(筆者がそうでした)。
何をしているのかを簡単に説明すると、連立なので複数の式があります。これらの式を足したり、引いたり、掛けたうえで加減したすることで変数を消していくというようなものです。最終的に一つの変数の値がわかれば、芋づる式に他の変数の値が求められるようになります。前半部分の変数を消していく作業は数式の左側の変数から順次やっていって、最終的に方程式の係数を表す行列である係数行列が上三角にするのが目的です。
具体的に例を使って説明してみます。

$3 x_0 + 5 x_1 - 7 x_2 = 0 \\ x_0 + x_1 - x_2 = 2 \\ 2 x_0 - 3 x_1 + x_2 = 5$
この方程式を解くことを考えます。
まず $x_0$ に着目します。
すると、2行目 - 1行目 * 1/3を行うことによって、2行目の $x_0$ は消え、3行目 - 1行目 * 2/3を行うことで3行目の $x_0$ は消えます。
このように1つずつ変数を消していく作業を行うのですが、次の $x_1$ については、3行目にだけ行ってあげ、上の1行目には行いません。これは後で簡単に値が求められるようになるので、わざわざ無駄な計算をする必要は無いというものです。ちなみに逐一上の行の変数も消していく手法はガウス・ジョルダン法といいます。
最終的に
$3 x_0 + 5 x_1 - 7 x_2 = 0 \\ -2/3 x_1 + 4/3 x_2 = 2 \\ -7 x_2 = 14$
となります。
ここまでの操作を前進消去(Forward elimination)といいます。
ここから1つずつxの値を求めていく作業になりまして、これを後進代入(Backward substitution)といいます。
下の方の数式から順番に値を求めていくだけという簡単な話で
$x_2 = -2, x_1 = 1, x_0 = 3$ と計算できます。

このような形で求めていくのがガウスの消去法で非常に簡単に理解できるものです。
ただ、係数行列の対角成分がゼロであるような場合にはうまく消去が出来ません。そういった際には対象とする数式を変更するというようなことを行います。ピボット選択と言われる処理になります。ここでは面倒なので割愛します。

ここからが本番でプログラムを書いてみます。
ちなみに行列の表記についてですが、私は2次元配列を使うのがあまり好きではないので、いつも1次元で表しています。この際にメモリの連続性をどう取るかという話が出てきます。行方向に連続に取るか、列方向に連続に取るかというものです。ここでは行方向に連続になるように配置を考えます。つまり、メモリ上には $a_{1,1}, a_{1,2}, a_{1,3},...,a_{1,n}, a_{2,1},...$ と並んでいる状態になります。
解法部分のプログラムが下のものになります。
ちなみにrealは浮動小数点数を表す型でこちらで作成したものです。また、行列はMxNのものを考えています。

void solve_gauss(real *dmat, real *x, real *y, int M, int N) {

  int i, j, k;
  real d;
  real *eq;
  
  eq = (real *)malloc(sizeof(real) * M * (N + 1));
  for (i = 0; i < M; i++) {
    for (j = 0; j < N; j++) {
      eq[i * (N + 1) + j] =dmat->val[i * N + j];
    }
    eq[i * (N + 1) + N] = y[i];
  }
 
  /* Forward elimination */
  for (i = 0; i < M; i++) {
    for (j = i + 1; j < M; j++) {
      d = eq[j * (N + 1) + i] / eq[i * (N + 1) + i];
      for (k = i; k <= M; k++) {
        eq[j * (N + 1) + k] -= d * eq[i * (N + 1) + k];
      }
    }
  }

  /* Backward substitution */
  for (i = M - 1; i >= 0; i--) {
    d = eq[i * (N + 1) + N];
    for (j = i + 1; j < M; j++) {
      d -= eq[i * (N + 1) + j] * eq[j * (N + 1) + M];
    }
    eq[i * (N + 1) + M] = d / eq[i * (N + 1) + i];
  }
  
  /* Copy to x vector */
  for (i = 0; i < M; i++) {
    x[i] = eq[i * (N + 1) + N];
  }
}

順を追って最初から説明していきますと、まずはじめに行列と右辺ベクトルから数式eqを作成します。その後、前進消去を行い、上三角行列を作ったら、後進代入を行っていきます。最終的に右辺ベクトルが合った場所に求めるxの解がありますので、それをxにコピーしてあげることで終了となります。

配列のアドレス計算に気をつけてあげれば、難しくはないプログラムです。

ガウスの消去法では前進消去のところに関しては、並列化が簡単に出来そうだなあと思いました。2番目のjループを並列化すればある程度性能は上がりそうですね。iが大きくなると並列度は下がりますが。

というわけでざっくりした感じでしたが、ガウスの消去法についてはここまでということで。

全然関係ないけど、数式うまく書けるようにしときたいなあ

2016-02-28

連立一次方程式を解く

数値計算の分野ではどうやら連立一次方程式を解くということが多く行われているようです。私は詳しいことはわからないですが、微分方程式とかで表されているようなものもうまくやってあげると連続であったものを離散化して扱うことが出来るようでして、つまるところプログラム出来るような形式になるようです。

連立一次方程式が現れて、「さあ、これを解こう！」という流れになるようですが、その際には色々な方法があるようです。

方程式については各変数の係数をまとめた係数行列をAとして

$Ax = b (A:m \times n, x:n \times 1, b: m \times 1)$

みたいな形で表されます。

行列Aの特性には単純にエルミート行列(Aを転置して共役をとったものとAが等しい)であるだとか、正定値行列だとか様々なものがありまして、特性に合わせてより良い方法というものがあります。また、そんなものは関係なくとりあえず解くという感じのもあります。

ここからはガウスの消去法あたりからぼちぼち書いていこうかなというところで、また次回