3

I have two lists, er and ns:

Dimensions@er
Dimensions@ns

gives:

{343, 3}
{343, 3996, 3}

I want to do a dot product such that in loops it is:

For [i = 1, i <= 343, i++
  For j = 1, j <= 3996, j++
    result[[i,j]] = er[[i]].ns[[i,j]]
  ]
]

So I was thinking:

MapThread[Function[{x, y}, ((x.#) & /@ y)], {er, ns}]

Because

(er[[1]].#)&/@ns[[1]]

works fine.

But this is not it. Also I would like to parallelize it, and I do not know how well MapThread is paralelizable.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
atapaka
  • 3,954
  • 13
  • 33
  • 1
    (Almost) never use For; use Do, Table, While etc. In this case I'd suggest ParallelDo (maybe with ParallelEvaluate; see here). – corey979 Oct 11 '16 at 18:19
  • 2
    @corey979 I do not use For, I simply mentioned it because it is absolutely obvious with that notation what I want to do. I am trying to do it with Map and others, it is written in the question though... – atapaka Oct 11 '16 at 18:21
  • Well, that wasn't entirely clear to me, and when I see a For I automatically discourage from using it. – corey979 Oct 11 '16 at 18:23
  • 1
    Why MapThread "is not it"? What's wrong with it? It's correct, and the fastest among the methods here. – corey979 Oct 11 '16 at 19:10
  • As a dot product: Transpose[er.Transpose[ns, {2, 3, 1}], {1, 1, 2}]. – J. M.'s missing motivation Mar 21 '18 at 02:47

2 Answers2

4

The straghtforward way is

ans1 = Table[er[[i]].ns[[i, j]], {i, 343}, {j, 3996}]

on my Macbook Pro this gives an AbsoluteTiming of 0.611 seconds. Using ParallelTable cuts that time in half if multiple kernels are already launched.

I always feel that using vectorized operations gives much bigger boosts to speed than parallellization does, so here is what I came up with:

ans2 = Block[{erCopies, erT, nsT},
 erCopies = ConstantArray[#, 3996] & /@ er;
 erT = Transpose[erCopies, {2, 3, 1}];
 nsT = Transpose[ns, {2, 3, 1}];
 Plus @@ (erT*nsT)
 ];

This uses more RAM than the straightforward way because it enlarges er to match the dimensions of ns, but is a bit faster for me, 0.177 seconds. Note that I could not afford to run these timings on a completely fresh kernel, so the timings may not be reproducible.

Marius Ladegård Meyer
  • 6,805
  • 1
  • 17
  • 26
1

First, some data:

er = RandomReal[1, {343, 3}];
ns = RandomReal[1, {343, 3996, 3}];

A compiled, auto-parallelized way, which is quite fast:

cf = Compile[{{v, _Real, 1}, {w, _Real, 1}},
  v.w,
  RuntimeAttributes -> {Listable}, Parallelization -> True
  ];

(result = cf[er, ns]) // Dimensions // AbsoluteTiming
(*  {0.04556, {343, 3996}}  *)

An array-based approach, which is somewhat faster than the OP's MapThread code:

(result2 = MapThread[Dot, {er, #}] & /@ Transpose[ns, {2, 1, 3}] // Transpose) //
      Dimensions // AbsoluteTiming
(*  {0.147916, {343, 3996}}  *)

result2 == result
(*  True  *)

OP's:

MapThread[Function[{x, y}, ((x.#) & /@ y)], {er, ns}] // Dimensions // AbsoluteTiming
(*  {0.376943, {343, 3996}}  *)

Marius's solutions:

Table[er[[i]].ns[[i, j]], {i, 343}, {j, 3996}] // Dimensions // AbsoluteTiming
(*  {0.855036, {343, 3996}}  *)

ans2 = Block[{erCopies, erT, nsT}, 
    erCopies = ConstantArray[#, 3996] & /@ er;
    erT = Transpose[erCopies, {2, 3, 1}];
    nsT = Transpose[ns, {2, 3, 1}];
    Plus @@ (erT*nsT)]; // AbsoluteTiming
(*  {0.149258, Null}  *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747