N-dim array operations
37 views
Skip to first unread message
Francesco Bonazzi
Feb 4, 2016, 3:08:46 PM2/4/16
to sympy
Hi there,
ISSUE 1:
I had a discussion with Kalevi Suominen about two operations I introduced in my PR:
https://github.com/sympy/sympy/pull/10511
Tensor product acts on two N-dim arrays A[m1, m2, ... ] and B[n1, n2, ... ] to give A[m1, m2, ... ]*B[n1, n2, ... ] (that is an N-dim array of rank the sum of the two original ranks).
Tensor contraction sums the indices given by the specified positions: A[m1, ... , mk, ... , mp , ... ] ===> Sum(A[m1, ... , i , ... , i , ... ], i)
The point to decide is wheter tensorproduct( ) and tensorcontraction( ) should be class methods or functions in the module namespace.
My point for the public namespace functions is that those operations shall also act on scalars. Scalars can be thought of as zero-dimensional N-dim arrays. Every SymPy Expr is a scalar (maybe with some minor exceptions). A function in public namespace can be called as tensorproduct(a, b) even if a is a scalar (say a Symbol or number), while a.tensorproduct(b) would fail because only N-dim arrays would be associated the method.
On the other hand, for type stability reasons one could wish to represent scalars still as array objects, just with no shape defined. They would have all tensorial methods defined. I don't like this idea because it looks more natural to just pick a SymPy object an multiply it by an array. Type stability could be the pro of this point.
Read the discussion for details.
Wolfram Mathematica does it my way:
https://reference.wolfram.com/language/ref/TensorProduct.html?q=TensorProduct
https://reference.wolfram.com/language/ref/TensorContract.html?q=TensorContract
(though it may be argued that Mathematica language has no syntax like a.method(b) )
What do you think? Anyone supporting my API choice?
Can we solve this issue quickly before the release process for version 1.0 starts?
ISSUE 2:
In Wolfram Mathematica it is possible to derive an object by an array of symbols:
https://reference.wolfram.com/language/ref/D.html?q=D
In their documentation:
D[f,{{x1,x2,…}}] for a scalar f gives the vector derivative
.
Basically, if f is an N-dim array (scalar, vector, matrix or higher rank), it raises the array rank by one by deriving elementwise f by each {x1, x2, ... }.
What about introducing this to SymPy? I could quickly write this routine for the newly created N-dim array module.
If these two issues are solved, it will get much more easy to use SymPy to deal with modern physics, such as General Relativity. I counted at least 5 PR in the last two years by minor contributors trying to introduce some sort of tensorial stuff applied to relativity.
ISSUE 1:
I had a discussion with Kalevi Suominen about two operations I introduced in my PR:
https://github.com/sympy/sympy/pull/10511
Tensor product acts on two N-dim arrays A[m1, m2, ... ] and B[n1, n2, ... ] to give A[m1, m2, ... ]*B[n1, n2, ... ] (that is an N-dim array of rank the sum of the two original ranks).
Tensor contraction sums the indices given by the specified positions: A[m1, ... , mk, ... , mp , ... ] ===> Sum(A[m1, ... , i , ... , i , ... ], i)
The point to decide is wheter tensorproduct( ) and tensorcontraction( ) should be class methods or functions in the module namespace.
My point for the public namespace functions is that those operations shall also act on scalars. Scalars can be thought of as zero-dimensional N-dim arrays. Every SymPy Expr is a scalar (maybe with some minor exceptions). A function in public namespace can be called as tensorproduct(a, b) even if a is a scalar (say a Symbol or number), while a.tensorproduct(b) would fail because only N-dim arrays would be associated the method.
On the other hand, for type stability reasons one could wish to represent scalars still as array objects, just with no shape defined. They would have all tensorial methods defined. I don't like this idea because it looks more natural to just pick a SymPy object an multiply it by an array. Type stability could be the pro of this point.
Read the discussion for details.
Wolfram Mathematica does it my way:
https://reference.wolfram.com/language/ref/TensorProduct.html?q=TensorProduct
https://reference.wolfram.com/language/ref/TensorContract.html?q=TensorContract
(though it may be argued that Mathematica language has no syntax like a.method(b) )
What do you think? Anyone supporting my API choice?
Can we solve this issue quickly before the release process for version 1.0 starts?
ISSUE 2:
In Wolfram Mathematica it is possible to derive an object by an array of symbols:
https://reference.wolfram.com/language/ref/D.html?q=D
In their documentation:
D[f,{{x1,x2,…}}] for a scalar f gives the vector derivative
. Basically, if f is an N-dim array (scalar, vector, matrix or higher rank), it raises the array rank by one by deriving elementwise f by each {x1, x2, ... }.
What about introducing this to SymPy? I could quickly write this routine for the newly created N-dim array module.
If these two issues are solved, it will get much more easy to use SymPy to deal with modern physics, such as General Relativity. I counted at least 5 PR in the last two years by minor contributors trying to introduce some sort of tensorial stuff applied to relativity.
Aaron Meurer
Feb 4, 2016, 3:26:10 PM2/4/16
to sy...@googlegroups.com
On Thu, Feb 4, 2016 at 3:08 PM, Francesco Bonazzi <franz....@gmail.com> wrote:
Hi there,
ISSUE 1:
I had a discussion with Kalevi Suominen about two operations I introduced in my PR:
https://github.com/sympy/sympy/pull/10511
Tensor product acts on two N-dim arrays A[m1, m2, ... ] and B[n1, n2, ... ] to give A[m1, m2, ... ]*B[n1, n2, ... ] (that is an N-dim array of rank the sum of the two original ranks).
Tensor contraction sums the indices given by the specified positions: A[m1, ... , mk, ... , mp , ... ] ===> Sum(A[m1, ... , i , ... , i , ... ], i)
The point to decide is wheter tensorproduct( ) and tensorcontraction( ) should be class methods or functions in the module namespace.
My point for the public namespace functions is that those operations shall also act on scalars. Scalars can be thought of as zero-dimensional N-dim arrays. Every SymPy Expr is a scalar (maybe with some minor exceptions). A function in public namespace can be called as tensorproduct(a, b) even if a is a scalar (say a Symbol or number), while a.tensorproduct(b) would fail because only N-dim arrays would be associated the method.
The function approach sounds fine, but to me, it seems that you would mainly want to override __mul__ and use that.
On the other hand, for type stability reasons one could wish to represent scalars still as array objects, just with no shape defined. They would have all tensorial methods defined. I don't like this idea because it looks more natural to just pick a SymPy object an multiply it by an array. Type stability could be the pro of this point.
Read the discussion for details.
Wolfram Mathematica does it my way:
https://reference.wolfram.com/language/ref/TensorProduct.html?q=TensorProduct
https://reference.wolfram.com/language/ref/TensorContract.html?q=TensorContract
(though it may be argued that Mathematica language has no syntax like a.method(b) )
What do you think? Anyone supporting my API choice?
Can we solve this issue quickly before the release process for version 1.0 starts?
Is this going to be an API change from what's in master. I plan to start the release branch *very soon*.
ISSUE 2:
In Wolfram Mathematica it is possible to derive an object by an array of symbols:
https://reference.wolfram.com/language/ref/D.html?q=D
In their documentation:
D[f,{{x1,x2,…}}] for a scalar f gives the vector derivative
Basically, if f is an N-dim array (scalar, vector, matrix or higher rank), it raises the array rank by one by deriving elementwise f by each {x1, x2, ... }.
What about introducing this to SymPy? I could quickly write this routine for the newly created N-dim array module.
If these two issues are solved, it will get much more easy to use SymPy to deal with modern physics, such as General Relativity. I counted at least 5 PR in the last two years by minor contributors trying to introduce some sort of tensorial stuff applied to relativity.
--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b9d5e371-f858-413e-b466-4a1bcc810732%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Francesco Bonazzi
Feb 4, 2016, 3:47:46 PM2/4/16
to sympy
On Thursday, 4 February 2016 21:26:10 UTC+1, Aaron Meurer wrote:
The function approach sounds fine, but to me, it seems that you would mainly want to override __mul__ and use that.
Not sure about __mul__, when you apply __mul__ to two matrices you are performing a tensor product followed by a contraction of the central indices, that is
- A[i, j], B[m, l] ===> A[i, j]*B[m, l] (got a rank-4 array, matrix is a rank-2 array)
- A[i, j]*B[m, l] ===> Sum(A[i, k]*B[k, l], k) (got back a rank-2 array by contracting indices j and m)
I don't want to override __mul__ because the behavior with matrices is different.
Can we solve this issue quickly before the release process for version 1.0 starts?Is this going to be an API change from what's in master. I plan to start the release branch *very soon*.
No API changes. Just a small bug fix in the PR (unrelated to new functionality).
ISSUE 2:
Had a quick look at that matrix cookbook in the link. The examples I verified appear to be the same of my proposal, and to Wolfram Mathematica's API, but better doublecheck.
Anyway, I would temporarily add it to the N-dim array only. Extending this functionality to matrices should be procrastinated.
Aaron Meurer
Feb 4, 2016, 3:53:33 PM2/4/16
to sy...@googlegroups.com
On Thu, Feb 4, 2016 at 3:47 PM, Francesco Bonazzi <franz....@gmail.com> wrote:
On Thursday, 4 February 2016 21:26:10 UTC+1, Aaron Meurer wrote:The function approach sounds fine, but to me, it seems that you would mainly want to override __mul__ and use that.
Not sure about __mul__, when you apply __mul__ to two matrices you are performing a tensor product followed by a contraction of the central indices, that is
- A[i, j], B[m, l] ===> A[i, j]*B[m, l] (got a rank-4 array, matrix is a rank-2 array)
- A[i, j]*B[m, l] ===> Sum(A[i, k]*B[k, l], k) (got back a rank-2 array by contracting indices j and m)
I don't want to override __mul__ because the behavior with matrices is different.
It's worth pointing out that Python 3.5 has the new __matmul__ operator (A@B). You're the best one to decide what these operators should do, and whether not defining them is too much of an inconvenience. I would warn against using non-standard operators (like the logic operators &, |, and ^) for multiplication, since they have different precedence rules with other operators (like +) than what you would generally expect.
Aaron Meurer
Can we solve this issue quickly before the release process for version 1.0 starts?Is this going to be an API change from what's in master. I plan to start the release branch *very soon*.
No API changes. Just a small bug fix in the PR (unrelated to new functionality).
ISSUE 2:Is this the same issue as https://github.com/sympy/sympy/issues/5858?
Had a quick look at that matrix cookbook in the link. The examples I verified appear to be the same of my proposal, and to Wolfram Mathematica's API, but better doublecheck.
Anyway, I would temporarily add it to the N-dim array only. Extending this functionality to matrices should be procrastinated.
--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/1e94aac1-6c07-4a24-9a3d-83e501bd4037%40googlegroups.com.
Francesco Bonazzi
Feb 4, 2016, 4:12:02 PM2/4/16
to sympy
It's worth pointing out that Python 3.5 has the new __matmul__ operator (A@B). You're the best one to decide what these operators should do, and whether not defining them is too much of an inconvenience. I would warn against using non-standard operators (like the logic operators &, |, and ^) for multiplication, since they have different precedence rules with other operators (like +) than what you would generally expect.
In the aforementioned PR there is no operator overloading at all, and probably there won't ever been.
Consider that multiplication among matrices or array-like operations is not uniquely defined, there are many kinds. To get a list from Mathematica's documentation:
- https://reference.wolfram.com/language/ref/TensorProduct.html
- https://reference.wolfram.com/language/ref/Outer.html
- https://reference.wolfram.com/language/ref/TensorWedge.html
- https://reference.wolfram.com/language/ref/KroneckerProduct.html
- https://reference.wolfram.com/language/ref/Inner.html
- https://reference.wolfram.com/language/ref/Dot.html
- https://reference.wolfram.com/language/ref/Cross.html
- https://reference.wolfram.com/language/ref/Times.html
In Mathematica only Times (elementwise multiplication) and Dot (matrix multiplication) have an operator form.
Maybe operators aren't even enough to cover all kinds of product.
I suggest that __matmul__ be the same as __mul__ for matrices in SymPy. Maybe __matmul__ could be the matrix multiplication for N-dim arrays.
By the way, I set __mul__ to raise an exception when multiplying two N-dim arrays.
Reply all
Reply to author
Forward
0 new messages