Using mtimesx-function multiple times without a loop-function

Edit: each slice of matrix P is of tridiagonal form.

Hi,

I need to compute the m-th power of every single 2D-layer of a 3D matrix P with size(P)=[50,50,1000] in an efficient/fast way.

I tried using mtimesx, like

 if true
   P_s=P;
   for i=1:m
      P_s=mtimesx(P_s,P_s);
   end
 end

However it takes too much time, since m is in the order of 1e5. Is there are way to implement m>2 without using a loop-function?

I'm glad for any advice, thanks

markus

[Matt Fig]

What is mtimesx? If you ask a question that refers to a non-standard function, it might help to clarify what you are talking about...

[Markus]

You're right, sorry. mtimesx is available at FEX (written by James Tursa) and basically a powerful tool for matrix multiplications in N dimensions. In my case I tried to use it, since it allows to multiply 3D-matrices slice by slice without calling each slice separately in a loop.

Answer 1

[Matt J]

Don't know what kind of speed you're looking for but this doesn't seem too bad. It works by converting P to a block diagonal sparse matrix whose blocks are the slices of the original P,

    P=rand(50,50,1000); m=1e5;

    Pcell=num2cell(P,[1,2]);
    Pcell=cellfun(@(c) sparse(c), Pcell,'uni',0);

    P_s=blkdiag(Pcell{:});

    tic;
     P_s=P_s^m;
    toc;
    %Elapsed time is 6.253806 seconds.

I leave it to you to convert P_s back to a 3D stack if that's what you need.

[James Tursa]

I don't know what OP's data looks like, but my guess is that using the mpower function directly is going to be more accurate, especially for large m, vice doing repeated matrix multiplications. Thus I would probably opt for doing the calculations directly in a loop using mpower, or doing the sparse approach ala Matt J.

[Jan]

@James: Would your FEX: mpower2 be helpful?

[James Tursa]

@Jan: Probably not. For reasonably low powers it can outperform the built-in mpower function for speed. So if you are doing a low power calculation many times (e.g. in a loop) it can be useful. But for higher powers it will eventually suffer from accuracy loss just like any other repeated-matrix-multiply method. So I would still opt for calling mpower in a loop or use the sparse approach ala Matt J for higher powers.

[Markus]

Thanks for the suggestions and helpful discussion so far. And sorry for late reply; couldn't access my code last weekend.

Comparing speed of slice-by-slice mpower-loop to sparse-matrix method by Matt J the loop seems to run faster on my system.

However, I just noticed that each slice of P is of tridiagonal form. Does this give more opportunities in speeding up the code?

Thanks a lot, markus