safran-lab / hodlr

A fast, accurate direct solver and determinant computation for dense linear systems

HODLRlib

HODLRlib is a flexible library for performing matrix operations like matrix-vector products, solving and determinant computation in near-linear complexity(for matrices that resemble a HODLR structure). The solver has also been extended to matrices not necessarily arising out of kernels and also to higher dimensions. Further, the solver has been optimized and the running time of the solver is now massively (a few orders of magnitude) faster than the running times reported in the original articles[1][2]. Low-rank approximation of the appropriate blocks are obtained using the rook pivoting algorithm. The domain is sub-divided based on a KDTree. The solver is fairly general, works with minimal restrictions and has been parallelized using OpenMP.

For more details on the usage of the library, visit the documentation page.

Features

MatVecs: Obtains at a cost of

Factorization: Factors the matrix into the desired form at a cost of

Cholesky-like Symmetric Factorization: Obtains at a cost of

Solve: Solves linear systems at an additional cost of

Determinant Computation: Additional Cost of

Version 3.1415

Date: January 6th, 2019

Copyleft 2019: Sivaram Ambikasaran

Developed by Sivaram Ambikasaran, Karan Raj Singh, Shyam Sundar Sankaran

License

This program is free software; you can redistribute it and/or modify it under the terms of MPL2 license. The Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of the MPL was not distributed with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

hodlr's Issues

running error: kepler_test

