Java中的CHOLMOD

Question

我已经问了类似的东西,但这次我会更具体.

我需要在for循环内执行通常大的正定对称矩阵(约1000x1000)的Cholesky分解.现在,要做到这一点,我一直试图:

1)Apache Math库

2)并行Colt库

3)JLapack库

在上述三种情况中的任何一种情况下,例如,与MATLAB相比,时间消耗非常长.

因此,我想知道在Java中是否存在用于Cholesky分解的高度优化的外部工具:例如,我一直在考虑CHOLMOD算法,其实际上是内部调用的MATLAB和其他工具.

我真的很感激能够对此事进行全面的反馈.

Answer 1

以下是Java的一些BLAS库的一个很好的总结:性能的java-matrix-math-libraries.您还可以在Java-Matrix-Benchmark上看到许多这些库中的许多库的基准.

但是,根据我的经验,大多数这些库似乎没有针对解决大型稀疏矩阵进行调整.在我的例子中,我所做的是通过JNI 实现使用Eigen的求解.

Eigen对其线性求解器进行了很好的讨论,包括一个关于CHOLMOD的线性求解器.

对于8860x8860的情况,使用Eigen解算器通过JNI的稀疏矩阵比平行小马快20倍,比我自己的密集求解器快10倍.更重要的是,它看起来像是缩放n^2而不是n^3它使用的内存比我的密集解算器少得多(我的内存扩展耗尽).

实际上有一个名为JEigen的 Eigen包装器,它使用JNI.但是,它没有实现稀疏矩阵求解,因此它不会包装所有内容.

我最初使用JNA但对开销不满意.维基百科就如何使用JNI有一个很好的例子.一旦编写了函数声明并用它们编译它们就javac可以用来javah为C++创建头文件.

例如

//Cholesky.java
package cfd.optimisation;
//ri, ci, v : matrix row indices, column indices, and values
//y = Ax where A is a nxn matrix with nnz non-zero values
public class Cholesky {
    private static native void solve_eigenLDLTx(int[] ri, int[] ci, double[] v, double[] x, double[] y, int n, int nnz);
}

使用javah生成带有声明的头文件cfd_optimization_Cholesky.h

JNIEXPORT void JNICALL Java_cfd_optimisation_Cholesky_solve_1eigenLDLTx
        (JNIEnv *, jclass, jintArray, jintArray, jdoubleArray, jdoubleArray, jdoubleArray, jint, jint); 

以下是我实现求解器的方法

JNIEXPORT void JNICALL Java_cfd_optimisation_Cholesky_solve_1eigenLDLTx(JNIEnv *env, jclass obj, jintArray arrri, jintArray arrci, jdoubleArray arrv, jdoubleArray arrx, jdoubleArray arry, jint jn, jint jnnz) {
    int n = jn;
    int *ri = (int*)env->GetPrimitiveArrayCritical(arrri, 0);
    int *ci = (int*)env->GetPrimitiveArrayCritical(arrci, 0);
    double *v = (double*)env->GetPrimitiveArrayCritical(arrv, 0);
    int nnz = jnnz;

    double *x = (double*)env->GetPrimitiveArrayCritical(arrx, 0);
    double *y = (double*)env->GetPrimitiveArrayCritical(arry, 0);

    Eigen::SparseMatrix<double> A = colt2eigen(ri, ci, v, nnz, n);
    //Eigen::MappedSparseMatrix<double> A(n, n, nnz, ri, ci, v);

    Eigen::VectorXd a(n), b(n);
    for (int i = 0; i < n; i++) a(i) = x[i];
    //a = Eigen::Map<Eigen::VectorXd>(x, n).cast<double>(); 
    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > solver;
    solver.setMode(Eigen::SimplicialCholeskyLDLT);
    b = solver.compute(A).solve(a);
    for (int i = 0; i < n; i++) y[i] = b(i);
    env->ReleasePrimitiveArrayCritical(arrri, ri, 0);
    env->ReleasePrimitiveArrayCritical(arrci, ci, 0);
    env->ReleasePrimitiveArrayCritical(arrv, v, 0);
    env->ReleasePrimitiveArrayCritical(arrx, x, 0);
    env->ReleasePrimitiveArrayCritical(arry, y, 0);
}

该函数colt2eigen从两个整数数组创建一个稀疏矩阵,该数组包含行和列索引以及值的双数组.

Eigen::SparseMatrix<double> colt2eigen(int *ri, int *ci, double* v, int nnz, int n) {
    std::vector<Eigen::Triplet<double>> tripletList;
    for (int i = 0; i < nnz; i++) { 
        tripletList.push_back(Eigen::Triplet<double>(ri[i], ci[i], v[i]));  
    }
    Eigen::SparseMatrix<double> m(n, n);
    m.setFromTriplets(tripletList.begin(), tripletList.end());
    return m;
}

其中一个棘手的部分是从Java和Colt获取这些数组.为此,我做到了这一点

//y = A x: x and y are double[] arrays and A is DoubleMatrix2D
int nnz = A.cardinality();
DoubleArrayList v = new DoubleArrayList(nnz);
IntArrayList ci = new IntArrayList(nnz);
IntArrayList ri = new IntArrayList(nnz);

A.forEachNonZero((row, column, value) -> {
    v.add(value); ci.add(column); ri.add(row); return value;}
);

Cholesky.solve_eigenLDLTx(ri.elements(), ci.elements(), v.elements(), x, y, n, nnz);