Make it more like StaticArrays.jl

[KeitaNakamura]

No description provided.

[codecov]

Codecov Report

Attention: Patch coverage is 93.33333% with 1 line in your changes missing coverage. Please review.

Project coverage is 93.73%. Comparing base (09b5860) to head (2a7a3f7).
Report is 2 commits behind head on main.

Files with missing lines	Patch %	Lines
src/Space.jl	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #222      +/-   ##
==========================================
+ Coverage   93.68%   93.73%   +0.05%     
==========================================
  Files          18       18              
  Lines        1583     1581       -2     
==========================================
- Hits         1483     1482       -1     
+ Misses        100       99       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[KeitaNakamura]

Hello, @ExpandingMan. Based on your comment in #210, I’m experimentally considering changing the package to be more consistent with StaticArrays.jl, making it easier to switch from StaticArrays.jl to Tensorial.jl. To do this, I plan to:

1. Revert the definition of multiplication between tensors
2. Revert the definition of the dot function
3. Define single contraction using a different operator

However, (3) is difficult because similar-looking operators, like • (\bullet), cannot currently be used in Julia. Do you have any suggestions for addressing this?

[ExpandingMan]

Cool! This is an exciting development, I think this is very worth doing. StaticArrays.jl has a lot of its own problems, so I think it's possible to go too far in emulating it, but on the other hand I think maintaining consistency with standard conventions of the Julia AbstractArray interface is a really big deal, so I would view this as being more about Base.AbstractArray than about StaticArrays per say. (Like I had commented before, I've been a bit afraid of this package because I did not have confidence that swapping between it and StaticArrays in existing code would be straightforward.)

Let me re-familiarize myself with the API here and see if I can make an educated suggestion.

[ExpandingMan]

I'm not entirely clear on all the changes that are implicit in this PR, but I'll give an outline of where the current version is incompatible with standard conventions of AbstractArray and where that might be relevant here. Again, I don't think StaticArrays in particular is relevant here, StaticArrays conforms to these same conventions but they are not set by it, rather they are set by the base language (which is the reason why this package breaking those conventions is potentially problematic).

First, and perhaps most importantly, * must perform the following operations, depending on arguments::

$$C_{ik} = A_{ij} B_{jk}$$ $$b_i = A_{ij} a_j$$

as well as for scalars (which Tensorial already does). As far as I know, Base and popular packages such as StaticArrays do not define any * operations for higher rank objects, however, extrapolating on this pattern I would argue that * should always contract the last index of the first argument with the first index of the second argument. This is currently the single thing that most badly breaks interoperability between Tensorial.jl and other array packages at present.

⋅ is more ambiguous, as far as I know the only methods are currently for vectors. I would say that it is implied that this is an inner product, so it should return a Number, but I don't think anything actually enforces that, and it's not clear to me whether anything would break if it were left alone (though if you take the previous suggestion and not this one, both operations will do the same thing).

As for single contractions, I'm a bit surprised that there aren't already methods for contracting individual indices apart from @einsum, for example I'd have expected something like contract(A, B, Val(1), Val(2)) to do A_{ij}B_{ki}. I do think it would be good to have a lower-level interface for some of the operations which (it appears to me) are only now possible with @einsum. This also makes me a bit puzzled by your comment since I had thought that if you would

revert the definition of multiplication

this would mean that * is now a single contraction as I've outlined above? Or am I reading that wrong?

[KeitaNakamura]

Thanks for your comments.

First, and perhaps most importantly, * must perform the following operations, depending on arguments::

$$C_{ik} = A_{ij} B_{jk}$$ $$b_i = A_{ij} a_j$$

as well as for scalars (which Tensorial already does). As far as I know, Base and popular packages such as StaticArrays do not define any * operations for higher rank objects, however, extrapolating on this pattern I would argue that * should always contract the last index of the first argument with the first index of the second argument. This is currently the single thing that most badly breaks interoperability between Tensorial.jl and other array packages at present.

...