I just begin to try the example for kepler_test.cpp, but it stopped halfway.
Here the error report:
Factoring the matrix...
kepler: /home/yu/下载/eigen-eigen-dd44029906a1/Eigen/src/Core/Redux.h:413: typename Eigen::internal::traits::Scalar Eigen::DenseBase::redux(const Func&) const [with BinaryOp = Eigen::internal::scalar_max_op<double, double>; Derived = Eigen::PartialReduxExpr<const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op, const Eigen::Matrix<double, -1, -1> >, Eigen::internal::member_sum, 0>; typename Eigen::internal::traits::Scalar = double]: Assertion `this->rows()>0 && this->cols()>0 && "you are using an empty matrix"' failed.
已放弃 (核心已转储)
Is that my Eigen installing problem?

HODLR doesn't work with latest versions of Eigen

I want to thank you for your work and for making the code publicly available. I am using your library heavily to evaluate determinants that are numerically out of reach for standard algorithms like Cholesky decomposition. Without your work my current studies would not be possible.

When working on my code, I recently updated my Eigen library. After the update, the HODLR library gave wrong results. I found that the problem is caused by Eigen versions starting from 3.2.8. Up to 3.2.7, however, everything works fine.

The reason is that in version 3.2.8 Eigen made a change in FullPivLU.h. Here is the relevant quote from the changelog:

Make FullPivLU::solve use rank() instead of nonzeroPivots().

Here is a diff of the change that causes the HODLR library to break:

*** FullPivLU.h_327     2018-03-23 16:02:57.268313250 +0100
--- FullPivLU.h_328     2018-03-23 16:03:22.409744343 +0100
*************** struct solve_retval<FullPivLU<_MatrixTyp
*** 688,694 ****
       */  
  
      const Index rows = dec().rows(), cols = dec().cols(),
!               nonzero_pivots = dec().nonzeroPivots();
      eigen_assert(rhs().rows() == rows);
      const Index smalldim = (std::min)(rows, cols);
  
--- 688,694 ----
       */  
  
      const Index rows = dec().rows(), cols = dec().cols(),
!               nonzero_pivots = dec().rank();
      eigen_assert(rhs().rows() == rows);
      const Index smalldim = (std::min)(rows, cols);

Note that a single line causes the problem: The change from nonzero_pivots = dec().nonzeroPivots(); to nonzero_pivots = dec().rank();.

The problem may be reproduced using the HODLR_Test program. Here is the output using Eigen 3.2.7

Number of particles is: 50000

Setting things up...
Time taken is: 0.001317

Assembling the matrix in HODLR form...
Time taken is: 0.140323

Exact matrix vector product...
Time taken is: 46.2392

Fast matrix matrix product...
Time taken is: 0.013503

Factoring the matrix...
Time taken is: 0.445196

Solving the system...
Time taken is: 0.016266

Error in computed solution: 2.79761e-09

Error in matrix matrix product: 9.35927e-11

Computing the log determinant...
Time taken is: 0.008733

Log determinant is: -15227.07956830352

and here is the output of Eigen 3.3 (note that the result is the same also with other versions greater than 3.2.7)

Number of particles is: 50000

Setting things up...
Time taken is: 0.001308

Assembling the matrix in HODLR form...
Time taken is: 0.139998

Exact matrix vector product...
Time taken is: 45.5071

Fast matrix matrix product...
Time taken is: 0.012069

Factoring the matrix...
Time taken is: 0.453063

Solving the system...
Time taken is: 0.015234

Error in computed solution: -nan

Error in matrix matrix product: 1.45292e-10

Computing the log determinant...
Time taken is: 0.007363

Log determinant is: nan

Unfortunately, I don't have enough oversight over the HODLR library or knowledge of C++ to further trace down the bug. If you decide only to support Eigen up to version 3.2.7, a comment in the readme would be helpful.

HODLR sometimes yields inaccurate results

I found that sometimes the HODLR library yields inaccurate results. This problem might be caused by an empty self-interaction matrix.

As the problem is somehow hard to reproduce, I wrote a test program HODLR_Test.zip.
I am using the latest version of the HODLR library available in this repository, except for one change:
I have commented out the call to srand in HODLR_NODE.hpp in order to be able to reproduce the problem. Also, I am using Eigen in version 3.2.7.

The code computes the determinant of a matrix. The first command line parameter gives the seed of the random number generator. Most of the times, the logdet of the matrix is accurate up to an error of 1e-10 which is a very reasonable result. A typical output looks like:

$ ./HODLR_Test 1522915804
initialize PRNG with seed 1522915804
exact: -117.00871690447
HODLR: -117.00871688608
relative error: 1.57196e-10

However, sometimes (although not very often) the relative error is much bigger:

$ ./HODLR_Test 1522913970
initialize PRNG with seed 1522913970
exact: -117.00871690447
HODLR: -117.00813598053
relative error: 4.96479e-06

Maybe the problem is caused by an empty self-interaction matrix The self-interaction matrix might be empty if the interaction between different blocks is extremely weak (this is the case with this matrix).

After quite some time of debugging, I figured out that one problem might be in
HODLR_Node.hpp. Line 154 Kinverse.compute(K) in compute_Inverse fails if K is
empty. The same problem occurs in the method apply_Inverse.

I have also tested the code using valgrind and the undefined-behavior sanitizer of LLVM: While valgrind doesn't report any errors, LLVM reports:

$ ./HODLR_Test 1522913970
initialize PRNG with seed 1522913970
/path/to/eigen-3.2.7/Eigen/src/Core/PlainObjectBase.h:156:16: runtime error: reference binding to null pointer of type 'double'
/path/to/eigen-3.2.7/Eigen/src/Core/MapBase.h:103:14: runtime error: reference binding to null pointer of type 'const Scalar' (aka 'const double')
/path/to/eigen-3.2.7/Eigen/src/Core/MapBase.h:204:14: runtime error: reference binding to null pointer of type 'double'
/path/to/eigen-3.2.7/Eigen/src/Core/util/BlasUtil.h:141:14: runtime error: reference binding to null pointer of type 'const double'
/path/to/eigen-3.2.7/Eigen/src/Core/PlainObjectBase.h:169:16: runtime error: reference binding to null pointer of type 'const double'
exact: -117.00871690447
HODLR: -117.00813598053
relative error: 4.96479e-06

I have compiled the code using the command:

clang++ -fsanitize=undefined,unsigned-integer-overflow -fno-omit-frame-pointer -w -O3 HODLR_Test.cpp -I /path/to/eigen-3.2.7/ -I ../header

At this point it is really hard for me to further track down the problem, and so I rely on your help. Thank you so much for your library and your hard work. :)

Parallelization of HODLR::solveNonSPD()

Hi,

The function HODLR::solveNonSPD(Mat) has not been parallelized. The following is a straightforward parallelization, which may be useful for some users.

Mat HODLR::solveNonSPD(Mat b)
{
Mat x = Mat::Zero(b.rows(),b.cols());

int r = b.cols();

// Solving over the leaf nodes:
#pragma omp parallel for
for(int k = 0; k < nodes_in_level[n_levels]; k++) 
{
    int start = tree[n_levels][k]->n_start;
    int size  = tree[n_levels][k]->n_size;

    x.block(start, 0, size, r) = this->solveLeafNonSPD(k, b.block(start, 0, size, r));
}

b = x;

// Solving over nonleaf levels:
for(int j = n_levels - 1; j >= 0; j--) 
{
    #pragma omp parallel for
    for (int k = 0; k < nodes_in_level[j]; k++) 
    {
        int start = tree[j][k]->n_start;
        int size  = tree[j][k]->n_size;
        
        x.block(start, 0, size, r) = this->solveNonLeafNonSPD(j, k, b.block(start, 0, size, r));
    }

    b = x;
}

return x;

}

Failed to install Python package - update request

I cannot install and use the Python wrapper for HODLR.
Could you please let me know if you are going to update the package?

install.sh needs to download a newer version of CMake for `add_compile_definitions` support

install.sh is currently downloading CMake 3.3.2, but add_compile_definitions (which is currently used by HODLR) requires CMake >= 3.12. As is, attempting to configure for the unit tests via

../deps/cmake/bin/cmake .. -DCMAKE_BUILD_TYPE=COVERAGE -DINPUT_FILE=../test/test_HODLR.cpp

yields an error

CMake Error at CMakeLists.txt:24 (add_compile_definitions):
  Unknown CMake command "add_compile_definitions".

I therefore recommend using the following replacement within install.sh (which I have verified to fix this issue):

    CMAKE_URL="https://cmake.org/files/v3.12/cmake-3.12.4-Linux-x86_64.tar.gz"

(As an aside, I recommend making install.sh executable so that ./install.sh is possible and adding an echo warning the user that the dependencies may take some time to download and build.)

Python installation failure

The python installation fails with the following error message:

hodlr/_hodlr.cpp: In function ‘int _hodlr_init(_hodlr*, PyObject*, PyObject*)’:
hodlr/_hodlr.cpp:93:42: error: no matching function for call to ‘HODLR_Tree<Gaussian_Matrix>::assemble_Matrix(Eigen::VectorXd&, double&)’
     self->solver->assemble_Matrix(dv, tol);
                                          ^
In file included from hodlr/_hodlr.cpp:6:0:
../header/HODLR_Tree.hpp:73:7: note: candidate: void HODLR_Tree<MatrixType>::assemble_Matrix(HODLR_Node<MatrixType>*&) [with MatrixType = Gaussian_Matrix]
  void assemble_Matrix(HODLR_Node<MatrixType>*& node) {
       ^
../header/HODLR_Tree.hpp:73:7: note:   candidate expects 1 argument, 2 provided
../header/HODLR_Tree.hpp:188:7: note: candidate: void HODLR_Tree<MatrixType>::assemble_Matrix(Eigen::VectorXd&, double, char) [with MatrixType = Gaussian_Matrix; Eigen::VectorXd = Eigen::Matrix<double, -1, 1>]
  void assemble_Matrix(VectorXd& diagonal, double lowRankTolerance, char s) {
       ^
../header/HODLR_Tree.hpp:188:7: note:   candidate expects 3 arguments, 2 provided

It looks like the problem is that char s is not being provided by _hodlr_init to HODLR_Tree<MatrixType>::assemble_Matrix.

Checking for failure

What happens when we get a system that can't be solved? Does it just fail silently and spit out garbage? If so, that's probably something that we should look into...

clang error [Fixed]

When installing HODLR on a fresh install of OSX Mavericks, I got this error:

clang: error: -O4 is equivalent to -O3

Google found openMVG/openMVG#102

The solution is to change line 81 of CMakeLists.txt from "-O4" to "-flto". Everything else stayed the same. The issue is the the newest Xcode does not attribute LTO directly with -O4 (I'm paraphrasing the link above).

I was hoping there could be a way for future builds to avoid this issue.

Thank you all the same for a great package

-Jonathan

Improvements to CMake code needed

The cmake code for this library causes the exact same (very expensive) compilation to happen many times over when building examples.
Building anything with Eigen is already expensive, so I think it would be best to avoid this.

If there is no good reason for this then I would be happy to make some changes to the CMake code and submit a PR.
Currently there isn't actually a HODLR library, and changing that could drastically improve compilation time.

I am also trying to use HODLR in my own project, and several improvements to the CMake code could greatly improve the user experience.
Again, I would be happy to make some improvements in a PR if that is OK.