2
$\begingroup$

I have just created a code snippet in Octave/Matlab that aims to create a plot which shows the accuracy of an initial probability vector $\vec{\pi}$ derived from the transition probability matrix $\textbf{P}$ of a Markov chain using its powers ($\textbf{P}^n$).

n = 1:10000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^10000,2);

for i=1:n(end)
  y(i) = (tProbs^n(i)*[1;1]-iProbs)(1);
end

Is there a way, using vectorisation on the line mentioned, to remove the main loop and hopefully improve running times? I have not been able to raise a matrix to the power of a vector $\vec{n}$ .

P.D
This does not solve my problem. I am looking for a vectorisation technique. Just in case there is any confusion.

I also hope that this was the right Stack Exchange in which to post this question. If it is completely out of place, I would appreciate it if a Mod. could migrate this question somewhere more appropriate.

$\endgroup$
4
$\begingroup$

Before even discussing any vectorization, there are gross inefficiencies in the code you provided. Just by employing basic techniques, I got 21-fold speedup:

n = 1000000;
y = zeros(1,n);

tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
iProbsSqr = iProbs.^2;
V0 = [1;1];

for i=1:n
  V0 = tProbs*V0;
  y(i) = sqrt(dot(V0.^2,iProbsSqr));
end

If the previous version takes 35 seconds on your machine, this should take about 1.5. Further improvement is not possible without increasing the memory usage significantly. Also, it would look very ugly and cryptic (see How to Write Unmaintainable Code by Roedy Green - here is a link: https://www.se.rit.edu/~tabeec/RIT_441/Resources_files/How%20To%20Write%20Unmaintainable%20Code.pdf ). If you still want to do vectorization, let me know; I think using repmat one could precompute all the matrix powers in a vectorized fashion, and then do a tensor reduction. This may work, potentially.

Edit: There is a minor modification which improves the efficiency in Octave (does not affect MATLAB). Octave has huge overhead with function calls and just calling dot rather than using '* costs a lot. sum and sqrt seem to be fine.

n = 1000000;
y = zeros(1,n);

tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);
iProbsSqr = iProbs.^2;
V0 = [1;1];

for i=1:n
  V0 = tProbs*V0;
  y(i) = sqrt((V0.^2)'*iProbsSqr);
end

Vectorization of this code is very difficult. It requires computing all the matrix powers up to n and that is not possible to do efficiently with built-in functions (well, the whole idea is just redundant computation but seeing past through that). If you can precompute all the necessary matrix powers, this becomes an exercise in matrix-vector and vector-vector multiplications.

You can try to precompile the for loop using oct-files (see https://octave.org/doc/v6.2.0/Getting-Started-with-Oct_002dFiles.html )

Edit 2: Unmaintainable piece of code delivered.

n = 1:1000000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);

[V,D] = eig(tProbs);
d = diag(D);
dp = d.^n;

aux1 = inv(V)*[1;1];
aux1 = dp.*aux1;
aux1 = V*aux1;
aux1 = aux1.^2;
aux1 = aux1'*(iProbs.^2);
aux1 = sqrt(aux1);

I would feel bad if I didn't explain why this code is unmaintainable. So here I go: we have to compute all the matrix powers simultaneously, but that is stupid, and it is a lot of redundant calculation. However, we are not using C here; redundant calculations are good tools to gain speedup.

But how do we compute all the matrix powers at once, you ask? We assume that the matrix is diagonalizable (phew, fortunately this time it is. Hope the next guy is smart enough to check that), because if a matrix $A$ is diagonalizable, i.e. if there is a diagonal matrix $D$ and an invertible matrix $V$, then $A^n = VD^nV^{-1}$.

But isn't this just pushing the problem around? Not really! We take the diagonal part out of the matrix $D$, put it into a vector, and we can take an elementwise power of that vector to the vector n to obtain a 2-by-howevermuch matrix dp. (Yay! The code is becoming harder to follow!). Now what?

Just compute the result using a string of operations: multiply $V^{-1}$ and $[1,1]^T$, elementwise multiply the result with the vector dp, multiply the result with $V$,... Hope you followed.

The resulting code runs in 0.1 seconds in Octave (close to what I get with a simpler and easier to maintain code using MATLAB) and the same code runs in 0.01 seconds in MATLAB. My preference would be to use the simpler version on MATLAB, but MATLAB is proprietary so feel free to use the complicated/unmaintainable version :) Don't forget to comment it. At least if you have to go back and read it yourself, you can understand it.

Edit 3: Modifications to accommodate the edit to the question. Here is the code:

n = 1:1000000;
tProbs = [0.1 0.5; 0.9 0.5];
iProbs = sum(tProbs^100000,2);

[V,D] = eig(tProbs);
d = diag(D);
dp = d.^n;

aux1 = inv(V)*[1;1];
aux1 = dp.*aux1;
aux1 = V*aux1;
aux1 = aux1-iProbs;
aux1 = aux1(1,:);

This is another reason why this code is unmaintainable. A simple operation change requires large amount of code restructuring.

$\endgroup$
8
  • $\begingroup$ Yeah the repmat thing did occur to me. Oh well I did make some pretty big mistakes at the beginning but I was kind of hurrying through that part in all fairness. I'll change it. Yes I'm interested in the vectorisation part of it's OK with you. $\endgroup$ Jun 11 at 8:11
  • $\begingroup$ BTW I just tried both versions and your code runs 40% faster (50s ->30s) using a different computer $\endgroup$ Jun 11 at 11:09
  • $\begingroup$ Are you using Octave or MATLAB? Using MATLAB 2020b, the original version takes 5 seconds to run on my 4-thread i5 and the "optimized" version runs in 0.20 seconds. I am guessing Octave because the original code takes about 15seconds on a single cpu high thread (64) count machine and mine takes 9 seconds on the same machine. $\endgroup$ Jun 11 at 14:17
  • $\begingroup$ I will test a few things for vectorization and update my answer if they work. If I don't update my answer today, probably my ideas did not really work. $\endgroup$ Jun 11 at 14:19
  • $\begingroup$ Yep I'm using Octave $\endgroup$ Jun 11 at 17:30

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.