operations on Mat.empty

[nilsbecker]

as a corner case, i happened on gemm Mat.empty Mat.empty which i expected to yield Mat.empty. however, instead this resulted in an exception from blas. just wanted to point this out in case it is not intentional. one could of course test explicitly for empty matrices. i guess it's a tradeoff between convenience and potentially some overhead if Mat.empty were handled by all functions where it makes sense.

[nilsbecker]

interestingly, Mat.add does work as expected, as do other unop and binop functions. i conclude that only gemm is so low-level that it does not contain any checking.
in that case, would it make sense to provide a matrix multiplication which is a binop? (i have not found one?)

[mmottl]

BLAS semantics is sometimes not exactly well-thought-out when it comes to corner cases. LACAML should not bomb with a Fortran exception here, of course (please keep reporting any such occurrences). But I don't think following the exact BLAS semantics makes sense either. BLAS allows M and N to be 0, but the result matrix must have a leading dimension of 1. This runs counter to the intuition that gemm Mat.empty Mat.empty should return Mat.empty.

I'm not sure what "sane" behavior should be. I'm actually leaning towards just raising exceptions for empty arguments when BLAS/LAPACK would bomb. Empty matrices can be strange anyway. E.g. gemm (Mat.create 10 0) (Mat.create 0 5) returns a 10x5 matrix even though neither argument actually contains data. What should the result matrix contain? Zeros? Uninitialized data? I think users should not be calling functions on empty data. This almost surely indicates a bug in the user code.

Note that Mat.add, etc., are not part of BLAS/LAPACK, but were implemented in C within LACAML. I don't think we should add a binop for matrix multiplication. There are important performance considerations here, e.g. gemm is noticeably faster when one (and only one) argument is transposed.

@Chris00 what do you think about empty matrix arguments? Any preference for how they should be handled?

[nilsbecker]

I think users should not be calling functions on empty data. This almost surely indicates a bug in the user code.

I think you're probably right in most cases. My case was like this: I construct a block matrix from four input matrices by first allocating it, then looping over the inputs and filling rows and column of the block matrix. I was pleased to see that this actually worked out of the box for Mat.empty inputs. What then did not work anymore was taking gemm a b where both a and b were empty, as one of the input blocks. I am now using an option type to handle that.

Its' interesting to think about what a dimension of 0 should behave like for best consistency. Thinking about matrices as maps on vector spaces, I guess the corresponding concept (say on the real numbers) would be the 0-dimensional space R^0 = {\vec{0}} consisting only of the zero vector, which would be Vec.create 0. R^0 can also be naturally identified with the 0 element in any of the other R^n. Then Mat.create 0 2 would be a map from R^2 to R^0, while Mat.make0 1 2 would be a map from R^2 to R^0 embedded in R^1, which is identified with Vec.make0 1. Conversely any map from R^0 to R^n must have the n-dim zero vector as its image since linear maps map zero to zero, and there is nothing else in the domain to map. If this is correct then gemv (Mat.create 10 0) (Vec.create 0) = Vec.make0 10 should hold. Similarly gemm (Mat.create n 0) (Mat.create 0 m) = Mat.make0 n m. I think this could be a consistent interpretation - maybe even the only one? Not sure it it's practical.

I wasn't aware that the binops were separately implemented, I thought they were wrappers. In that case maybe a gemm binop in C is not the way to go.

[Chris00]

For gemm (Mat.create 10 0) (Mat.create 0 5), mathematically speaking the result should be filled with 0 because that's what empty sums gives:

 ∑  aᵢⱼ bⱼₖ = 0 
j∈∅

The problem is that there is no such thing as the empty matrix—dimensions matter! Users would be equally surprised if multiplying a m×n matrix by a n×p one did not give one of size m×p.

I only see two ways: either empty matrices are considered a programming error—but that means more work to circumvent that limitation when one actually needs empty matrices (and Mat.empty would be seen as a matrix to use to please the type system that would always raise an exception when used in an operation)—or they are handled in a correct fashion from a mathematical viewpoint—in which case gemm Mat.empty Mat.empty should indeed return Mat.empty. I am not sure I used empty matrices for anything useful (besides having a consistent interface or returning a matrix I do not care about).

I do not agree that 0 ∈ ℝ⁰ should be identified with 0 ∈ ℝᵈ for d > 0 (which would imply that Vec.create0 0 can be used in place of Vec.create0 n) but, in any case, there is no need to do that to deal with empty matrices (note the plural). (Consistently with what I said above, Mat.make0 1 2 is the matrix of a linear map from ℝ² to ℝ¹, not from ℝ² to ℝ⁰.) Note that, in this context, we do not even identify ℝ¹ with ℝ (although they enjoy a very natural isomorphism) because, in OCaml, the representation for the two is completely different. One can see this as similar to the non-identification of ℤ (int) to a subset of ℝ (float) although in maths, we do not hesitate to write ℤ ⊆ ℝ.

[mmottl]

Thanks for your take on this problem, Chris. I agree that zeroing out the result of multiplying two empty matrices with non-zero dimensions seems like the mathematically right thing to do. In that respect BLAS is also not quite consistent, e.g. gemm ~k:0 (Mat.make 10 5) (Mat.make 5 3) does not zero out the result though K = 0 probably should.

