## 3 Computational Results
### Computational Cost Optimizations  
One feature of the Markov chain approximation is that the computational cost remains relatively constant with different discretization schemes, phase functions, or different ODs; given that the size of transition matrix is fixed (which is determined by the number of meshes in the system). In the previous works, the computational cost for anisotropic multiple scattering was reported as \(\sim\)8 min. Such a computational cost is obviously superior compared to many Monte Carlo simulations, especially with larger OD values, because with more scattering events, the computational cost for a single photon increases considerably. Despite that, it is still desirable to further decrease the computational cost of the original Markov chain program, so that more computationally efficient reconstruction algorithms can be made possible. To analyze the computational costs, the approximation process was broken down into three parts. The three parts were: (1) Solving Equations (3) and (4), which builds the connection between the propagation angle before scattering (\(\theta_{1}\)) and after scattering (\(\theta_{2}\)); (2) Solving Equations (1) and (2), which incorporates the OD and part 1 to form the transition matrix; and (3) Perform a matrix inversion and obtain the results. For Markov chain methods, part (1) and part (3) consume most of the computational cost in the simulation, with each process consuming a computational time of more than 100 s.  
To reduce the computational cost, for part (1), it is worth noting that \(P\left(\theta_{i},\theta_{j},n\right)\) is only dependent on the phase function in this specific layer where the scattering happens, which means \(P\left(\theta_{i},\theta_{j},n\right)\) can be calculated in advance and a 180 by 180 matrix can be stored to represent each phase function. When the transition matrix is needed during the calculation, the simulation can load the matrices directly with negligible computational costs. This is also an added benefit of the Markov chain approximations when compared to the Monte Carlo method, because the trajectory of each photon must be calculated in every evaluation in Monte Carlo simulations. Figure 10: Absolute relative error of Markov chain approximations compared to Monte Carlo simulations with zenith angle discretization scheme 3. Phase function 3 was selected. Panel (**a**) OD = 2, Panel (**b**) OD = 5. The angular resolution was 0.25\({}^{\circ}\) for \(\theta\) < 15\({}^{\circ}\) and \(\theta\) > 165\({}^{\circ}\); 2.5\({}^{\circ}\) for \(\theta\) > 15\({}^{\circ}\) and \(\theta\) < 165\({}^{\circ}\). In part (3), the computational cost for matrix inversion is also significant, mostly due to the computational complexity of inversing a significantly large transition matrix (18,000 by 18,000). From the definition of the fundamental matrix:  
\[N=\left(I-T\right)^{-1}=I+T+T^{2}+T^{3}+\ldots+T^{n},n\,\sim\infty \tag{11}\]  
where \(N\) is the fundamental matrix, and \(T\) is the transition matrix obtaining the transition probabilities. Note that the angular distribution can be expressed as (based on Markov chain Theory):  
\[ \begin{array}{l}Q_{total}=P^{\prime}\,\times N\times R=P^{\prime}\,\times \left(I-T\right)^{-1}\times R\\ =P^{\prime}\,\times\left(I+T+T^{2}+T^{3}+\ldots+T^{n}\right)\times R\end{array}  \tag{12}\]  
where \(Q_{total}\) is the sought angular distribution, \(R\) is the absorption matrix, and \(P^{\prime}\) is the initial scattering distribution. One way to reduce the computational cost of the Markov chain inversion is to use the expanded form of \(N\) to solve for \(N\), instead of obtaining the inverse of \(\left(I-T\right)\) directly. If it is defined that:  
\[M_{t}=P^{\prime}\,\times\,T^{t-1},Q_{t}=M_{t}\times R \tag{13}\]  
\(Q_{t}\) stands for the angular distribution of photons that experienced exactly \(t\) times of scattering events (which is defined as the scattering order of this transmitted photon), before exiting the turbid slab. Then an iterative relationship can be easily derived:  
\[