Robust Multiple Importance Sampling with Tsallis φ-Divergences

Abstract

Multiple Importance Sampling (MIS) combines the probability density functions (pdf) of several sampling techniques. The combination weights depend on the proportion of samples used for the particular techniques. Weights can be found by optimization of the variance, but this approach is costly and numerically unstable. We show in this paper that MIS can be represented as a divergence problem between the integrand and the pdf, which leads to simpler computations and more robust solutions. The proposed idea is validated with 1D numerical examples and with the illumination problem of computer graphics.

1. Introduction and Previous Work

Multiple Importance Sampling (MIS) [1,2] has been proven efficient in Monte Carlo integration. It is able to preserve the advantages of the combined techniques and requires only the calculation of the pdfs of all methods when a sample is generated with one particular method. The weighting scheme applied in MIS depends on the pdfs of the individual techniques and also on the number of samples generated with each of them.

MIS has been applied in a number of rendering algorithms, and its variance is extensively studied [1]. Several estimators have been proposed that are better than balance heuristics with equal sample budgets [3,4,5,6]. Lu et al. [7] proposed an adaptive algorithm for environment map illumination, which used the Taylor series approximation of the variance around an equal sample budget case. In [8,9], different equal sample number strategies were analyzed. Sbert et al. [10] considered the cost associated with the sampling strategies and obtained an adaptive solution by optimizing the variance using the Newton–Raphson method [11]. In [12], the Kullback–Leibler divergence was optimized instead of the variance. Several authors have shown that the variance of an importance sampling estimator is equal to a chi-square divergence [13,14,15,16], which, in turn, can be approximated by the Kullback–Leibler divergence up to the second order [17]. The optimal sample budget has also been targeted with neural networks [15,18]. Recently, a theoretical formula has been elaborated for the weighting functions [19]. In [16], the balance heuristic estimator was generalized by decoupling the weights from the sampling rates, and implicit solutions for the optimal case were given.

These techniques offer lower variance and, therefore, theoretically outperform MIS with equal number of samples. However, equations determining the optimal weighting and sample budget require the knowledge of the integrand and numerical solution methods. In computer graphics, for example, this integrand is not analytically available, so previous discrete samples should be used for the approximation, which introduces errors in the computation. Thus, it is not guaranteed that a theoretically superior estimator also performs better in practice.

This paper proposes an adaptive approach to automatically determine the sampling budgets of the combined methods based on the statistics of previous samples. To improve the robustness, instead of directly optimizing the variance, we make the combined pdf mimic the integrand by minimizing the divergence between them. From the possible alternatives, the Tsallis $\phi$-divergence is used since this leads to a general and stable approach. We show that finding the optimum in all cases, including also Hellinger, chi-square, Kullbach–Leibler divergences, and the variance, boils down to a non-linear equation stating that the $\gamma$-moments of the different combined techniques must be equal. This equation has unique root that can be found by Newton–Raphson iteration.

The organization of the paper is the following. In Section 2, the multi-sample version of the MIS estimator is reviewed where the number of samples allocated to each combined technique is fixed deterministically. Section 3 summarizes the one-sample MIS estimator where the particular sampling method is also selected randomly. In Section 4, the adaptation of the MIS weighting functions is discussed and we argue that the adaptation should be controlled not by the estimated variance but by an appropriately chosen divergence. Section 5 reformulates the MIS problem as the optimization of divergences and introduces the $\gamma$-moments. In Section 6, the details of the numerical optimization are discussed. Finally, the results are demonstrated with numerical examples and the direct lighting solution of computer graphics.

2. Multi-Sample Balance Heuristic Mis

Suppose we wish to estimate the value of integral $\int f\left(x\right)\mathrm{d}x$ and have m proposal pdfs $p_i\left(x\right)$ to generate the jth sample $X_{ij}$ of method i in the domain of the integral. With method i, $n_i$ independent samples are drawn, so the total number of samples is ${\sum }_{i=1}^m{n}_i=N$. When the sample numbers $n_i$ are fixed deterministically, the approach is called multi-sample.

The multi-sample MIS estimator [1] has the following expression:

$$F=\sum _{i=1}^m\frac{1}{n_i}\sum _{j=1}^{n_i}{w}_i\left(X_{ij}\right)\frac{f\left(X_{ij}\right)}{p_i\left(X_{ij}\right)},$$

(1)

where the weights satisfy the normalization condition:

$$\sum _{i=1}^m{w}_i\left(x\right)=1.$$

Integral estimator F is unbiased, as its expected value $\mu$ is equal to the integral:

$$\mu =E\left[F\right]=\sum _{i=1}^m\frac{1}{n_i}\sum _{j=1}^{n_i}\int \frac{w_i\left(x\right)f\left(x\right)}{p_i\left(x\right)}{p}_i\left(x\right)\mathrm{d}x=\sum _{i=1}^m{\mu }_i=\int f\left(x\right)\mathrm{d}x$$

