# Access Large Array with Preassigned Index

I wish to use a preassigned variable to index a large array many times. For example:

a = {i,j+1,k-1};
b = {i+1,j,k-2};


I know I can do this several ways but for conciseness in formula I want to use the double bracket assignment such as:

Ex[[ Sequence@@a ]] + Ey[[ Sequence@@b]]


In a large number of calculations however the sequence operation adds about 50% time overhead. Is there a way without the time cost overhead? thanks

The form of these two expressions are constructed elsewhere in the program

ind = {{i, j + 1, k}, {i - 1, j, k - 1}};


and need to be passed to a function, e.g.

f[in_] := Module[{a = in[[1]], b = in[[2]]},
Table[Ex[[Sequence @@ a]] + Ey[[Sequence @@ b]], {m, 500}, {i, 2, 8},
{j, 2, 8}, {k, 2, 8}]];


If I did not have the requirement to construct the index form somewhere else I could simply put:

g[]:= Table[
Ex[[i, j + 1, k]] + Ey[[i - 1, j, k - 1]], {m, 500}, {i, 2, 8}, {j,
2,8}, {k, 2, 8}]

Ex = RandomReal[{0, 10}, {10, 10, 10, 10}];
Ey = RandomReal[{0, 10}, {10, 10, 10, 10}];

AbsoluteTiming[f[ind];]


{0.651337, Null}

AbsoluteTiming[g[];]


{0.0892097, Null}

• You say "assignment" but I see only an addition. Is this a typo? – Henrik Schumacher Dec 5 '18 at 22:37

You may use Extract:

m = 100;
n = 10000;
Ex = RandomReal[{-1, 1}, {m, m, m}];
Ey = RandomReal[{-1, 1}, {m, m, m}];
a = RandomInteger[{1, m}, {n, 3}];
b = RandomInteger[{1, m}, {n, 3}];

r1 = Table[Ex[[Sequence @@ a[[i]]]] + Ey[[Sequence @@ b[[i]]]], {i, 1, Length[a]}]; // AbsoluteTiming // First
r2 = Extract[Ex, a] + Extract[Ey, b]; // AbsoluteTiming // First
r1 == r2


0.103876

0.000698

True

# Edit

In the meantime, I got to the conclusion that the reason for f being slower than g is that Sequence is not compilable. So let us perform the compilation by hand:

f2[in_] := Quiet@Block[{X, XX, Y, YY},
With[{XX = Part[X, Sequence @@ in[[1]]], YY = Part[Y, Sequence @@ in[[2]]]},
Compile[{{X, _Real, 4}, {Y, _Real, 4}},
Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
CompilationTarget -> "WVM"
];
cf[Ex, Ey]
]
];


Now we have:

aa = f[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
bb = g[]; // RepeatedTiming // First
cc = f2[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == bb == cc


0.692

0.090

0.0775

True

A further speed-up can be obtained by replacing Part by CompileGetElement:

f3[in_] := Block[{X, XX, Y, YY},
With[{XX = CompileGetElement[X, Sequence @@ in[[1]]], YY = CompileGetElement[Y, Sequence @@ in[[2]]]},
Compile[{{X, _Real, 4}, {Y, _Real, 4}},
Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
CompilationTarget -> "WVM",
RuntimeOptions -> "Speed"
];
cf[Ex, Ey]
]
];


Runtime test:

dd = f3[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == dd


0.028

True

Finally, let me say that handing over the code for the indices by the symbolic code {{i, j + 1, k}, {i - 1, j, k - 1}} is very fragile. A more robust way would be to submit functions that generate the symbolic code:

f4[in_] := Quiet@Block[{X, XX, Y, YY, i, j, k},
With[{XX = CompileGetElement[X, Sequence @@ in[[1]][i, j, k]], YY = CompileGetElement[Y, Sequence @@ in[[2]][i, j, k]]},
Compile[{{X, _Real, 4}, {Y, _Real, 4}},
Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
CompilationTarget -> "WVM",
RuntimeOptions -> "Speed"
];
cf[Ex, Ey]
]
];


Test:

ee = f4[{{i, j, k} \[Function] {i, j + 1, k}, {i, j, k} \[Function] {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == ee
`

0.028

True

• I have not explained myself well. The indices to be used are created – user8281 Dec 5 '18 at 23:12
• I have not explained myself well: – user8281 Dec 5 '18 at 23:12
• Then please clarify in your post. It may help to post the code that you actually used. – Henrik Schumacher Dec 5 '18 at 23:16
• Thank you for your detailed explanation. It is very helpful – user8281 Dec 9 '18 at 19:15
• You're welcome. – Henrik Schumacher Dec 9 '18 at 19:19