Fair enough. However you can see how this inconsistency will have people scratching their heads, as can be seen in the OP. The sizes will contradict the types…
Please do propose any documentation improvements you think would make this clearer.
Not really. Vectors and column matrices are behaviorally identified in the same way that adjoint vectors and row matrices are. There are also many situations where trailing singleton dimensions can behave as if they don’t exist or as if more of them exist than actually do. What is the actual problem here?
size(a)
(1, 1)
size(b)
(1, 2)
Clearly the objects are compatible for matrix-matrix multiplication (by size), however the types say they are not compatible.
EDIT: Sorry.
a = transpose([-6.5])
b = [0, .95]'
size(a)
(1, 1)
size(b)
(1, 2)
I’ve lost track of which a
and b
definitions you’re referring to.
I agree that should work; this just seems like a missing method bug to me. It’s a bit questionable from the perspective of variance analysis, if we were to keep track of up and down dimensions and/or covariant and contravariant dimensions, but we don’t do that in general, so I don’t think we should do it here either, so this should “just work”. It’s a little weird since b'*a'
doesn’t work, but I think that’s the consequence of how we’ve chosen to identify vector transposes with row matrices.
I might disagree. I don’t think a*b
or b'*a'
should work: vector-vector or row vector-row vector products don’t make sense. What bothers me is the disconnect between the sizes and the types. The sizes provide hints inconsistent with the behaviour based on types.
I am not sure why people use transpose
in general, but since the definition is for AbstractMatrix{<:Number}
(see, permutedims
for non-numeric), I would hope one could use these for matrix multiplication with a close relation to the mathematical concept, especially with the size
suggesting these work as such). A second best, would be to have to call Matrix(transpose(::AbstractMatrix{<:Real}))
to get that behavior. For Factorization
there is a good reason behind it, but here I fail to see why it should just work.
I think we should tread carefully here.
Right now if you do
a = rand(4,4)
b = rand(4)
a*b
You are doing {a^{i}}_{j}b^{j}. Note that the order of oprations determines the contractions. This is why you currently cannot do b*a
which would give b^{i}{a^{i}}_{j}. Meanwhile b'*a
does b_{i}{a^{i}}_{j}.
In this particular case we have two objects with the index structure {a^{i}}_{j}, {b^{i}}_{\alpha} where the latin and greek indices have different numbers of dimensions. So b'*a'
would give you {b_{\alpha}}^{i}{a_{i}}^{j}, so yes, that should work.
Right now I think everything is correct in the general metric signature case. It would be nice if we can make sure that we do not specialize to the Euclidean metric signature case as this might make it harder to implement other metric signatures.