17
解決したいAx = bここで、Aは非常に大きな正方形の正定値対称ブロック行列であり、xとbはベクトルです。私が大きく言うと、nxnの行列はnと大きく、300,000と大きくなります。ブロック疎行列で大きな線形システムを解く
私が解決したいと考えているはるかに小さいが代表的な行列の例です。
そして、ここでそれが密集matriciesのブロックで構成されていることを示すにはズームイン同じ行列です。
私は以前（here、hereを参照してください、とまで遡っhereなど）n<10000のためにうまく働いた固有のコレスキーソルバを使用しますがn=300000としてコレスキーソルバが遅すぎます。しかし、私はブロックマトリックスを持っているという事実を利用しませんでした。明らかに、疎ブロック行列（例えば、block cholesky factorization）を解くためのアルゴリズムが存在する。
Eigenが、私が使用できる疎密ブロック行列に対して、分解または反復法を使って最適化されたアルゴリズムを持っているかどうかを具体的に知りたいですか？
また、私のマトリックスを解くのに理想的かもしれない他のアルゴリズムを提案できますか？私は、少なくともパーミュテーション行列を見つけることはNP困難であることを理解する限り、多くの異なる発見的方法が存在し、人々が異なる行列構造（例えばbanded matrix）の直感を構築することができ、それら。私はまだこの直感を持っていない。
私は現在、conjugate gradient methodを使用しています。私はこれをC++で実装しましたが、行列がブロック行列であるという事実を利用していません。
//solve A*rk = xk
//Eigen::SparseMatrix<double> A;
//Eigen::VectorXd rk(n);
Eigen::VectorXd pk = rk;
double rsold = rk.dot(rk);
int maxiter = rk.size();
for (int k = 0; k < maxiter; k++) {
Eigen::VectorXd Ap = A*pk;
double ak = rsold /pk.dot(Ap);
xk += ak*pk, rk += -ak*Ap;
double rsnew = rk.dot(rk);
double xlength = sqrt(xk.dot(xk));
//relaxing tolerance when x is large
if (sqrt(rsnew) < std::max(std::min(tolerance * 10000, tolerance * 10 * xlength), tolerance)) {
rsold = rsnew;
break;
}
double bk = rsnew/rsold;
pk *= bk;
pk += rk;
rsold = rsnew;
}
本当にこれを実装しますか？おそらく、あなたはこの種のものをやっている利用可能なライブラリを試してみるべきでしょうか？例えば。 [elemental]（http://libelemental.org/）。 – sascha
@サシャ、私はライブラリを使用して幸せですが、私はすでにEigen（カスタム共役勾配法を除いて）でそれをやっていますが、Eigenがこの場合十分でない場合は他のものを検討します。 –
固有音は、私にとって基本的な線形代数のように聞こえますが、（より複雑な構造を使用する）最適化アルゴリズム自体はありません。しかし、私はそれを使ったことはありません。私も* elemental *での経験は全くありませんが、あなたがやっていることは何とかできていると思います。悲しいことに、あなたに与えることのできるヒントはもうありません。おそらく、このドキュメントの[this]（http://libelemental.org/documentation/dev/optimization/models.html）の部分を見てください。がんばろう！ – sascha