Convergence issues in derivatives of Monte Carlo null-collision integral formulations: a solution

TL;DR

This paper investigates convergence issues in Monte Carlo derivatives of null-collision integral formulations, especially in radiative transfer, and proposes a solution to improve the accuracy of sensitivity evaluations.

Contribution

It identifies convergence problems in Monte Carlo sensitivity analysis for null-collision algorithms and offers a theoretical analysis and an alternative method to address these issues.

Findings

Null-collision algorithms cause convergence difficulties in sensitivity evaluations.

Theoretical analysis reveals why simultaneous Monte Carlo estimators may diverge.

An alternative solution improves the accuracy of derivative estimations.

Abstract

When a Monte Carlo algorithm is used to evaluate a physical observable A, it is possible to slightly modify the algorithm so that it evaluates simultaneously A and the derivatives $\partial$ $\varsigma$ A of A with respect to each problem-parameter $\varsigma$. The principle is the following: Monte Carlo considers A as the expectation of a random variable, this expectation is an integral, this integral can be derivated as function of the problem-parameter to give a new integral, and this new integral can in turn be evaluated using Monte Carlo. The two Monte Carlo computations (of A and $\partial$ $\varsigma$ A) are simultaneous when they make use of the same random samples, i.e. when the two integrals have the exact same structure. It was proven theoretically that this was always possible, but nothing insures that the two estimators have the same convergence properties: even when a large enough sample-size is used so that A is evaluated very accurately, the evaluation of $\partial$ $\varsigma$ A using the same sample can remain inaccurate. We discuss here such a pathological example: null-collision algorithms are very successful when dealing with radiative transfer in heterogeneous media, but they are sources of convergence difficulties as soon as sensitivity-evaluations are considered. We analyse theoretically these convergence difficulties and propose an alternative solution.

Figures

Equations

$$\left\{\begin{array}[]{rl}&\partial_{t}f+c\bm{\vec{\omega}}.\bm{\vec{\nabla}}f=-(k_{a}+k_{s})cf+k_{a}cf^{eq}+\displaystyle\int_{4\pi}k_{s}cf^{\prime}p_{\mathcal{S}}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime},~{}~{}~{}~{}\forall\bm{\vec{x}}\in\Omega,\forall\bm{\vec{\omega}}\in\mathbb{S}^{2}\\ &f(\bm{\vec{y}},\bm{\vec{\omega}}_{+})=f_{\partial\Omega}(\bm{\vec{y}},\bm{\vec{\omega}}_{+}),~{}~{}~{}~{}\forall\bm{\vec{y}}\in\partial\Omega,~{}\forall\bm{\vec{\omega}}_{+}\in\mathbb{S}_{+}^{2}\\ &f(\bm{\vec{x}},\bm{\vec{\omega}},0,\varsigma)=f_{0}(\bm{\vec{x}},\bm{\vec{\omega}}),~{}~{}~{}~{}\forall\bm{\vec{x}}\in\Omega,\forall\bm{\vec{\omega}}\in\mathbb{S}^{2}\end{array}\right.$$

$$\left\{\begin{array}[]{rl}&\partial_{t}f+c\bm{\vec{\omega}}.\bm{\vec{\nabla}}f=-\hat{k}cf+k_{a}cf^{eq}+\displaystyle\int_{4\pi}k_{s}cf^{\prime}p_{\mathcal{S}}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}+\displaystyle\int_{4\pi}k_{n}cf^{\prime}\delta(\bm{\vec{\omega}}-\bm{\vec{\omega}^{\prime}})\,\mathrm{d}\omega^{\prime},~{}~{}~{}~{}\forall\bm{\vec{x}}\in\Omega,\forall\bm{\vec{\omega}}\in\mathbb{S}^{2}\\ &f(\bm{\vec{y}},\bm{\vec{\omega}}_{+})=f_{\partial\Omega}(\bm{\vec{y}},\bm{\vec{\omega}}_{+}),~{}~{}~{}~{}\forall\bm{\vec{y}}\in\partial\Omega,\forall\bm{\vec{\omega}}_{+}\in\mathbb{S}_{+}^{2}\\ &f(\bm{\vec{x}},\bm{\vec{\omega}},0,\varsigma)=f_{0}(\bm{\vec{x}},\bm{\vec{\omega}}),~{}~{}~{}~{}\forall\bm{\vec{x}}\in\Omega,\forall\bm{\vec{\omega}}\in\mathbb{S}^{2}\end{array}\right.$$

$$\begin{array}[]{rl}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=&\exp\left(-\displaystyle\int_{0}^{\lambda_{\partial\Omega}}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)f_{\partial\Omega}(\bm{\vec{y}},\bm{\vec{\omega}})\\ +&\displaystyle\int_{0}^{\lambda_{\partial\Omega}}\exp\left(-\displaystyle\int_{0}^{\lambda}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)\left(\begin{array}[]{rl}&k_{a}(\bm{\vec{x}^{\prime}},\varsigma)f^{eq}(\bm{\vec{x}^{\prime}})\\ +&k_{s}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&k_{n}(\bm{\vec{x}^{\prime}},\varsigma)f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\mathrm{d}\lambda\end{array}$$

$$\exp \left(-{\int }_0^{{\lambda }_{\partial \Omega }}\hat{k}\left(\overrightarrow{\stackrel{~}{\mathbf{x}}}\right)\mathrm{d}\tilde{\lambda }\right)={\int }_{{\lambda }_{\partial \Omega }}^{+\infty }\hat{k}({\vec{\mathbf{x}}}^{\prime })\exp \left(-{\int }_0^{\lambda }\hat{k}\left(\overrightarrow{\stackrel{~}{\mathbf{x}}}\right)\mathrm{d}\tilde{\lambda }\right)\mathrm{d}\lambda$$

$$\begin{array}[]{rl}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=&\displaystyle\int_{\lambda_{\partial\Omega}}^{+\infty}\hat{k}(\bm{\vec{x}^{\prime}})\exp\left(-\displaystyle\int_{0}^{\lambda}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)f_{\partial\Omega}(\bm{\vec{y}},\bm{\vec{\omega}})\mathrm{d}\lambda\\ +&\displaystyle\int_{0}^{\lambda_{\partial\Omega}}\hat{k}(\bm{\vec{x}^{\prime}})\exp\left(-\displaystyle\int_{0}^{\lambda}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)\left(\begin{array}[]{rl}&\frac{k_{a}(\bm{\vec{x}^{\prime}},\varsigma)}{\hat{k}(\bm{\vec{x}^{\prime}})}f^{eq}(\bm{\vec{x}^{\prime}})\\ +&\frac{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{\hat{k}(\bm{\vec{x}^{\prime}})}\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\frac{k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}{\hat{k}(\bm{\vec{x}^{\prime}})}\displaystyle f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\mathrm{d}\lambda\end{array}$$

$$\begin{array}[]{rl}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=&\\ \displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda&\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})w_{\partial\Omega}\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)w_{A}\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=\displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})w^{\varsigma}_{\partial\Omega}\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)w^{\varsigma}_{A}\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\left(\begin{array}[]{rl}&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)\left(\begin{array}[]{rl}&\frac{\partial_{\varsigma}k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\left\{f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma);\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)\right\}=\displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\left\{w_{\partial\Omega};w^{\varsigma}_{\partial\Omega}\right\}\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left\{w_{A};w^{\varsigma}_{A}\right\}\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\left\{f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma);\left(\begin{array}[]{rl}&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\right\}\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)\left\{f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma);\left(\begin{array}[]{rl}&\frac{\partial_{\varsigma}k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\right\}\end{array}\right)\end{array}\right)\end{array}$$

$$\sigma_{\partial_{\varsigma}T}=\frac{\sqrt{L^{2}\mathrm{e}^{-k_{a}L}\Big{(}\frac{k_{n}+1/L}{k_{n}}\Big{)}-L^{2}\mathrm{e}^{-2k_{a}L}}}{\sqrt{N}}$$