(2)

where ${\mu }_i=\int w_i\left(x\right)f\left(x\right)\mathrm{d}x$. The variance of the integral estimator is

$$V\left[F\right]=\sum _{i=1}^m\int \frac{w_i^2\left(x\right)f^2\left(x\right)}{n_i{p}_i\left(x\right)}\mathrm{d}x-\sum _{i=1}^m\frac{{\mu }_i^2}{n_i}.$$

(3)

The variance of the estimator depends on the combination weights $w_i\left(x\right)$. The first step of their definition is to propose an algebraic form. The balance heuristic states that the weight of method i at domain point x should be proportional to the density of samples generated by method i in this point:

$$w_i\left(x\right)=\frac{{\alpha }_i{p}_i\left(x\right)}{p(\alpha ,x)}$$

(4)

where ${\alpha }_i=n_i/N$ is the fraction of the samples allocated to method i, and $p(\alpha ,x)$ is the mixture pdf:

$$p(\alpha ,x)=\sum _{k=1}^m{\alpha }_k{p}_k\left(x\right).$$

(5)

Substituting this weighting function into the MIS estimator formulas, we obtain the multi-sample balance heuristics estimator:

$$F=\frac{1}{N}\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{f\left(X_{i,j}\right)}{p(\alpha ,X_{i,j})}.$$

(6)

The variance of the balance heuristics estimator is

$$V\left[F\right]=\frac{1}{N}\left(\int \frac{f^2\left(x\right)}{p(\alpha ,x)}\mathrm{d}x-\sum _{i=1}^m{\alpha }_i{{\mu }^{\prime }}_i^2\right),$$

(7)

where ${{\mu }^{\prime }}_i={\mu }_i/{\alpha }_i$.

3. One-Sample Balance Heuristic MIS

In the one-sample version of MIS, the numbers of samples used by particular techniques are not determined before starting the sampling process, but for each sample, the sampling method itself is selected randomly with the probability of fraction parameter ${\alpha }_i$. As this approach introduces additional randomization, its variance is larger than that of the multi-sample version, but can also be used in cases when the number of total samples is less than the number of combined methods.

The one-sample balance heuristic estimator is given by

$$\mathcal{F}=\frac{f\left(X\right)}{p(\alpha ,X)}$$

(8)

with variance

$$V\left[\mathcal{F}\right]=\int \frac{f^2\left(x\right)}{p(\alpha ,x)}\mathrm{d}x-{\mu }^2.$$

(9)

4. Adaptive MIS

Having the algebraic form of the weight function, the task is to find the optimal fractions ${\alpha }_i$ with the constraint that the sum of sample numbers must be equal to the total sample budget, i.e., ${\sum }_{i=1}^m{\alpha }_i=1$. One possibility is to use some heuristics before starting the sampling process. A simple example is the equal sample count MIS that sets all fractions to be equal.

Adaptive techniques, on the other hand, try to minimize the error by controlling the fractions on-the-fly as new samples introduce additional information. During this, the following issues need to be considered:

• The integrals in the variance cannot be computed analytically, but must be estimated from the available finite number of samples drawn from the sampling distribution. This uncertainty may significantly affect the goodness of the final results.
• The optimization process should be fast and should not introduce too high overhead.

Unfortunately, the direct optimization for the variance does not meet these requirements. The integral of the variance can be estimated from the available samples $X_i$ drawn from mixture distribution $p(\alpha ,x)$ as:

$$\int \frac{f^2\left(x\right)}{p(\alpha ,x)}\mathrm{d}x=E_p\left[\frac{f^2\left(x\right)}{p^2(\alpha ,x)}\right]\approx \frac{1}{N}\sum _{i=1}^N{\left(\frac{f\left(X_i\right)}{p(\alpha ,X_i)}\right)}^2$$

(10)

where $E_p\left[\xi \right]$ is the expectation of random variable $\xi$ of pdf p. In this estimation, the ratios of the integrand and the pdf are squared causing high fluctuation and making the optimization process unstable.

Thus, instead of the variance, we need other optimization targets with more robust estimators in order to find the fractions ${\alpha }_k$ during adaptation.

5. MIS as Divergence between Distributions

MIS consists in looking for a distribution $p(\alpha ,x)$ that mimics $f\left(x\right)$ as much as possible. If we have non-negative integrand $f\left(x\right)\ge 0$ with integral $\mu =\int f\left(x\right)dx$, then the integrand scaled down by the integral $g\left(x\right)=f\left(x\right)/\mu$ is also a distribution. The MIS problem can be stated as finding mixture pdf $p(\alpha ,x)$ that minimizes the divergence between two distributions.

