連立一次方程式で遊ぶ

プログラミング

はじめに

連立一次方程式(Ax=b)を直接法(ガウスジョルダン法)、反復法(最急降下法、共役勾配法)でちょっと解いてみた。

コード

#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>

static const double EPS = 1e-8;
typedef std::vector<double> vec; //ベクトル
typedef std::vector<vec> mat; //行列

struct mat_utility { //ちゃんとチェックしてない
  //余因子
  double cofact(const mat& A, int i, int j){
    mat B(A.size()-1, vec(A.size()-1, 0));
    int k = 0, l = 0;
    for(int ii=0; ii<A.size(); ii++) if(ii != i){
	l = 0;
	for(int jj=0; jj<A[ii].size(); jj++) if(jj != j){
	    B[k][l] = A[ii][jj];
	    l++;
	  }
	k++;
      }
    return det(B) * ((i+j)%2==0?1:-1);
  }
  
  //行列式
  double det(const mat& A){
    //1x1行列の場合
    if(A.size() == 1) return A[0][0];
    //2x2行列の場合
    if(A.size() == 2) return A[0][0] * A[1][1] - A[0][1] * A[1][0];
    //それ以外の場合(余因子展開)
    double d = 0;
    for(int i=0; i<A.size(); i++){
      d += A[i][0] * cofact(A, i, 0);
    }
    return d;
  }

  //正定値判定(先頭主部分行列の行列式がすべて正かどうかで判定)
  // A : 対称行列
  bool is_positive_definite(mat A){
    for(int i=0; i<A.size(); i++){
      double d = det(A);
      if(d < EPS) return false;
      A.pop_back();
      for(int j=0; j<A.size(); j++){
	A[j].pop_back();
      }
    }
    return true;
  }

  //対称行列判定
  bool is_symmetric(const mat& A){
    for(int i=0; i<A.size(); i++){
      if(A.size() != A[i].size()) return false;
      for(int j=i+1; j<A[i].size(); j++){
	if(A[i][j] != A[j][i]) return false;
      }
    }
    return true;
  }
};

//ガウスジョルダン法
class gauss_jordan {
public:
  vec calc(const mat& A, const vec& b){
    int n = A.size();
    mat B(n, vec(n+1));
    //拡大係数行列
    for(int i=0; i<n; i++){
      for(int j=0; j<n; j++){
	B[i][j] = A[i][j];
      }
    }
    for(int i=0; i<n; i++){
      B[i][n] = b[i];
    }

    for(int i=0; i<n; i++){
      //注目している変数の係数の絶対値が大きい式をi番目に持ってくる
      int pivot = i;
      for(int j=i; j<n; j++){
	if(fabs(B[j][i]) > fabs(B[pivot][i])){
	  pivot = j;
	}
      }
      swap(B[i], B[pivot]);

      //解なし、一意でない(pivotが小さすぎる)
      if(fabs(B[i][i]) < EPS) return vec();

      //注目している変数の係数を1にする
      for(int j=i+1; j<=n; j++){
	B[i][j] /= B[i][i];
      }

      //j番目の式からi番目の変数を消去する
      for(int j=0; j<n; j++){
	if(i != j){
	  for(int k=i+1; k<=n; k++){
	    B[j][k] -= B[j][i] * B[i][k];
	  }
	}
      }
    }

    //結果
    vec x(n);
    for(int i=0; i<n; i++){
      x[i] = B[i][n];
    }
    return x;
  }
};

//最急降下法
class steepest_decent_method {
  vec prod(const mat& A, const vec& b){
    vec ret(b.size(), 0);
    for(int i=0; i<b.size(); i++){
      for(int j=0; j<b.size(); j++){
	ret[i] += A[i][j] * b[j];
      }
    }
    return ret;
  }
public:
  vec calc(const mat& A, const vec& b, int iter){
    int n = A.size();
    vec x_k(n, 0);
    vec r_k(n, 0);
    for(int k=0; k<iter; k++){
      //rの更新
      vec Ax_k = prod(A, x_k);
      bool flg = true;
      for(int i=0; i<n; i++){
	r_k[i] = b[i] - Ax_k[i];
	if(r_k[i] > EPS) flg = false;
      }
      if(flg) break;

      //alphaの計算
      double r2 = 0, ra2 = 0;
      for(int i=0; i<n; i++) r2 += r_k[i] * r_k[i];
      vec Ar_k = prod(A, r_k);
      for(int i=0; i<n; i++) ra2 += r_k[i] * Ar_k[i];
      if(ra2 < 0) return vec();

      //xの更新
      for(int i=0; i<n; i++){
	x_k[i] += (r2/ra2) * r_k[i];
      }
    }

    return x_k;
  }
};