$$\begin{array}[]{rl}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)&=\displaystyle\int_{0}^{+\infty}p_{\Lambda}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})w_{\partial\Omega}\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)w_{A}\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=\displaystyle\int_{0}^{+\infty}p_{\Lambda}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\left(-w_{\partial\Omega}\displaystyle\int_{0}^{\lambda_{\partial\Omega}}\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)\mathrm{d}l+w^{\varsigma}_{\partial\Omega}\right)\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left(-w_{A}\displaystyle\int_{0}^{\lambda}\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)\mathrm{d}l+w^{\varsigma}_{A}\right)\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\\ &\hskip 56.9055pt\left(\begin{array}[]{rl}-&\displaystyle\int_{0}^{\lambda}\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)\mathrm{d}lf(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=\displaystyle\int_{0}^{+\infty}p_{\Lambda}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\displaystyle\int_{0}^{\lambda_{\partial\Omega}}p_{L}(l|\lambda_{\partial\Omega})\mathrm{d}l\left(-w_{\partial\Omega}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda_{\partial\Omega})}+w^{\varsigma}_{\partial\Omega}\right)\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\displaystyle\int_{0}^{\lambda}p_{L}(l|\lambda)\mathrm{d}l\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left(-w_{A}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}+w^{\varsigma}_{A}\right)\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\\ &\hskip 56.9055pt\begin{array}[]{rl}&\left(\begin{array}[]{rl}-&\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\end{array}\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=\displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\displaystyle\int_{0}^{\lambda_{\partial\Omega}}p_{L}(l|\lambda_{\partial\Omega})\mathrm{d}l\left(-w_{\partial\Omega}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda_{\partial\Omega})}+w^{\varsigma}_{\partial\Omega}\right)\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\displaystyle\int_{0}^{\lambda}p_{L}(l|\lambda)\mathrm{d}l\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left(-w_{A}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}+w^{\varsigma}_{A}\right)\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\\ &\hskip 56.9055pt\begin{array}[]{rl}&\left(\begin{array}[]{rl}-&\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\end{array}\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)\left(\begin{array}[]{rl}-&\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\end{array}\right)\end{array}\right)\end{array}$$

$$\begin{array}[]{rl}&\left\{f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)~{};~{}\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)\right\}=\displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda\\ &\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\displaystyle\int_{0}^{\lambda_{\partial\Omega}}p_{L}(l|\lambda_{\partial\Omega})\mathrm{d}l\left\{w_{\partial\Omega}~{};~{}-w_{\partial\Omega}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda_{\partial\Omega})}+w^{\varsigma}_{\partial\Omega}\right\}\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\displaystyle\int_{0}^{\lambda}p_{L}(l|\lambda)\mathrm{d}l\left(\begin{array}[]{rl}&P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left\{w_{A}~{};~{}-w_{A}\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}+w^{\varsigma}_{A}\right\}\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}\\ &\hskip 56.9055pt\begin{array}[]{rl}&\left\{f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)~{};~{}\left(\begin{array}[]{rl}-&\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\frac{\partial_{\varsigma}k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{s}(\bm{\vec{x}^{\prime}},\varsigma)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right)\right\}\end{array}\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)\left\{f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)~{};~{}\left(\begin{array}[]{rl}-&\frac{\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)}{p_{L}(l|\lambda)}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right)\right\}\end{array}\right)\end{array}\right)\end{array}$$

$$\frac{{\sigma }_{{\partial }_{\varsigma }T}}{{\partial }_{\varsigma }T}=\frac{{\sigma }_T}{T}=\frac{\sqrt{1-e^{-k_{a}L}}}{\sqrt{N}}$$

$$\begin{split}\partial_{\varsigma}T(L,\varsigma)=&\displaystyle\int_{0}^{\infty}\hat{k}e^{-\hat{k}\lambda}\,\mathrm{d}\lambda\left\{\begin{array}[]{rl}&\mathcal{H}(\lambda-L)~{}0\\ +&\mathcal{H}(L-\lambda)\left\{\begin{array}[]{rl}&P_{A}(\varsigma)~{}0\\ +&P_{N}(\varsigma)\Big{(}\frac{-1}{k_{n}(\varsigma)}T(L-\lambda,\varsigma)+\partial_{\varsigma}T(L-\lambda,\varsigma)\Big{)}\Big{)}\end{array}\right\}\end{array}\right\}\end{split}$$

$$\begin{array}[]{rl}Var_{Alg.\ref{alg:AlgoDf1}}(\partial_{\varsigma}T(L,\varsigma))&=\sum_{n=1}^{\infty}p_{n}w_{n}^{2}-\Big{(}\sum_{n=1}^{\infty}p_{n}w_{n}\Big{)}^{2}\\ &=L^{2}\mathrm{e}^{-k_{a}L}\Big{(}\frac{k_{n}+1/L}{k_{n}}\Big{)}-L^{2}\mathrm{e}^{-2k_{a}L}\end{array}$$

$$\begin{split}\partial_{\varsigma}T(L,\varsigma)=&\displaystyle\int_{0}^{\infty}\hat{k}e^{-\hat{k}\lambda}\,\mathrm{d}\lambda\left\{\begin{array}[]{rl}&\mathcal{H}(\lambda-L)~{}(-L)\\ +&\mathcal{H}(L-\lambda)\left\{\begin{array}[]{rl}&P_{A}(\varsigma)~{}0\\ +&P_{N}(\varsigma)~{}\partial_{\varsigma}T(L-\lambda,\varsigma)\end{array}\right\}\end{array}\right\}\end{split}$$

$$Va{r}_{Alg.\text{ref}alg:AlgoDf2}({\partial }_{\varsigma }T(L,\varsigma ))=L^2{\mathrm{e}}^{-k_{a}L}\left(1-{\mathrm{e}}^{-k_{a}L}\right)$$

$$Va{r}_{Alg.\text{ref}alg:Algof}({\partial }_{\varsigma }T(L,\varsigma ))={\mathrm{e}}^{-k_{a}L}\left(1-{\mathrm{e}}^{-k_{a}L}\right)$$

$${\partial }_{\varsigma }f(\overrightarrow{{\mathbf{x}}_0},\overrightarrow{{\mathbf{\omega }}_0},\varsigma )={\partial }_{\varsigma }f(\overrightarrow{{\mathbf{x}}_1},\overrightarrow{{\mathbf{\omega }}_0},\varsigma )+\frac{{\partial }_{\varsigma }{k}_n(\overrightarrow{{\mathbf{x}}_1},\varsigma )}{k_n(\overrightarrow{{\mathbf{x}}_1},\varsigma )}f(\overrightarrow{{\mathbf{x}}_1},\overrightarrow{{\mathbf{\omega }}_0},\varsigma )$$

$$\begin{array}[]{rl}\partial_{\varsigma}f(\vec{\bm{x_{0}}},\vec{\bm{\omega_{0}}},\varsigma)=&\partial_{\varsigma}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)+\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)+\frac{\partial_{\varsigma}k_{n}(\vec{\bm{x_{1}}},\varsigma)}{k_{n}(\vec{\bm{x_{1}}},\varsigma)}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)\\ =&\partial_{\varsigma}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)+\left(\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}+\frac{\partial_{\varsigma}k_{n}(\vec{\bm{x_{1}}},\varsigma)}{k_{n}(\vec{\bm{x_{1}}},\varsigma)}\right)f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)\end{array}$$

