Regression tests for matrix calculations

[Friedrich Boeckh]

Hi, Dmitry Gronchev and Mike Anderson,

I am refering to my blog post, on which I received comment from you this week.

As I understood that you find that kind of analysis useful, I thought about how to adapt it for your purposes. From what I understand, the tests would need to be integrated "somehow" into the following parts:

matrix-api/src/test/clojure/core/matrix
matrix-api/src/main/clojure/clojure/core/matrix/compliance_tester.clj
matrix-api/src/dev/clojure/clojure/core/matrix/docgen/bench.clj#L89

My analysis of your code structure makes me think that our test strategies differ significantly and may be therefore complementary.
For example, I purposely did not want to have a random matrix generator like in bench.clj, but I wanted to have full control about the test matrix, e.g. to be able import a result from NumPy for functional tests.

I had chosen the Midje framework mainly because of the filter mechanism for test cases, the REPL integration by 'midje.repl' and the (controlable) higher level of verbosity, which helps me during development, but I do not want to burden you with a discussion on this.

Under the bottom line, I think it is best to keep my project separate and update it accordingly to your amazing progress of developing the matrix library. Please let me know, if I could contribute something specific to support you.

Thank you for your great work on core.matrix!

[Mike Anderson]

Glad you like core.matrix! 

And it's great to have extra test coverage, both in the functional and performance sense. Thank you for your efforts here.

If you'd like to move some tests into core.matrix itself I'm totally fine with that. I can give you commit access providing you have a Clojure CA. On the other hand, if you want to keep the project separate that's fine too - in some ways this has advantages in terms of being able to mix and match different core.matrix versions / various implementations for testing.

[Friedrich Boeckh]

