+/ . * is so much faster than (+/@:(*"1 _)) that even if you were given
a notation for selecting an axis, the implementation would do the
transpose so that it could use the faster product.
+/@:*"1 is pretty fast for small arguments but for big ones its cache
characteristics mean you end up limited by slower memory.
I have never done multiple-axis contractions, but I think the right
approach would be to transpose the axes to the end and combine them into
a single axis so that you could use +/ . * .
Your trick using i. to generate a transpose is a good way to make the
results intuitive.
Henry Rich
On 2/6/2020 5:32 PM, Marshall Lochbaum wrote:
Here is a way to permute the axes of an array that you might find easier
to work with:
A =. i. 2 3 4 5
$ ('ijkl'i.'ljik') |: A
5 3 2 4
The idea is that we start with A, which has axes ijkl, and want to
shuffle them around to get ljik. So we look up each of the letters ljik
in ijkl, which gives us a permutation 3 1 0 2 which should be applied to
the axes (you could also just write out the permutation by hand). This
is then used as the left argument to transpose.
Although a matrix product will execute faster, if I want to write tensor
code that is straightforward, and extends to multiple contracted axes, I
would transpose the contracted axes to the end and combine them with
(+/@:*"1). This pairs row vectors from the two arrays and takes their
dot product. For three contracted axes I might write (+/@:,@:*"3), which
differs in that it pairs up rank-3 subarrays, and ravels before summing,
since sum only operates on the first axis.
The idea with (|:~ 1 -.~ i.@(#@$)) is to get the list of indices for
every axis
i.@(#@$) A
0 1 2 3
and then remove a 1
1 -.~ i.@(#@$) A
0 2 3
to obtain a left argument to transpose. Transpose adds missing axes to
the front, so this will move axis 1 to the beginning, for an array of
any shape.
I think this is a weak point for J's array programming style. Just like
tacit code works well for small manipulations but gets confusing and
difficult as the number of arguments increases, J primitives like rank,
transpose, and inner product are excellent for a few axes, but don't
work so well when the number of axes increases and the ways they are
manipulated become more complicated. In J's defense, most programming
languages are pretty awful at this (how many times are you going to
write .map before you just give up?). But it's something that bothers
me. Sure, I can write a bunch of nested rank operators to do what I
want, but why is it so hard to just say which axes correspond to which?
Marshall
On Thu, Feb 06, 2020 at 09:44:32PM +0000, [email protected] wrote:
Henry Rich <[email protected]> writes:
It depends on what the '*' means in the definition of the product.
C_ijlmnop = sum_k A_ijkl * B_mknop
I really mean the usual multiplication (of numbers).
Have a look at
a =. i. 2 3 4 5
b =. i. 2 4 3 5 6
$ (2 |: a) +/ . * ((|:~ 1 -.~ i.@(#@$)) b)
2 3 5 2 3 5 6
Well, I hoped for seeing something simpler. :)
If I understand your example a bit, you sum along the 4-long axis, and
you make this axis the last in the case of 'a', ie. by (2 |: a), then
use the so called dot product (+/ . *), and then something I cannot
easily decipher, but I guess it must, among other things, make the
4-long axis (number 1) the first.
Does one really has to shuffle with the axes?
Does one really need to use the for me almost incomprehensible
conjunction . (which no one dares to explain simply enough...)?
Why is it so simple to say what I want mathematically, and so awkward
to do it programmatically?
Thanks!
Ruda
On 2/6/2020 2:46 PM, [email protected] wrote:
having two multidimensional matrices A and B,
with some indices, say, A_ijkl and B_mknop,
how can I obtain a matrix C, where
C_ijlmnop = sum_k A_ijkl * B_mknop
ie., C has all indices of A and B but for the index k,
which was summed over.
Thanks for your suggestions.
Best regards
Ruda
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