5.1. $\phi$-Divergences

A $\phi$-divergence ${\mathcal{D}}_{\phi }(q,r)$[20,21] is a measure of the dissimilarity between two probability distributions $q\left(x\right)$ and $r\left(x\right)$. It is defined by a convex function $\phi \left(t\right)$ satisfying $\phi \left(1\right)=0$:

$${\mathcal{D}}_{\phi }(q,r)=\int r\left(x\right)\phi \left(\frac{q\left(x\right)}{r\left(x\right)}\right)\mathrm{d}x.$$

(11)

$\phi$-divergences are always positive (by Jensen inequality) and are zero iff $q\equiv r$ (for strict convexity) [21].

For instance, for $\phi \left(t\right)=tlogt$, we obtain the Kullback–Leibler divergence:

$${\mathcal{D}}_{\phi }(q,r)=\int q\left(x\right)log\frac{q\left(x\right)}{r\left(x\right)}\mathrm{d}x=KL(q,r).$$

(12)

Note that $\phi$-divergence is not symmetric, thus we have in principle two options to measure the dissimilarity of probability densities q and r, using either ${\mathcal{D}}_{\phi }(q,r)$ or ${\mathcal{D}}_{\phi }(r,q)$. Any of these two options can be used in practice, although they are not independent. If we define ${\phi }^{*}\left(t\right)=t\phi \left(\frac{1}{t}\right)$ then ${\mathcal{D}}_{{\phi }^{*}}(r,q)={\mathcal{D}}_{\phi }(q,r)$[22]. For $\phi \left(t\right)=tlogt$, we have ${\phi }^{*}\left(t\right)=-logt$.

The objective of importance sampling is to make probability density $p(\alpha ,x)$ similar to scaled integrand $g\left(x\right)=f\left(x\right)/\mu$, i.e., to minimize $\phi$-divergence ${\mathcal{D}}_{\phi }(g,p\left(\alpha \right))$. A $\phi$-divergence is a convex function in both arguments and $p\left(\alpha \right)$ is a linear function of $\alpha$, thus the functional is convex in $\alpha$. This will ensure the existence of a minimum. With choosing $\phi \left(t\right)=t^2-1$, we can get the variance back as a special case of this more general approach:

$$V\left[\mathcal{F}\right]=\int \frac{f^2\left(x\right)}{p(\alpha ,x)}\mathrm{d}x-{\mu }^2={\mu }^2\int p(\alpha ,x)\left({\left(\frac{g\left(x\right)}{p(\alpha ,x)}\right)}^2-1\right)\mathrm{d}x={\mu }^2{\mathcal{D}}_{\phi }(g\left(x\right),p(\alpha ,x)).$$

(13)

The optimal MIS problem is reduced to find the minimum of ${\mathcal{D}}_{\phi }(g,p\left(\alpha \right))$ with the constraints ${\sum }_{i=1}^m{\alpha }_i=1$ and ${\alpha }_i\ge 0$. Using Lagrange multipliers, the solution is

$$\frac{\partial {\mathcal{D}}_{\phi }(g,p\left(\alpha \right))}{\partial {\alpha }_i}-\lambda =0,$$

(14)

which means that the derivatives of the divergence with respect to fractions ${\alpha }_i$ must be the same for each sampling method i.

The derivative can be expressed as

$$\frac{\partial {\mathcal{D}}_{\phi }(g,p\left(\alpha \right))}{\partial {\alpha }_i}=\int p_i\left(x\right)\left[\phi \left(\frac{g\left(x\right)}{p(\alpha ,x)}\right)-\frac{g\left(x\right)}{p(\alpha ,x)}{\phi }^{\prime }\left(\frac{g\left(x\right)}{p(\alpha ,x)}\right)\right]\mathrm{d}x.$$

(15)

This integral is the expectation when samples are generated with pdf $p_i\left(x\right)$, thus we can write

$$E_{p_i}\left[\phi \left(\frac{g}{p\left(\alpha \right)}\right)-\frac{g}{p\left(\alpha \right)}{\phi }^{\prime }\left(\frac{g}{p\left(\alpha \right)}\right)\right]=E_{p\left(\alpha \right)}\left[\phi \left(\frac{g}{p\left(\alpha \right)}\right)-\frac{g}{p\left(\alpha \right)}{\phi }^{\prime }\left(\frac{g}{p\left(\alpha \right)}\right)\right]=\lambda .$$

(16)

Equation (16) only holds when ${\alpha }_i>0$, i.e., in the interior of the simplex domain, not on the border.

Swapping the two distributions in the $\phi$-divergence, we have another measure for the dissimilarity between $p\left(\alpha \right)$ and g, ${\mathcal{D}}_{\phi }(p\left(\alpha \right),g)$. Using Lagrange multipliers again, this leads to the following equation for unknown fractions ${\alpha }_i$

