java - Java での大きな疎行列のコレスキー分解

Question

Java で大規模な疎行列のコレスキー分解を実行したいと考えています。現在、Parallel Coltライブラリ クラスSparseDoubleCholeskyDecompositionを使用していますが、JNI を使用して Java で使用する高密度行列用に C で記述したコードを使用するよりもはるかに遅いです。

たとえば、SparseDoubleCholeskyDecomposition を使用した 0.25% の非ゼロ密度の 5570x5570 行列の場合、因数分解に26.6 秒かかり、高密度ストレージを使用する同じ行列の自分のコードは1.12 秒しかかかりません。ただし、密度を 0.025% に設定すると、colt ライブラリは 0.13 秒しかかかりません。

また、圧縮された行ストレージと OpenMPを使用して、C で独自の疎行列コレスキー分解を作成しました。これは、SparseDoubleCholeskyDecomposition よりもかなり高速ですが、密なストレージを使用するアルゴリズムよりも遅いです。おそらくさらに最適化することもできますが、いずれにしてもコレスキー因子が密集しているため、速度もストレージの問題も解決されません。

しかし、ミリ秒単位の時間が必要であり、最終的には 10000x10000 を超えるマトリックスにスケーリングする必要があるため、密なマトリックス用の独自の密なコードでさえ遅くなり、メモリを使いすぎます。

疎行列のコレスキー分解は、はるかに高速に実行でき、メモリ使用量も少ないと信じる理由があります。たぶん、より良い Java BLAS ライブラリが必要ですか? 疎行列のコレスキー分解を効率的に行うJava用のライブラリはありますか? Parallel Colt を最適に使用していない可能性があります。

import cern.colt.matrix.tdouble.algo.decomposition.SparseDoubleCholeskyDecomposition;
import cern.colt.matrix.tdouble.impl.*;
import java.util.Random;

public class Cholesky {
    static SparseRCDoubleMatrix2D create_sparse(double[] a, int n, double r) {
        SparseRCDoubleMatrix2D result = new SparseRCDoubleMatrix2D(n, n);
        Random rand = new Random();
        for(int i=0; i<n; i++) {
            for(int j=i; j<n; j++) {
                double c = rand.nextDouble();
                double element = c<r ? rand.nextDouble() : 0;
                a[i*n+j] = element;
                a[j*n+i] = element;
                if (element != 0) {
                    result.setQuick(i, j, element);
                    result.setQuick(j, i, element);
                }
            }
        }
        for(int i=0; i<n; i++) {
            a[i * n + i] += n;
            result.setQuick(i,i, result.getQuick(i,i) + n);
        }
        return result;
    }

    public static void main(String[] args) {
        int n = 5570;
        //int n = 2048;
        double[] a = new double[n*n];
        SparseRCDoubleMatrix2D sparseMatrix = create_sparse(a, n, 0.0025);

        long startTime, endTime, duration;

        startTime = System.nanoTime();
        SparseDoubleCholeskyDecomposition sparseDoubleCholeskyDecomposition = new SparseDoubleCholeskyDecomposition(sparseMatrix, 0);
        endTime = System.nanoTime();
        duration = (endTime - startTime);
        System.out.printf("colt time construct %.2f s\n", 1E-9*duration);
    }
}