Ax=b 線形代数システム向けの C++ メモリ効率の高いソリューション

https://stackoverflow.com/questions/1242190

12-09-2019
|

質問

単純な線形システムを解くために、Boost UBlas の数値ライブラリバインディングを使用しています。比較的小さな「M」のマトリックスa（MXM）の処理に限定されていることを除いて、次のものは正常に機能します。

実際には、次元 m= 10^6 (最大 10^7) のはるかに大きな行列があります。
メモリを効率的に使用して Ax=b を解決するための既存の C++ アプローチはありますか。

#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/io.hpp>
#include<boost/numeric/bindings/traits/ublas_matrix.hpp>
#include<boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_vector2.hpp>

// compileable with this command


//g++ -I/home/foolb/.boost/include/boost-1_38 -I/home/foolb/.boostnumbind/include/boost-numeric-bindings solve_Axb_byhand.cc -o solve_Axb_byhand -llapack


namespace ublas = boost::numeric::ublas;
namespace lapack= boost::numeric::bindings::lapack;


int main()
{
    ublas::matrix<float,ublas::column_major> A(3,3);
    ublas::vector<float> b(3);


    for(unsigned i=0;i < A.size1();i++)
        for(unsigned j =0;j < A.size2();j++)
        {
            std::cout << "enter element "<<i << j << std::endl;
            std::cin >> A(i,j);
        }

    std::cout << A << std::endl;

    b(0) = 21; b(1) = 1; b(2) = 17;

    lapack::gesv(A,b);

    std::cout << b << std::endl;


    return 0;
}

解決

短い答え：ブーストは使わないでください LAPACK バインディング、これらはスパースマトリックスではなく、密なマトリックス用に設計されています。 UMFPACK その代わり。

長い答え: UMFPACK は、A が大きくて疎な場合に Ax=b を解くための最適なライブラリの 1 つです。

以下はサンプルコードです（ベース） umfpack_simple.c) 単純な A そして b そして解決します Ax = b.

#include <stdlib.h>
#include <stdio.h>
#include "umfpack.h"

int    *Ap; 
int    *Ai;
double *Ax; 
double *b; 
double *x; 

/* Generates a sparse matrix problem: 
   A is n x n tridiagonal matrix
   A(i,i-1) = -1;
   A(i,i) = 3; 
   A(i,i+1) = -1; 
*/
void generate_sparse_matrix_problem(int n){
  int i;  /* row index */ 
  int nz; /* nonzero index */
  int nnz = 2 + 3*(n-2) + 2; /* number of nonzeros*/
  int *Ti; /* row indices */ 
  int *Tj; /* col indices */ 
  double *Tx; /* values */ 

  /* Allocate memory for triplet form */
  Ti = malloc(sizeof(int)*nnz);
  Tj = malloc(sizeof(int)*nnz);
  Tx = malloc(sizeof(double)*nnz);

  /* Allocate memory for compressed sparse column form */
  Ap = malloc(sizeof(int)*(n+1));
  Ai = malloc(sizeof(int)*nnz);
  Ax = malloc(sizeof(double)*nnz);

  /* Allocate memory for rhs and solution vector */
  x = malloc(sizeof(double)*n);
  b = malloc(sizeof(double)*n);

  /* Construct the matrix A*/
  nz = 0;
  for (i = 0; i < n; i++){
    if (i > 0){
      Ti[nz] = i;
      Tj[nz] = i-1;
      Tx[nz] = -1;
      nz++;
    }

    Ti[nz] = i;
    Tj[nz] = i;
    Tx[nz] = 3;
    nz++;

    if (i < n-1){
      Ti[nz] = i;
      Tj[nz] = i+1;
      Tx[nz] = -1;
      nz++;
    }
    b[i] = 0;
  }
  b[0] = 21; b[1] = 1; b[2] = 17;
  /* Convert Triplet to Compressed Sparse Column format */
  (void) umfpack_di_triplet_to_col(n,n,nnz,Ti,Tj,Tx,Ap,Ai,Ax,NULL);

  /* free triplet format */ 
  free(Ti); free(Tj); free(Tx);
}


int main (void)
{
    double *null = (double *) NULL ;
    int i, n;
    void *Symbolic, *Numeric ;
    n = 500000;
    generate_sparse_matrix_problem(n);
    (void) umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    (void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null);
    umfpack_di_free_symbolic (&Symbolic);
    (void) umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null);
    umfpack_di_free_numeric (&Numeric);
    for (i = 0 ; i < 10 ; i++) printf ("x [%d] = %g\n", i, x [i]);
    free(b); free(x); free(Ax); free(Ai); free(Ap);
    return (0);
}