$$\frac{\partial {\mathcal{D}}_{\phi }(p\left(\alpha \right),g)}{\partial {\alpha }_i}=\int p_i{\phi }^{\prime }\left(\frac{p\left(\alpha \right)}{g}\right)\mathrm{d}x={\lambda }^{*}.$$

(17)

Using expected values, this option is expressed as

$$E_{p_i}\left[{\phi }^{\prime }\left(\frac{p\left(\alpha \right)}{g}\right)\right]=E_{p\left(\alpha \right)}\left[{\phi }^{\prime }\left(\frac{p\left(\alpha \right)}{g}\right)\right]=\lambda$$

(18)

The problem with the optimization of the variance was that in the resulting equation the ratios of the sampled integrand and the pdf are squared (see Equation (10)), making the equation sensitive to non-optimal pdfs. With the proper definition of $\phi$, we can solve this issue. For example, if $\phi \left(t\right)=-logt$, then ${\phi }^{\prime }=-1/t$, and the condition of optimality is

$$E_{p_i}\left[\frac{g}{p\left(\alpha \right)}\right]=E_p\left[\frac{g}{p\left(\alpha \right)}\right]=-{\lambda }^{*}.$$

(19)

Remembering that $g\left(x\right)=f\left(x\right)/\mu$, the last equation can also be written as

$$E_{p_i}\left[\frac{f}{p\left(\alpha \right)}\right]=E_p\left[\frac{f}{p\left(\alpha \right)}\right],$$

(20)

which contains the ratios without squaring. This gives back the solution in [23], which is also conjectured in [16]. To further exploit this idea, we take the parametric family of Tsallis divergences. With its parameter the compromise between the robustness of the optimization equation and the similarity to the variance can be controlled.

5.2. Tsallis Divergences

Tsallis divergence [24] is associated with the Tsallis entropy [25]. Let us consider the one-parameter Tsallis divergence between distributions $q\left(x\right)$ and $r\left(x\right)$:

$$\begin{array}{ccc}\hfill {\mathcal{D}}_{\gamma }^T(q,r)& =& \frac{1}{\gamma -1}\left(\int \frac{q^{\gamma }\left(x\right)}{r^{\gamma -1}\left(x\right)}\mathrm{d}x-1\right)\hfill \\ \hfill & =& \frac{1}{\gamma -1}\int \left(\frac{q^{\gamma }\left(x\right)}{r^{\gamma -1}\left(x\right)}-r\left(x\right)\right)\mathrm{d}x=\frac{1}{\gamma -1}\int r\left(x\right)\left(\frac{q^{\gamma }\left(x\right)}{r^{\gamma }\left(x\right)}-1\right)\mathrm{d}x.\hfill \end{array}$$

(21)

Tsallis divergence ${\mathcal{D}}_{\gamma }^T(q,r)$ is a $\phi$-divergence defined by

$$\phi \left(t\right)=\frac{t^{\gamma }-1}{\gamma -1}.$$

Parameter $\gamma$ should satisfy $\gamma >0$ and $\gamma \ne 1$, which guarantees that $\phi \left(t\right)$ exists and is convex.

It can be shown that ${lim}_{\gamma \to 1}{\mathcal{D}}_{\gamma }^T(q,r)=KL(q,r)$ is the Kullbach–Leibler divergence. When $\gamma =2$, Tsallis divergence turns to be the chi-square divergence:

$$\begin{array}{ccc}\hfill {\chi }^2(q,r)& =& \int \frac{{(q\left(x\right)-r\left(x\right))}^2}{r\left(x\right)}\mathrm{d}x\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & =\int \frac{q^2\left(x\right)}{r\left(x\right)}\mathrm{d}x-2\int q\left(x\right)\mathrm{d}x+\int r\left(x\right)\mathrm{d}x\hfill \\ \hfill & =\int \frac{q^2\left(x\right)}{r\left(x\right)}\mathrm{d}x-1={\mathcal{D}}_2^T(q,r).\hfill \end{array}$$

(23)

For $\gamma =1/2$, Tsallis divergence ${\mathcal{D}}_{1/2}^T(q,r)$ becomes the squared Hellinger distance $H^2(q,r)$:

$$\begin{array}{ccc}\hfill {\mathcal{D}}_{1/2}^T(q,r)& =& -2\int (q{\left(x\right)}^{1/2}r{\left(x\right)}^{1/2}-q\left(x\right))\mathrm{d}x\hfill \\ \hfill & =& 2\left(1-\int q{\left(x\right)}^{1/2}r{\left(x\right)}^{1/2}\mathrm{d}x\right)=2H^2(q,r).\hfill \end{array}$$

(24)

Note that $H(q,r)$ is a true metric.