$$\begin{array}[]{rl}\partial_{\varsigma}f(\vec{\bm{x_{0}}},\vec{\bm{\omega_{0}}},\varsigma)=&\partial_{\varsigma}f(\vec{\bm{x_{3}}},\vec{\bm{\omega_{3}}},\varsigma)+\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{3}}},\varsigma)}{k_{s}(\vec{\bm{x_{3}}},\varsigma)}f(\vec{\bm{x_{3}}},\vec{\bm{\omega_{3}}},\varsigma)+\left(\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}+\frac{\partial_{\varsigma}k_{n}(\vec{\bm{x_{1}}},\varsigma)}{k_{n}(\vec{\bm{x_{1}}},\varsigma)}\right)f(\vec{\bm{x_{3}}},\vec{\bm{\omega_{3}}},\varsigma)\\ =&\partial_{\varsigma}f(\vec{\bm{x_{3}}},\vec{\bm{\omega_{3}}},\varsigma)+\left(\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{3}}},\varsigma)}{k_{s}(\vec{\bm{x_{3}}},\varsigma)}+\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}+\frac{\partial_{\varsigma}k_{n}(\vec{\bm{x_{1}}},\varsigma)}{k_{n}(\vec{\bm{x_{1}}},\varsigma)}\right)f(\vec{\bm{x_{3}}},\vec{\bm{\omega_{3}}},\varsigma)\end{array}$$

$$\begin{array}[]{rl}\partial_{\varsigma}f(\vec{\bm{x_{0}}},\vec{\bm{\omega_{0}}},\varsigma)=&w^{\varsigma}_{\partial\Omega}+\left(\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{3}}},\varsigma)}{k_{s}(\vec{\bm{x_{3}}},\varsigma)}+\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}+\frac{\partial_{\varsigma}k_{n}(\vec{\bm{x_{1}}},\varsigma)}{k_{n}(\vec{\bm{x_{1}}},\varsigma)}\right)w_{\partial\Omega}\end{array}$$

$${\int }_0^{\lambda }{\partial }_{\varsigma }{k}_e(\overrightarrow{{\mathbf{x}}_{\mathbf{l}}},\varsigma )\mathrm{d}l={\int }_0^{\lambda }{p}_L(l)\mathrm{d}l{\partial }_{\varsigma }{k}_e(\overrightarrow{{\mathbf{x}}_{\mathbf{l}}},\varsigma )/p_L(l)$$

$$\begin{array}[]{rl}\displaystyle\int_{0}^{\lambda}\partial_{\varsigma}k_{e}(\vec{\bm{x_{l}}},\varsigma)\mathrm{d}l=&\displaystyle\int_{0}^{\beta_{1}}p_{L}(l_{1})\mathrm{d}l_{1}~{}\partial_{\varsigma}k_{e}(\vec{\bm{\zeta_{1}}},\varsigma)/p_{L}(l_{1})\\ +&\displaystyle\int_{0}^{\beta_{2}}p_{L}(l_{2})\mathrm{d}l_{2}~{}\partial_{\varsigma}k_{e}(\vec{\bm{\zeta_{2}}},\varsigma)/p_{L}(l_{2})\\ +&...\\ +&\displaystyle\int_{0}^{\beta_{n}}p_{L}(l_{n})\mathrm{d}l_{n}~{}\partial_{\varsigma}k_{e}(\vec{\bm{\zeta_{n}}},\varsigma)/p_{L}(l_{n})\end{array}$$

$${\partial }_{\varsigma }\hat{k}=\frac{{\partial }_{\varsigma }{k}_e\hat{k}}{k_e}$$

$$\begin{array}[]{rl}\partial_{\varsigma}f(\bm{\vec{x}},\bm{\vec{\omega}},\varsigma)=&\\ \displaystyle\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda&\left(\begin{array}[]{rl}&\mathcal{H}(\lambda-\lambda_{\partial\Omega})\left[-w_{\partial\Omega}\displaystyle\int_{0}^{\lambda_{\partial\Omega}}\partial_{\varsigma}\hat{k}(\bm{\vec{x}_{l}},\varsigma)\mathrm{d}l+w_{\partial\Omega}^{\varsigma}\right]\\ +&\mathcal{H}(\lambda_{\partial\Omega}-\lambda)\left(\begin{aligned} &P_{A}(\bm{\vec{x}^{\prime}},\varsigma)\left[\begin{array}[]{rl}&\psi(\bm{\vec{x}},\lambda,\varsigma)w_{A}\\ +&\left(\frac{\partial_{\varsigma}k_{a}}{k_{a}}-\frac{\partial_{\varsigma}\hat{k}}{\hat{k}}\right)_{|\bm{\vec{x}^{\prime}},\varsigma}w_{A}\\ +&w_{A}^{\varsigma}\end{array}\right]\\ +&P_{S}(\bm{\vec{x}^{\prime}},\varsigma)\displaystyle\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}\left[\begin{array}[]{rl}&\psi(\bm{\vec{x}},\lambda,\varsigma)f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\left(\frac{\partial_{\varsigma}k_{s}}{k_{s}}-\frac{\partial_{\varsigma}\hat{k}}{\hat{k}}\right)_{|\bm{\vec{x}^{\prime}},\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}^{\prime}},\varsigma)\end{array}\right]\\ +&P_{N}(\bm{\vec{x}^{\prime}},\varsigma)\left[\begin{array}[]{rl}&\psi(\bm{\vec{x}},\lambda,\varsigma)f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\\ +&\partial_{\varsigma}f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)\end{array}\right]\end{aligned}\right)\end{array}\right)\end{array}$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

When numerically simulating linear-transport physics using Monte Carlo algorithms, one of the most recurrent difficulties is the handling of highly non-homogenous or fast-variating media. This difficulty was encountered since the beginning of neutron-transport and plasma-physics modelling. But a quite elegant trick was soon identified as a way to bypass this difficulty : virtual collisionners can be added where the true collisionners are scarce so that the total collisionner-density is homogeneous. Of course, in order to ensure that the physical problem is unchanged, when a particle interacts with a virtual collisionner, it simply continues its path as if no collision had occurred [1, 2, 3, 4]. This is the meaning of the denomination null-collision algorithm or fictitious-collision algorithm111Similar keywords are pseudo-collision, null-events, fictitious-events, null-collisions, Woodcock tracking and maximum cross-section. The first practical benefit is that the next collision event can be sampled as if the medium was homogenous. Then the choice is made to select a true-collision or a virtual-collision as function of their local respective-amounts and this is how the spatial information is recovered. But several other benefits were recently foreseen in [1] and practically tested in [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16], mainly for radiative-transfer applications. The main idea is that null-collision algorithms transform the non-linearity of Beer-extinction into a linear-recursive problem that Monte Carlo handles without approximation[14]. This was for instance used in [5] to deal with absorption-spectra of molecular gases combining very numerous transitions: the summation over all transitions could be treated by the Monte Carlo algorithm itself, which was previously assumed impossible because this summation was inside the exponential of Beer-extinction. Similarly, the vanishing of the exponential allowed the extension of implicit Monte Carlo algorithms for inversion of absorption and scattering coefficients from intensity measurements [6]. Outside radiative transfer, a very similar idea was used to solve Electromagnetic Maxwell equations for energy propagation in particle-ensembles of statistically-distributed shapes despite of the nonlinearity associated to the square of the electric field[13]. Again similar is the algorithm proposed in [14] solving Boltzmann equation for micro-fluidics applications despite of the nonlinearity of the collision operator.

Back to radiative-transfer applications, the ideas suggested in [1] have motivated significant developments in the computer-graphics community for the cinema industry. Here the benefit of using null-collisions is that it extends to participating media (aerosols or clouds) the orthogonality between data-description and data-treatment that was at the heart of the most recent use of Monte Carlo for rendering complex scenes [7, 9, 8, 10]. The algorithm is indeed processed without any knowledge of the exact spatial-information, and it is only when a collision occurs that access to the field is required: the interaction between the radiative-transfer algorithm and the field-data is strictly restricted to this very moment. This allows the implementation of numerous acceleration techniques with little changes by comparison with those developed for handling complex surfaces. One of these techniques consists in the setting of an acceleration grid, adjusting the amount of virtual collisionners so that the total collisionner-density is both homogeneous in part and close enough to the real density-field. This avoids the sampling of too many useless virtual-collisions. This is one of the starting points of the present paper: null-collision algorithms allow the use of any amount of virtual-collisionners but numerical efficiency justifies that one tries to reduce them to the minimum.

