Expose Eigen's Direct and Iterative Sparse System Solvers #18
Conversation
| spmat row_major_coeffs = c_to_eigen(c_coefficients).pruned(); | ||
| spmat row_major_observ = c_to_eigen(c_observations).pruned(); | ||
|
|
||
| coefficients = new col_major_spmat(row_major_coeffs); |
There was a problem hiding this comment.
I am a bit confused by this. Unless I am mistaken, we always use row major matrices in this module, but then we pretend that they are col major here. I understand this for symmetric or self-adjoint matrices, but many of the exposed solvers work on non-symmetric matrices.
Moreover, wouldn't we return a result in the wrong shape?
There was a problem hiding this comment.
The eigen solvers complained when I tried to pass in row-major matrices, so I added these conversions to make it work. I didn't realize this was just recasting the existing memory layout as column-major. My use case is symmetric, so I think I didn't notice that this is actually a bug. This sounds like it will be an issue. Do you know a way to run the solvers on row-major matrices?
There was a problem hiding this comment.
I suppose, in the absence of a way to call the solvers with row-major matrices we could (1) only expose those solvers that operate on symmetric matrices or (2) do a real memory re-map to column major format, though (2) might be too expensive to be worthwhile?
There was a problem hiding this comment.
I don’t know eigen well enough, I’ll look into it sometimes in the next days. @ryanrhymes you wrote the original bindings, maybe you know? Let’s wait a few days and see if we can find a workaround for this. If I don’t find it, I’ll ask for help on eigen repository.
We could make a test with a rectangular matrix and see what happens if we want to double check my guess.
There was a problem hiding this comment.
I don't know it very well either honestly, sorry. Owl internally always uses row-major format, I originally suspected the conversion took too long but become unsure when Tyler said most of the time spent in solve function. I think an easier way is to run the same calculation in other libs like matab/python/julia so that we know how to benchmark.
There was a problem hiding this comment.
What does the solver compain actually? I am reading Eigen's doc, e.g. https://eigen.tuxfamily.org/dox/classEigen_1_1SparseQR.html Although the first parameter is matrix type, but the doc says it must col-major unfortunately. SparseLU has the same requirement, but not for SparseCholesky.
So I think the conversion from row-major to col-major is necessary for QR and LU.
There was a problem hiding this comment.
I began testing with LU and QR, I didn't try to user row-major matrices with other solvers after it failed with LU/QR. The error I saw was a compilation error saying "these functions may only be called with colum-major inputs."
There was a problem hiding this comment.
@steelet41 I think we only need an additional step in the templates.
According to the documentation (see https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html), ```
Matrices and arrays using one storage order can be assigned to matrices and arrays using the other storage order, as happens in the above program when Arowmajor is initialized using Acolmajor. Eigen will reorder the entries automatically. More generally, row-major and column-major matrices can be mixed in an expression as we want.
We could check if a matrix is symmetric, and in that case do nothing and proceed with the actual code, and otherwise have an intermediate cast to col-major as done in the documentation snippet, use that for the sparse solver and then cast the results back to row-major where needed.
There was a problem hiding this comment.
Okay, that's good news. I have actually settled on a different solution for my application (the python solver was just way faster...? a mystery I haven't solved) and I'm leaving the country for two weeks tomorrow, so I won't have a chance to work on this until next year at the earliest.
There was a problem hiding this comment.
There is no rush! Thank you very much. Is it the one in numpy/scipy another one?
| /******************** pointer conversion ********************/ | ||
|
|
||
| typedef SparseMatrix<spmat_c_elt, Eigen::RowMajor, INDEX> spmat_c; | ||
| typedef SparseMatrix<spmat_c_elt, Eigen::ColMajor, int> col_major_spmat_c; |
There was a problem hiding this comment.
Is there a reason for not using INDEX?
There was a problem hiding this comment.
Yeah, that was a weird one. Using INDEX broke compilation for me. Maybe it was some issue local to my environment? Can you try changing it to see if you can reproduce that issue?
|
I used scipy.sparse.linalg.spsolve which solves my system in about 30
seconds, compared to many minutes with eigen.
…On Thu, Dec 12, 2019 at 11:22 AM Marcello Seri ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In eigen_cpp/lib/sparse_solver.cpp
<#18 (comment)>:
> +::eigen_to_c(spmat& ref)
+{
+ return reinterpret_cast<c_spmat*>(&ref);
+}
+
+template <class Solver,
+ typename spmat,
+ typename spmat_elt,
+ typename c_spmat>
+SparseSolver<Solver, spmat, spmat_elt, c_spmat>
+::SparseSolver(c_spmat *c_coefficients, c_spmat *c_observations)
+{
+ spmat row_major_coeffs = c_to_eigen(c_coefficients).pruned();
+ spmat row_major_observ = c_to_eigen(c_observations).pruned();
+
+ coefficients = new col_major_spmat(row_major_coeffs);
There is no rush! Thank you very much. Is it the one in numpy/scipy
another one?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#18?email_source=notifications&email_token=ADXJVQYORQBIAP6TSNNW6M3QYJQK7A5CNFSM4I7OI372YY3PNVWWK3TUL52HS4DFWFIHK3DMKJSXC5LFON2FEZLWNFSXPKTDN5WW2ZLOORPWSZGOCPAA76Q#discussion_r357239640>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ADXJVQ6DB23TCSXQCMOQXVLQYJQK7ANCNFSM4I7OI37Q>
.
|
|
I see. They use https://portal.nersc.gov/project/sparse/superlu/ as solver, it may well be that it is much more optimised. Afaik eigen does transform the sparse matrix to a dense matrix internally in some cases |
|
Perf really depends on underlying solvers and eigen is definitely not the best one. Eigen was originally designed for providing tensor operations for TF's DNN applications, and can be optimised on various archs by exploiting c++ meta-proramming ... not really for solving sparse problems. Tho, Google itself now changes to new solutions. For solving sparse problems, the state-of-the-art is suiteparse, and we are working on this hard, and soon we will release it https://github.com/owlbarn/owl_suitesparse |
I've done the bare minimum here. Unfortunately in my use case (A = 20k x 15k, b = 20k x 7k) none of the solvers complete in less than 10 minutes. These seems egregious, so I suspect I must be doing something inefficient, but I haven't been able to figure out what exactly. It does appear that the time is being spent in the solve() function, and not in the preparation or unpacking functions.
Thanks for your help in getting this to work, and also for any advice you can offer.