Fix mmul for complex matrices

mikera commented 8 years ago

Current mmul is broken with core.matrix 0.50.0. Needs a fix to ensure that this works with complex numbers.

In the meantime, inner-product works as a replacement

gituser768 commented 8 years ago

With core.matrix 0.34.0, (m/mmul (complex-array [[1 1] [1 1]]) (complex-array [[1 1] [1 1]])) gives me:

ClassCastException org.apache.commons.math3.complex.Complex cannot be cast to java.lang.Number clojure.lang.Numbers.multiply (Numbers.java:148)

As well as on 0.50.0. Is this the use case you're addressing?

I would expect to receive this exception since PMatrixMultiply is not implemented in core.matrix.complex; correct?

mikera commented 8 years ago

Yes, that is the problem. Two options to fix it, we should probably consider both: a) Implement PMatrixMultiply for complex arrays in core.matrix.complex b) Alter the default implementation in core.matrix so that it can potentially work with non-real scalars

gituser768 commented 8 years ago

It looks like a the default implementation can be slightly modified to accommodate a more general set of scalars. Changing the end of the extend-protocol mp/PMatrixMultiply matrix-multiply definition in core.matrix to:

(mp/set-2d! new-m i j (new-add (mp/get-2d new-m i j)
                                        (new-mul (mp/get-2d m i k)
                                           (mp/get-2d a k j)))))
             new-m))))
...

where the preceeding let statement contains:

new-add (mp/generic-add new-m)
new-mul (mp/generic-mul new-m)

And then calling:

(def A (cm/complex-array (m/identity-matrix 2)))
(mp/matrix-multiply A A)

Gives:

IllegalArgumentException No implementation of method: :generic-add of protocol: #'clojure.core.matrix.protocols/PGenericOperations found for class: clojure.core.matrix.impl.ndarray_object.NDArray clojure.core/-cache-protocol-fn (core_deftype.clj:554)

Which can be solved by assigning the generic operations for NDArray. I'll look more into this tomorrow. What are your thoughts?

gituser768 commented 8 years ago

Additionally, the definition of mp/generic-add and mp/generic-mul for NDArray should also be generic arithmetic functions which dispatch on the arguments. I expect that doing this would enable quite a few of the default implementation functionality for complex-array.

algo.generic seems to be what we need. Perhaps not all the functionality of algo.generic, but at least a similar approach.

Alternatively, if core.matrix.complex included a mutable implementation, NDArray would not be used in matrix-multiply. Accordingly, the generic-add and generic-mul of complex-array (already defined as c/+ and c/*) would be used, giving us mmul on complex-array.

mikera commented 8 years ago

I feel a bit uncomfortable doing generic-add and generic-mul at the per-element level. Two big issues:

Adds a lot of overhead at the per-element level when it isn't needed (e.g. for normal numbers)
Requires implementations to do all the generic operations

I'd rather see a default implementation for mmul that expresses the inner product in terms of lower dimensional operations, ultimately handled by stuff like dot-product or scale at the bottom level. so for matrices it would just be dot products of rows with columns etc. This is most likely to work well with a wide range of different underlying implementations (including complex).

GWP commented 8 years ago

Hi Mikera and dmh43,

I'm very much interested in using this library for quantum computer simulation, but that's impossible if there is no way to multiply complex matrices. If this issue could be resolved, that would be great. I would attempt it myself, but I don't know where to start.

GWP commented 8 years ago

I haven't been able to get inner-product working for complex matrices either; it ends up throwing the same error, that Complex cannot be cast to Number.

gituser768 commented 8 years ago

@mikera, looks to me like matrix-multiply should use m/scale instead of * but there isn't an equivalent for + that I know of which will dispatch according to the type without doing much unnecessary boxing and unboxing of the values. For example, using matrix-add. Any other suggestions?