Open rebcabin opened 6 years ago
It's a reasonable feature request, providing we can be very clear on the precise definition of "solve".
The challenge is that I think it is hard to provide a good default implementation, this is a tricky operation to implement. If anyone has good ideas, happy to discuss.
Other idea would be to create a protocol and make it an "optional" operation for implementations to support. But then we would need to figure out how to test / validate effectively.
I've had good luck with Cholesky solvers in Eigen/C++ and LAPACKE/C (see around line 191 in this gist). I think LAPACK doesn't even have a matrix inverse routine! If my good luck were to fail, then I would go to Golub & Van Loan for advice. I like the protocol idea (with a reasonable default implementation like Cholesky so that the new user doesn't have a barrier to getting started), and I wonder whether we could just borrow test suites from the LAPACK world?
Apologies if this is the not the proper place to propose new features. Please let me know if this is improper and I'll transport my request to the proper place. In case it is proper:
To avoid computing inverses and because linear solvers are usually expected to be stable and robust, I normally translate matrix equations like this
into equations like this
The current implementation of
solve
incore.matrix
only works whenB
is a 1D vector, not when it is a 2D column vector, let alone an nxm matrix. I therefore propose something like the following