However, we show here that reducing the amount of virtual-collisionners to a minimum leads to convergence difficulties when evaluating sensitivities. Sensitivity evaluation is a very general feature of Monte Carlo techniques: when a Monte Carlo algorithm is used to evaluate a physical observable $A$, it is always possible to modify the algorithm in such a way that it evaluates both $A$ and the derivatives $\partial_{\varsigma}A$ of $A$ with respect to each problem-parameter $\varsigma$, and most commonly the corresponding implementation is quite straightforward [17, 18, 19, 20, 21, 22, 23]222This is not at all straightforward for domain-deformation sensitivities[19, 20], but we here stick to pure parametric sensitivities. But evaluating sensitivities using null-collision algorithms is pathological: the better we adjust the acceleration grid, the worse the statistical convergence rate. In Sec. 2 we will illustrate this pathological behaviour evaluating the transmissivity of a beam through a non-homogeneous column. Then we propose an alternative approach in Sec. 3 where the design of the sensitivity-evaluation algorithm starts from the standard integral solution of the Boltzmann equation, i.e. without virtual-collisionners. The resulting sampling requirements are then addressed with the null-collision approach viewed as a simple rejection-sampling approach. This introduces the cost of sampling an additional random variable, but at this cost the convergence difficulties vanish. We illustrate the numerical behaviour of this modified algorithm is Sec. 4 using a benchmark inspired of [1].

2 Convergence difficulties when evaluating sensitivities

In this section we design a Monte-Carlo algorithm using the standard null-collision approach for the evaluation of the distribution function $f$, i.e. the solution of Boltzmann equation, and apply the technique of [19, 17] for simultaneous evaluation of a sensitivity $\partial_{\varsigma}f$ with respect to $\varsigma$, where $\varsigma$ is a parameter appearing in the absorption and scattering coefficients. The Boltzmann equation that we use is linear with constant speed particles. It matches the monochromatic radiative-transfer equation exactly and all application examples will be restricted to radiative transfer. We choose to make use of the notation $f$ instead of the more radiative-transfer oriented notation $I=h\nu cf$ (the specific intensity) in order to simplify the access for readers of the plasma and neutronics communities. The monochromatic radiative-transfer equation becomes

[TABLE]

where $f\equiv f(\bm{\vec{x}},\bm{\vec{\omega}},t,\varsigma)$ with $\bm{\vec{x}}$ the location, $\bm{\vec{\omega}}$ the propagation direction and $t$ the time. For incoming scattering in any direction $\bm{\vec{\omega}^{\prime}}$ of the unit sphere $\mathbb{S}^{2}$, we write $f^{\prime}\equiv f(\bm{\vec{x}},\bm{\vec{\omega}^{\prime}},t,\varsigma)$ and $p_{\mathcal{S}}$ is the single scattering phase function, i.e. $p_{\mathcal{S}}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\mathrm{d}\omega^{\prime}$ is the probability density that the scattering direction is $\bm{\vec{\omega}}$ for this incoming direction $\bm{\vec{\omega}^{\prime}}$. The constant particle-speed is $c$ and the coefficients $k_{a}\equiv k_{a}(\bm{\vec{x}},t,\varsigma)$, $k_{s}\equiv k_{s}(\bm{\vec{x}},t,\varsigma)$ and $k_{e}=k_{a}+k_{s}$ are the absorption coefficient, the scattering coefficient and the extinction coefficient respectively. $f^{eq}\equiv f^{eq}(\bm{\vec{x}},t)$ is the equilibrium distribution (following the Planck function). $\Omega$ is the geometrical domain and $\partial\Omega$ its boundary at which the distribution function $f_{\partial\Omega}$ is known for all locations $\bm{\vec{y}}$ and all directions $\bm{\vec{\omega}}_{+}$ of the incoming hemisphere $\mathbb{S}_{+}^{2}$. $f_{0}$ is the initial condition.

Introducing null-collisions

In order to design a null collision algorithm (NCA)[1] we add a field of virtual collisionners such that the total extinction coefficient is practicable, in the sense that we can sample the corresponding beer extinction:

[TABLE]

where $k_{n}\equiv k_{n}(\bm{\vec{x}},t,\varsigma)$ is the null-collision coefficient, $\hat{k}=k_{a}+k_{s}+k_{n}$ is the total extinction-coefficient and $\delta$ is the Dirac distribution. Equation(2) is strictly equivalent to Eq.(1) because of the Dirac distribution that insures $\int_{4\pi}k_{n}cf^{\prime}\delta(\bm{\vec{\omega}}-\bm{\vec{\omega}^{\prime}})\,\mathrm{d}\omega^{\prime}=k_{n}cf$.

When numerically adressing the solution $f(\bm{\vec{x}}_{0},\bm{\vec{\omega}}_{0})$ of this transport equation at $(\bm{\vec{x}}_{0},\bm{\vec{\omega}}_{0})$ (also solution of Eq.(1)) using the Monte Carlo method, one of the most standard approach consists in a simple statistical reading that allows to view $f(\bm{\vec{x}}_{0},\bm{\vec{\omega}}_{0})$ as an average over radiative paths that are tracked backward from the observation location $(\bm{\vec{x}}_{0},\bm{\vec{\omega}}_{0})$ to the sources[1]. In this reading, the pure transport term $\partial_{t}f+c\bm{\vec{\omega}}.\bm{\vec{\nabla}}f$ corresponds to the spatial and temporal propagation of $f$ in direction $\bm{\vec{\omega}}$ at constant-speed $c$. The collisional term $-\hat{k}cf$ corresponds to either an absorption or a scattering event (including the null-collision events that are forward scattering events). When combining it with the transport term this leads to collision locations that are distributed exponentially along the line of sight (Beer law). Tracking the path backward, this means that the preceeding collision at $\bm{\vec{x}}_{1}$ is at a distance $\lambda_{0}$ that is a realisation of a random variable $\Lambda_{0}$ of probability density $p_{\Lambda_{0}}(\lambda_{0})=\exp(-\hat{k}\lambda_{0})$ (see Fig. 1). Once $\bm{\vec{x}}_{1}$ is sampled, the collision type is sampled in turn to decide wether an absorption, a true scattering or a null collision occurs. In the backward tracking picture, this corresponds respectively to the three remaining terms

with $k_{a}cf^{eq}$ an absorption event is translated into thermal emission and the algorithm stops with the Monte Carlo weight $f^{eq}(\bm{\vec{x}}_{1})$ (the source at $\bm{\vec{x}}_{1}$),

• 2.

with $\int_{4\pi}k_{s}cf^{\prime}p_{\mathcal{S}}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}$ a scattering event is translated into the sampling of a “previous” direction $\bm{\vec{\omega}}_{1}$ and the algorithm continues recursively as if evaluating $f(\bm{\vec{x}}_{1},\bm{\vec{\omega}}_{1})$,

• 3.

with $\int_{4\pi}k_{n}cf^{\prime}\delta(\bm{\vec{\omega}}-\bm{\vec{\omega}^{\prime}})\,\mathrm{d}\omega^{\prime}$ and its Dirac function, a null collision event is translated into a pure forward scattering event, i.e. the “previous” direction $\bm{\vec{\omega}}_{1}$ is equal to $\bm{\vec{\omega}}_{0}$.

Of course the statistical translation includes the boundary conditions: when backward reaching the boundary at a location $\bm{\vec{x}}_{i}$ and direction $\bm{\vec{\omega}}_{i}$, the algorithm stops with the Monte Carlo weight $f(\bm{\vec{x}}_{i},\bm{\vec{\omega}}_{i})$ (the incoming source at the boundary). The corresponding Monte Carlo algorithm is detailed in Alg. 1 and illustrated in Fig. 1.

Integral formulation

This null-collision algorithm belongs to the family of analog Monte Carlo algorithms, i.e. algorithms that can be designed without any formal development because they only numerically-implement the well established statistical pictures of radiation physics. However, in the present context it is very much useful to also choose a viewpoint under which the same algorithm appears as a statistical estimate of the integral solution of Eq.(2). For sake of clarity we only write this integral solution at the stationary limit:

[TABLE]

where $\bm{\vec{\tilde{x}}}=\bm{\vec{x}}-\tilde{\lambda}\bm{\vec{\omega}}$, $\bm{\vec{x}^{\prime}}=\bm{\vec{x}}-\lambda\bm{\vec{\omega}}$, $\bm{\vec{y}}=\bm{\vec{x}}-\lambda_{\partial\Omega}\bm{\vec{\omega}}$, with $\lambda_{\partial\Omega}$ the distance to the first boundary-intersection starting at $\bm{\vec{x}}$ in the direction $-\bm{\vec{\omega}}$, i.e. $\lambda_{\partial\Omega}=\mathrm{min}\{\left\lVert{\bm{\vec{x}}-\bm{z}}\right\rVert$; $\bm{z}\in\mathrm{Vect}^{-}(\bm{\vec{x}},\bm{\vec{\omega}})\cap\partial\Omega\}$ where $\mathrm{Vect}^{-}(\bm{\vec{x}},\bm{\vec{\omega}})=\{\bm{\vec{x}}-\lambda^{\prime}\bm{\vec{\omega}}$;  $\lambda^{\prime}\in\mathbb{R}_{+}\}$. This standard Fredholm equation, typical of the formal solution of linear-transport physics, can be transformed using the following property

[TABLE]

to give

[TABLE]

Then

$p_{\hat{\Lambda}}(\lambda)=\hat{k}(\bm{\vec{x}^{\prime}})\exp\left(-\int_{0}^{\lambda}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)$ can be viewed as the probability density function of the free path $\hat{\Lambda}$ (the distance until next collision),

• 2.

$P_{A}=\frac{k_{a}}{\hat{k}}$, $P_{S}=\frac{k_{s}}{\hat{k}}$ and $P_{N}=\frac{k_{n}}{\hat{k}}$ can be viewed as the probabilities that the collision is an absorption, a scattering event or a null-collision respectively,

• 3.