5.3. Optimal MIS Solution with Tsallis Divergence

Optimal MIS needs to find fractions ${\alpha }_i$ that minimize ${\mathcal{D}}_{\gamma }^T(g\left(x\right),p(\alpha ,x))$ with the constraint ${\sum }_{i=1}^m{\alpha }_i=1$. Substituting function $\phi \left(t\right)$ of Tsallis divergence into Equation (16), we obtain that for all i, the quantities, called γ-moments,

$$\begin{array}{c}\hfill M_{\gamma ,i}=\int {\left(\frac{f\left(x\right)}{p(\alpha ,x)}\right)}^{\gamma }{p}_i\left(x\right)\mathrm{d}x=E_{p_i}\left[{\mathcal{F}}^{\gamma }\right],\end{array}$$

(25)

have to be equal. Equation (25) guarantees a global minimum, as the function ${\mathcal{D}}_{\gamma }^T(g\left(x\right),p(\alpha ,x))$ is convex for all $\gamma >0$. Considering the case when $g\left(x\right)$ and $p(\alpha ,x)$ are swapped in the divergence, i.e., substituting function $\phi \left(t\right)$ of Tsallis into Equation (17), the requirement of the equality of the $\gamma$-moments is established again, but now $\gamma$ should be replaced by $1-\gamma$. Examining Equation (20) we can realize that optimizing for $\gamma =1$ is the same as minimizing the cross entropy or the Kullbach–Leibler divergence. Thus, the criterion of the equality of the $\gamma$-moments covers all cases.

5.4. Other $\phi$-Divergences: Total Variation Divergence, ${\chi }_P^k$ and ${\left|\chi \right|}_P^k$ Divergences

The total variation can also be interpreted as a $\phi$-divergence defined by $\phi \left(t\right)=\frac{1}{2}|t-1|$:

$$TV(q,r)=\frac{1}{2}\int |q\left(x\right)-r\left(x\right)|\mathrm{d}x$$

This is the only $\phi$-divergence that is a true metric [26]. The condition of optimality is

$$\int \mathrm{sign}(g\left(x\right)-p(\alpha ,x))p_i\left(x\right)\mathrm{d}x=\int \mathrm{sign}(g\left(x\right)-p(\alpha ,x))p(\alpha ,x)\mathrm{d}x$$

(26)

for all i, although we can hardly use Equation (26) to approximate a solution.

$\phi$-divergences also include the Pearson–Vajda ${\chi }_P^k$ and ${\left|\chi \right|}_P^k$ divergences:

$${\chi }_P^k(q,r)=\int \frac{{(q\left(x\right)-r\left(x\right))}^k}{q^{k-1}\left(x\right)}\mathrm{d}x,{\left|\chi \right|}_P^k(q,r)=\int \frac{{|q\left(x\right)-r\left(x\right)|}^k}{q^{k-1}\left(x\right)}\mathrm{d}x$$

defined by the functions $\phi \left(t\right)={(t-1)}^k$ and $\phi \left(t\right)={|t-1|}^k$, respectively [17]. For $k=2$, we have ${\chi }_P^2(q,r)={\left|\chi \right|}_P^2(q,r)={\chi }^2(r,q)$, and for $k=1$, identities ${\chi }_P^1(q,r)=0$ and ${\left|\chi \right|}_P^1(q,r)=2TV(q,r)$ hold.

6. Finding the Optimal Solution

For the sake of simplicity, we assume that two sampling methods of pdfs $p_1\left(x\right)$ and $p_2\left(x\right)$, associated with fractions $\alpha$ and $1-\alpha$ are combined. The mixture density is

$$p(\alpha ,x)=\alpha p_1\left(x\right)+(1-\alpha )p_2\left(x\right).$$

From the optimality condition of Equation (25), we obtain

$$\int {\left(\frac{f\left(x\right)}{p(\alpha ,x)}\right)}^{\gamma }(p_1\left(x\right)-p_2\left(x\right))\mathrm{d}x=\zeta \left(\alpha \right)=0.$$

(27)

This is a non-linear equation for unknown fraction $\alpha$, which is solved with Newton–Raphson iteration. The derivative of $\zeta \left(\alpha \right)$ is

$${\zeta }^{\prime }\left(\alpha \right)=\frac{\partial \zeta \left(\alpha \right)}{\partial \alpha }=-\gamma \int \frac{f^{\gamma }\left(x\right)}{p^{\gamma +1}(\alpha ,x)}{(p_1\left(x\right)-p_2\left(x\right))}^2\mathrm{d}x,$$

(28)

as $\partial p(\alpha ,x)/\partial \alpha =p_1\left(x\right)-p_2\left(x\right)$.

The value of $\zeta \left(\alpha \right)$ and its derivative can be estimated in the following way. Let us re-write $\zeta \left(\alpha \right)$ as an expectation:

$$\zeta \left(\alpha \right)=\int \frac{f^{\gamma }\left(x\right)}{p^{\gamma +1}(\alpha ,x)}(p_1\left(x\right)-p_2\left(x\right))p(\alpha ,x)\mathrm{d}x=E_{p\left(\alpha \right)}\left[\frac{f^{\gamma }\left(X\right)}{p^{\gamma +1}(\alpha ,X)}(p_1\left(X\right)-p_2\left(X\right))\right].$$

(29)

Taking N samples $\{X_1,X_2,\dots ,X_N\}$ according to pdf $p(\alpha ,x)$, we have the estimator

$$\widehat{\zeta \left(\alpha \right)}=\frac{1}{N}\sum _{i=1}^N\frac{f^{\gamma }\left(X_i\right)}{p^{\gamma +1}(\alpha ,X_i)}(p_1\left(X_i\right)-p_2\left(X_i\right)),$$

(30)

and in the same way,

$$\widehat{{\zeta }^{\prime }\left(\alpha \right)}=-\frac{\gamma }{N}\sum _{i=1}^N\frac{f^{\gamma }\left(X_i\right)}{p^{\gamma +2}(\alpha ,X_i)}{(p_1\left(X_i\right)-p_2\left(X_i\right))}^2.$$

(31)

The integrals are estimated with an iterative algorithm updating $\alpha$ after each iteration. Starting with ${\alpha }^{\left(0\right)}=1/2$, in iteration k we draw N samples according to $p({\alpha }^{(k-1)},x)$, compute $\widehat{\zeta \left({\alpha }^{(k-1)}\right)}$ and $\widehat{{\zeta }^{\prime }\left({\alpha }^{(k-1)}\right)}$, and use the Newton–Raphson formula to obtain updated fraction ${\alpha }^{\left(k\right)}$:

$${\alpha }^{\left(k\right)}={\alpha }^{(k-1)}-\frac{\widehat{\zeta \left({\alpha }^{(k-1)}\right)}}{\widehat{{\zeta }^{\prime }\left({\alpha }^{(k-1)}\right)}}.$$

(32)

7. Numerical Examples

We present here three 1D examples. We used the Mathematica package for the computations.

7.1. Example 1

Suppose we want to evaluate the integral (see Figure 1)

$$\mu ={\int }_{0.01}^{3.5\pi }(\sqrt{x}+sinx)\mathrm{d}x\approx 25.3065$$

(33)

by MIS combining sampling methods of pdfs $Gauss(2,1)$ and $Gauss(8,2)$, where $Gauss(m,\sigma )$ stands for the normal distribution of mean m and standard deviation $\sigma$. For this example, equal sample number MIS has variance $V\left[F\right]=24.1152$, and the minimum variance is $V\left[F\right]=13.4788$.

In Figure 2, we show the variances $V\left[F\right]$ and $V\left[\mathcal{F}\right]$ for the optimal solution of Tsallis divergences from $\gamma =0.1$ to $\gamma =8$ in steps of 0.1. Dots depict the values of $V\left[F\right]$ for the optimal $\alpha$ fractions for each parameter $\gamma$ using 5 Newton–Raphson iterations with 100 total samples in each iteration.

Figure 2. Example 1: Values of $V\left[F\right]$ (dashed line) and $V\left[\mathcal{F}\right]$ (continuous line) for the solution of Equation (25) for each value of parameter $\gamma$ from 0.1 to 8, in the horizontal axis. Variance $V\left[\mathcal{F}\right]$ (continuous line) is minimal when $\gamma =2$, as $V\left[\mathcal{F}\right]$ corresponds to ${\chi }^2$ divergence (Equation (13)). Dots show the variance $V\left[F\right]$ for the optimal $\alpha$ values for each parameter $\gamma$ using 5 Newton–Raphson iterations with 100 total samples in each iteration. Horizontal line corresponds to the minimum variance $V\left[F\right]=13.4788$. A zoom-out is shown in Figure 3.

Figure 2. Example 1: Values of $V\left[F\right]$ (dashed line) and $V\left[\mathcal{F}\right]$ (continuous line) for the solution of Equation (25) for each value of parameter $\gamma$ from 0.1 to 8, in the horizontal axis. Variance $V\left[\mathcal{F}\right]$ (continuous line) is minimal when $\gamma =2$, as $V\left[\mathcal{F}\right]$ corresponds to ${\chi }^2$ divergence (Equation (13)). Dots show the variance $V\left[F\right]$ for the optimal $\alpha$ values for each parameter $\gamma$ using 5 Newton–Raphson iterations with 100 total samples in each iteration. Horizontal line corresponds to the minimum variance $V\left[F\right]=13.4788$. A zoom-out is shown in Figure 3.

