Calculating probability to get from s1:s400 to sn with large Markov transition matrix in R takes forever

Question

I successfully implemented the following task, but I only got it to work with a small test dataset. For my real dataset, it just keeps calculating forever.

I have a 400x400 transition probability matrix in R. A user hits "Conversion" if she converts on the Markov walk. The absorbing state is "Null" for all users. "Start" is my beginning state.

Two things I need to calculate:

Hit state s_j on a random walk beginning at "Start"
Hit "Conversion" on a random walk beginning at each of the 397 other states

First one is easy in R:

v <- numeric(length = ncol(transitionMatrix1))
v[1] <- 1
i <- 2
R0 <- v%*%(transitionMatrix1 %^% 1)
R <- R0
repeat {
  R1 <- v%*%(transitionMatrix1 %^% i)
  R <- rbind(R, R1)
  if (rowSums( R[nrow(R)-1,] - R1 ) == 0) {
    #if (rowSums( R[nrow(R)-1,] - R1 ) < epsilon) {
    break
  }
  else {
    i <- i+1
  }
}
visit1 <- colSums(R)

I successfully implemented 2., but I only got it to work with a small matrix. It takes forever with a big one:

w <- i
C1 <- matrix( nrow = w, ncol = ncol(transitionMatrix1))
for (i in 1:ncol(transitionMatrix1)) {
  x <- numeric(length = ncol(transitionMatrix1))
  x[i] <- 1
  for (j in 1:w) {
    C1[j,i] <- x%*%(transitionMatrix1 %^% j)[,ncol(transitionMatrix1)-1]
  }
}
convert1 <- colSums(C1)

I should probably not use loops. Unfortunately, I did not succeed in vectorizing said operations.

What helps me was to keep TWO matrix in memory at the same time - original one and transposed one, and use one or another accordingly — Severin Pappadeux, Apr 25 '16 at 16:16
While I don't completely understand what you are trying to calculate, I think you can get the expressions you want without repeatedly multiplying matrices. For a quick overview see https://en.wikipedia.org/wiki/Absorbing_Markov_chain and for a complete analysis see http://www.math.pku.edu.cn/teachers/yaoy/Fall2011/Kemeny-Snell1976.pdf . You need to calculate the fundamental matrix and then you can obtain all sorts of performance measures. — Forzaa, Apr 26 '16 at 11:25

Sebastian Lechner · Accepted Answer · 2016-05-18T11:09:00.227

Thanks @Forzaa, with the help of your links, I was able to find a solution in MatLab (https://de.mathworks.com/matlabcentral/newsreader/view_thread/48984), which I converted to R:

 f.absorb <- function(P) {
  # input: P an absorbing Markov transition matrix in 'from rows to column' form. Rows are
  # probability vectors and must sum to 1.0. At least one main diagonal
  # element=1 output: 
  # tau: Time before absorption from a transient state. 
  # B: P(ending up in absorbing state from a transient state) 
  # tau2: variance matrix for tau; 
  # N: fundamental matrix 
  # N2: covariance matrix for N; 
  # CP: transition matrix in canonical form. 
  # Q: Submatrix of transient states, 
  # R: from transient to absorbing state submatrix, 
  # H: hij the probability process will ever go to transient state j starting 
  # in transient state i (not counting the initial state (K & S, p. 61)) 
  # MR: expected number of times that a Markov process remains in a transient 
  # state once entered (including entering step) 
  # VR: Variance for MR (K & S, Theorem 3.5.6, p.
  # 61) written by E. Gallagher, Environmental Sciences Program UMASS/Boston
  # Internet: Gallagher@umbsky.cc.umb.edu written 9/26/93, revised: 10/26/94
  # refs: Kemeny, J. G. and J. L. Snell. 1976. Finite Markov chains.
  # Springer-Verlag, New York, New York, U.S.A. Roberts, F. 1976. Discrete
  # Mathematical Models with applications to social, biological, and
  # environmental problems. Prentice-Hall
  # R adaption from MatLab Code found at 
  # https://de.mathworks.com/matlabcentral/newsreader/view_thread/48984
  pdim = nrow(P)
  if (rowSums ((colSums(t(P))-matrix(1,1,pdim)) * (colSums(t(P))-matrix(1,1,pdim))) > 0.1) {
    stop('Rows of P must sum to 1.0')
  }
  dP=diag(P);
  ri=which(dP==1.0);
  if (is.null(ri)) {
    stop('No absorbing states (reenter P or use ergodic.m)')
  }
  rdim=length(ri);
  qi=which(dP!=1.0);
  qdim=length(qi);
  I=diag(qdim)
  Q=P[qi,qi];
  N=solve(I-Q); # the fundamental matrix
  tau=data.matrix(colSums(t(N)));
  CP=P[c(t(ri),t(qi)),c(t(ri),t(qi))]; # the canonical form of P
  R=CP[(rdim+1):pdim,1:rdim];
  B=N%*%R;
  Ndg=diag(diag(N));
  N2=N%*%(2*Ndg-I)-N*N
  tau2=(2*N-I)%*%tau-tau*tau;
  H=solve(Ndg,N-I);
  dQ=diag(Q);
  oneq=matrix(1,qdim,1);
  MR=oneq/(oneq-dQ);
  VR=(dQ/(oneq-dQ)) * (dQ/(oneq-dQ));
  list(tau=tau, N=N, B=B, CP=CP, Q=Q, R=R, N2=N2, H=H, MR=MR, VR=VR)
}

The values I was searching for can be found in the matrix H

Calculating probability to get from s1:s400 to sn with large Markov transition matrix in R takes forever

1 Answers1