Feature Request: add Matrix-Matrix solve

[rebcabin]

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

X = A^{-1} B

into equations like this

X = solve(A, B)

The current implementation of solve in core.matrix only works when B is a 1D vector, not when it is a 2D column vector, let alone an nxm matrix. I therefore propose something like the following

(require '[clojure.core.matrix :as ccm])
(ccm/set-current-implementation :vectorz)
(require '[clojure.core.matrix.linear :as ccml])

(defn solve-matrix
  "The 'solve' routine in clojure.core.matrix only works on Matrix times Vector.
  We need it to work on Matrix times Matrix. The equation to solve is

  Ann * Xnm = Bnm

  Think of the right-hand side matrix Bnm as a matrix of columns. Iterate over
  its transpose, treating each column as a row, then converting that row to a
  vector, to get the transpose of the solution X."
  [Ann Bnm]
  (ccm/transpose (mapv (partial ccml/solve Ann) (ccm/transpose Bnm))))

[mikera]

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.

[mikera]

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.

[rebcabin]

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?