7
$\begingroup$

I have a list of numerical complex 2x2 matrices

list = {M[4,3],T[3],M[3,2],T[2],M[2,1],T[1]}

(in reality this list contains over 1000 matrices and I have to do many runs, that's why fast would be good. I hope the structure is clear)

I would like to create a list of 2 dimensional vectors structured as follows:

vectorlist = {{1,B0}, M[2,1].T[1].{1,B0}, M[3,2].T[2].M[2,1].T[1].{1,B0}, M[4,3].T[3].M[3,2].T[2].M[2,1].T[1].{1,B0}}

Note that I do not need this structure, just the result of this multiplications. And since the j vector is M[j,j-1].T[j-1] dot the j-1 vector, it should be possible to exploit this to save time rather than calculating each vector on its own. But I am not sure how to do that.

Any additional ideas how to do this faster are very appreciated.

$\endgroup$
9
$\begingroup$

First a CompiledFunction; it expects the T-matrices and the M matrices is separate lists.

cf = Compile[{{T, _Complex, 3}, {M, _Complex, 3}, {B0, _Complex}},
   Block[{u, n},
    n = Min[Length[T], Length[M]];
    u = {1. + 0. I, B0};
    Join[
     {u},
     Table[
      u = Compile`GetElement[M, i].(Compile`GetElement[T, i].u),
      {i, 1, n}]
     ]
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Now you can do quite many such operations per second:

n = 1000000;
T = 0.5 RandomComplex[{-1 - I, 1 + I}, {n, 2, 2}];
M = 0.5 RandomComplex[{-1 - I, 1 + I}, {n, 2, 2}];
cf[T, M, 2.]; // AbsoluteTiming // First

0.236832

Actually, I made the experience that Dot in compiled code for small matrices/vectors tends to be a bit slow. Here is a handwritten expansion of Dot, that crucially depends on the undocumented Compile`GetElement.

cg = Compile[{{T, _Complex, 3}, {M, _Complex, 3}, {B0, _Complex}},
   Block[{a, b, n, A, an, bn},
    n = Min[Length[T], Length[M]];
    a = 1. + 0. I;
    b = B0;
    A = Table[0. + 0. I, {n + 1}, {2}];
    A[[1, 1]] = a;
    A[[1, 2]] = b;
    Do[
     an = Compile`GetElement[T, i, 1, 1] a + Compile`GetElement[T, i, 1, 2] b;
     bn = Compile`GetElement[T, i, 2, 1] a + Compile`GetElement[T, i, 2, 2] b;
     a = Compile`GetElement[M, i, 1, 1] an + Compile`GetElement[M, i, 1, 2] bn;
     b = Compile`GetElement[M, i, 2, 1] an + Compile`GetElement[M, i, 2, 2] bn;
     A[[i + 1, 1]] = a;
     A[[i + 1, 2]] = b;
     , {i, 1, n}];
    A
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

This is almost five times as fast on my machine:

A1 = cf[T, M, 2.]; // RepeatedTiming // First
A2 = cg[T, M, 2.]; // RepeatedTiming // First
A1 == A2

0.24

0.050

True

Beware that calling cg with T and M being lists of matrices of smaller size than $2 \times 2$ will cause Mathematica's kernel to crash without warning (because Compile`GetElement[M, i, 2, 2] is a forbidden memeory access which the OS will prevent by killing the Mathematica kernel). At least that is what seems to happen on macos.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much! Unfortunately, I get two error messages: "CreateLibrary::nocomp: A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers." and "Compile::nogen: A library could not be generated from the compiled function." From the Documentation I learn that I need to install a C Compiler on my computer, right? Do you have a recommendation for a free C Compiler? $\endgroup$ – Luke Mar 6 at 17:55
  • 1
    $\begingroup$ On Linux: gcc. On macos: Just install XCode. On Windows: I have heard that Visual Studio were the compiler of choice. But I am not sure. If you do not care that much for speed, just change CompilationTarget -> "C" to CompilationTarget -> "WVM" (Wolfram Virtual Machine). That helps already a lot, but is not as fast as "C". $\endgroup$ – Henrik Schumacher Mar 6 at 18:15
  • $\begingroup$ thanks a lot! Just a side note: In your code the order of the T and M matrices are reversed with respect to the order of my initial list, right? That's not a problem at all, just to avoid confusion :) $\endgroup$ – Luke Mar 6 at 18:57
  • $\begingroup$ Hm. I did not really care about it. But the order of the list T should be the same as yours: The first T is multiplied first. Same with M. It is in the nature of successive multiplication from the left, that the last matrix appears as first in the product when seen from the left... $\endgroup$ – Henrik Schumacher Mar 6 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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