Going through the test suite again, I had the idea to generate a Householder marix (as named P in section 'Definition' of http://en.wikipedia.org/wiki/Householder_reflection) from an arbitrary random vector of, say, size 50.

This approach would provide nicely a non-trivial closed test loop of several functions (outer-product/mmul, normalise/dot product, identity, sub, inverse, transpose, mmul), which can be verified by the specific properties of P being symmetric, unitary and involutary. It would be especially beneficial in adding the first test case for inverse. I suggest to add these tests in matrix-api/src/main/clojure/clojure/core/matrix/compliance_tester.clj.

In a second step, I could enhance the benchmark tests in matrix-api/src/test/clojure/clojure/core/matrix. Finally, the CFD Solver, preferably lesson 12, could be added as a demo case.

If this is ok, I´d start with task 1 and in parallel file in an application for a CA (no idea yet, who this works).

[Mike Anderson]

Good idea: it's a nice way to test various functions.

The householder matrix is pretty useful in its own right, I think it would be worth implementing a "householder" API function in clojure.core.matrix itself.

The only mild complexity with this kind of test is that you need to restrict it to numerical matrix types that support floating-point values. An array of longs wouldn't pass, for example :-). Perhaps this also needs an API predicate to test for this property, but it's not totally clear what it should test for in the general sense - it gets tricky when you think about complex arrays, binary matrices etc.

[Friedrich Boeckh]

You may want to check the additional test here: https://github.com/fbmnds/matrix-api/blob/master/src/main/clojure/clojure/core/matrix/compliance_tester.clj#L374. It is as easy as expected.

The test is safeguarded by 'test-numeric-instance' in line 403. Although redundant, as the test matrix is generated independently, it is logically correct at this position, I´d argue. A nice extention would be to use the test for complex matrices (with Hermite transformation in lieu of 'transpose').

It should be equally straightforward to add a 'householder' function to the API. It just operates on any given numerical array/vector. I would need to look a little bit deeper in the implementation layer, though.

[Mike Anderson]

Looks good! Happy to take a PR on this assuming you're going through the CA process. Needs to be a PR against develop though rather than master (we're roughly following the git-flow conventions where master is for releases: see http://nvie.com/posts/a-successful-git-branching-model/)

Ideally the API householder function should be backed by protocol: there are various special cases where a fast underlying implementation can potentially do something smarter than (I - 2vvt)

Guess we want a hermitian-transpose function too, right? And complex-conjugate while we're going down that path? At least the implementation is nice and easy for real vectors/matrices  :-)

[Friedrich Boeckh]

Concerning a householder function, I would assume the implementation is optimal, if the outer product vv^t is implemented optimal. I only see the corner case of 1x1 vectors at the moment. I think, the operation should be defined for all numeric types. The result type is either a rea or complex matrix.

For complex numbers, this would involve to switch to hermitian-transpose in order to be consistent in an algebraic sense, by the way. For practical reasons and technical consistency, the implementation of hermitian-transpose would however better rely on a second protocol in analogy to PTranspose. PHermitian would be definied for all number types and coerce with PTranspose for non-complex numbers.

I hope to be able to provide a proposal later this week. 

[Friedrich Boeckh]

Hi Mike,

this is a sketch of my thoughts so far:

Given an implementation of complex numbers, introducing a Hermitian transposition could follow the implementation steps of the ordinary matrix transposition:

1. Introduce a protocol PHermitian: https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix/protocols.clj#L375
2. Introduce a API function hermitian-transpose:
https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix.clj#L563
3. Provide a default implementation:
https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix/impl/default.clj#L323

Seems like to be easy task from this perspective.

For introducing complex numbers, I see 2 alternatives:

1. Provide an extension to java.lang.Number such that complex numbers could be treated as a first class number, recognized by 'number?', all math ops properly defined, etc.

2. Provide a new complex type "between" java.lang.Number and java.lang.Object, where the later is used to identify the different matrix implementations in the current code base.

There would also be the question of choosing a concrete implementation. Should an existing implementation, e.g. JBLAS, be used via Java interop, or should it be reengineered, build from scratch?

These questions made me deliberately think about the scope and priorities of core.matrix. Complex numbers have a low, nice-to-have priority with that respect in my opinion because matrix equations on complex numbers can typically be expressed in terms of separated real and imaginary parts.

I would prefer to see enforced progress of integrating other libraries. Beside of the addressed top priority of NDArray, I´d argue for JBLAS/improved Clatrix and Colt. While I personally appreciate your work on Vectorz and regard it to be sufficient for practical purposes, this approach would be the best answer to X-is-better-than-Y discussions in the community which probably keeps people away from using Clojure for interesting topics. Updating the benchmark overview (http://mikera.github.io/matrix-api/bench.html) in sync with every progress would be top priority from this perspective.

[Mike Anderson]

On Thursday, 12 September 2013 17:52:53 UTC+8, Friedrich Boeckh wrote:

Hi Mike,

this is a sketch of my thoughts so far:

Given an implementation of complex numbers, introducing a Hermitian transposition could follow the implementation steps of the ordinary matrix transposition:

1. Introduce a protocol PHermitian: https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix/protocols.clj#L375
2. Introduce a API function hermitian-transpose:
https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix.clj#L563
3. Provide a default implementation:
https://github.com/mikera/matrix-api/blob/master/src/main/clojure/clojure/core/matrix/impl/default.clj#L323

Yes that sounds exactly right. 

At the same time, it probably makes sense to add a "complex-conjugate" API function and protocol. This would obviously just be identity for real values or matrices, but would be needed for complex implementations and it also would make it easy to implement the default hermitian transpose as (comp transpose complex-conjugate)

 

Seems like to be easy task from this perspective.

For introducing complex numbers, I see 2 alternatives:

1. Provide an extension to java.lang.Number such that complex numbers could be treated as a first class number, recognized by 'number?', all math ops properly defined, etc.

Don't think this works: I think the API for java.lang.Number assumes real values. We'd be violating Liskov Substitution principle if we tried to squeeze in complex values, I think?

 

2. Provide a new complex type "between" java.lang.Number and java.lang.Object, where the later is used to identify the different matrix implementations in the current code base.

There would also be the question of choosing a concrete implementation. Should an existing implementation, e.g. JBLAS, be used via Java interop, or should it be reengineered, build from scratch?

I'm hoping that we can define the numeric protocols in a way that extend naturally to complex number implementations.

 

These questions made me deliberately think about the scope and priorities of core.matrix. Complex numbers have a low, nice-to-have priority with that respect in my opinion because matrix equations on complex numbers can typically be expressed in terms of separated real and imaginary parts.

Agreed. It's mostly a nice to have right now. But I think is also good to bear in mind because it forces us to get the API design logically clear. Also it means that Clojure+core.matrix may ultimately become much more appealing to segments of the scientific community, which I think is a good thing.

 

I would prefer to see enforced progress of integrating other libraries. Beside of the addressed top priority of NDArray, I´d argue for JBLAS/improved Clatrix and Colt. While I personally appreciate your work on Vectorz and regard it to be sufficient for practical purposes, this approach would be the best answer to X-is-better-than-Y discussions in the community which probably keeps people away from using Clojure for interesting topics. Updating the benchmark overview (http://mikera.github.io/matrix-api/bench.html) in sync with every progress would be top priority from this perspective.

Well - the whole point of core.matrix is to enable many implementations! Also core.matrix means that you often can mix and match implementations fairly easily. I'd love to see more work on all these. 

It's mostly just dependent on people volunteering their time and energy to make it happen.

[Friedrich Boeckh]

Fair point. I´ll check my options with regard to JBLAS.
I´d keep the habit to outline first the approach to you before starting to hack, so I hope to come back soon.