行列演算 - mist-team/mist GitHub Wiki

はじめに

MISTでは数値計算を行うための行列演算クラスを提供しています．

ヘッダの準備

行列演算を行う場合は，次のヘッダを必要なところでインクルードしてください．

#include <mist/matrix.h>

このままでは，行列の四則演算程度しか行うことができません． 行列式，LU分解，固有値計算等を行う場合には matrix.h とあわせて次のヘッダもインクルードしてください．

#include <mist/numeric.h>

基本的なこと

行列を格納するオブジェクトは次のように用意します.

mist::matrix< float >  m;              // 0 行 0 列の行列を確保.
mist::matrix< double > m( 4, 4 );      // 4 行 4 列の行列を確保, 要素は 0.0 で埋められる.
mist::matrix< double > m( 4, 4, 1.0 ); // 要素は 1.0 で埋められる.

また，複素数を要素とする行列を確保する場合は次のように書きます．

mist::matrix< std::complex< float > >  m;              // 0 行 0 列の行列を確保.
mist::matrix< std::complex< double > > m( 4, 4 );      // 4 行 4 列の行列を確保, 要素は ( 0.0, 0.0 ) で埋められる.
mist::matrix< std::complex< double > > m( 4, 4, 1.0 ); // 要素は ( 1.0, 0.0 ) で埋められる.

その際に，次のヘッダもインクルードしてあげてください．

#include <complex>

例えば次のように行列の要素にアクセスします.

mist::matrix< double >::size_type r, c;
mist::matrix< double > ma( 2, 3, 4.0 ), mb( 2, 3 );

for( c = 0 ; c < ma.cols( ) ; c++ )
{
  for( r = 0 ; r < ma.rows( ) ; r++ )
  {
    mb( r, c ) = ma( r, c ) * 2.0;
  }
}

行列とスカラーの乗算も用意されているので, 次のようにも書けます.

mb = ma * 2.0;

行列同士の演算も可能です.

mist::matrix< double > A( 2, 2 ), x( 2, 1 ), b;
A( 0, 0 ) =  1.0; A( 0, 1 ) = 2.0;
A( 1, 0 ) = -3.0; A( 1, 1 ) = 1.0;
x( 0, 0 ) =  2.0;
x( 1, 0 ) =  7.0;

b = A * x; // b の必要とする領域は自動的に割り当てられます.

numeric.h

上で書いたように, 行列式, LU分解, 固有値計算等を行う場合には numeric.h をインクルードします. さらに CLAPACK という行列計算ライブラリが必要になるので, ここから下を読み進める前に, 次のチュートリアルを参考にして予めライブラリを準備しておいてください. LU 分解を用いた逆行列の計算は次のように書けます.

mist::matrix< double > m;
mist::matrix< double > m_inv;
/*
 (逆行列を持つ) m を設定. 
*/
m_inv = mist::inverse( m );

ここで, 与える行列がある条件を満たしている場合には, より高速・高精度に計算が可能です. 行列が対称行列であるなら,

m_inv = mist::inverse( m, mist::matrix_style::sy ); // m が対称行列

と書けます.

/* ... */
  ge, // 一般行列
  gb, // 一般帯行列
  gt, // 一般3重対角行列
  sy, // 対称行列
  sb, // 対称帯行列
  st, // 対称3重対角行列
  he, // エルミート行列
  hb, // エルミート帯行列
  ht, // エルミート3重対角行列
/* ... */

となっています. その他の数値計算に関してはドキュメントを参照してください. 基本的には同じです.

サンプルプログラム

以下は主成分分析を行なう簡単なプログラムです. プログラムを入力する前に, 固有値分解と特異値分解での固有値の並び順を同じにするために, mist_conf.h のオプションを変更しておきます.

#define _DESCENDING_ORDER_EIGEN_VALUE_          0

の行を探します. "固有値・固有ベクトルを計算した時に, 降順に並べる時は 1, 昇順に並べる時は 0 にする" とのことなので, ここは特異値分解での並びに合わせて, 降順に並べるように変更します. 解くべき問題によっては, そのままの方が都合の良いことも多いのでここは好みの問題になります. 今回は主成分分析を行ないたいので, 1 を設定しておきましょう.

#define _DESCENDING_ORDER_EIGEN_VALUE_          1

#include <iostream>

// MIST の一番メインとなるヘッダファイル
#include <mist/mist.h>

// 行列演算
#include <mist/matrix.h>
#include <mist/numeric.h>

// 擬似乱数生成
#include <mist/random.h>
#include <time.h>

//
// 主成分分析サンプルプログラム
//
int main( int argc, char *argv[] )
{
  mist::matrix< double >::size_type r, c;

  //
  // データの作成
  //
  mist::matrix< double > x( 1000, 2 );
  mist::gauss::random rand( time( NULL ) );
  for( r = 0 ; r < x.rows( ) ; r++ )
  {
    x( r, 0 ) = rand.generate( ) * 10.0 + 175.0;                       // 身長
    x( r, 1 ) = ( x( r, 0 ) - 100.0 ) * 0.9 + rand.generate( ) * 10.0; // 体重
  }

  //
  // データの正規化
  //
  mist::matrix< double > x_( x.rows( ), x.cols( ) );
  for( c = 0 ; c < x.cols( ) ; c++ )
  {
    double ave = 0.0;
    for( r = 0 ; r < x.rows( ) ; r++ )
    {
      ave += x( r, c );
    }
    ave /= x.rows( );

    double dev = 0.0;
    for( r = 0 ; r < x.rows( ) ; r++ )
    {
      double s = x( r, c ) - ave;
      dev += s * s;
    }
    dev = sqrt( dev / x.rows( ) );

    for( r = 0 ; r < x.rows( ) ; r++ )
    {
      x_( r, c ) = ( x( r, c ) - ave ) / dev;
    }
  }

  //
  // 共分散行列の作成 (データを正規化しているので相関行列に等しい)
  //
  mist::matrix< double > v;
  v = ( x_.t( ) * x_ ) / ( x.rows( ) - 1.0 );

  //
  // 固有値計算
  //
  mist::matrix< double > eval, evec;
  mist::eigen( v, eval, evec, mist::matrix_style::sy ); // 共分散行列は対称行列

  // 行列の << による出力をサポート
  std::cout << eval << std::endl; // mist_conf.h 設定の結果, 固有値の大きいほうから並ぶ
  std::cout << evec << std::endl; // mist_conf.h 設定の結果, 固有値の大きいほうから並ぶ
                                      // Row 方向に固有ベクトルの要素が並ぶ
                                      // Column 方向に固有ベクトルが並ぶ

  //
  // 主成分
  //
  mist::matrix< double > z;
  z = x_ * evec;

  //
  // 特異値分解 (おまけ)
  //
  mist::matrix< double > u, s, vt;
  mist::svd( x_, u, s, vt, mist::matrix_style::ge );

  std::cout << vt << std::endl; // 固有値の大きいほうから並ぶ
                                     // Column 方向に固有ベクトルの要素が並ぶ
                                     // Row 方向に固有ベクトルが並ぶ
                                    // ただし, 固有値計算の固有値ベクトルと符号が
                                    // 反転していることもあります

  std::cout << vt.t( ) << std::endl; // vt.t( ) にすると, 固有値計算と同じく並ぶ

  // evec と vt.t( ) で同じ値が出力されているのを確認してみてください.

  return( 0 );
}
⚠️ **GitHub.com Fallback** ⚠️