目次


C、C++ 用行列ライブラリー

Meschach、Cooperware matrix、Blitzの評価と比較

Comments

以下の記事は、C/C++ をある程度理解し、C/C++自体 に行列処理の機能がないことを知っているものとして書かれています。本稿の読者は、計量経済学のデータを分析する方かもしれませんし、熱帯雨林のモデリングを行っている方かもしれません。私自身は、ニューラル・ネットの研究を行っており、行列を1個か2個使うことで、ニューラル・ネットを非常に簡単に実装することができます。C/C++ には行列に類するコンテナーがいろいろと用意されていますが (配列の他、標準ライブラリーのベクトル、リスト、マップなど)、行列そのものを表すコンテナーがあると、いま現実に直面している作業が極めて簡単に処理できるようになります。というわけで、行列をいちから構築するのではなく、行列ライブラリーに少し手を加えるだけで済むような3つのオープン・ソースの製品について検討してみたいと思います。とくに、これまでとは違うやり方でライブラリーを応用したいと考えている方にとって、この研究は役に立つはずです。

Meschach: Cでの選択肢

Cでコーディングされるプロジェクト用としては、Meschach (ミーシャークと読む) が行列やベクトルの演算を行うルーチンを提供します。Linuxを始め、他のほとんどのオペレーティング・システムでコンパイルできるという点が特長で、慣例的な謝辞 (customary acknowledgment) を提示し、エラーを報告するという条件の下、著作権のある形で、オープンに利用することができます。Meschachは、とくに、密または疎な連立一次方程式 (systems of dense or sparse linear equations) を解いたり、固有値や固有ベクトルを計算したり、最小二乗 (least squares) 問題を解くために設計されたものです。倍精度の値や複素数を扱うことのできる約400個の関数が用意されています。また、説明用の小規模なケース・スタディーに沿って、これらの関数の使い方を紹介するチュートリアルも付属しています。David StewartとZbigniew Leykが、過剰決定された方程式 (over-determined equations) に対する一般化された最小二乗方程式の解法 (least square equation solver) とか、疎なマトリックスが関係する問題などの話題を用いて、Meschachを紹介しています。さらに、このチュートリアルでは、三次元行列とかエラー報告などのもう少し高度な話題も扱っています。

オブジェクトやクラスの機能は、しばしば、コードと関係付けられ、Cの構造体は、少しわかり難いと思われる場合もあるため、Cのライブラリーは、将来性のあるソリューションかどうかという点で、採用されないことがよくあります。しかし弁護するとすれば、このライブラリーは非常にうまく構成されており、最初から不採用と決めてかかることはありません。私の場合、Meschachをダウンロードして15分で、行列を作成し、データを埋め込み、表示を行うことができました (Hello World! プログラムを作成することと同等のこと)。参考文献に、"Meschach: Matrix Computations in C" という、手頃な値段の製本されたマニュアルもあります (稿末の参考文献参照)。とくに、tortureというテスト・プログラムには、参考になるヒントが数多く含まれています。

行列は、簡単にファイルや標準出力に送ることができます。Meschachは、高速フーリエ変換 (Fast Fourier Transforms) を計算したり、列や行を抽出したり、対称行列の固有値を計算することができます。行列にランダムな整数や複素数を埋め込むこともできます。驚くことに、このライブラリーには、行列を加算する機能まで用意されています。Meschachは、足がかりとなるプログラムを簡単に記述できるようにする機能として、[0,1) の範囲内でランダムな倍精度の値を返す関数も備えています。Meschachには、行列を短精度の1.0で埋める関数はありますが、残念ながら、任意の倍精度値やランダムな倍精度値で行列を埋めるための関数はありません。しかし、そうした関数を追加するのは簡単です。

リスト1. 任意の倍精度値またはランダムな倍精度値で行列を埋める方法
MAT *m_fill( MAT *A, double x)
/* MAT *m_random_fill( MAT *A ) */
{
   int i, j;
   for ( i = 0; i < A->m; i++ ) for ( j = 0; j < A->n; j++ )
      { A->me[i][j] = x ; }
      /* { A->me[i][j] = m_random ; } */
   return A;
}