関数 generate_sparse_matrix_problem マトリックスを作成します A そして右側 b. 。マトリックスは最初にトリプレット形式で構築されます。ベクトルTI、TJ、およびTXはAを完全に説明していますトリプレットフォームは簡単に作成できますが、効率的なスパースマトリックスメソッドには、圧縮されたスパースカラム形式が必要です。変換はで実行されます umfpack_di_triplet_to_col.

記号因数分解は次のように実行されます。 umfpack_di_symbolic. 。のまばらなlu分解 A で実行されます umfpack_di_numeric. 。下三角および上三角の解決は次のように実行されます。 umfpack_di_solve.

と n 私のマシンでは 500,000 なので、プログラム全体の実行に約 1 秒かかります。Valgrind は、369,239,649 バイト (352 MB 強) が割り当てられたと報告しています。

これに注意してくださいページトリプレット（座標）および圧縮形式のスパースマトリックスに対するBoostのサポートについて説明します。必要に応じて、これらのブーストオブジェクトを単純な配列に変換するルーチンを作成できます UMFPACK 入力として必要です。

他のヒント

あなたの巨大な行列は、私は、彼らはそのサイズであることを望む PARDISO <を見ている、スパースあると仮定すると、 />スパース線形ソルバーであるプロジェクトは、これはあなたが言ったように大きな行列を扱うようにしたい場合はあなたが必要とするものです。唯一の非ゼロ値の効率的な記憶を可能にし、密行列の同じシステムを解くことよりもはるかに高速である。

マトリックスが密であると仮定します。スパースであれば、すでに述べたように、多数の特殊なアルゴリズムを見つけることができます。デウスアデューロそしてダフィーモ.

自由に使える (十分な大きさの) クラスターがない場合は、アウトオブコアアルゴリズムを検討する必要があります。スカラパックの一部としてコア外のソルバーがいくつかあります。プロトタイプパッケージ, 、ドキュメントを参照してくださいここそしてグーグル詳細については。Web で「アウトオブコア LU / (マトリックス) ソルバー / パッケージ」を検索すると、さらに豊富なアルゴリズムやツールへのリンクが表示されます。私はそれらの専門家ではありません。

ただし、この問題に対して、ほとんどの人はクラスターを使用するでしょう。ほぼすべてのクラスターで見つかるパッケージは、やはり ScaLAPACK です。さらに、通常、典型的なクラスターには他にも多数のパッケージがあるため、問題に合ったものを選択できます (例ここそしてここ).

コーディングを開始する前に、問題を解決するのにどれくらい時間がかかるかを簡単に確認したいと思うでしょう。一般的なソルバーは約 O(3*N^3) フロップを要します (N は行列の次元です)。N = 100000 の場合、3000000 Gflops になることになります。インメモリソルバーがコアあたり 10 Gflops/秒を実行すると仮定すると、単一コアで 3 日半かかることになります。アルゴリズムは適切に拡張できるため、コアの数を増やすと時間はほぼ直線的に短縮されるはずです。その上に I/O が追加されます。

C++ の実装についてはわかりませんが、メモリが問題になる場合は、扱っている行列の種類に応じていくつかの対処法があります。

行列がスパースまたはバンド状である場合は、スパースまたはバンド幅ソルバーを使用できます。これらはバンドの外側にゼロ要素を保存しません。
波面ソルバーを使用すると、行列をディスクに保存し、分解のために行列波面のみを取り込むことができます。
行列を完全に解くことを避け、反復法を使用することもできます。
モンテカルロ法による解法を試すことができます。

解決のために自由に利用可能なソフトウェアのリストを見てください線形代数の問題のジャック・ドンガラとハテムLtaiefでコンパイル、します。

私は、あなたが見ている問題の大きさのために、あなたはおそらく、反復アルゴリズムが必要だと思います。あなたはスパース形式で行列Aを保存したくない場合は、マトリックスフリーの実装を使用することができます。通常、行列Aの個々のエントリにアクセスする必要はありません反復アルゴリズムは、彼らが唯一の行列 - ベクトル積のAV（そして時にはA ^ T Vをベクトルと転置行列の積）を計算する必要があります。ライブラリはうまく設計されている場合、あなたはそれを行列 - ベクトル積を行う方法を知っているクラスを渡すのであれば、それは十分なはずです。

受け入れ答えが示すように、

UMFPACKがあります。あなたはBOOSTを使用している場合しかし、あなたはまだBOOSTでコンパクトな行列を使用してシステムを解決するためにUMFPACKを使用することができます。それが容易になりますどのバインディングがあります：

http://mathema.tician.de/software/boost-numeric-bindings の

その約2年後には日付が、そのわずか（数他の人と一緒に）結合ます。

関連する質問を参照してください。

の UMFPACKとBOOSTのuBLASスパース行列

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow