Transformations of Markov Random Processes and Modifications of the Monte Carlo Method for Solving the Radiative Transfer Equation

G. P. Tarasov

(Presented by the Department of Engineering and Computational Mathematics)

1. Introduction

1.1. If the transfer of radiation is described by specifying the collision density \(\Psi(x)\), satisfying the integral equation

\[ \Psi(x)=K\Psi(x)+\Psi^1(x), \tag{1.1} \]

then for a linear functional \(I_\xi=(\Psi,\varphi)\) of the solution of (1.1) the following assertion is well known (countably additive functions \(\Psi(dx)\) and their densities \(\Psi(x)\) are not distinguished).

Proposition. If \(\Psi^1\in M^+(X,A)\), \(K\in N^+(M)\), \(\|K\|<1\), then

\[ I_\xi=(G\Psi^1,\varphi), \tag{1.2} \]

where \(M^+(X,A)\) is the nonnegative cone in the \(B\)-space \(M(X,A)\) of countably additive functions on \(X\); \(N^+(M)\) is the nonnegative cone in the space of endomorphisms of \(M\); \(G\) is the Green operator of equation (1.1), belongs to the weak closure \(G(N)\) of the Neumann algebra in \(N\) generated by the operator \(K\), i.e.,

\[ G=\lim_{m\to\infty}\sum_{r=0}^{m}K^{(r)}. \tag{1.3} \]

1.2. An extensive literature \([2, 3, 5, 7]\) is devoted to the computation of the functional (1.2) by the Monte Carlo method (MCM). The modifications of the MCM that have been developed, pursuing the goal of improving the efficiency of the calculations, are based on various conceptual premises \([2]\). In the present note it is proposed to consider these modifications of the MCM from the standpoint of transformations over the Markov random process (MRP) of walks. Being statistical in its nature, the MCM has three standard aspects: 1) organization of walks (organization of observations), 2) construction of estimating functionals from the trajectories of walks (choice of observed quantities), 3) construction of an estimate \(\hat I_\xi\) of the desired functional (processing of the results of observations). The last point is uninteresting, since as \(\hat I_\xi\) one almost always takes the sample mean (and sample-

root variance); the first two points are of fundamental importance. First we discuss the construction of the MRP; then modifications are considered. The proofs are brief or omitted altogether. They may be found in the cited literature.

2. Construction of the MRP and general propositions of the MMC

2.1. Consider the construction of an MRP [1] for a jump-type random-walk process in transfer [4]. Let \((U,\sigma,P)\) be the basic probability space, \(X=X^* \cup O,\ X^*=V\times\Omega\times E,\ O=O_1\cup O_2\), (\(X^*\) is the set of existence, \(O_1\) is the natural absorbing set, and \(O_2\) is the artificial absorbing set), \((X,A)\) the phase space, \(Z=0,\infty\) (\(m\in Z\) is the collision order), and \(x[k,u]\) the trajectory of the process (for each \(u\in U\) it takes values in \((X,A)\) for \(k\in Z\)). The function \(n:U\to Z\), defined by the relation

\[ n(u)=\sup(k:x[k,u]\in X^*), \tag{2.1} \]

determines the termination of the trajectory (which for \(k>n\) is completed by setting \(x[k,u]\in O\)). Let
\[ \tau_m=\{x_m^{-1}\Gamma\subset U:\Gamma\in A,\ m\le n\},\qquad \tau=\{x_k^{-1}\Gamma\subset U:\Gamma\in A\} \]
be the algebras of events observed over \(m\) collisions, and let \(\sigma\) be the algebra generated by the process. The relation

\[ p(x_{k-1}\to\Gamma_k)=P_{x_{k-1}}\{x^{-1}\Gamma_k,\tau_{k-1}\} = P_{x_{k-1}}\{x_k\in\Gamma_k,k<n\} \tag{2.2} \]

defines the transition probability from the state \(x_{k-1}\) to the set \(\Gamma_k\in A\) in one step; it satisfies the recurrence

\[ p^{(m)}(x\to\Gamma)=\int_X p^{(m-1)}(x\to dy)p(y\to\Gamma), \tag{2.3} \]

where

\[ p(x_{k-1}\to\Gamma_k)= \begin{cases} \displaystyle \int_{\Gamma_k} c_k e^{-\tau_k}\,d\tau_k\,h(x_{k-1}'\to dx_k'), & \text{if } \begin{cases} x_{k-1}\notin O,\\ O\cap\Gamma_k=\varnothing, \end{cases}\\[1.2em] (1-c_k), & \text{if } \begin{cases} x_{k-1}\notin O,\\ \Gamma_k=O, \end{cases}\\[0.6em] 0, & \text{if } x_{k-1}\in O; \end{cases} \]

\[ h(x_{k-1}'\to dx_k')= \begin{cases} \delta(x_{k-1}'-x_k')\,dx_k', & \text{if } k=1,\\ \text{indicatrix}, & \text{if } k>1. \end{cases} \]

2.2. Let \((X,B)\) be the sample space of the process, \(D=\{D_m,\ m=1,\infty\}\) a partition of \(X\) into disjoint subsets (\(D_m\) contains only \(m\)-link trajectories), \(C=\{C_m,\ m=1,\infty\}\), \(C_m=x\setminus\bigcup_{k=1}^m D_m\), the primitive class of cylindrical sets (\(C_m\) contains trajectories having no fewer than \(m\) links), and let \(\{\Gamma_k\in A,\ k=\overline{1,m}\}\) be the base of the cylindrical set \(B_m\); then

\[ p^{[m]}(B_m)=P(B_m\cap C_m)=\int \pi(dx_0)\left[\prod_{k=1}^{m}p(x_{k-1}\to dx_k)\right] = X\times\prod_{r=1}^{m}\Gamma_r. \tag{2.4} \]

defines a probability distribution on cylindrical sets, and

\[ P^{(m)}(\widetilde{B}_m)=P(B_m\cap D_m)=\int \pi(dx_0)\left[\prod_{k=1}^{m}p(x_{k-1}\to dx_k)\right]g(x_m)=X\times \prod_{k=1}^{m}\Gamma_r . \tag{2.5} \]

defines probabilities on disjunctive partitions. Here \(\pi(dx)\) is the initial distribution \((\pi(X^*)=1)\), \(g(x_m)=p(x_m\to O)\) is the probability of termination at \(x_m\).

The distributions (2.4) and (2.5) generate one and the same probability measure \(P\) on \(B\), which makes it possible to consider the directly specified probability space \((X,B,P)\). We shall denote the constructed Markov stochastic process by \(m=(U,X,p,n)\). If \(p\) is naturally determined by the kernel of the operator \(K\) (1.1), then \(m\) is called consistent with the transfer equation (1.1). If \(\widetilde{m}=(\widetilde{U},\widetilde{X},\widetilde{p},\widetilde{n})\) is another Markov stochastic process, then the value of the Radon–Nikodym derivative is defined on the trajectory \(x[k]\in C_m\) by

\[ \frac{dP}{d\widetilde{P}}(x[k])= \frac{\pi(dx_0)}{\widetilde{\pi}(dx_0)} \left[ \prod_{k=1}^{m} \frac{p(x_{k-1}dx_k)}{\widetilde{p}(x_{k-1}\to dx_k)} \right]. \tag{2.6} \]

If \(x[k]\in D_m\), then

\[ \frac{dP}{d\widetilde{P}}(x[k])= \frac{\pi(dx_0)}{\widetilde{\pi}(dx_0)} \left[ \prod_{k=1}^{m} \frac{p(x_{k-1}\to dx_k)}{\widetilde{p}(x_{k-1}\to dx_k)} \right] \frac{g(x_m)}{\widetilde{g}(x_m)} . \tag{2.7} \]

2.3. As estimating functionals of trajectories of the Markov stochastic process, additive functionals [1] are used, i.e., mappings \(F_\varphi:X\to R_1\) possessing the property

\[ F_\varphi(x[k]^l_i)+F_\varphi(x[k]^j_{l+1})=F_\varphi(x[k]^j_i), \]

where \(x[k]^j_i\) is the segment of the trajectory from the \(i\)-th to the \(j\)-th collision (inclusive). A large class of additive functionals \(F_\varphi(x[k])\), sufficient for practical purposes, can be represented by the formula (for \(x[k]\in D_n\))

\[ \eta=F_\varphi(x[k,u])=\sum_{l=1}^{n(u)}\int_X T(x_l(u),dx)\varphi(x). \tag{2.8} \]

For example, \((\delta_{ij}\) is the Kronecker symbol).

1) The kernel

\[ T(x_i,dx)=\delta(x_i-x)\,dx \tag{2.9} \]

defines the simple estimate \(I_1\) [2],

\[ F_\varphi(x[k,u])=\sum_{i=1}^{n}\varphi(x_i). \]

2) The kernel

\[ T(x_i,dx)=\sum_{r=1}^{k} a_{ir}p^{(r)}(x_{i-1}\to dx), \tag{2.10} \]

where

\[ a_{ir}= \begin{cases} \delta_{li}, & \text{if } r<k-1,\\ 1, & \text{if } r=k, \end{cases} \]

defines the estimator \(F_\varphi(x[n,u])\), which performs an analytic calculation of the contributions up to \((k-1)\)-fold scattered radiation and a statistical estimation of the contributions of \((k\) and higher)-fold scattered radiation.

3) The kernel

\[ T(x_i,dx)=\hat{\delta}_{in}\left[\frac{\delta(x_i-x)}{g(x)}\right]dx \tag{2.11} \]

defines the simple estimator \(I_2\) [2], \(F_\varphi(x[k,u])=\dfrac{\varphi(x_n)}{g(x_n)}\). The following assertions hold.

2.4. Theorem. If \(\mathbf{m}\) is an MSP consistent with the transport equation (1.1), then the estimating functionals (2.9), (2.10), (2.11) provide an unbiased estimate \(I_\zeta\) (1.2).

Proof [2, 4]. In accordance with formula (2.5),

\[ M\eta=\int F_\zeta(x[k])P(dx[k])= \sum_{m=1}^{\infty}\int_{D_m} F_\zeta(x[k])P^{(m)}(dx[k])=I_\zeta . \tag{2.12} \]

2.5. Theorem. If the conditions

\[ \text{a) } g(x)>0,\quad x\in X^*, \]

\[ \text{b) } \left|\int_X T(x,dy)\,\zeta(y)\right|<\infty\quad(\bmod P), \]

are satisfied, then the estimate

\[ \hat I_\zeta=\frac{1}{N}\sum_{i=1}^{N}\eta_i \]

is consistent.

Proof. If a) holds, the chain of jumps \(\mathbf{m}\) \((\bmod P)\) is finite, and hence, taking b) into account, the estimating functional \(F_\zeta(x[k])\) is bounded \((\bmod P)\), so that

\[ D\eta=\int_X |F_\zeta(x[k])|^2P(dx[k])=I_\zeta^2<\infty . \tag{2.13} \]

(2.12) and (2.13) satisfy the conditions of the law of large numbers, which proves the theorem.

2.6. Theorem. Let \(\mathbf{m}\) and \(\hat{\mathbf{m}}\) be two MSPs, with \(\mathbf{m}\) consistent with (1.1). Then \(\mathbf{m}\) can be used for an unbiased estimation of \(I_\zeta\), if the following conditions are satisfied:

a) there exists a measurable mapping

\[ S:(X,B,P)\longrightarrow (X,\bar B,\bar P), \]

such that the measure \(PS^{-1}(\cdot)\) induced by it is absolutely continuous with respect to the measure \(\bar P\), i.e., there exists

\[ \frac{dPS^{-1}}{d\bar P}(x[k]). \]

b) there exists a measurable function \(F_\zeta:X\to R_1\) such that the original estimating functional (2.8) has the representation