上のリストからわかるように、Meschachのコードを拡張するのは簡単です。理想を言えば、コードにもっと懇切丁寧にコメントを入れてあれば、もっと良いのですが。とはいっても、Cで行列を必要とする計算処理を行うことになったら、このライブラリーは非常に便利です。

Cooperware matrix: 基本的な機能

オブジェクト指向のコードをC++ で記述したいと考え、速度よりも概念的な明晰さを重視したいという方は、Harry KuiperのCooperware Matrix (CwMtx) が非常に役に立つことでしょう。本稿で取り上げた3つの行列ライブラリーの中では、このライブラリーの概念的アーキテクチャーが最も直観的に分かりやすいものでした。行列の作成には、以下のような分かりやすい記述が利用できます。

リスト2. 行列の作成
CWMatrix A( rows, columns ) ;

今回取り上げた3つのライブラリーの中では、詳細については以下の速度の項目で説明しますが、3つの評価用のタスクに関するかぎり、CwMtxが成績がよくありませんでした。しかし、性能よりも明快さのほうが重要な場合 (たとえば、データが正しく計算されているかどうかを確認したいというような場合) には、CwMtxが好都合です。まずは正しく、次に速くです。

CwMtxの行列には、ベクトルや正方行列も含まれます。ベクトルには空間ベクトルや四元数 (quaternion) も含まれます。行列を別の行列に写像 (map) したり、行列にいくつかの要素を埋め込んだり、行列を転置したり、行列に通常の数学的な演算を施すといったことができます。Kuiperは、もともと、分散した対話式のステート・マシン群を基に構築されるシステムをシミュレートするために、CwMtxを使用しました。当然必要な行列クラスに加えて、四元数クラスも存在します。四元数とは何かというと、q = r + xi + yi + ziで、rが実数、iが -1の平方根、x、y、zが複素数のときのqのことです。四元数を導入すると、三次元での回転の概念を四次元に拡張することができます (四元数についての参考文献へのリンクは、参考文献を参照してください)。

CwMtxには、ランダム数のジェネレーターは用意されていませんし、行列をランダムな要素で埋めるクラス関数もありませんが、GNU LGPLライセンスの下、無料で配布されているライブラリーですので、必要なら、そういう関数を作成してもかまいません。行列をランダムな要素で埋めることに関して言えば、以下の方法が便利であり、簡単です。

リスト3. 行列をランダムな要素で埋める方法
#include <stdlib.h>
...
void random_fill( CWMatrix &M )
{
 int SIZE = M.GetRows() ;
 for ( int r=0; r<SIZE; r++ ) for ( int c=0; c<SIZE; c++ ) { M[r][c]= drand48(); }
}

ドキュメントは、わずかですが、明快でよく構成されています。クラス階層、コンストラクター、メンバー関数のオプションなど、簡単に見つけ出すことができます。チュートリアルはありませんが、それが必要だとは感じないでしょう。ドキュメントとテスト・プログラムの間に、そういうものは必要ないでしょう。

Blitz: Fortran並みの速さ ?

Blitzも、GNU GPLの下で配布されるC++ ライブラリーで、これを拡張して自由にオブジェクトを作成することができます。KAI、Intel、gcc、Metroworks、Cray 3.0.0.0といったC++ コンパイラーをサポートし、整数、浮動小数点数、複素数、およびwell-behavedユーザー定義の型を要素としてとることのできるn次元の配列クラスを用意しています。コンストラクターは、以下の例からわかるように、CwMtxのコンストラクターよりも複雑です。

リスト4. Blitzのコンストラクター
Array<double,2> A(4,7) ;

これによって、倍精度の値からなる4行7列、階数 (rank) 2の配列が作成されます。しかし、これでは少しわかりにくいので、Blitzは、1個の配列を行列とみなすようにさせます。また、行列の関数の数多くを実装していません。たとえば、行列の固有値を返す関数は用意されていません。配列や行列をランダムな倍精度の値で埋める関数もありません。しかし、Blitzには、基本的な特長が2つあります。