Figure 3. Example 1. Zoom-out of Figure 2. Observe that $V\left[\mathcal{F}\right]<V\left[F\right]$ except for $\gamma =1$ where $V\left[\mathcal{F}\right]=V\left[F\right]$. For the equal sample budget case, $V\left[F\right]=24.1152$, while the minimum is $V\left[F\right]=13.4788$.

Figure 3. Example 1. Zoom-out of Figure 2. Observe that $V\left[\mathcal{F}\right]<V\left[F\right]$ except for $\gamma =1$ where $V\left[\mathcal{F}\right]=V\left[F\right]$. For the equal sample budget case, $V\left[F\right]=24.1152$, while the minimum is $V\left[F\right]=13.4788$.

7.2. Example 2

Let us consider integral (see Figure 4)

$$\mu ={\int }_{-4}^{4}Gauss(-1.5,1)+2Gauss(1.5,0.75)\mathrm{d}x\approx 2.9929$$

(34)

by MIS using functions $Gauss(-1.5,1)$ and $Gauss(1.5,0.75)$. For this example, equal sample number MIS has a variance of $V\left[F\right]=0.1134$, and the minimum variance is $V\left[F\right]=0$. Figure 5 and Figure 6 present the optimal solutions and the results after 5 Newton–Raphson iterations with 100 total samples in each iteration.

7.3. Example 3

Finally, consider the approximation of the following integral (see Figure 7)

$$\begin{array}{c}\hfill \mu ={\int }_{0.01}^{\pi /2}\left(\sqrt{x}+sinx\right)\mathrm{d}x\approx 2.3118.\end{array}$$

(35)

by MIS using functions $2-x$, and ${sin}^2\left(x\right)$. For this example, equal sample budget MIS has a variance of $V\left[F\right]=0.2772$, and the minimum variance is $V\left[F\right]=0.09032$. In Figure 8 and Figure 9 we show the variances $V\left[F\right]$ and $V\left[\mathcal{F}\right]$ for the optimal solution for parameters from $\gamma =0.1$ to $\gamma =8$. Dots are the results of $V\left[F\right]$ using 5 Newton–Raphson iterations with 100 total samples in each iteration.

8. Combination of Light Source Sampling and BRDF Sampling in Computer Graphics

Here we consider a classical problem of computer graphics, the combination of light source sampling and BRDF (Bidirectional Reflectance Distribution Function) sampling. To compute the reflected radiance $L^r(\overrightarrow{p},\overrightarrow{\omega })$ of a surface point $\overrightarrow{p}$ at direction $\overrightarrow{\omega }$, the following integral should be evaluated in the hemispherical domain of incident directions $\mathsf{\Omega }$:

$$L^r(\overrightarrow{p},\overrightarrow{\omega })={\int }_{\mathsf{\Omega }}L({\overrightarrow{p}}^{\prime },{\overrightarrow{\omega }}^{\prime })f_r({\overrightarrow{\omega }}^{\prime },\overrightarrow{p},\overrightarrow{\omega })cos{\theta }^{\prime }\mathrm{d}{\omega }^{\prime }$$

(36)

where $L({\overrightarrow{p}}^{\prime },{\overrightarrow{\omega }}^{\prime })$ is the radiance of point ${\overrightarrow{p}}^{\prime }$ visible from point $\overrightarrow{p}$ in incident direction $-{\overrightarrow{\omega }}^{\prime }$, $f_r({\overrightarrow{\omega }}^{\prime },\overrightarrow{p},\overrightarrow{\omega })$ is the BRDF expressing the portion of the light beam that is reflected from direction ${\overrightarrow{\omega }}^{\prime }$ to $\overrightarrow{\omega }$ at point $\overrightarrow{p}$ and ${\theta }^{\prime }$ is the angle between the surface normal at point $\overrightarrow{p}$ and incident direction $-{\overrightarrow{\omega }}^{\prime }$. While solving such problems, we have two sampling methods, $p_1\left({\overrightarrow{\omega }}^{\prime }\right)$ approximately mimicking the $f_r(\overrightarrow{\omega },\overrightarrow{p},{\overrightarrow{\omega }}^{\prime })cos{\theta }^{\prime }$ factor and $p_2\left({\overrightarrow{\omega }}^{\prime }\right)$ mimicking the incident radiance $L({\overrightarrow{p}}^{\prime },{\overrightarrow{\omega }}^{\prime })$.

MIS would use a combined pdf in the following form

$$p(\alpha ,{\overrightarrow{\omega }}^{\prime })=\alpha p_1\left({\overrightarrow{\omega }}^{\prime }\right)+(1-\alpha )p_2\left({\overrightarrow{\omega }}^{\prime }\right).$$

(37)

where the task is the optimal determination of weight $\alpha$.

