This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

Bug 1179

Summary: Hypersparse matrix class proposal
Product: Eigen Reporter: Dmitry Zhdanov <dm.zhdanov>
Component: SparseAssignee: Dmitry Zhdanov <dm.zhdanov>
Status: REVIEWNEEDED ---    
Severity: Feature Request CC: Angelos.Mantzaflaris, chtz, dm.zhdanov, gael.guennebaud, r.moraes
Priority: Normal    
Version: 3.4 (development)   
Hardware: All   
OS: All   
Whiteboard:
Bug Depends on: 1104    
Bug Blocks: 1608    
Attachments:
Description Flags
test program
none
HyperSparseMatrix class
none
sparse_accumulator container
none
test program
none
HyperSparseMatrix class
none
HyperSparseMatrix class
none
sparse_accumulator container
none
test program
none
HyperSparseMatrix class
none
HyperSparseMatrix class
dm.zhdanov: review? (gael.guennebaud)
test program none

Description Dmitry Zhdanov 2016-03-15 19:30:13 UTC
I implemented a class for hypersparce (NxM) matrices such that nnz<<max(N,M), where nnz is number of nonzero elements. This class allows to dramatically reduce the storage requirements (O(nnz) instead of O(min(N,M))) and affords an opportunity to essentially speed up basic matrix operations with such matrices.  

Technically, the class is based on doubly compressed sparse column\row (DCSC(R)) data structure (see e.g. A. Buluc¸ and J. R. Gilbert "On the representation and multiplication of hypersparse matrices", in IPDPS’08: Proceedings of the 2008 IEEE International Symposium on Parallel and Distributed Processing, IEEE Computer Society, Washington, DC, 2008, pp. 1–11).

My implementation seems to be compatible with the Eigen 3.3 and is based on step-by-step reproduction of SparseMatrix.h functionality (except for diagonal traits which I did not test). So, by their construction, the DCSC(R) matrices are well integrated into Eigen and can be used virtually everywhere instead of ordinary Eigen sparse matrices. However, the computational performance is not optimal (i.e. it is similar to one for usual sparse matrices) because of some generic limitations of generic Eigen sparse iterators. To avoid these limitations wrote some ad-hoc utilities which might be interesting on their own (for example, parallel vectorized algorithm for (hyper)sparse-(hyper)sparse matrix multiplication which resembles one used in Matlab). I am actively using DCSC class in my quantum chemistry research where it is really indispensable.

So, there are 3 questions to Eigen developers.
1) Would you be interested to see the hypersparse matrices as a part of Eigen project?
2) If so, would you be willing to add some features to generic SparseMatrix implementation to allow for the performance improvements for DCSC(R) matrices?
3) Are you interested in including implementations of high-performance (hyper)sparse-(hyper)sparse matrix multiplications?
Comment 1 Gael Guennebaud 2016-03-21 09:55:40 UTC
Hi,

sure! this would be a nice addition to Eigen, and it perfectly fits within Eigen's scope. I guess, that the main missing bit would be a general iterator over *all* non-zeros, or perhaps an "outer-iterator" to leverage a bit more the structure?
Comment 2 Dmitry Zhdanov 2016-03-21 19:11:10 UTC
Hi Gael.

I am very glad that you are interested in this proposal!

Regarding your question: yes, exactly. In the specific case of (DCSC(R)) matrix structure I am thinking about outer iterator over all non-zero rows\columns. But in general, it would make sense to template the outer iterator in order to allow for different access schemes (including the iteration over all non-zeros as it was suggested in some previous feature requests) in order to fit the needs of different algorithms and to facilitate painless introduction of new sparse formats in the future.

For the second part of proposal (sparse-sparse matrix multiplications) I do not have the immediate suggestion about the best interface within the Eigen ideology. The conceptual issue is that such multiplication needs for temporary structures which size might be large compared to the sizes of sparse matrices. So, creation and destruction of these structures for each multiplication looks as unnecessary waste of resources. On other side, involvement of the external structures conflicts with preserving the nice form of multiplication C = A * B;. So far, I implemented this functionality as an external class sparse_evaluator in which the multiplication operations are realized as functions, e.g.
sparse_evaluator<std::complex<double> /*Scalar*/, int /*Index*/> sev(maxOuterSize);
sev.multiply(C,A,B); //A, B, C is any combination of sparse of hypersparse matrices, C also can be dense.