特長の1つは、その間口の広さにあります。以下の例からわかるように、Blitzがもともと備えている機能を利用することで、ランダムな倍精度値を埋め込む関数を簡単に実装し、構築できるという点です。

リスト5. Blitzでランダムな倍精度値を埋め込むための関数
template <class Array, class Uniform>
void fill( Array &a, int size, Uniform &u )
{
   for ( int i=0; i<size; i++ )
      for ( int j=0; j<size; j++ )
         { a(i,j) = u.random() ; }
}

BlitzのUniformクラスには、[0,1) の範囲内でランダムな倍精度の値を返すメンバー関数が用意されています。また、配列の要素にアクセスするためのメソッドが3種類用意されています。すなわち、標準的な添え字によるもの、サブアレイ (subarray) を作成するもの、および配列セグメントを次元の低い形にして提示するスライシング (slicing) の3種類のメソッドです。さらにBlitzには、標準的な電卓型の関数も用意されており、配列を標準出力に表示したり、ファイルとの間で入出力を行うこともできます。ラプラシアン、グラジエント、ヤコビアンといった演算子は、Blitzの組み込み関数の3つの代表例です。

もう1つの特長は、Blitzの速度にあります。使用するコンパイラーにもよりますが、科学技術計算に関して高い性能を示すことで有名なFortranとほぼ同等の性能をC++ で引き出すことができます。以下の表に性能比を示しますが、このデータの分析およびそれに基づいて示唆される性能については速度 の項目を参照してください。

表1. いろいろなプラットフォームでのBlitzの性能
プラットフォームコンパイラーキャッシュ外キャッシュ内
HPC-160KAI C++100.2%97.5%
Pentium IIegcs98.4%79.6%
Cray T3EKAI C++95.7%98.1%
Origin 2000KAI C++88.1%79.8%

Blitzには、HTMLとPostscriptのマニュアルが添付されていますが、残念ながらチュートリアルはありません。ただし、説明用のコードがかなり用意されていて、そこからBlitzの構文の微妙なところを学びとることができます。また、クラス・リファレンスに、通常の情報が収められています。役に立つメーリング・リストも多数あり、アーカイブされていて、検索できるようになっています (参考文献参照)。

速度

ライブラリーは、その関数リソース、ドキュメント、チュートリアルの質、拡張の容易さなどによって評価することができます。また、性能や速度に関しても評価することができます。しかし、(われわれの場合がそうであるように) どれも同じ言語で記述されているわけではありませんし、それぞれの本来的な機能があまり共通していないため、ライブラリーを比較することが難しい場合もあります。われわれの場合、ライブラリー間に重複する部分がかなりありますので、比較的単純ながら差異を明らかにしてくれる3つのタスクに基づいて比較することが可能です。以下の3つの例で、それを紹介し、説明します。

リスト6. タスク1
For ( d=2; d<7; d++)
   Construct 3 dxd matrices: A, B, C
   Start Clock: Do the following one million times:
Fill A and B with 1.0s
      Let C = A + B
   Stop clock: Report elapsed time in seconds.

今回取り上げたライブラリーを使って、このアルゴリズムを実装し実行したところ、以下の結果が得られました。

表2. タスク1の結果
ライブラリー2x23x34x55x56x6
Blitz0.400.480.620.750.91
CwMtx2.643.574.585.606.60
Meschach0.170.270.450.600.79

CwMtxは、直観的に分かりやすいアーキテクチャーなのですが、残念ながら、この場合、性能はよくありません。またBlitzは、Meschachよりも劣っていますが、同じオブジェクト指向の競合製品であるCwMtxよりははるかに高い性能を示している点は注目に値します。