In order to test the proposed method, we render the classic scene of Veach [1] with combined light source and BRDF sampling. The shiny rectangles have max-Phong BRDF [27] with shininess parameters 50, 100, 500, and 1000, respectively. The four spherical light sources emit the same power.

For each pixel, we used 50 samples in total organized in 5 iterations. The process starts with 5 BRDF and 5 light source samples per pixel, and the per-pixel $\alpha$ weights are updated at the end of each iteration. Figure 10 shows the rendered images together with the $\alpha$-maps, and we compare the original sampling techniques, equal count MIS, variance minimization, and Tsallis divergence. The reference obtained by high number of samples is shown by Figure 11.

Table 1. RMSE of different MIS methods.

Method	RMSE
BRDF sampling	613
Light source sampling	295
Equal count MIS	244
Tsallis ($\gamma =0.01$)	216
Tsallis ($\gamma =0.1$)	213
Tsallis ($\gamma =0.5$) (Hellinger)	215
Tsallis ($\gamma =1$) (Kullback-Leibler)	218
Tsallis ($\gamma =1.5$)	218
Tsallis ($\gamma =2$) (chi-square, variance)	221
Tsallis ($\gamma =2.5$)	254

depicts the Root Mean Square Error (RMSE) values. Figure 12 shows the convergence plots.

Finally, we note that $\phi$-divergences were used for global illumination as oracles for adaptive sampling in [28].

9. Discussion

The goal of adaptive MIS is to find the weights of given pdfs that eventually lead to minimal integration error. The definition of the error may be application dependent, in this paper we chose the RMSE, which describes the variance of the estimator. We argued that instead of the direct optimization of the estimated variance, it is better to use other measures for the optimization target since the variance estimation may strongly fluctuate if the sample number is moderate, which prevents the adaptation process from converging to the real optimum. The alternative criterion is the Tsallis divergence, which includes important particular divergence types and by setting its parameter, more robust targets can be defined.

In order to demonstrate the proposed method, we took three numerical examples and a fundamental problem of computer graphics. We can observe in all examples that the adaptive MIS can significantly outperform the equal sample count MIS, thus examining adaptation strategies has theoretical and practical relevance. In the first numerical example, none of the combinations of the two sampling pdfs could well mimic the integrand, thus the computation errors are quite high and the influence of parameter $\gamma$ is small. In the second numerical example, the integrand is proportional to an appropriate linear combination of the two sampling pdfs, thus zero variance estimator is possible. If parameter $\gamma$ is not too small, the method could indeed compute the integral with high accuracy. This means that in easy integration problems, Tsallis divergence can be safely used. The third numerical example is a more difficult, realistic case when sufficiently high number of samples are used in each iteration step to estimate the optimization target. As expected, the optimal value of parameter $\gamma$ is close to 2 since this is the point where Tsallis divergence becomes the variance $V\left[\mathcal{F}\right]$.

In the image synthesis example, the integrand is not a weighted sum but a product where the sampling pdfs approximately mimic the factors. In this case, no zero variance MIS is possible. Unlike in the numerical examples, we took small number of samples for each pixel according to practical scenarios. The small sample number made the optimization target equation unreliable unless parameter $\gamma$ is reduced from 2 that corresponds to the variance. Looking at Figure 10, we can observe that when $\gamma$ becomes smaller, the computed weight maps become slightly less accurate in average but also less noisy, which eventually reduces the computation error.

We have also seen in the examples that the differences between the variances corresponding to the optimal values for the different parameters $\gamma$ of the Tsallis divergence are small. This is in accordance with any $\phi$-divergence being a second order approximation of ${\chi }^2$ divergence [17], thus obtaining any minimum within a wide range of $\gamma$ parameters would be a good solution. Aside from $\gamma =2$, which corresponds to $V\left[\mathcal{F}\right]$, the most interesting Tsallis divergence from a practical perspective is the Kullback-Leibler divergence, corresponding to $\gamma =1$.

10. Conclusions and Future Work

We have shown in this paper that finding the optimal weights for MIS can be presented as finding the minimum $\phi$-divergence between the normalized integrand and the importance pdf. Being convex, a $\phi$-divergence allows convex optimization. We have singled out the Tsallis divergences with parameter $\gamma$ as they have the further advantage that the equations for the optimal value consist of making the $\gamma$-moments equal. We have tested our results with numerical examples and solving an illumination problem.

In our future work we will investigate the optimal efficiencies, i.e., taking into account the cost of sampling. This has been done for the variance $V\left[\mathcal{F}\right]$ (i.e., $\gamma =2$) in [4,6] where the implicit equation for the optimal solution was given. However, the efficiency equation does not allow convex optimization. Using the $\gamma$-moments, or a different family of $\phi$-divergences, could help in obtaining a robust solution for the efficiency.