\[ F_\zeta(x[k])=F_\zeta S(x[k]). \]

Proof. According to the “change of variables rule” [6], we have

\[ \int_x F_\varphi S(x[k])P(dx[k])= \int_x \left\{F_\varphi(\tilde{x}[k])\frac{dPS^{-1}}{d\tilde{P}}(x[k])\right\}\tilde{P}(dx[k]). \]

From the condition of unbiasedness of the estimator it is necessary, for all \(\varphi\in(M+)^*\),

\[ I_\varphi=\int_x F_\varphi(x[k])P(dx[k]) =\int_x F_\varphi S(x[k])P(dx[k]). \]

This proves the necessity of a), b).

2.7. Corollary. Let \(\mathbf{m}\) be an MSP consistent with (1.1), and let \(\tilde{\mathbf{m}}\) be an MSP obtained by transforming \(\mathbf{m}\) without changing the “trajectory map” \((x=\tilde{x})\). Then \(\tilde{\mathbf{m}}\) can be used for unbiased estimation of \(I_\varphi\), if the measure \(P\) is absolutely continuous with respect to \(\tilde{P}\), i.e., the derivative (2.6) exists.

Proof.

\[ I_\varphi=\int_x F_\varphi(x[k])P(dx[k]) =\int_x\left\{F_\varphi(x[k])\frac{dP}{d\tilde{P}}(x[k])\right\}\tilde{P}(dx[k]). \]

3. Relation of MMC modifications to transformations of MSP

Transformations over \(\mathbf{m}=(U,X;p,n)\) that do not take it out of the class of MSP are connected with transformations over the objects \(U, X, p, n\) [1]. To these transformations of the random-walk process \(\mathbf{m}\), consistent with the transport equation (1.1), one can reduce the known modifications of the MMC for computing \(I_\varphi\). We dwell on several of them.

3.1. The method of analytical averaging (of absorption, termination) reduces to transforming \(\mathbf{m}\) by cleaning the space of elementary events [1] \(U\) down to the set \(U\), corresponding to the part of the process \(\mathbf{m}\) on the marked set \(\tilde{X}\subset X\). The new MSP, \(\tilde{\mathbf{m}}=(U,\tilde{X},p,n)\), has transition function

\[ \tilde{p}(x\to\Gamma)=\frac{p(x\to\Gamma)}{p(x\to\tilde{X})},\quad x\in\tilde{X},\ \Gamma\in A\cap\tilde{X}. \tag{3.1} \]

For example, if \(\tilde{X}=\tilde{V}\times\Omega\times E\), where \(\tilde{V}\) is a convex set in \(R_3\), then analytical averaging of absorption and leakage gives

\[ p(x_k\to\tilde{X})= c_k\left(1-\exp\left[-\int_0^{L_k}\mu(\vec r_k+t\vec\Omega_k,E_k)\,dt\right]\right), \]

where

\(L_k\) is the distance from \(\vec r_k\) to \(d\tilde{V}\) in the direction \(\vec\Omega_k\).

If \(x[k]\in C_m\), then

\[ \frac{dPT^{-1}}{dP}=\prod_{i=1}^{m}p(x_i\to\tilde{X}). \tag{3.2} \]

3.2. Modification of the MMC with the aid of branching reduces to transforming \(\mathbf{m}\) by splitting elementary events \(u\in U\) [1] and randomly replacing time.

If \(\gamma^{-1}(u)=\{\tilde u:\gamma(\tilde u)=u\}\), then \(P_x(\gamma^{-1}A)=P_x(A)\) and, correspondingly, \((x=\tilde x)\); for \(x[k]\subset C_m\),

\[ \frac{dP}{d\tilde P}(\tilde x[k])=\prod_{i=1}^{m}\frac{1}{r(x_i)}, \tag{3.3} \]

where

\[ r(x_i):X\longrightarrow Z \]

indicates into how many parts the trajectory-branch splits at the point \(x_i\).

