LU分解

LU分解 (LU decomposition) の用途： 連立方程式を解くなど．正方行列でなくともよい．
行列 (Matrix/ComplexMatrix) A に対するLU分解は

  LU  fact (const Matrix& A);
  ComplexLU  fact (const ComplexMatrix& A);

によって得られる． octave では，行列Aは

  PA = LU

のように分解される．ここで L, U はLU分解された行列， P は置換行列 (Permutation matrix) である． P は直行行列なので，

  A = P'LU

のようにして A を復元できる（ P' は P の転置行列）．

余談だが， octave で

  [l,u]=lu(A)

とした場合には， l=P'L, u=U が代入され，

  [l,u,p]=lu(A)

とした場合には， l=L, u=U, p=P が代入される．

ここで LU クラス, ComplexLU クラス のオブジェクトにはLU分解の結果が代入される．

LU オブジェクトは以下の public メンバ関数を持つ:

  // L を返す(Matrix) :
  Matrix  L (void) const;
  // U を返す(Matrix) :
  Matrix  U (void) const;
  // P を返す(Matrix) :
  Matrix  P (void) const;

サンプル:

  Matrix A(3,3);
  stringstream ss (
    "  3 1 2"
    "  5 0 0"
    "  4 8 1"); ss >> A;
  cout<<"A="<<endl<<A;
  // LU オブジェクトの生成 :
  LU  fact (A);
  cout<<"fact.L()= "<<endl<< fact.L();
  cout<<"fact.U()= "<<endl<< fact.U();
  cout<<"fact.P()= "<<endl<< fact.P();
  // P'LU を計算（Aに一致するか?） :
  cout<<"fact.P().transpose()*fact.L()*fact.U()= "<<endl
    << fact.P().transpose()*fact.L()*fact.U();
  // 行列Aを変更
  A.resize(3,6);
  ss.clear(); ss.str(
    "  4   0   4   6   5   4"
    "  9   6   2   4   1   0"
    "  1   6   9   0   7   5"); ss >> A;
  cout<<"A="<<endl<<A;
  // 再度LU分解 :
  fact= LU(A);
  cout<<"fact.L()= "<<endl<< fact.L();
  cout<<"fact.U()= "<<endl<< fact.U();
  cout<<"fact.P()= "<<endl<< fact.P();
  // P'LU を計算（Aに一致するか?） :
  cout<<"fact.P().transpose()*fact.L()*fact.U()= "<<endl
    << fact.P().transpose()*fact.L()*fact.U();

結果:

  A=
    3 1 2
    5 0 0
    4 8 1
  fact.L()=
    1 0 0
    0.8 1 0
    0.6 0.125 1
  fact.U()=
    5 0 0
    0 8 1
    0 0 1.875
  fact.P()=
    0 1 0
    0 0 1
    1 0 0
  fact.P().transpose()*fact.L()*fact.U()=
    3 1 2
    5 0 0
    4 8 1
  A=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5
  fact.L()=
    1 0 0
    0.111111 1 0
    0.444444 -0.5 1
  fact.U()=
    9 6 2 4 1 0
    0 5.33333 8.77778 -0.444444 6.88889 5
    0 0 7.5 4 8 6.5
  fact.P()=
    0 1 0
    0 0 1
    1 0 0
  fact.P().transpose()*fact.L()*fact.U()=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5