and the two integrals over $[0,\lambda_{\partial\Omega}[$ and $[\lambda_{\partial\Omega},+\infty[$ can be gathered into a single integral over $[0,+\infty[$ using the Heaviside function $\mathcal{H}$

to give

[TABLE]

with $w_{\partial\Omega}=f_{\partial\Omega}(\bm{\vec{y}},\bm{\vec{\omega}})$ and $w_{A}=f^{eq}(\bm{\vec{x}^{\prime}})$. This last equation is the integral formulation that we needed in order to construct Alg. 1 at the stationary limit: Alg. 1 is indeed nothing more than the algorithmic-reading of Eq. 5 (and reciprocally Eq. 5 is nothing more than the integral translation of Alg. 1, [18]):

$\int_{0}^{+\infty}p_{\hat{\Lambda}}(\lambda)\mathrm{d}\lambda$ stands for the sampling of the distance of the collision (according to the $\hat{k}$-field),

• 2.

$\mathcal{H}(\lambda-\lambda_{\partial\Omega})$ stands for the case where the sampled collision is outside the boundary, then the algorithm stops at the boundary with the Monte Carlo weight $w_{\partial\Omega}$ (the value of $f$ corresponding to the incoming radiation),

• 3.

$\mathcal{H}(\lambda_{\partial\Omega}-\lambda)$ stands for the case where the sampled collision is at a location $\bm{\vec{x}^{\prime}}$ inside the volume, and then three collision types are possible:

• (a)

$P_{A}(\bm{\vec{x}^{\prime}},\varsigma)$ stands for the case where the collision is an absorption, then the algorithm stops at the boundary with the Monte Carlo weight $w_{A}$ (the value of $f^{eq}$ at the collision location),

• (b)

$P_{S}(\bm{\vec{x}^{\prime}},\varsigma)$ stands for the case where the collision is a scattering event, then $\int_{4\pi}p_{S}(-\bm{\vec{\omega}^{\prime}}|-\bm{\vec{\omega}})\,\mathrm{d}\omega^{\prime}$ stands for the sampling of a new direction $\bm{\vec{\omega}^{\prime}}$ according to the phase function and the algorithm continues recursively with the estimation of $f$ at $\bm{\vec{x}^{\prime}}$ in direction $\bm{\vec{\omega}^{\prime}}$,

• (c)

$P_{N}(\bm{\vec{x}^{\prime}},\varsigma)$ stands for the case where the collision is null, then the algorithm continues recursively with the estimation of $f$ at $\bm{\vec{x}^{\prime}}$ in the unchanged direction $\bm{\vec{\omega}}$

Straightforward application of sensitivity-evaluation techniques

Now that we have constructed the integral formulation of Alg. 1 we can apply the sensitivity-evaluation technique introduced in [17, 19, 20]. It consists in derivating Eq. 5 with respect to $\varsigma$ and multiplying and dividing by each of the probabilities and probability density functions that depend on $\varsigma$. This leads to an integral formulation of the sensitivity that has the very same structure as that of Eq. 5:

[TABLE]

with $w^{\varsigma}_{\partial\Omega}=0$ and $w^{\varsigma}_{A}=\frac{\partial_{\varsigma}k_{a}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{a}(\bm{\vec{x}^{\prime}},\varsigma)}f^{eq}(\bm{\vec{x}^{\prime}})$. Because of their identical structure, we can gather Eq. 5 and 6 into one using the vectorial notation $\left\{w;w^{\varsigma}\right\}$:

[TABLE]

The algorithmic-reading of 7 leads to Alg. 2 that evaluates simultaneously $f$ and $\partial_{\varsigma}f$. The recursive nature of this algorithm comes from the fact that the final brackets in the scattering and null-collision terms contain $f$ and $\partial_{\varsigma}f$ at the same location in the same direction. The fact that their sensitivity part includes a summation is translated into an algorithm incrementing the Monte Carlo weight as explained in Appendix B.

Simulation examples

At this stage, we designed a null-collision algorithm, constructed the corresponding integral formulation and applied the proposition of [17, 19, 20] in a straightforward manner so that the algorithm also evaluates sensitivities. We now test this simulation strategy by evaluating the transmissivity of a non-diffusive heterogeneous column and also evaluating the sensitivity of this transmissivity w.r.t. $\varsigma$, a parameter influencing the absorption coefficient. Hereafter this configuration is called heterogeneous-slab (see Fig. 2): $\Omega$ is a column of length $L$ with $\bm{\vec{e_{x}}}$ the normal incoming at location $y=0$. The equilibrium distribution is null (cold medium, $f^{eq}\equiv 0$). The boundary conditions are $f^{\partial\Omega}(0,\bm{\vec{e_{x}}})=0$ and $f_{\partial\Omega}(L,-\bm{\vec{e_{x}}})=f_{inc}$. The absorption and scattering coefficients are $k_{a}(x,\varsigma)=\left(\varsigma-\gamma\right)\frac{\mathrm{atan}(-\alpha\left(x-\beta\right)+\frac{\pi}{2})}{\pi/2}+\gamma$ and $k_{s}\equiv 0$. Alg 2 is used to evaluate both $f(0,\bm{\vec{e_{x}}},\varsigma)$ and $\partial_{\varsigma}f(0,\bm{\vec{e_{x}}},\varsigma)$ that correspond to the transmissivity $T$ and its derivative $\partial_{\varsigma}T$ respectively: $T=f(0,\bm{\vec{e_{x}}},\varsigma)/f_{inc}$ and $\partial_{\varsigma}T=\partial_{\varsigma}f(0,\bm{\vec{e_{x}}},\varsigma)/f_{inc}$. We chose this particular profile of $k_{a}$, because it is possible to calculate $T$ and $\partial_{\varsigma}T$ analytically (see the caption of Fig. 2). Example Monte Carlo results, using $N=10000$ samples, are compared to the analytical solution in Tables 5a and 5b. The statistical uncertainty is noted $\sigma$ (the standard deviation of the Monte Carlo estimator). In Table 5b we also provide the number of samples $N_{1\%}$ required to achieve a $1\%$ accuracy. The simulations were made using five different $\hat{k}$-profiles (each overestimating $k_{a}$ at all locations), with acceleration-grids, $\hat{k}$ being uniform within each mesh (see Fig. 2):

for $\hat{k}_{20\%}$ no grid is used: the profile of $\hat{k}$ is uniform, equal to $1.2$ time the maximum $k_{a}$-value.

• 2.

for $\hat{k}_{1}$ no grid is used: the profile of $\hat{k}$ is again uniform, exactly equal to the maximum $k_{a}$-value.

• 3.

for $\hat{k}_{10}$ the grid is constructed in such a way that across each mesh the variations of $k_{a}$ are $1/10$ of the maximum $k_{a}$-value, and the profile of $\hat{k}$ is uniform within each mesh, exactly equal to the maximum $k_{a}$-value inside the mesh.

• 4.

for $\hat{k}_{100}$ and $\hat{k}_{1000}$ the grid is constructed the same way with $1/100$ and $1/1000$ variation respectively.

The transmissivity results of Table 5a confirm that the estimation of $T$ is insensitive to the adjustment of the $\hat{k}$-field (only the computation time is affected). But the sensitivity results of Table 5b clearly indicate the opposite: the statistical convergence is worse when $\hat{k}$ is close to $k$ and the number of samples required to reach a given accuracy level can be risen up to infinity when matching $\hat{k}$ to $k$ exactly. This is the pathological behavior that we announced in introduction: sensitivities cannot be evaluated accurately when using acceleration grids reducing the number of virtual collisions.

The variance of the sensitivity estimate

For a better understanding of this behavior, we studied a homogeneous-slab for which the variance of the Monte Carlo estimate can be calculated analytically. This case is identical to the previous one (transmissivity of a purely absorbing column) but now $k=k_{a}$ is uniform: $k_{a}(\varsigma)\equiv\varsigma$, $T=exp(-\varsigma L)$ and $\partial_{\varsigma}T=-LT$. Of course, there is no need to make use of a null-collision algorithm as soon as $k$ is uniform. We only do it for theoretical reasons (with $\hat{k}>k$ uniform). This allows us to fully identify the reasons why the variance of the sensitivity estimate rises when reducing $k_{n}=\hat{k}-k$. This may sound trivial as soon as when encountering a null-collision event, the Monte Carlo weight of the sensitivity algorithm includes a factor $\frac{\partial_{\varsigma}k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}{k_{n}(\bm{\vec{x}^{\prime}},\varsigma)}=1/k_{n}$ (see Eq. 7), but reducing $k_{n}$ also reduces the number of such null-collision occurrences. This may lead to a compensation, maintaining the variance at a finite value. The developments of appendix A.1 indicate the opposite: the statistical uncertainty is indeed

[TABLE]

Figure 3 illustrates the meaning of this dependance of $\sigma_{\partial_{\varsigma}T}$ with the problem parameters. In this idealised case, looking at the behavior of such an algorithm applying sensitivity-evaluation techniques in a straightforward manner, the difficulty is well identified: when $\frac{k_{n}}{\hat{k}}$ approaches zero, the number of samples required for a $1\%$ accurate evaluation of the sensitivity tends to infinity (see Fig. 3c). This figure also displays the behavior of an algorithm implementing the very same sensitivity-evaluation technique, but without the use of null-collisions (which is possible here in this idealised uniform case). Without null-collisions, the relative value of the standard deviation of the sensitivity-estimate (Fig. 3b) is identical to that of the main quantity (the transmissivity-estimate, Fig. 3a). This is an ideal behaviour: the sensitivity is estimated with the same relative accuracy as that of the main quantity. Altogether in this simple example, we see that evaluating sensitivities can be perfectly costless before using null-collisions and may become pathological when null-collisions are introduced.

Note that in the general case, even without null-collisions, evaluating sensitivities can be truely difficult. Understanding the relative variance of sensitivity estimates and comparing them to the relative variance of the algorithm estimating the main quantity was indeed one of the main concerns of the initial work of De Lataillade[17]. Essentially, serious difficulties arise as soon as the scattering optical-thickness is high. The objective of the present paper is not at all to address this specific issue: at the end of the following section, when an alternative solution will be proposed for evaluating sensitivities in null collisions algorithms, the problems associated to highly scattering media will remain unsolved.

3 An alternative approach

The preceding section identifies convergence difficulties when evaluating sensitivities using null-collisions. Theses difficulties are not associated to the standard sensitivity-evaluation algorithm itself: considering slab transmission, we have seen that when we do not make use of null-collisions, the sensitivity-evaluation algorithm converges as well as the algorithm evaluating the main quantity. So the observed difficulties are only the consequences of introducing virtual-collisionners. At this stage, null-collision algorithms appear therefore as perfect tools for handling heterogeneous fields, but are incompatible with the simultaneous evaluation of sensitivities.

We have seen that this problem is related to the term $\frac{1}{k_{n}}$ appearing in the Monte Carlo weight of the sensitivity algorithm. At which stage did this term appear and can we bypass this step? Clearly, $\frac{1}{k_{n}}$ appeared when derivating with $\varsigma$ the null-collision probability $P_{N}(\varsigma)=1-P_{A}(\varsigma)-P_{S}(\varsigma)$, with $P_{A}(\varsigma)=k_{a}(\varsigma)/\hat{k}$ and $P_{S}=k_{s}(\varsigma)/\hat{k}$. A first way to suppress this $\frac{1}{k_{n}}$ term consists in making $\hat{k}$ dependent on $\varsigma$. This is always possible because $\hat{k}$ is a free parameter and we can therefore adjust it to the variation of $k_{a}(\varsigma)+k_{s}(\varsigma)$ so that $P_{N}$ does not depend on $\varsigma$ anymore. We first tested this solution and it proved itself already quite practical: the corresponding details are provided in Appendix D. But we finally retained another algorithm, starting from the integral solution of the original Boltzmann Eq. 1, i.e. prior to the introduction of virtual-collisionners. The idea consists in first designing an algorithm evaluating simultaneously $f$ and $\partial_{\varsigma}f$ as if the heterogeneity of the field could be handled without difficulty and only introduce null-collisions in a second stage. For this, we can simply rewrite Eq. 5 with $k_{n}=0$ (no virtual collisionners):

[TABLE]

The only differences with Eq. 5 are that

$P_{N}=0$,

• 2.

the random variable $\hat{\Lambda}$ of probability density $p_{\hat{\Lambda}}(\lambda)=\hat{k}(\bm{\vec{x}^{\prime}})\exp\left(-\int_{0}^{\lambda}\hat{k}\left(\bm{\vec{\tilde{x}}}\right)\mathrm{d}\tilde{\lambda}\right)$ (the free path in the $\hat{k}$-field) is replaced with the random variable $\Lambda$ of probability density $p_{\Lambda}(\lambda)=k_{e}(\bm{\vec{x}^{\prime}},\varsigma)\exp\left(-\int_{0}^{\lambda}k_{e}(\bm{\vec{\tilde{x}}},\varsigma)\mathrm{d}\tilde{\lambda}\right)$ (the free path in the original $k_{e}$-field).

This equation can then be derivated with respect to $\varsigma$ and multiplied/divided by each of the probabilities and probability density functions that depend on $\varsigma$ (exactly the same way Eq. 7 was constructed from Eq. 5) to give

[TABLE]

A main point of the present paper is that this integral equation, although it was derived the same way as Eq. 7, cannot be interpreted in algorithmic terms: the integral pattern $\int\partial_{\varsigma}k_{e}\mathrm{d}l$ is not yet transformed into a statistical expectation. An additional random generation will be required. At this stage let us introduce an arbitrary random variable $L$ of probability density function $p_{L}$ and write $\int\partial_{\varsigma}k_{e}\mathrm{d}l=\int p_{L}(l)\mathrm{d}l\frac{\partial_{\varsigma}k_{e}}{p_{L}(l)}$. Reporting this into Eq. 10 and using $\int p_{L}(l)\mathrm{d}l=1$ leads to

[TABLE]

At this stage, null-collisions have not been introduced. Therefore, the algorithmic reading of Eq. 11 would not be practical as soon as the $k_{e}$-field is heterogeneous: the difficulty would come from the sampling of $\Lambda$. The objective of introducing null collisions will therefore be to replace $\Lambda$ with another path-length $\hat{\Lambda}$, shorter in average but easy to sample, and compensate the too many collisions by the fact that some of them are null. However, this not as trivial as in the algorithm for the main quantity because of the new random variable $L$ that we needed to introduce when transforming Eq. 10 into Eq. 11 (transforming it into an expectation). Indeed $\partial_{\varsigma}k_{e}$ needs to be integrated along the whole path, now including null collisions. In statistical terms, this means that the recursivity of the path-sampling algorithm is only insured if the Monte Carlo weight associated to null collisions includes the term $-\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)/p_{L}(l|\lambda)f(\bm{\vec{x}^{\prime}},\bm{\vec{\omega}},\varsigma)$, exactly like for true scattering events:

[TABLE]

Equations 5 and 12 have now a similar structure: all the samples used to evaluate $f$ can also be used for the evaluation of $\partial_{\varsigma}f$. But in order to complete the evaluation of sensitivity, we must add one sample (of $L$) per collision. Thanks to this similar structure, we can gather them into a single vectorial writting (exactly the same way Eq. 7 was constructed from Eq. 5 and 6):

[TABLE]

The algorithmic-reading of 13 leads to Alg. 3 that is an alternative to Alg. 2 for evaluating simultaneously $f$ and $\partial_{\varsigma}f$. As explained in the algorithmic reading of Eq.7, the recursive nature of this algorithm comes from the fact that the final brackets in the scattering and null-collision terms contain $f$ and $\partial_{\varsigma}f$ at the same location in the same direction. The fact that their sensitivity part includes a summation is translated into an algorithm incrementing the Monte Carlo weight as explained in B. It is interesting to note that the integral $\int p_{L}(l)\mathrm{d}l~{}\partial_{\varsigma}k_{e}(\bm{x_{l}},\varsigma)/p_{L}(l)$ can be more or less difficult to evaluate depending on the profile of $\partial_{\varsigma}k_{e}$. But this can be easily handled using importance sampling based on the $\hat{k}$-adjustment grid, as explained in C.

4 Simulations using the alternative approach

4.1 Transmissivity of a purely absorbing column

Applying the alternative approach of Sec. 3 to the evaluation of column-transmissivities leads to the content of Fig. 3d and Table 5c. For the homogeneous-slab configuration, Fig. 3 shows that not only the pathological behavior of Sec. 2 is removed, but the sensitivity is estimated with a statistical uncertainty that is perfect for a simultaneous evaluation: its dependence on the parameters of the problem is identical to that of the main quantity. As above, in this very simple case, this uncertainly can be expressed analytically (see A.2) and indeed

[TABLE]

For the heterogeneous-slab configuration, the uncertainty cannot be predicted theoretically, but the conclusions of Fig. 5 are identical to those of Fig. 3: in terms of relative accuracy, the convergence rate is equal to that of the algorithm evaluating the main quantity. It is therefore strictly independent of the adjustment of $\hat{k}$ to $k$. The use of an acceleration grid does only what we expect: it reduces the number of null-collisions but does not impact the variance anymore.

4.2 Full radiative transfer in a 3D configuration

In [1], a cubic benchmark configuration was used to test null-collision algorithms when dealing with three-dimension highly-heterogeneous fields for all ranges of optical thickness and single-scattering albedo. We here make use of the same configuration, named heterogeneous-cube hereafter, in order to test our alternative approach with 3D radiation (see Fig. 4):

radiation is monochromatic;

• 2.

the cube is of side $2L$, with $0K$ black faces;

• 3.

the inside-temperature field is such that $f^{eq}$ varies from $f^{eq}=f^{eq}_{max}$ (at the center of the face at $x=-L$) to $f^{eq}=0$ (at $x=L$ and $(y=\pm L,z=\pm L)$) and mimics the shape of flame: $f^{eq}(x,y,z)=\eta(x,y,z)f^{eq}_{max}$ (see the $\eta$ profile in Fig. 4);

• 4.

the fields of absorption and scattering coefficients follows the same spatial dependence: $k_{a}(x,y,z)=\eta(x,y,z)k_{a,max}$ and $k_{a}(x,y,z)=\eta(x,y,z)k_{s,max}$;

• 5.

The single-scattering phase function is that of Henyey-Greenstein with a uniform value of the asymmetry parameter g;

• 6.

$\hat{k}$ is adjusted to $k$ using a regular cubic-grid ($\hat{k}$ uniform within each mesh): the only parameter for $\hat{k}$ is therefore the number of mesh per direction.

The evaluated quantity $A(x,y,z)$ is the stationary net-power density and the free physical parameters are $k_{a,max}L$, $k_{s,max}L$ and $g$. In Table 6 we reproduce the computations of Table 1 in [1], i.e. testing wide ranges of optical thicknesses but fixing $g=0$ (isotropic scattering). In the same table we also provide two sensitivities, $\partial_{k_{a,max}}$ and $\partial_{k_{s,max}}$, that we evaluated simultaneously with $A$. As in [1], although they are not displayed, we checked that simulation results with non-isotropic scattering lead to the exact same conclusions.

These conclusions are very similar to those reached on the slab-transmissivity example: Table 7 highlights the same features as Fig. 5, the standard deviation of Alg. 3 being independent of $\hat{k}$, whereas it increases when $\hat{k}$ gets close to $k_{e}$ for Alg. 2.

5 Conclusion

The simulation examples of the preceding section indicate that the solution proposed in Sec. 3 is practical. Sensitivities can be accurately evaluated even when using null-collision algorithms. We still face convergence difficulties for highly-scattering media, but this is an open question identified in all linear-transport physics, independently of the use of null-collisions.

The way we bypassed the convergence difficulties specific to null-collisions is in rupture with the principle of evaluating sensitivities $\partial_{\varsigma}A$ simultaneously with the main quantity $A$: in all previous works, the two evaluations were truly simultaneous in the sense that the very same samples were used in the $A$ and $\partial_{\varsigma}A$ Monte Carlo algorithms. Here, all the samples used to evaluate $A$ are also used for the evaluation of $\partial_{\varsigma}A$, but in order to complete the evaluation of $\partial_{\varsigma}A$, some additional random variables must be sampled. The sensitivity evaluation has therefore a specific computation-cost. In all our test-cases, this cost was very small compared to the total computation cost, and considering the wide range of tested optical depths we may state that this additional cost has no practical significance for radiative transfer applications. Other linear-transport physicas wold have to be investigated specifically.

In the present text, we only gave very little indications about absolute computer-requirements. We provided some computation-time examples in Fig. 6, without the use of any acceleration grid. Most of our analysis was made with another measure: the average number of random generation per path, which should be proportional to the computation time, but only once all optimisations will be implemented. As illustrated in [7, 9, 8, 10], such data-access considerations in relation with null-collision algorithms (i.e the time associated to a close adjustment of $\hat{k}$ to $k$) is an ongoing computer-science research. The question of implementing the corresponding tools (developped by the computer-graphics community) together with the present sensitivity evaluation algorithm will be discussed in a separate paper.

Appendix A Analytical variances for homogeneous-slab

$\Omega$ is a slab of length $L$ with $\bm{\vec{e_{x}}}$ the incoming normal at location $\bm{\vec{y}}=\bm{0}$ and $f_{\partial\Omega}(L,-\bm{\vec{e_{x}}},.)=f_{inc}$ as boundary condition (see figure 2). The equilibrium distribution $f^{eq}$ is null (no emssion) and the absorption and scattering coefficients are $k_{a}\equiv\varsigma$ and $k_{s}\equiv 0$. The transmissivity is $T=f(\bm{0},\bm{\vec{e_{x}}},.,\varsigma)/f_{inc}$.

A.1 For the standard approachs

Using Eq.(7),

[TABLE]

To get the analytical variance for the standard approach (Alg. 2) in homogeneous-slab, we started by designing a countable categorisation of paths: $\mathcal{C}_{n}$ ($n\in\mathbb{N}^{*}$) is the set of paths that crossed the column and met $n$ null-collisions and $\mathcal{C}$ all other paths (path without null-collisions). The probability $p_{n}$ that a path belongs to the $\mathcal{C}_{n}$ category is

p_{n}=\frac{\hat{k}^{n}}{n!}L^{n}\mathrm{e}^{-\hat{k}L}\Big{(}1-\frac{k_{a}(\varsigma)}{\hat{k}}\Big{)}

The Monte Carlo weight $w_{n}$ obtained at the end of a path belonging to the $\mathcal{C}_{n}$ is

$w_{n}=\frac{-n}{k_{n}}$

The weight of a path belonging to $\mathcal{C}$ is zero. Therefore, the variance of the Alg. 2 is

[TABLE]

A.2 For the alternative approach

Using Eq.(13),

[TABLE]

The variance of a Bernoulli law of $p$ parameter is $p(1-p)$, and therefore the variance of Alg. 3 is

[TABLE]

Similarly for Alg. 1,

[TABLE]

Appendix B Incrementation of the Monte Carlo weights

The standard approach of Alg. 2 is designed from the algorithmic-reading of Eq. 7 and the recursive nature of this algorithm comes from the fact that the final brackets in the scattering and null-collision terms contain $f$ and $\partial_{\varsigma}f$ at the same location in the same direction. As the integral formulation of $f$ and $\partial_{\varsigma}f$ have identical structures, we can just use one unique sampled-path to evaluate simultaneously $f$ and $\partial_{\varsigma}f$. For a better understanding of the meaning of this unique sampled-path, we express here the Monte Carlo weight for one path example: the path displayed in Fig. 1, i .e. one null-collision at location $\vec{\bm{x_{1}}}$, two scatterings at location $\vec{\bm{x_{2}}}$ and $\vec{\bm{x_{3}}}$, and one boundary-collision at location $\vec{\bm{y_{4}}}$. Thanks to the algorithmic-reading of Eq. 7, after one null-collision at location $\vec{\bm{x_{1}}}$. The notations are simplified by writing $"f"$ instead of $w_{"f"}$. We have

[TABLE]

At this stage, $\partial_{\varsigma}f$ and $f$ need to be evaluated at the same location ($\vec{\bm{x_{1}}}$) in the same direction ($\vec{\bm{\omega_{0}}}$). According to Eq. 7, it’s possible to use the same path to evaluate $\partial_{\varsigma}f$ and $f$. Therefore, after one scattering-collision at location $\vec{\bm{x_{2}}}$ in direction $\vec{\bm{\omega_{2}}}$, we have $\partial_{\varsigma}f(\vec{\bm{x_{1}}},\vec{\bm{\omega_{0}}},\varsigma)=\partial_{\varsigma}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)+\frac{\partial_{\varsigma}k_{s}(\vec{\bm{x_{2}}},\varsigma)}{k_{s}(\vec{\bm{x_{2}}},\varsigma)}f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)$ and $f(\vec{\bm{x_{1}}},\vec{\bm{\omega_{0}}},\varsigma)=f(\vec{\bm{x_{2}}},\vec{\bm{\omega_{2}}},\varsigma)$,

[TABLE]

Similarly, after one scattering at location $\vec{\bm{x_{3}}}$ in direction $\vec{\bm{\omega_{3}}}$,

[TABLE]

and after the boundary-collision at location $\vec{\bm{y_{4}}}$ in direction $\vec{\bm{\omega_{3}}}$,

[TABLE]

where $w^{\varsigma}_{\partial\Omega}=0$ and $w_{\partial\Omega}=f_{\partial\Omega}(\vec{\bm{y_{4}}},\vec{\bm{\omega_{3}}},\varsigma)$. The summation of $f$ and $\partial_{\varsigma}f$ is therefore translated into an algorithm incrementing the Monte Carlo weight. This is how the Monte Carlo weight is incremented in the fully recursive algorithm of Alg. 3.

Appendix C The additional sampling

To evaluate $\partial_{\varsigma}T$ with the alternative approach of Alg. 3, we must evaluate $\displaystyle\int_{0}^{\lambda}\partial_{\varsigma}k_{e}(\vec{\bm{x_{l}}},\varsigma)\mathrm{d}l$. We make use of one additional sampling:

[TABLE]

However, depending on the the profile of $\partial_{\varsigma}k_{e}$, this integral can be source of variance if the spatial variations of $p_{L}(l)$ do not sufficiently match those $\partial_{\varsigma}k_{e}(\vec{\bm{x_{l}}}$. We did not yet adress this importance-sampling issue and bypassed it using stratified sampling: we increased the number of ”collisions” by including the grid-intersections to write

[TABLE]

where $n$ is the number of grid-intersections along the path $[\vec{\bm{x_{0}}},\vec{\bm{x_{\lambda}}}]$. For all $i\in[1,n]$, $\vec{\bm{z_{i}}}$ is the location of the $ith$ ”collision”, $\beta_{i}$ is the length between two successive ”collisions” and $\vec{\bm{\zeta_{i}}}=\vec{\bm{z_{i}}}-l_{i}\vec{\bm{\omega}}$ (see Fig. 8).

This stratified sampling is only meaningful if we assume that the grid is sufficiently refined to capture the spatial variation of $\partial_{\varsigma}k_{e}(\vec{\bm{x_{l}}})$. But initially, the grid was only designed to reduce the amount of null-collisions. Therefore, another grid-refinement criterium should be introduced if this strategy is retained. But again, stratified sampling is not the only possible solution and the grid may be used differently: make pre-computations in each mesh and sample $L$ only once along the whole path using importance sampling, which would suppress the increase of computation costs when refining the grid in Fig. 5 ($R$ increasing in Fig. 5c).

Appendix D Another possible approach: $\hat{k}$ varying with $k_{e}$

In Eq. (5), we chose to use a constant upper bound of $k_{e}$ noted $\hat{k}$. Nevertheless, nothing prevents us from varying $\hat{k}$ according to the parameter $\varsigma$, while preserving the fact that $\hat{k}(.,\varsigma)$ is independent of space in order to still enable an easy sampling of $\hat{\Lambda}$. Thus, we have complete freedom in the choice of these variations, and we will choose them to reduce the variance as much as possible.

In addition, we have shown that the explosion is due to the fact that Monte Carlo’s weight includes a factor $1/k_{n}$. This factor appeared when applying the sensitivity evaluation technique to Eq. (5), especially because of the term \partial_{\varsigma}\Big{(}\frac{k_{n}}{\hat{k}}\Big{)}.

Therefore, and now that $\hat{k}$ depends on $\varsigma$, we may impose \partial_{\varsigma}\Big{(}\frac{k_{n}}{\hat{k}}\Big{)}=0, which implies the following relation:

[TABLE]

Thus, by derivating Eq. (5), we obtain:

[TABLE]

with $\psi(\bm{\vec{x}},\lambda,\varsigma)=\left(\frac{\partial_{\varsigma}\hat{k}(\bm{\vec{x}^{\prime}},\varsigma)}{\hat{k}(\bm{\vec{x}^{\prime}},\varsigma)}-\displaystyle\int_{0}^{\lambda}\partial_{\varsigma}\hat{k}(\bm{\vec{x}_{l}},\varsigma)\mathrm{d}l\right)$

We have implemented the algorithm resulting from this integral formulation and compared it with the approach proposed in the core of the present paper. Whatever the case, it was always less efficient.

Acknowledgments

We acknowledge support from the Agence Nationale de la Recherche (ANR, grant HIGH-TUNE ANR-16-CE01-0010, http://www.umr-cnrm.fr/high-tune), from the french Programme National de Teledetection Spatiale (PNTS-2016-05), from Region Occitanie (Projet CLE-2016 EDStar) and from the French Minister of Higher Education, Research and Innovation for the PhD scholarship of the first author.