Add a generic implementation of matrix determinant

[mikera]

We don't currently have a default implementation for "determinant".

We should ideally have something here that works, even if it is a bit naive / slow.

[BlanceXR]

How can I use the determinant and inverse function for now? I'm kind of new to clojure, not sure how to call the non-default implementation... When I call them directly, I got: IllegalArgumentException No implementation of method: :determinant of protocol: #'clojure.core.matrix.protocols/PMatrixOps found for class: clojure.lang.PersistentVector clojure.core/-cache-protocol-fn

[mikera]

Currently you need to use a core.matrix implementation that has implemented the determinant and inverse. e.g. you can use vectorz-clj : https://github.com/mikera/vectorz-clj

(set-current-implementation :vectorz)

(det (matrix [[1 2] [3 4]]))
=> -2.0

(inverse (matrix [[1 2] [3 4]]))
=> #<Matrix22 [[-2.0,1.0],[1.5,-0.5]]>

Note that using (set-current-implementation :vectorz) will cause functions like array and matrix to produce Vectorz matrices by default. You can also get this behaviour explicitly with something like (matrix :vectorz [[1 2] [3 4]])

Once someone addresses this issue and writes a default implementation in core.matrix itself, these operations should work on any matrices.

Note that as a user, this should be completely transparent: det and inverse will just magically start working for all matrices in a future release.

[BlanceXR]

It works! Thanks!

[mikera]

@si14 - should be possible to implement this with your LU decomposition algorithm, which provides a handy way to compute the determinant efficiently:

See - http://en.wikipedia.org/wiki/LU_decomposition#Computing_the_determinant

[si14]

@mikera I've added a similar determinant function to ndarray-double in develop branch. An additional benefit is that the particular LU-decomposition that I use produces L with singular diagonal, so only multiplication of U's diagonal is needed. I think I'll also add a default implementation based on coercing to ndarray-double.

[si14]

Done in 533e8cb

[mikera]

Look good! I think coercing to ndarray-double is a decent strategy for the default implementation.

[lilian0312]

This looks strange:

=> (def mmm [[0.0 0.0 0.0 0.0 1.0] [0.0 0.0 1.0 0.0 0.0] [0.0 1.0 0.0 1.0 0.0] [0.0 0.0 0.0 1.0 0.0] [0.0 0.0 1.0 0.0 1.0]]) mmm => (det mmm) IllegalArgumentException lu-decompose can't decompose singular matrix clojure.core.matrix.impl.ndarray/lu-decompose!-double/fn--9591 (ndarray.clj:1169)

Looks like a bug to me. I am using the default matrix implementation on an ubuntu machine.

[mikera]

Hi @lilian0312 - agreed, that looks like a singular matrix so the result should be zero, right?

As an immediate workaround you can use vectorz-clj which gets the result right.

@si14 - are you able to fix this one in NDArray? If not we may need to plug in an alternative determinant calculation for the time being.

[dmarjenburgh]

This one is open for a while. One easy fix would be to catch the exception from lu-decompose! and return zero. The only problem I see is that the cause for the exception might not be certain. We could throw ex-data instead of iae's, at least within the ndarray namespace and wrap them with iae's for the users of the namespace.

Related, the ndarray invert fn throws an exception while the protocol docstring of inverse says it should be nil for singular matrices. The solution would be the same.

Are you interested in a PR for this if @si14 is not working on it anymore?

[mikera]

Yeah would be good to have a PR if this is still an issue.

I believe the decompose function shouldn't be throwing an error anyway, it should return nil I think if the matrix is singular.