[request] Matrix Inverse

[fosskers]

massiv has transpose already - it would be handy for mathematical libraries/applications to have inverse as well.

[fosskers]

Matrix multiplation would be cool too.

[lehins]

Matrix multiplication is already implemented as (|*|), but currently it isn't particularly fast, so it does need some love.

Edit |*| did receive a bit of love ❤️ in caf38ee and got much faster.

Matrix inverse would certainly be a good thing to have for and array library :) So I'll work on it when I get some free time.

[ocramz]

The consensus in the numerics community is to not compute inverse matrices directly, as this is . If the matrix is dense (i.e. mostly nonzero), one could implement a LU factorization once and solve the two associated linear systems in linear time, otherwise there's a wide choice of iterative methods. In both cases, the solution to the linear system might not exist (if for example the matrix has numerically parallel columns or rows, i.e. it's "rank-deficient"). What sort of problems are you looking at @fosskers ?

[ocramz]

I should also mention that linear has 2,3 and 4-dimensional matrix inverses available : https://github.com/ekmett/linear/blob/master/src/Linear/Matrix.hs#L344

[lancelet]

I think matrix inverse can sometimes be useful for didactic purposes, but it would be amazing to have practical linear solvers in pure Haskell.

Is there any way we can feasibly abstract over matrix types for these kinds of algorithms? I'm not sure what it would look like exactly (maybe a standard set of operations in ST/PrimMonad; something like what vector uses? maybe something higher-level? maybe a backpack thing?), but I think it would be great to have some kind of "functional blas" that provides re-usable low-level operations that abstract themselves enough that they could be used on top of vector, massiv, some future sparse matrix format, etc.

(Maybe that's too-ambitious a goal without some concrete implementations first...?)

[lehins]

@lancelet I don't think it would be a good goal to have a general abstract interface that would work for all array libraries in Haskell. Certainly not a backpack thing :)
I do agree though that it would be great to have a sort of "functional blas", but it is much more likely to be successful if it is implemented for a particular library, and I obviously would be happy if that library would be massiv. As far as vector package, there are O(1) conversion functions from/to massiv arrays, which means that any operations implemented in massiv can be easily applied to vector as well.

All of the functionality is already available for adding matrix inverse, it's just a matter of implementing it. Unfortunately, there is no support for sparse matrices yet, but I do have plans on adding it. Created #50 so not to forget and discuss a possible implementation.

Currently I am busy with translating hip to use massiv, so this ticket isn't on top of my list at the moment, but I'll be happy to accept a PR.