3.3. Modification of the MMC by means of “roulette” reduces to the transformation of \(m\) by constructing an \(\alpha\)-subprocess [1]. This is connected with shortening the lifetime \(n(u)\) to \(r(u)\). Here

\[ \tilde p(x_{k-1}\to\Gamma_k)=P_{x_{k-1}}(x_k\in\Gamma_k,\ k<\tau). \tag{3.4} \]

and, correspondingly, on \(\tilde x[k]\subset C_m\),

\[ \frac{dP}{d\tilde P}(\tilde x[k])=\prod_{i=1}^{m}\frac{1}{q(x_i)}, \tag{3.5} \]

where

\(q(x_i)\) is the probability (planned) of terminating the trajectory at collision \(x_i\).

3.4. Modification of the MMC by means of sampling according to “importance” reduces to the transformation of \(m\) by means of a \(\xi\)-transformation of the transition probabilities [1]. Here

\[ \tilde P_{x_{k-1}}(x_k\in\Gamma_k) = \int_{x^{-1}\Gamma_k}\xi(u)P_{x_{k-1}}(du). \tag{3.6} \]

There exists a \(T\)-measurable function \(\xi(u)\) such that the corresponding measure \(P\) on the sampling space \((\tilde X,\tilde B)\) yields zero variance for the estimate \(\eta\).

In this case, on \(x[k]\in C_m\),

\[ d\tilde P(d\tilde x[k]) = \frac{F_\varphi(x[k])}{I_\varphi}P(dx[k]). \tag{3.7} \]

3.5. Remark. The function \(\xi(u)\) is not measurable with respect to any of the \(\sigma\)-subalgebras \(T_k\), which greatly complicates the practical realization of such an algorithm. Various variants of this modification are discussed in [2, 5, 7].

3.6. Modification of the MMC by means of “dependent trials” with unchanged \(X\) reduces to the introduction of a family \(\{\tilde m\}\) with equivalent measures \(\{\tilde P\}\). One MSP \(\tilde m\in\{\tilde m\}\) is taken as the leading one. In this case, on \(x[k]\subset C_m\),

\[ \frac{dP}{d\tilde P}(\tilde x[k]) = \frac{\pi(dx_0)}{\tilde\pi(dx_0)} \prod_{i=1}^{m} \frac{p(x_{i-1}\to dx_i)}{\tilde p(x_{i-1}\to dx_i)}. \tag{3.8} \]

3.7. The modification of the MMC by means of a “correlated sample” [5] reduces to transforming \(\mathbf{m}\) by forming parts of the process [1] on sets \(X_i \subset X\).

Another class of modifications is associated with the use of various estimating functionals of trajectories of \(\mathbf{m}\), consistent with equation (1.1).

3.8. The modification of the MMC by means of the “collision method” [5] reduces to applying to \(\mathbf{m}\) the estimating functional \(F_{\zeta}(x[k])\) in the form (2.13) for \(k=1\).

3.9. Local and double-local estimates of the MMC reduce to applying to \(\mathbf{m}\) the estimating functional \(F_{\zeta}(x[k])\) in the form (2.13), respectively for \(k=2\) and \(k=3\).

In conclusion, let us note that the theory of transformations of MRP can be used to justify and construct modifications of the MMC with the aim of increasing the efficiency of computing functionals of the solution of the transport equation (1.1).

REFERENCES

1. E. B. Dynkin. Markov Processes. Moscow, Fizmatgiz, 1963.

2. S. M. Ermakov. The Monte Carlo Method and Related Questions. Moscow, “Nauka,” 1971.

3. The Monte Carlo Method in the Problem of Radiation Transport. Edited by G. I. Marchuk. Moscow, Atomizdat, 1967.

4. A. I. Khisamutdinov. Importance sampling in the theory of radiation transport. Journal of Computational Mathematics and Mathematical Physics, 1970, 10, No. 2, pp. 999–1005.

5. U. Fano, L. Spencer, M. Berger. Gamma-Ray Transport. Moscow, Atomizdat, 1963.

6. P. Halmos. Measure Theory. Moscow, IL, 1953.