## 3 Computational Results
### Computational Cost Optimizations
\begin{array}{l}M_{t+1}=P^{\prime}\,\times\,T^{t}=M_{t}\,\times\,T\\ Q_{t+1}=M_{t}\,\times\,T\,\times R\\ Q_{total,t+1}=Q_{total,t}+Q_{t+1}\end{array}  \tag{14}\]  
Additionally, the iteration process is stopped when the number of photons transmitted with a scattering order of \(t\) is sufficiently small compared to the total number of photons transmitted. In this study, the ratio of \(Q_{t}\) and \(Q_{total}\) was obtained during each iteration, and the program would stop if the maximum value of the ratio vector was smaller than \(10^{-4}\). Figure 11 shows the computational cost breakdown for solving for the angular distribution, with and without the iterative approach. Simulations with various OD were performed and the Phase Function 2 was selected. As can be seen in Figure 11, when utilizing the iteration process, the computational cost for the inversion process is reduced from 147 s to 4 s (OD = 2) or 7 s (OD = 5). By incorporating this approach with the pre-calculated phase function \(P(\theta_{i},\theta_{j},n)\), the total computational cost for executing Markov chain approximations can be reduced from 463 s to 15 s (for OD = 2) and 18 s (for OD = 5). Figure 11: Histogram of computational costs for the original and improved Markov chain approximations. Figure 12 examines the stopping scattering order (SO) and Markov Chan computational costs as compared to the Monte Carlo simulations. The same CPU was utilized for Markov chain approximation and Monte Carlo simulations and 2.5 billion photons were sent in each case. Since MATLAB matrix manipulations involve parallel computing, the computational cost of the Monte Carlo simulation was scaled accordingly. From Figure 11, it is evident that with a higher OD, the computational cost of the modified Markov chain method was also higher because the stopping scattering order \(t\) was also higher. The stopping scattering order was roughly equal to the OD multiplied by 12, representing a quasi-linear relationship. The computational time for the modified Markov inversion process ranged from 4 s for OD = 2 cases to 41 s for OD = 30 cases, which was fast enough even when global optimization algorithms were utilized. In comparison, the Monte Carlo simulation cost 3000 s for OD = 2 cases and 61,000 s for OD = 30 cases, and such computational costs would inevitably discourage any meaningful inversion process. It is worth noting that varying the discretization could slightly change the stopping scattering order, but the relative difference was within 10% as the authors examined. Therefore, it was established that the computational cost estimation in Figure 12 is valid for various discretization schemes. Figure 13 depicts the absolute relative error obtained with the modified Markov inversion process. Uniform zenith angle discretization was used and the OD discretization scheme 1 was chosen. As such, the conditions in Figure 13 are the same as Figure 5, except for where Markov inversion process was chosen. As can be seen in all six cases, from Panel (a) to Panel (f), the absolute relative error with the modified Markov inversion process was at the same level as the original Markov inversion with a higher computational cost, therefore verifying the application of the modified Markov inversion process, as well as the recommended value of \(10^{-4}\) as the stopping criterion. Note that this recommended stopping criterion is subject to change depending on the desired accuracy, and the program can be stopped even earlier if high accuracy is not required. Furthermore, with the modified method, transmitted angular distribution at each scattering order can also be derived without additional computational costs, thus this method is more efficient than directly performing matrix inversions in solving multiple scattering problems. Finally, it is worth mentioning that the modified Markov inversion process can also lead to the use of a greater transition matrix, thus enabling finer resolution for Markov chain process.