MeschachとBlitzにはランダムな倍精度値を供給する関数 (Random Number Generators: ランダム値ジェネレーター) が用意されているのに対して、CwMtxには独自のランダム生成機能がないことを思い出してください。行列をベースにして行われるシミュレーションの中には、ランダムな値を使うことが重要な役割を果たすものもありますので、ランダム値生成を行う必要のある状況で、これらのライブラリーがどんな性能を示すのかを調べることは有益なことです。

リスト7. タスク2
For ( d=2; d<7; d++)
   Construct 3 dxd matrices: A, B, C
   Start Clock: Do the following one million times:
Fill A and B with random doubles (using library RNG, if any)
      Let C = A + B
   Stop clock: Report elapsed time in seconds.

今回取り上げたライブラリーの性能は、以下のとおりでした。

表3. タスク2の結果
ライブラリー2x23x34x55x56x6
Blitz0.871.712.834.346.13
CwMtx3.675.878.5911.9315.61
Meschach0.420.801.321.862.52

Blitzの性能は、ここでもMeschachより劣っていますが、オブジェクト指向の競合製品CwMtxには大きな差をつけています。これはRNG (ランダム値ジェネレーター) の性能に因るものであると思われるといけませんので、3つ目の評価タスクを見てみることにしましょう。

リスト8. タスク3
For ( d=2; d<7; d++)
   Construct 3 dxd matrices: A, B, C
   Start Clock: Do the following one million times:
Fill A and B with random doubles (using shared RNG)
      Let C = A + B
   Stop clock: Report elapsed time in seconds.

以下のように、性能の順位は、前の2つのタスクの場合と変わりません。

表4. タスク3の結果
ライブラリー2x23x34x55x56x6
Blitz1.312.624.506.859.71
CwMtx3.675.878.5911.9315.61
Meschach1.172.454.286.639.40

予想どおりだと思いますが、CwMtxの性能の順位は変わっていません。また、BlitzとMeschachの相対的な順位も変わりはありません。生の速度が決定的な意味をもつのであれば、これら3つのライブラリーの間の順位は明らかとなりました。

Red Hat Linux 7.1でのインストールとコンパイル

これら3つのライブラリーをインストールおよびコンパイルする際の便を図るために、以下のメモを付記しておきます。ダウンロードのリンク先は、参考文献に示してあります。

リスト9. Blitz
tar zxf blitz-0.5beta3.tar.gz
cd blitz-0.5beta3
./configure --with-cxx=gcc
make all
cp -a blitz/ /usr/local/include/( or whatever you wish )
cp -a lib/ /usr/local/include/blitz/
Compile with: g++ -O2 pgm.cpp -o pgm
リスト10. Meschach
unzip -q mesch12b.zip -d mesch12b
cd mesch12b
./configure
make basic
mkdir /usr/local/include/meschach( or whatever you wish ) 
cp *.h meschach.a /usr/local/include/meschach/
Compile with: gcc -O2 pgm.c /usr/local/include/meschach/meschach.a -o pgm
リスト11. CwMtx
tar zxf cwmtx-0.3.0.tar.gz
cd cwmtx-0.3.0
Open Makefile, and set INSTALL_DIRto /usr/local/include/cwmtx( or whatever you wish )
Execute: make install
Compile with: g++ -O2 pgm.cpp -o pgm

まとめ

本稿では、3つの行列ライブラリーの特徴を紹介し、それぞれのライブラリーが本来備えている機能を詳しく吟味しました。また、これらのライブラリーが抱えている機能的な短所を確認するとともに、それらの短所を克服する方法を研究しました。本稿では、いくつか単純なテストを紹介し、生の数量的データを提示しました。これらのデータは、みなさんがコーディングを選択するときの手助けとなるかもしれませんが、最終的な決定はみなさんに委ねられており、個々の状況におけるライブラリーの速度だけでなく、みなさんのプロジェクトのいろいろな側面にしたがって決定すべきものです。


ダウンロード可能なリソース


関連トピック


コメント

コメントを登録するにはサインインあるいは登録してください。

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=60
Zone=Linux
ArticleID=226820
ArticleTitle=C、C++ 用行列ライブラリー
publish-date=07012002