5

我已经问过类似的问题，但这次我会更具体。

我需要在一个for循环中执行一个通常很大的正定对称矩阵（约1000x1000）的 Cholesky 分解。现在，为此，我一直在尝试：

1) Apache 数学库

2) 并行 Colt 库

3) JLapack 库

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

因此，我想知道在 Java 中是否有任何高度优化的外部 Cholesky 分解工具：例如，我一直在考虑CHOLMOD算法，它实际上在内部MATLAB和其他工具中调用。

我非常感谢您能就此事获得全面的反馈。

4

2 回答 2

5

以下是 Java 的一些 BLAS 库的一个很好的总结：performance-of-java-matrix-math-libraries。您还可以在Java-Matrix-Benchmark中查看其中许多库的基准测试。

但是，根据我的经验，这些库中的大多数似乎都没有针对解决大型稀疏矩阵进行调整。就我而言，我所做的是通过 JNI使用Eigen实现求解。

Eigen 对其线性求解器进行了很好的讨论，包括 CHOLMOD 中的一个。

对于我通过 JNI 使用 Eigen 求解器的 8860x8860 稀疏矩阵的情况，它比并行 colt 快 20 倍，比我自己的密集求解器快 10 倍。更重要的是，它似乎是按比例缩放的，而且它使用的内存n^2比n^3我的密集求解器要少得多（我的内存用完了）。

实际上有一个用于 Eigen 的 Java 包装器，称为JEigen，它使用 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);
于 2015-05-29T09:30:48.127 回答
1

我没有使用过这些工具中的任何一个，但我怀疑您会因为 Java 在某些版本/某些平台上不使用本机处理器浮点平方根指令这一事实而受到打击。

请参阅：在哪里可以找到 Java 平方根函数的源代码？

如果不需要绝对精度，您可以尝试切换上述实现之一以使用平方根的近似值（请参阅Java 中的 Fast sqrt 以牺牲准确性为代价获取建议），这应该会更快一些。

于 2013-06-11T14:42:06.357 回答