This also makes me a bit puzzled by your comment since I had thought that if you would

revert the definition of multiplication
this would mean that * is now a single contraction as I've outlined above? Or am I reading that wrong?

Yes, this PR simply restores the multiplication as follows:

julia> A = rand(Mat{3,2})
3×2 Mat{3, 2, Float64, 6}:
 0.265676  0.383535
 0.525552  0.664998
 0.27495   0.77021

julia> B = rand(Mat{2,3})
2×3 Mat{2, 3, Float64, 6}:
 0.66506   0.275736  0.206113
 0.145932  0.279541  0.616196

julia> A * B
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 0.23266   0.18047   0.291092
 0.446568  0.330808  0.518092
 0.295257  0.291119  0.531271

julia> a = rand(Vec{2})
2-element Vec{2, Float64}:
 0.24376793178856293
 0.8046192769558261

julia> A * a
3-element Vec{3, Float64}:
 0.3733630609973071
 0.6631827209074246
 0.6867495652822495

I also think that defining * as a single contraction makes sense, but in this case, a * a (where a is a vector) would need to return a scalar, which could cause some confusion. Additionally, operations like a * a' and a' * a might be difficult to interpret as tensor calculations. So, it might be better to keep * as regular multiplication between matrices/vectors, similar to how Array handles it.

⋅ is more ambiguous, as far as I know the only methods are currently for vectors. I would say that it is implied that this is an inner product, so it should return a Number, but I don't think anything actually enforces that, and it's not clear to me whether anything would break if it were left alone (though if you take the previous suggestion and not this one, both operations will do the same thing).

In my field, ⋅ is commonly understood as a single contraction without causing confusion. For instance, both Tensors.jl and Gridap.jl define single contraction using ⋅. So, if people in other fields don’t find this problematic, I would prefer to keep it as is.

As for single contractions, I'm a bit surprised that there aren't already methods for contracting individual indices apart from @einsum, for example I'd have expected something like contract(A, B, Val(1), Val(2)) to do A_{ij}B_{ki}. I do think it would be good to have a lower-level interface for some of the operations which (it appears to me) are only now possible with @einsum.

The lower-level interface has already been implemented to achieve the functionality of @einsum, so I will export it. Thanks.

[ExpandingMan]

I also think that defining * as a single contraction makes sense, but in this case, a * a (where a is a vector) would need to return a scalar

That's a good point. Anyway, I think it's fine, it's clearly absolutely required for matrix-vector and matrix-matrix to satisfy the interface, but if you leave the rest undefined, it's fine, I don't think any other package defines any other cases.

In my field, ⋅ is commonly understood as a single contraction without causing confusion.

I seem to remember having this conversation when I opened the initial issue, and at first I thought it was strange when you said that, but then I remembered this convention sometimes being used in literature I'm familiar with. I'm most familiar with tensors from general relativity and quantum field theory, which tend to use either component or abstract index notation nearly 100% of the time, so when I see $\cdot$ I tend to think it's ambiguous. Regardless, as I said the only thing that would be a clear requirement is for vector inner products, this is already used in many packages.

If you'd like, I can text out the changes you made here and verify that StaticArray (or indeed Array) can now be swapped for Tensorial tensors without breaking, but it might take me a day or two to get to it.

[ExpandingMan]

I've gone in and checked out this branch. I have done various tests of compatibility with AbstractArray, basically, for all operations which are defined for AbstractArray, I'm checking op(A,B) == op(Array(A), Array(B)). This looks good! I believe in this PR everything is compatible.

This package is really great, by the way, I'm excited for this. StaticArrays uses a horrible hack called MArray which makes it a lot more flexible than this in many cases, but it is also the cause of a huge amount of pain and headaches, and I think there are a huge number of obvious use cases where it simply is not needed. So not only does this package have lots of really great tensor functionality that's missing from StaticArrays but it is likely to be significantly easier to optimize.

[KeitaNakamura]

Thank you very much for checking the compatibility with AbstractArray. I'm glad to hear that the compatibility seems good.