//共役勾配法
class conjugate_gradient_method {
  vec prod(const mat& A, const vec& b){
    vec ret(b.size(), 0);
    for(int i=0; i<b.size(); i++){
      for(int j=0; j<b.size(); j++){
	ret[i] += A[i][j] * b[j];
      }
    }
    return ret;
  }
public:
  vec calc(const mat& A, const vec& b, int iter){
    int n = A.size();
    vec x_k(n, 0);
    vec r_k1(n, 0), r_k2(n, 0), p_k(n, 0);
    for(int k=0; k<iter; k++){
      //pの更新
      if(k==0){
	vec Ax_k = prod(A, x_k);
	bool flg = true;
	for(int i=0; i<n; i++){
	  r_k1[i] = p_k[i] = b[i] - Ax_k[i];
	  if(r_k1[i] > EPS) flg = false;
	}
	if(flg) break;
      }else{
	double r12 = 0, r22 = 0;
	for(int i=0; i<n; i++) r12 += r_k1[i] * r_k1[i];
	for(int i=0; i<n; i++) r22 += r_k2[i] * r_k2[i];
	for(int i=0; i<n; i++){
	  p_k[i] = r_k1[i] + (r12/r22) * p_k[i];
	}
      }
      //alphaの更新
      double r2 = 0, pa2 = 0;
      for(int i=0; i<n; i++) r2 += r_k1[i] * r_k1[i];
      vec Ap_k = prod(A, p_k);
      for(int i=0; i<n; i++) pa2 += p_k[i] * Ap_k[i];
      if(pa2 < 0) return vec();

      //r,xの更新
      r_k2 = r_k1;

      bool flg = true;
      for(int i=0; i<n; i++){
	x_k[i] += (r2/pa2) * p_k[i];
	r_k1[i] = r_k1[i] - (r2/pa2) * Ap_k[i];
	if(r_k1[i] > EPS) flg = false;
      }
      if(flg) break;
    }

    return x_k;
  }
};

int main(){
  // Ax = b の形式の線形連立一次方程式を解く (ただし、Aは正定値対称行列)
  // 以下、次の方程式を求める。
  // 10x + 2y + 3z = 6
  // 2x + 9y + 5z = 12
  // 3z + 5y + 13z = 21

  mat A(3, vec(3));
  A[0][0] = 10; A[0][1] = 2; A[0][2] = 3;
  A[1][0] = 2; A[1][1] = 9; A[1][2] = 5;
  A[2][0] = 3; A[2][1] = 5; A[2][2] = 13;

  vec b(3);
  b[0] = 6;
  b[1] = 12;
  b[2] = 21;

  //対称性、正定性の確認
  mat_utility mat_util;
  if(!mat_util.is_symmetric(A) || !mat_util.is_positive_definite(A)){
    std::cerr << "matrix A is not symmetric positive definite matrix." << std::endl;
    return 1;
  }

  {//ガウスジョルダン法
    std::cout << "gauss_jordan" << std::endl;
    gauss_jordan gj;
    vec x = gj.calc(A, b);
    if(x.size() > 0){
      std::cout << x[0] << std::endl;
      std::cout << x[1] << std::endl;
      std::cout << x[2] << std::endl;
    }
    std::cout << "==========" << std::endl;
  }

  {//最急降下法
    std::cout << "steepest_decent_method" << std::endl;
    steepest_decent_method sdm;
    vec x = sdm.calc(A, b, 1000);
    if(x.size() > 0){
      std::cout << x[0] << std::endl;
      std::cout << x[1] << std::endl;
      std::cout << x[2] << std::endl;
    }
    std::cout << "==========" << std::endl;
  }

  {//共役勾配法
    std::cout << "conjugate_gradient_method" << std::endl;
    conjugate_gradient_method cgm;
    vec x = cgm.calc(A, b, 1000);
    if(x.size() > 0){
      std::cout << x[0] << std::endl;
      std::cout << x[1] << std::endl;
      std::cout << x[2] << std::endl;
    }
    std::cout << "==========" << std::endl;
  }
  return 0;
}

結果

gauss_jordan
0.0743802
0.545455
1.38843
==========
steepest_decent_method
0.0743802
0.545455
1.38843
==========
conjugate_gradient_method
0.0743802
0.545455
1.38843
==========

みんな解は同じみたい。
最急降下法は25回で収束したけど、共役勾配法は2回で収束してた。

参考文献