The question really is whether we want to give users the freedom to hang themselves. It is probably true that most such occurrences in real code would constitute a bug. My only use case for Mat.empty so far was as a placeholder in e.g. still uninitialized references, etc.

I think for the while being just raising exceptions is the best solution. Users can always work around this by checking dimensions themselves. Though this is slightly less convenient, it also forces the user to be explicit about the assumption that empty matrices do not constitute an error. Should we ever decide to be more permissive, switching from exceptions to acceptance would not cause issues. The solution is also much easier to implement.

Just let me know if you find any objections to this approach.

[nilsbecker]

I do not agree that 0 ∈ ℝ⁰ should be identified with 0 ∈ ℝᵈ for d > 0 (which would imply that Vec.create0 0 can be used in place of Vec.create0 n)

I agree of course that Vec.create0 0 should not be equivalent to Vec.create0 n. This idea of identifying the zero elements was probably nonsense, and in any case is unnecessary. The meaningful point is this: Mat.empty is a (the) linear map from R^0 to R^0, that Mat.create 0 n is a map from R^n to R^0, and Mat.create n 0 is a map from R^0 to R^n whose entire image is 0 in R^n - this would require that gemv (Mat.create n 0) (Vec.create 0) = Vec.make0 n -- this empty sum gives 0 - but in the right space.

so i still think that this would be a mathematically consistent story - but of course that does not mean it is the best way to handle empty matrices in lacaml.

[Chris00]

On 2017-01-17, Markus Mottl wrote: I think for the while being just raising exceptions is the best solution. Users can always work around this by checking dimensions themselves. Though this is slightly less convenient, it also forces the user to be explicit about the assumption that empty matrices do not constitute an error. Should we ever decide to be more permissive, switching from exceptions to acceptance would not cause issues. The solution is also much easier to implement.

I agree that empty vectors & matrices are usually due to programming bugs. I’d like to suggest: - We may want to consult the list to have the opinions about the usage of others (I guess it will likely comfort us in our opinion but maybe better make sure). - Mat.create and Vec.create should not allow to create matrices and vectors with (at least one) null dimension. Only Mat.empty and Vec.empty vcan by used. (Of course, one can always to it through Bigarray but Lacaml should be consistent.) - Add Mat.is_empty and Vec.is_empty checking that at least one dimension is 0.

[mmottl]

• I'm not sure it would help much to ask on the list. This would likely only affect a tiny number of people. Given that nobody has ever complained until now in years, they probably wouldn't care anyway. I'm sure everybody will agree that the program shouldn't exit though. Note that the "mathematically right way", which will almost surely never be used in real code, would also require considerably more implementation effort.

• There might be bizarre future cases where empty matrices make sense so I wouldn't prohibit their creation. We cannot prohibit such a creation anyway if it happens outside of Lacaml so we'd gain little.

• I have added functions for checking for dimensions of size 0, but instead of is_empty I called themhas_zero_dim. Otherwise people might confuse this function with a comparison to the value Vec/Mat.empty. The name makes it clearer what the checked property is.

I'm currently adding checks for empty matrices, but only in places where the function would otherwise exit the program. That's an easy, incremental improvement that won't prevent other solutions in the future. I hope to push the changes tomorrow.

[Chris00]

On 2017-01-18, Markus Mottl wrote: * There might be bizarre future cases where empty matrices make sense so I wouldn't prohibit their creation. We cannot prohibit such a creation anyway if it happens outside of Lacaml so we'd gain little.

We would gain an additional check consistent with our take on the matter.

[mmottl]

It won't prevent us from having to raise exceptions from BLAS/LAPACK functions to avoid exiting the program when a user passes empty matrices. I'll implement this inevitable last step. We can decide on any additional checks later. Note that Mat.empty already exists anyway and is useful as a placeholder.

[Chris00]

On 2017-01-18, Markus Mottl wrote: It won't prevent us from having to raise exceptions from BLAS/LAPACK functions to avoid exiting the program when a user passes empty matrices.

Sure. It was more in the spirit “if it has to fail, let it fail early”.

I'll implement this inevitable last step. We can decide on any additional checks later. Note that Mat.empty already exists anyway and is useful as a placeholder.

Exactly. Hence my suggestion not to allow empty dimensions for (Mat|Vec).create

[mmottl]

The "empty" checks are in the distribution now. I'll close this issue for now. Should anybody absolutely need the mathematically saner but not quite BLAS-compatible semantics, they'd likely need to submit significant changes. Probably not worth the trouble.

[nilsbecker]

I just saw this new feature in ocaml 4.05 Bigarray: zero dimensions in bigarrays.
ocaml/ocaml#787
i thought i'd post this here for reference -- it sounds like there are possible interactions...

[mmottl]

I don't think this would affect Lacaml at this point, because zero-dimensional tensors (i.e. really inefficient but potentially useful representations for scalars) are not used anywhere.

[nilsbecker]

yes, after posting i realized that they mean rank 0 arrays, not dimension zero in a matrix. sorry. good that it's not a problem.