1
$\begingroup$

I want to calculate the determinant along the last slice of a 3-dimensional array. So for I do this by slow the Table command. I know that for time reasons I should use Map or Apply, however couldn't successful solve the problem.

m = 2; n = 3; o = 10;
SeedRandom[0];
x = RandomReal[{-1, 1}, {m, n, o}];
Table[Sqrt[Det[x[[;;,;;,i]].Transpose[x[[;;,;;,i]]]]],{i,1,o}]
Map[Sqrt[Det[# . Transpose[#]]] &, {x}, {2}]
(*{1.04663,0.437045,0.479911,0.260814,0.205563,0.171896,1.20112,1.00502,0.855893,0.125758}*)
(*{{5.38795,4.19589}}*)
$\endgroup$
1
  • $\begingroup$ Using TreeForm might help with understanding the difference between the Map[...,{2}] and the transpose solution given below. $\endgroup$ Oct 21, 2022 at 10:17

2 Answers 2

3
$\begingroup$

You must rearrange the matrix so that the blocks you want the Det from, appear at level 1:

Consider the first block or matrix:

x[[All, All, 1]] // MatrixForm

enter image description here

We can rearrange x so that these matrixes appear at level 1:

Transpose[x, {2, 3, 1}][[1]] // MatrixForm

enter image description here

With this we can calculate the Det's using Map:

Map[Sqrt[Det[# . Transpose[#]]] &, Transpose[x, {2, 3, 1}], {1}]

(* {1.04663, 0.437045, 0.479911, 0.260814, 0.205563, 0.171896, 1.20112, \
1.00502, 0.855893, 0.125758} *)
$\endgroup$
1
  • $\begingroup$ Map is here only 3-4 times faster than Table (tested for 1 million samples). However transposing the whole array before we use Map or taking the Sqrt after Mapdoes not improve the speed. Is there a faster way to calculate this? $\endgroup$ Oct 21, 2022 at 2:03
4
$\begingroup$

Update

I was curious to see what the speed would be if one used the neural network framework to do the computation. It turns out that it is faster than using Map and the computation is particularly fast with a GPU. However, error accumulates with the number of components

The network will be made of a single function layer reproducing the code of @DanielHuber but within the framework of neural networks (without much of a network as there will not be any weights or complicated structure).

The tensor/array on which the function will operate has the same structure as in the question but with more components :

m = 2; n = 3; o = 10^6;
  SeedRandom[0];
  x = RandomReal[{-1, 1}, {m, n, o}];

The function using a neural network layer and higher order operator.

f = FunctionLayer[Sqrt@Det[# . Transpose[#]] &];
mapf = NetMapOperator[f];

Comparison between the two methods

mapf[Transpose[x, {2, 3, 1}]]; // AbsoluteTiming

Map[Sqrt[Det[# . Transpose[#]]] &, Transpose[x, {2, 3, 1}]]; // AbsoluteTiming

(* {0.22, Null} *)

(* {1.1, Null} *)

Using a GPU (NVIDIA RTX 3050 Laptop) as the target device:

(* {0.067, Null} *)

However, the error between the two methods increases with the number of components on which Map operates (tensor length/dimension):

Table[m = 2; n = 3;
SeedRandom[0];
x = RandomReal[{-1, 1}, {m, n, o}];

{o, Max@
Abs[Map[Sqrt[Det[# . Transpose[#]]] &, Transpose[x, {2, 3, 1}]]/
   mapf[Transpose[x, {2, 3, 1}]] - 1]}, {o, 
Floor[10^(Subdivide[1, 6, 20])]}] // 
ListLogLogPlot[#, PlotRange -> All] &

Maximum relative error (vertical axis) vs tensor dimension on which Map operates (horizontal axis)

error


Previous version

(still interesting for the syntax)

This is a natural task for ArrayReduce introduced in 2020 (version 12.1). However, it is experimental and significantly slower than just using Table (for example using Norm with the example provided by OP ArrayReduce was 16 times slower than Table).

Thus using ArrayReduce at this moment is more for the convenience of writing code quickly than having high performance.

Discussion on how to use ArrayReduce

(code section below)

In the original Table version in the question, applying a function h to the two first levels of an array Array[a, {2, 2, 3}] would be:

Table[h@Array[a, {2, 2, 3}][[All, All, s]], {s, 3}]

To visualize the output of that code one can include a composition with MatrixForm:

Note : This is only to visualize the output, using MatrixForm for purposes other than visualization, in particular within calculations, can lead to problems. One may search for such examples here on stack exchange.

Table[h@MatrixForm@Array[a, {2, 2, 3}][[All, All, s]], {s, 3}]

Out: $$\left\{h\left(\left( \begin{array}{cc} a(1,1,1) & a(1,2,1) \\ a(2,1,1) & a(2,2,1) \\ \end{array} \right)\right),h\left(\left( \begin{array}{cc} a(1,1,2) & a(1,2,2) \\ a(2,1,2) & a(2,2,2) \\ \end{array} \right)\right),h\left(\left( \begin{array}{cc} a(1,1,3) & a(1,2,3) \\ a(2,1,3) & a(2,2,3) \\ \end{array} \right)\right)\right\}$$

How can one do this, functional style ?

ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 3}], {{1}, {2}}]

ArrayReduce computed the function h on the levels 1 and 2 where the input of the function h is a matrix made from these levels.

Why all those {} ? Why not just {1,2} instead of {{1},{2}} ?

Well, let's check

ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 3}], {1, 2}]

$$\{h(\{a(1,1,1),a(1,2,1),a(2,1,1),a(2,2,1)\}),h(\{a(1,1,2),a(1,2,2),a(2,1,2),a(2,2,2)\}),h(\{a(1,1,3),a(1,2,3),a(2,1,3),a(2,2,3)\})\}$$

The output is the same as if we had kept the {{1},{2}} but we flattened the input matrix before feeding it to h.

Maybe it might be hard to figure out how this generalizes so let's add another dimension.

ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 2, 2}], {1, 2, 3}]

$$\{h(\{a(1,1,1,1),a(1,1,2,1),a(1,2,1,1),a(1,2,2,1),a(2,1,1,1),a(2,1,2,1),a(2,2,1,1),a(2,2,2,1)\}),h(\{a(1,1,1,2),a(1,1,2,2),a(1,2,1,2),a(1,2,2,2),a(2,1,1,2),a(2,1,2,2),a(2,2,1,2),a(2,2,2,2)\})\}$$

Same ol' same ol'

take levels 1,2,3 -> flatten -> feed it to h

Ok,let's add some brackets

ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 2, 2}], {{1}, {2}, {3}}]

$$\left\{h\left(\left( \begin{array}{cc} \{a(1,1,1,1),a(1,1,2,1)\} & \{a(1,2,1,1),a(1,2,2,1)\} \\ \{a(2,1,1,1),a(2,1,2,1)\} & \{a(2,2,1,1),a(2,2,2,1)\} \\ \end{array} \right)\right),h\left(\left( \begin{array}{cc} \{a(1,1,1,2),a(1,1,2,2)\} & \{a(1,2,1,2),a(1,2,2,2)\} \\ \{a(2,1,1,2),a(2,1,2,2)\} & \{a(2,2,1,2),a(2,2,2,2)\} \\ \end{array} \right)\right)\right\}$$

That might be hard to read. Without the MatrixForm

(* {h[{{{a[1, 1, 1, 1], a[1, 1, 2, 1]}, {a[1, 2, 1, 1], a[1, 2, 2, 1]}}, {{a[2, 1, 1, 1], a[2, 1, 2, 1]}, {a[2, 2, 1, 1], a[2, 2, 2, 1]}}}], h[{{{a[1, 1, 1, 2], a[1, 1, 2, 2]}, {a[1, 2, 1, 2], a[1, 2, 2, 2]}}, {{a[2, 1, 1, 2], a[2, 1, 2, 2]}, {a[2, 2, 1, 2], a[2, 2, 2, 2]}}}]} *)

Not much better if not worse but notice the h[{{{ which means that h is fed a 3-tensor since there are 3 brackets. So with {{1},{2},{3}} h is fed the 3-tensor corresponding to the first three levels. Ok, last thing before we apply this to the question asked.

Given the discussion above what to expect from ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 2, 2}], {{1, 2}, {3}}] ? From the first example we see that {{something},{something}} means h is fed a 2-tensor or a matrix and {{something},{something},{something}} means it is fed a 3-tensor. So whatever happens, we know that h will be fed a matrix. Next we said that using {1,2} instead of {{1},{2}} means that levels 1 and 2 are flattened. So the result we find is :

ArrayReduce[h@*MatrixForm, Array[a, {2, 2, 2, 2}], {{1, 2}, {3}}]

$$\left\{h\left(\left( \begin{array}{cc} a(1,1,1,1) & a(1,1,2,1) \\ a(1,2,1,1) & a(1,2,2,1) \\ a(2,1,1,1) & a(2,1,2,1) \\ a(2,2,1,1) & a(2,2,2,1) \\ \end{array} \right)\right),h\left(\left( \begin{array}{cc} a(1,1,1,2) & a(1,1,2,2) \\ a(1,2,1,2) & a(1,2,2,2) \\ a(2,1,1,2) & a(2,1,2,2) \\ a(2,2,1,2) & a(2,2,2,2) \\ \end{array} \right)\right)\right\}$$

Ok, let's apply this to the original question.

Applying ArrayReduce to the question.

First, we might feel that all the brackets are tedious and we are not planning on flattening matrices anytime soon or at least not with ArrayReduce. So we can make a custom array reduce:

arrayreduce = ArrayReduce[#1, #2, List /@ #3] &

Now we can obtain the result by OP:

arrayreduce[Sqrt@Det[# . Transpose@#] &, x, {1, 2}] // AbsoluteTiming

Out: (* {0.000292, {1.04663, 0.437045, 0.479911, 0.260814, 0.205563, 0.171896, 1.20112, 1.00502, 0.855893, 0.125758}} *)

Verify that it worked:

Table[Sqrt[Det[x[[;; , ;; , i]] . Transpose[x[[;; , ;; , i]]]]],
{i, 1, o}] == arrayreduce[Sqrt@Det[# . Transpose@#] &,
 x,
{1, 2}]

Compare the timing:

ArrayReduce: $2.36\times 10^{-4} s$

Table: $7.4\times 10^{-5} s$

$\endgroup$
5
  • $\begingroup$ For a speed comparison one should use many samples. For 1 million cases ArrayReduce has the same performance as Map, i.e. about 3-4 times faster than Table. $\endgroup$ Oct 21, 2022 at 13:57
  • $\begingroup$ I just checked and for m = 10^2; n = 10^2; o = 10^3;SeedRandom[0]; x = RandomReal[{-1, 1}, {m, n, o}]; comparing Map, ArrayReduce and Table with RepeatedTiming, I found that ArrayReduce was a little bit slower than Table and that Map was 1.4 times faster than Table. $\endgroup$ Oct 22, 2022 at 0:01
  • $\begingroup$ I mean you should increase only o to 1 million and not change m and n. $\endgroup$ Oct 22, 2022 at 0:27
  • $\begingroup$ @granularbastard maybe Python would be faster not sure. If your computer has Python and Mathematica can find it, you can call numpy, which is a package of Python, from Mathematica. That said it might not be worth the trouble. I am not sure. $\endgroup$ Oct 22, 2022 at 0:28
  • $\begingroup$ Depending on the version of Mathematica you have, you can write code from a different language by typing > at the beginning of a cell and choosing the language. This might require setting up the language before. $\endgroup$ Oct 22, 2022 at 0:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.