I will add the code for (DCSC(R)) matrices as attachment to this proposal in few days.
Comment 3 r.moraes 2016-03-24 09:43:36 UTC
I am very much interested in this feature. I am developing an application that would very much benefit from it. Is it possible to include it, perhaps, in the current development version (3.3). I would be more than happy to perform tests on the implementation.
Comment 4 Gael Guennebaud 2016-03-24 12:27:03 UTC
Dmitry,

regarding the temporary buffers issue: this is a very general issue we have in several other places. My idea was to implement our own memory allocator for them: the allocation/deallocation pattern is like a stack, so extremely fast, and all we need compared to the default stack is to make it expandable through a linked list.
Comment 5 Dmitry Zhdanov 2016-03-26 01:46:58 UTC
(In reply to Gael Guennebaud from comment #4)
> Dmitry,
> 
> regarding the temporary buffers issue: this is a very general issue we have
> in several other places. My idea was to implement our own memory allocator
> for them: the allocation/deallocation pattern is like a stack, so extremely
> fast, and all we need compared to the default stack is to make it expandable
> through a linked list.

Gael,

your idea about stack sounds as a great solution to the problem though there might be some complications with multithreaded applications.

I am now annotating the code. What kind of credits or copyright statements can I put in the header? Can I name you as a contributor? (I took about 30% of code from SparseMatrix.h with minor adaptations.) From above comment it looks like there is already a person interested in starting using it. At this stage I am keeping everything in separate namespace.
Comment 6 Dmitry Zhdanov 2016-03-26 02:02:56 UTC
(In reply to r.moraes from comment #3)
> I am very much interested in this feature. I am developing an application
> that would very much benefit from it. Is it possible to include it, perhaps,
> in the current development version (3.3). I would be more than happy to
> perform tests on the implementation.

I am also glad you are interested in my work. I will attach code to this proposal after finishing annotating it. I heavily used the code for various manipulations with complex matrices, but there might be some issues due to minor changes I am introducing right now for publication and which are obviously not tested well. I also did not extensively check different combinations of sparse and dense Index types (there was an issue with specifying non-default index type for dense matrices).
Comment 7 Dmitry Zhdanov 2016-03-29 08:55:04 UTC
Created attachment 676 [details]
test program

test program which also shows the sample usage of the HyperSparseMatrix class (tested with MSVC 2015 only).
Comment 8 Dmitry Zhdanov 2016-03-29 08:57:25 UTC
Created attachment 677 [details]
HyperSparseMatrix class

implementation of Doubly Compressed Sparse matrices. Most of the code was tested with MSVC 2013, 2015 and the latest g++ compilers but the very last edits were tested with MSVC 2015 only.

Few known issues are marked by ??? signs in the code (and also in the test program dcs_test.cpp).
Comment 9 Dmitry Zhdanov 2016-03-29 08:59:25 UTC
Created attachment 678 [details]
sparse_accumulator container

sparse_accumulator utility container used by dcs.h It is useful for broad variety of linear algebra operations including sparse matrix multiplications, transpositions etc., so it is worth to keep it separated dcs.h
Comment 10 Dmitry Zhdanov 2016-03-29 09:13:54 UTC
Hello everyone interested. I attached the code of the proposed class. Gael, I copied the copyright preamble from the SparseMatrix.h file (Mozilla license bla-bla-bla) and stated there that the code is grounded on your work.

Few known issues are marked by ??? signs in the code. Please, note that I did not test conservativeResize() function extensively. 

This test implementation is immersed into the Eigen namespace, but one can easily do it by renaming namespace dcs by Eigen everywhere in dcs.h.

Looking forward for your feedback.
Comment 11 r.moraes 2016-03-29 09:17:29 UTC
Thank you very much, Dmitry, for sharing this implementation!

I'll be testing your implementation as soon as possible and be adding any feedback I have.
Comment 12 r.moraes 2016-03-29 11:19:28 UTC
Unfortunately I only have VS2010 installed at this point in time. If I try try to compile the code in this VS version I get the following error:

1>------ Build started: Project: EigenTests, Configuration: Release x64 ------
1>  dcs_test.cpp
1><PATH>\dcs.h(108): error C2143: syntax error : missing ';' before '<'
1><PATH>\dcs.h(108): error C2059: syntax error : '<'
1><PATH>\dcs.h(108): error C2065: '_Scalar' : undeclared identifier
1><PATH>\dcs.h(108): error C2065: '_Options' : undeclared identifier
1><PATH>\dcs.h(108): error C2065: '_Index' : undeclared identifier
1><PATH>\dcs.h(109): error C2065: 'SparseCompressedBase' : undeclared identifier
1><PATH>\dcs.h(109): error C2065: '_Scalar' : undeclared identifier
1><PATH>\dcs.h(109): error C2065: '_Options' : undeclared identifier
1><PATH>\dcs.h(109): error C2065: '_Index' : undeclared identifier
1><PATH>\dcs.h(109): error C2955: 'dcs::HyperSparseMatrix' : use of class template requires template argument list
1><PATH>\dcs.h(47) : see declaration of 'dcs::HyperSparseMatrix'
1><PATH>\dcs.h(109): fatal error C1903: unable to recover from previous error(s); stopping compilation
========== Build: 0 succeeded, 1 failed, 0 up-to-date, 0 skipped ==========

Perhaps something related to how the dcs::HyperSparseMatrix forward declaration is handled by the VS2010 compiler?

Anyways, I will try to have a newer VS version installed and test it again.
Comment 13 Dmitry Zhdanov 2016-03-29 13:56:40 UTC
(In reply to r.moraes from comment #12)
> Unfortunately I only have VS2010 installed at this point in time. If I try
> try to compile the code in this VS version I get the following error:
> 
> 1>------ Build started: Project: EigenTests, Configuration: Release x64
> ------
> 1>  dcs_test.cpp
> 1><PATH>\dcs.h(108): error C2143: syntax error : missing ';' before '<'
> 1><PATH>\dcs.h(108): error C2059: syntax error : '<'
> 1><PATH>\dcs.h(108): error C2065: '_Scalar' : undeclared identifier
> 1><PATH>\dcs.h(108): error C2065: '_Options' : undeclared identifier
> 1><PATH>\dcs.h(108): error C2065: '_Index' : undeclared identifier
> 1><PATH>\dcs.h(109): error C2065: 'SparseCompressedBase' : undeclared
> identifier
> 1><PATH>\dcs.h(109): error C2065: '_Scalar' : undeclared identifier
> 1><PATH>\dcs.h(109): error C2065: '_Options' : undeclared identifier
> 1><PATH>\dcs.h(109): error C2065: '_Index' : undeclared identifier
> 1><PATH>\dcs.h(109): error C2955: 'dcs::HyperSparseMatrix' : use of class
> template requires template argument list
> 1><PATH>\dcs.h(47) : see declaration of 'dcs::HyperSparseMatrix'
> 1><PATH>\dcs.h(109): fatal error C1903: unable to recover from previous
> error(s); stopping compilation
> ========== Build: 0 succeeded, 1 failed, 0 up-to-date, 0 skipped ==========
> 
> Perhaps something related to how the dcs::HyperSparseMatrix forward
> declaration is handled by the VS2010 compiler?
> 
> Anyways, I will try to have a newer VS version installed and test it again.

Heh... What happens if you move the entire traits declaration block to the very end of file and delete the forward declaration line?
Comment 14 Dmitry Zhdanov 2016-03-29 20:50:42 UTC
Created attachment 679 [details]
test program
Comment 15 Dmitry Zhdanov 2016-03-29 20:51:57 UTC
Created attachment 680 [details]
HyperSparseMatrix class
Comment 16 Dmitry Zhdanov 2016-03-29 22:52:07 UTC
Created attachment 681 [details]
HyperSparseMatrix class

minor compatibility fixes:
Now tested on compilers gcc-4.6.3 and gcc-5.1.0
Comment 17 Dmitry Zhdanov 2016-03-29 22:53:02 UTC
Created attachment 682 [details]
sparse_accumulator container

minor compatibility fixes:
Now tested on compilers gcc-4.6.3 and gcc-5.1.0
Comment 18 Dmitry Zhdanov 2016-03-29 22:53:49 UTC
Created attachment 683 [details]
test program

minor compatibility fixes:
Now tested on compilers gcc-4.6.3 and gcc-5.1.0
Comment 19 Dmitry Zhdanov 2016-03-30 03:36:48 UTC
Created attachment 684 [details]
HyperSparseMatrix class

One more little touch in order to support Intel compiler.
So, currently the confirmed compilers are:
- MSVC++ 14.0 (Visual Studio 2015)
- Intel 2011.3 and higher
- gcc-4.6.3 and higher
Comment 20 Dmitry Zhdanov 2016-04-01 19:00:10 UTC
Created attachment 690 [details]
HyperSparseMatrix class

Two bugs fixed on pruning empty matrices and constructing matrix from empty triplets list.
Comment 21 Dmitry Zhdanov 2016-04-01 19:04:21 UTC
Created attachment 691 [details]
test program

The tests are now performed on matrices having random number of nonzero elements (including empty matrices).
Comment 22 Dmitry Zhdanov 2016-04-07 06:01:30 UTC
(In reply to Gael Guennebaud from comment #4)
> Dmitry,
> 
> regarding the temporary buffers issue: this is a very general issue we have
> in several other places. My idea was to implement our own memory allocator
> for them: the allocation/deallocation pattern is like a stack, so extremely
> fast, and all we need compared to the default stack is to make it expandable
> through a linked list.

Gael, 
1) Should I start integrating the module into Eigen via Mercurial? What would be the best place for sparse_accumulator container? (Keep as separate file? Add to SparseMatrixBase.h?) (I would suggest to put it in internal namespace and, maybe, rename to match Eigen naming conventions.).

2) I thought about your idea on allocating temporary buffers. Now I think that it is not as good as it seemed initially. Any buffer must be initialized prior the first use. Hence, one has to either a) reinitialize the part of memory stack after each reallocation or b) ensure that each buffer will take care about preparing the allocated memory for the next use on destruction. In the case of sparse multiplications each solution may lead to notable write overheads (though, variant b) is substantially better than a)).
Comment 23 r.moraes 2016-04-13 09:30:07 UTC
(In reply to Dmitry Zhdanov from comment #13)
> (In reply to r.moraes from comment #12)
> > Unfortunately I only have VS2010 installed at this point in time. If I try
> > try to compile the code in this VS version I get the following error:
> > 
> > 1>------ Build started: Project: EigenTests, Configuration: Release x64
> > ------
> > 1>  dcs_test.cpp
> > 1><PATH>\dcs.h(108): error C2143: syntax error : missing ';' before '<'
> > 1><PATH>\dcs.h(108): error C2059: syntax error : '<'
> > 1><PATH>\dcs.h(108): error C2065: '_Scalar' : undeclared identifier
> > 1><PATH>\dcs.h(108): error C2065: '_Options' : undeclared identifier
> > 1><PATH>\dcs.h(108): error C2065: '_Index' : undeclared identifier
> > 1><PATH>\dcs.h(109): error C2065: 'SparseCompressedBase' : undeclared
> > identifier
> > 1><PATH>\dcs.h(109): error C2065: '_Scalar' : undeclared identifier
> > 1><PATH>\dcs.h(109): error C2065: '_Options' : undeclared identifier
> > 1><PATH>\dcs.h(109): error C2065: '_Index' : undeclared identifier
> > 1><PATH>\dcs.h(109): error C2955: 'dcs::HyperSparseMatrix' : use of class
> > template requires template argument list
> > 1><PATH>\dcs.h(47) : see declaration of 'dcs::HyperSparseMatrix'
> > 1><PATH>\dcs.h(109): fatal error C1903: unable to recover from previous
> > error(s); stopping compilation
> > ========== Build: 0 succeeded, 1 failed, 0 up-to-date, 0 skipped ==========
> > 
> > Perhaps something related to how the dcs::HyperSparseMatrix forward
> > declaration is handled by the VS2010 compiler?
> > 
> > Anyways, I will try to have a newer VS version installed and test it again.
> 
> Heh... What happens if you move the entire traits declaration block to the
> very end of file and delete the forward declaration line?

Hi Dmitry,

Thanks for all the updates. I've been testing your implementation and it might help me with my application.

I was able to compile and run the test in VS2010. Everything works fine.

I used to require a memory allocation of around 1.1 GB even if I created empty matrices (please refer to https://forum.kde.org/viewtopic.php?f=74&t=131714). However, with your implementation I now require 400 MB. The allocation of data structures seems to be faster as well.

One note: apparentely there is no method reserve(sparsityPattern) where sparsityPattern is a vector that indicates the number of non-zero entries. Would it make sense to have such method in this new matrix class?

Thanks!
Comment 24 Gael Guennebaud 2016-04-13 13:37:55 UTC
Dmitry, since this feature requires some changes in the current SparseCore module to update iteration on the outer-dimension to use iterators instead of manual for loops (where it makes sense), I guess that the best would be to create a fork on bitbucket dedicated to this feature. This will also ease testing and code review.


Regarding the implementation of HyperSparseMatrix (HSM) itself, I'm wondering whether we could share more code with SparseMatrix (SM) by implementing HSM on top of SM: HSM would store a SM and simply add one indirection maps from the compact SM to the HSM. One could also access to the internal SM for random column access to the non empty columns if needed. I haven't though much about it though...


Regarding temp buffer, hard to comment without seeing the multiplication code.
Comment 25 Dmitry Zhdanov 2016-04-13 15:51:47 UTC
(In reply to Gael Guennebaud from comment #24)
> Dmitry, since this feature requires some changes in the current SparseCore
> module to update iteration on the outer-dimension to use iterators instead
> of manual for loops (where it makes sense), I guess that the best would be
> to create a fork on bitbucket dedicated to this feature. This will also ease
> testing and code review.
> 
> 
> Regarding the implementation of HyperSparseMatrix (HSM) itself, I'm
> wondering whether we could share more code with SparseMatrix (SM) by
> implementing HSM on top of SM: HSM would store a SM and simply add one
> indirection maps from the compact SM to the HSM. One could also access to
> the internal SM for random column access to the non empty columns if needed.
> I haven't though much about it though...
> 
> 
> Regarding temp buffer, hard to comment without seeing the multiplication
> code.

Gael, on iterators issue: maybe it would be a good idea to open the unifying feature request on this issue in order to more clearly understand the necessary functionality? The problem is obviously broader than it is described in bug 1104. I thought about this problem only in context of matrix multiplication.  

On HSM implementation: Using the SparseMatrix as a private inner container for HyperSparseMatrix is simpler in some sense but will substantially reduce integrity of code and would complicate e.g. matrix multiplication implementation (though, this severity of this problem will be greatly reduced after introduction of improved iterators). For these reasons, my initial idea was to implement HSM as a class derived from SparseMatrix. However, if you compare dcs.h and SparseMatrix.h line by line you will see that almost every function (except for constructors and couple of little routines like default_prunning_func) is changed and despite the number of unchanged lines of code is relatively large (see e.g. insertCompressed) they can not be directly utilized from current SparseMatrix implementation. I do not think that touching SparseMatrix.h right now is a right way to go and would suggest to introduce HyperSparseMatrix as an independent feature. As a proposal for a long future, I would suggest to add one more intermediate base class which would implement some common low-level functionality and to derive both SparseMatrix and HyperSparseMatrix from it. 

On containers issue: I can extract and attach the matrix-multiplication-related parts of code from my project. However, it might be difficult to understand the code because it will include many unrelated features (including transposeInPlace implementation which I have mentioned in relation to bug 226) and will not be documented. I can supply some working examples though. Would you be willing to take a look at the code in such a form?
Comment 26 Dmitry Zhdanov 2016-04-13 16:14:49 UTC
> Hi Dmitry,
> 
> Thanks for all the updates. I've been testing your implementation and it
> might help me with my application.
> 
> I was able to compile and run the test in VS2010. Everything works fine.
> 
> I used to require a memory allocation of around 1.1 GB even if I created
> empty matrices (please refer to
> https://forum.kde.org/viewtopic.php?f=74&t=131714). However, with your
> implementation I now require 400 MB. The allocation of data structures seems
> to be faster as well.
> 
> One note: apparentely there is no method reserve(sparsityPattern) where
> sparsityPattern is a vector that indicates the number of non-zero entries.
> Would it make sense to have such method in this new matrix class?
> 
> Thanks!

The sparsityPattern feature is based on SparseMatrix uncompressed mode. I intentionally dropped support for this mode because it might lead to notable memory overheads. As a workaround, I would suggest to reserve the memory for total known number of nonzeros and then sequentially fill the matrix contents row-by-row (for RowMajor alignment) using low-level API functions (see line 364 and below in dcs.h).
Comment 27 Gael Guennebaud 2016-11-30 16:32:14 UTC
*** Bug 1353 has been marked as a duplicate of this bug. ***
Comment 28 Nobody 2019-12-04 15:32:54 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1179.