Regression-based Monte Carlo Integration

Corentin Salaün

Max-Planck-Institut für Informatik

Adrien Gruson

McGill University & École de Technologie Supérieure

Binh-Son Hua

VinAI Research

Toshiya Hachisuka

University of Waterloo

Gurprit Singh

Max-Planck-Institut für Informatik

Given an integrand (a), we first sample (b) $f(x)$ as in Monte Carlo (MC) integration.
(c-d) Traditional MC estimator can be interpreted as fitting a  constant model function to the sample values, with the integral of this constant function equals to $\Integral$.
(e) We, instead, propose to use a non-constant model function such as a polynomial, which
is then fitted to the sampled values.
(f) The resulting estimator is based on control variates; we add the analytical integral of the model function to MC integration of the difference between the original integrand and the model function.
The bottom row shows renderings and the corresponding error images to demonstrate the impact of our regression approach against the traditional MC integration.
The insets on the right compare our method with different orders ($O$x) of polynomials.
Our method has significant error reduction at equal time.

Monte Carlo integration is typically
interpreted as an estimator of the expected value using stochastic samples.
There exists an alternative interpretation in calculus where Monte Carlo integration can be seen as estimating a constant function—from the stochastic evaluations of the integrand—that integrates to the original integral.
The integral mean value theorem states that this constant function should be the mean (or expectation) of the integrand.
Since both interpretations result in the same estimator, little attention has been devoted to the calculus-oriented interpretation.
We show that the calculus-oriented interpretation actually implies the possibility of using a more complex function than a constant one to construct a more efficient estimator for Monte Carlo integration.
We build a new estimator based on this interpretation and relate our estimator to control variates with least-squares regression on the stochastic samples of the integrand.
Unlike prior work, our resulting estimator is provably better than or equal to the conventional Monte Carlo estimator.
To demonstrate the strength of our approach, we introduce a practical estimator that can act as a simple drop-in replacement for conventional Monte Carlo integration. We experimentally validate our framework on various light transport integrals.
The code is available at <https://github.com/iribis/regressionmc>.

<concept_desc>Computing methodologies Rendering</concept_desc>

[500]Computing methodologies Rendering

§ INTRODUCTION

MC integration is a basic tool for many numerical integration problems.
In textbooks, MC integration is typically interpreted as a statistical process where it samples the integrand and estimates the expected value by taking the average of these samples.
While less common, there is another interpretation of MC integration based on the mean value theorem in calculus.
The mean value theorem for definite integrals tells us that there exists a constant that integrates exactly to the same solution as the given definite integral (see <ref>).
This constant, by definition, is equal to the mean of the integrand, which can be again estimated by taking the average of random samples, and thus both interpretations result in the same estimator.
This alternative interpretation, however, raises interesting questions: Can we use a non-constant function to construct a MC estimator? If so, what happens?
These questions do not appear in the statistics-oriented interpretation as the expected value is fundamentally tied to the idea of samples' average.

We propose a new MC estimator that takes the same input as MC integration (see <ref>b) with no extra assumptions.

Instead of merely taking the average to find a constant function (<ref>c), we propose to fit a non-constant model function (<ref>e) to the samples via least-squares regression and utilize its analytical integral to define an estimator (<ref>f).

We show that our approach can be formulated as the method of control variates, where control variates are obtained via least-squares regression.

While the use of regression in control variates is not new, we show, for the first time how regression, control variates, and MC integration are closely coupled.

MC integration is a special case under our formulation where a constant function is the solution of regression.
Our approach is provably shown to be better than or equivalent to MC integration; in the worst case, it will revert back to conventional MC integration as long as a constant function is included in regression.
This theoretical guarantee has not been explored in prior work on regression-based control variates.

To summarize our contributions:

* We introduce a novel formulation of MC integration which connects least-squares regression, control variates, and conventional MC integration.

* We prove that our regression-based formulation leads to a class of estimators that are better than or equal to conventional MC integration in terms of variance.

* We develop a practical algorithm based on our formulation and polynomials, which can work as drop-in improvement over conventional MC integration.

* We demonstrate the effectiveness of our algorithms over light transport integral estimations by performing experiments over multiple dimensions (2D to 15D).

§ RELATED WORK

§.§ Numerical integration

Integration is a classical subject in numerical analysis <cit.>.
The calculus-oriented approach has lead to deterministic quadrature techniques such as a Riemann sum.

The quadrature techniques sample the integrand at predetermined points and approximate the given integrand using an analytically integrable function (e.g., quadrangles, trapezoids, and certain basis functions).

The calculus-oriented approach, in general, suffers from the exponential growth of the number of samples for higher dimensional functions.
A statistics-oriented approach, MC integration, performs a similar summation but instead evaluates the integrand at stochastically sampled locations.
This stochastic nature of sample generation extends well to higher dimensions <cit.> without the exponential growth of the number of samples, which has led to its wide adoption in the rendering community.
We take a step back from this trend and show how a calculus-oriented reformulation of MC integration can lead to more efficient numerical integration methods.

§.§ Control variates

Control variates (CVs) <cit.> are variance reduction techniques which utilize MC integration of the difference between a given integrand $\Integrand$ and another function $\CVfunction$ to perform numerical integration.
This difference estimate is added to the analytical integral of the function $\CVfunction$ to estimate the integral of $\Integrand$.
The improvement over conventional MC integration is dependent on how well $\CVfunction$ approximates $\Integrand$, and the improvement is not guaranteed in general.

<cit.> discussed a practical approach to perform importance sampling with control variates for variance reduction. They showed that adaptive control variates using a mixture density cannot be much worse (can only be slightly worse) than conventional MC using the best mixture component. Our approach, on the other hand, is provably never worse than conventional MC using the same PDF (shown in <Ref> and <Ref>). <cit.> used the method of Owen and Zhou in rendering with an additional approximation.

Multi-level Monte Carlo <cit.> shows how low resolution renderings can be used as CVs for high-resolution synthesis.
<cit.> shows a connection between CVs and the reconstruction step in gradient-domain rendering <cit.>.
Recently, <cit.> showed a connection between weights in multiple importance sampling <cit.> and CVs.
We show connections among control variates, least-squares regression, and MC integration to have a provable variance reduction.
We will discuss more related work <cit.> later once we introduced our approach.

Construction of control variates in rendering
Constructing control variates for complex integrands in light transport simulation is a challenging task.
In rendering, the simplest form of control variate uses the constant ambient term to approximate incoming illumination <cit.>.
This, however, does not work well in complex lighting conditions.

<cit.> used visibility correlations to design control variates for direct illumination. <cit.> used bayesian regression to construct a light clustering. For further improvement, a control variate function is derived from the respective estimated contribution of each cluster at any shading point.
<cit.> proposed a method of control variate extension for heterogenous participating media rendering. It is based on the decomposition of the medium into a control and a residual components. These two parts can be sampled separately.
The control variate has also been used by <cit.> to extend their method of integration of spherical harmonic expansions over polygonal domains to non-analytic integrands.
Neural control variates <cit.> fit a neural network to define control variates.
They also proposed to simultaneously train another network to importance sample the difference function.

The resulting estimator shows remarkable improvements, but with a significant computational overhead and changes to the sampling process that have to be accelerated by utilizing a high-performance multi-GPU workstation.
Recently, <cit.> propose a neural network Q-NET that uses a neural proxy as a control variate.
Our control variates are constructed from the samples within each pixel using least squares regression.

The expected improvement from these control variates is limited by construction due to the lack of complexity captured by these control variates.
It is thus reasonable to assume that control variates usually do not approximate the integrand well in light transport simulation.
Recognizing this limitation, our approach is designed to be provably better than MC integration as long as least-squares regression can find any non-constant solution.

A major deviation from prior work is that we show how regression is closely connected to control variates and MC integration. Many prior work on control variates instead focus on introducing a sophisticated model function for control variates.
This conceptual difference allows us to realize that we can use a significantly simpler model function, such as polynomials, and still provably outperforms MC integration.

We demonstrate that our estimators are already practical without any special optimization for regression (e.g., use of high performance hardware like GPUs) and without any change to the sampling process.

§.§ Regression in rendering

Regression has been applied for adaptive sampling and reconstruction in rendering <cit.>.
<cit.> proposed to locally approximate the unknown image function by locally regressing a low-order polynomial on auxiliary buffers (normals, depth, textures).
To alleviate the computational overhead of the method, <cit.> used a iterative linear model to simultaneously reconstruct multiple pixels. These methods control the filtering bandwidths locally to increase numerical accuracy.

<cit.> analyze various zero-order and first-order regression models and use auxiliary data just to fit the data and considers only the pixel color to compute the regression weights.
Our framework is largely agnostic to how regression is implemented, and we demonstrate multiple implementations as examples.
The regression part of our framework can thus utilize findings in those previous.

§ TWO INTERPRETATIONS OF MC INTEGRATION

Let us motivate our main question by recapitulating the two different interpretations of MC integration.
Suppose that we want to find a definite integral $\Integral$ of a non-negative function $\Integrand(x)$ over the domain $\Omega$ where

\begin{align}
\Integral = \int_{\Omega} \Integrand(x) \diff x.
\end{align}

In many practical problems, it is difficult to find the analytical solution to the above.
In light transport simulation, the function $\Integrand(x)$ is defined as the measurement contribution function for a light transport path $x$ <cit.>.
The definite integral in this case is the solution to light transport simulation, which is generally unavailable in an analytical form.

Rather than trying to find an analytical solution, MC integration uses a numerical estimator $\FEstimator{}$ where

\begin{align}
\label{eq:mc_estimator}
\Integral \approx \FEstimator = \frac{1}{N} \sum_{i=1}^{N} \frac{\Integrand(x_i)}{\PDF(x_i)}.
\end{align}

The variables $x_i$ are $N$ independent random samples distributed according to the probability density function $\PDF(x)$ which satisfies $\PDF(x) \neq 0$ whenever $\Integrand(x) \neq 0$ in the domain $\Domain$.

This estimator $\FEstimator{}$ is said to be unbiased in the sense that its expected value is exactly equal to the integral $\Integral$:

\begin{align}
\ExpectedOp{\FEstimator} = \frac{1}{N} \sum_{i=1}^{N} \ExpectedOp{\frac{\Integrand(x_i)}{\PDF(x_i)}} =  \frac{1}{N} \sum_{i=1}^{N} \int_{\Omega}  \frac{\Integrand(x)}{{p(x)}} {p(x)} \diff x = \Integral.
\end{align}

Change of variables
The same estimator can also be redefined in a unit hypercube $U$ of uniformly distributed random numbers $u_i$ that are used for generating samples $x_i$.
Suppose that we have a mapping $x = \Phi(u)$ where the Jacobian $|\diff x/ \diff u| = 1 / \PDF(x)$.
By change of variables, we can rewrite the integral into

\begin{align}
\Integral =  \int_{\Omega}  \Integrand(x) \diff x = \int_{U}  \Integrand(\Phi(u)) \left| \frac{\diff x}{\diff u} \right| \diff u = \int_{U} \IntegrandPSS(u) \diff u,
\end{align}

where $\IntegrandPSS(u) = \frac{\Integrand(\Phi(u))}{\PDF(\Phi(u))}$, and rewrite the estimator $\FEstimator$ into

\begin{align}
\Integral \approx \FEstimator= \frac{1}{N} \sum_{i=1}^{N} \frac{\IntegrandPSS(u_i)}{\PDF(u_i)} = \frac{1}{N} \sum_{i=1}^{N} \IntegrandPSS(u_i).
\end{align}

In other words, by properly redefining the integrand as above, one can always consider having uniformly distributed samples $u_i \sim \PDF(u) = 1$ in this unit hypercube $U$, instead of samples $x_i \sim \PDF(x)$ in the domain $\Omega$ which are generally not uniformly distributed.
This definition is known as as primary sample space <cit.> in rendering.
We use this formulation in the rest of the paper. One can generally convert Equation 1 to Equation 4 as long as change of variables from $x$ to $u$ is possible.

§.§ Statistic-oriented interpretation

The process of MC integration can be seen as an estimation of the expected value of $\IntegrandPSS(u)$ since we have

\begin{align}
\Integral = \int_{U} \IntegrandPSS(u) p(u) \diff u
= \ExpectedOp{\IntegrandPSS(u)} \approx \frac{1}{N} \sum_{i=1}^{N} \IntegrandPSS(u_i) = \FEstimator
\end{align}

where $u$ is distributed according to $\PDF(u)=1$.
Under this interpretation, one can also derive an expression of the expected squared errors of $\FEstimator$ as

\begin{align}
\ExpectedOp{\left( \FEstimator{} - \Integral \right)^2} = \VarianceOp{\FEstimator{}} =  \frac{1}{N} \VarianceOp{\IntegrandPSS(u)}
\end{align}

where we used the unbiasedness of the estimator $\ExpectedOp{\FEstimator{}} = \Integral$.
It is a well known expression that shows that the expected squared error of MC integration is proportional to the variance of $\IntegrandPSS(u)$ and inversely proportional to the number of samples $N$.

§.§ Calculus-oriented interpretation

One can also interpret MC integration as a two-step process; replacement of the integrand $\IntegrandPSS$ by a model function $\CVfunctionPSS$ and succeeding analytical integration of $\CVfunctionPSS$.

Let us consider a constant function $\CVfunctionPSS(u) = c$ as the model function.
To keep its integral $G$ equal to the original integral $F$, we would like to find a value for $c$ such that

\begin{align}
\Integral = \int_{U} \IntegrandPSS(u) \diff u = \int_{U} \CVfunctionPSS(u) \diff u = G.
\end{align}

The integral mean value theorem in calculus shows that there exist such a constant $c$ for a given definite integral.
Since the function $\CVfunctionPSS(u)$ can be analytically integrated, we have

\begin{align}
F = G = \int_{U} \CVfunctionPSS(u) \diff u = c |U|  \quad \therefore c  = \frac{1}{|U|} \int_{U} \CVfunctionPSS(u) \diff u
\end{align}

where $|U| = 1$ is the volume of the unit hypercube $U$.
Therefore, the constant $c$ can be derived as the average of $\IntegrandPSS(u)$ since

\begin{align}
c  = \frac{1}{|U|} \int_{U} \CVfunctionPSS(u) \diff u = \frac{1}{|U|} \int_{U} \IntegrandPSS(u) \diff u
\end{align}

where the last expression is the definition of the average of $\IntegrandPSS(u)$ over the unit hypercube $U$.
Note that there is no approximation introduced until this point.
Since $c$ is defined as the average of $\IntegrandPSS(u)$, we can estimate $c$ by taking the average of samples;

\begin{align}
c  \approx \frac{1}{N} \sum_{i=1}^{N} \IntegrandPSS(u_i) = \FEstimator.
\end{align}

This interpretation thus leads to the same estimator as the statistic-oriented interpretation since

\begin{align}
\Integral = \int_{U} \CVfunctionPSS(u) \diff u = \int_{U} c \diff u \approx \int_{U} \left( \frac{\FEstimator}{|U|}  \right) \diff u = \FEstimator. %= \left( \frac{\FEstimator}{|U|}  \right) |U| =
\end{align}

Under this interpretation, MC integration can be seen as the process of first estimating a function $\CVfunctionPSS(u) = c$, and then analytically integrating (the approximation of) the function $\CVfunctionPSS(u)$.

§.§ Problem statement

One might see that the calculus-oriented interpretation is  rather redundant at first, since the end result is the same. However, this interpretation actually points out the existence of $\CVfunctionPSS$ in the definition of the estimator, and also the fact that MC integration is derived by using a constant function $\CVfunctionPSS = c$.
This existence of $\CVfunctionPSS$ is not apparent in the statistics-oriented interpretation.

Our main question is whether making $\CVfunctionPSS$ non-constant is possible
and whether it is beneficial to consider such $\CVfunctionPSS$ to begin with.
We show that such a $\CVfunctionPSS$ can be defined based on control variates and least-squares regression of $\CVfunctionPSS$ to $\IntegrandPSS$.
Our formulation shows that conventional MC integration is a special case where regression happens to return a constant function for $\CVfunctionPSS$.
We also show that considering a more complex $\CVfunctionPSS$ is provably better than or equal to simply having a constant $\CVfunctionPSS$ in terms of its expected error.

§ REGRESSION-BASED MC INTEGRATION

We introduce a novel formulation of MC integration, regression-based MC integration, which allows us to consider any analytically integrable functions for $\CVfunctionPSS$, in place of the constant function.
Fig. <ref> illustrates our approach for a 1D function.

We start by considering an arbitrary model function $\CVfunctionPSS(u)$ which can be analytically integrated.
Unlike the calculus-oriented interpretation, our function $\CVfunctionPSS(u)$ can be non-constant and the integral $\CVIntegral$ does not need to be equal to the original integral $\Integral$.
We can use such $\CVfunctionPSS(u)$ by adding the difference from the original integral as

\begin{align}
\Integral = \CVIntegral + (\Integral - \CVIntegral)
\end{align}

where $(\Integral - \CVIntegral)$ is estimated by MC integration;

\begin{align}
\Integral - \CVIntegral
\approx \frac{1}{N} \sum_{i=1}^{N} \left( \IntegrandPSS(u_i) - \CVfunctionPSS(u_i) \right) = \DEstimator.
\end{align}

The final estimator is defined as a sum of two terms where

\begin{align}
\label{eq:control_variate_estimator}
\Integral \approx \CVEstimator = \CVIntegral + \DEstimator.
\end{align}

The first term $\CVIntegral$ is analytically evaluated for this given $\CVfunctionPSS(u)$ in a deterministic manner.
This estimator is not new and is known as the method of control variates.
The function $\CVfunctionPSS(u)$ is called the control variate which is typically given by a user.
In our method, we find $\CVfunctionPSS(u)$ by (least-squares) regression.

Control variates are unbiased regardless of  $\CVfunctionPSS(u)$.
Since $\CVIntegral$ is deterministic, the expected squared error of $\CVEstimator$ is given by the difference estimator

\begin{align}
\label{eq:control_variate_error}
\ExpectedOp{\left( \CVEstimator - \Integral \right)^2}  =  \VarianceOp{\DEstimator} = \frac{1}{N} \VarianceOp{\IntegrandPSS(u) - \CVfunctionPSS(u)}.
\end{align}

Control variates thus lead to a smaller expected square error compared to conventional MC integration only when $\operatorname{Var}[\IntegrandPSS(u) - \CVfunctionPSS(u)] < \operatorname{Var}[\IntegrandPSS(u)]$.
Error reduction is thus not guaranteed and it is widely recognized so <cit.>.

Our key contribution is to show that least-squares regression not only helps us automatically finding $\CVfunctionPSS(u)$ but also provably reduces the expected squared error  when a constant function is included in regression.

As we discuss later, a constant function has been usually ignored in (adaptive) control variates. Our work is the first to show that this combination of the least-square regression and inclusion of a constant function results in the provable improvement with a proper reduction to MC integration.

Fig. <ref> showcases our method for simple 1D integration problems.
The labels $O1$, $O3$, and $O5$ in our method mean different orders of polynomials used for $\hat{g}$ as we explain later.
The convergence plots show that, after a certain number of samples, our method is no worse than MC integration regardless of the order of polynomials.
It is also worth noting that regression might not well approximate the integrand (e.g., step and complex cases).
The method of control variates alone with such a bad approximation could increase the error.
Even in such cases, our method still outperforms MC integration once regression converges after some numbers of samples.
When the integrand can be modeled exactly by $\hat{g}$, as shown in the last row where the integrand is a 5th-order polynomial, our method dramatically reduces numerical error compared to MC integration.
We now explain how our estimator ${\langle F \rangle}^*$ achieves this reduction of errors by showing its connections to control variates and MC integration.

§.§ Least-squares regression of $\CVfunctionPSS(u)$

We model $\CVfunctionPSS(u, \theta)$ as a parametric function with the parameters $\theta = (\Coefficient_0, \cdots, \Coefficient_M)$.
We define the residual for $\theta$ as

\begin{align}
\label{eq:residual_cont}
\RIntegral(\theta) = \int_U \left( \IntegrandPSS(u) - \CVfunctionPSS(u, \theta) \right)^2 \diff u.
\end{align}

Let us then consider regression of $\CVfunctionPSS(u)$ to $\IntegrandPSS(u)$ such that the residual $\RIntegral(\theta)$ is minimized.
Since analytical solution to the integral in $ \RIntegral(\theta)$ is usually not available, we first consider using $N$ random samples $u_i$ to approximate $\RIntegral(\theta)$

\begin{align}
\label{eq:residual}
\RIntegral(\theta) \approx \frac{1}{N} \sum_{i=1}^{N} \left( \IntegrandPSS(u_i) - \CVfunctionPSS(u_i, \theta) \right)^2 = \REstimator(\theta),
\end{align}

which is equivalent to MC integration of the function $(\IntegrandPSS(u) - \CVfunctionPSS(u, \theta))^2$ using the estimator $\REstimator(\theta)$.
This unbiased approximation thus converges to $\RIntegral(\theta)$ as we increase the number of samples.

This expression of $\REstimator(\theta)$ can also be seen as the squared $\Lnorm{2}$ norm of the differences between $\IntegrandPSS(u_i)$ and $\CVfunctionPSS(u_i, \theta)$ on the samples $u_i$.
The solution to this regression can be obtained via least-squares regression of $\CVfunctionPSS(u, \theta)$ to $N$ pairs of $(u_i, \IntegrandPSS(u_i))$.

Regression-based Monte Carlo integration in 1D: a step function, a shifted Gaussian, a high-frequency function, and a polynomial.
Our method is no worse than MC integration. It is not required that the model function well approximates the integrand for our method to work (e.g., the step function and high-frequency function case where our model function is very different from the integrand). When the model function is a good fit (e.g., the shifted Gaussian, polynomial function case), the error can be reduced dramatically.

§.§ Connection to control variates

The original definition of control variates considers a given function $h(u)$ and a scalar parameter $\alpha$ such that

\begin{align}
\label{eq:control_variate_form}
\Integral = \int_{U} \IntegrandPSS(u)  \diff u = \alpha H + \int_U ( \IntegrandPSS(u) - \alpha h(u)) \diff u
\end{align}

where $H = \int_U h(u) \diff u$.
Given this definition, it is known that

\begin{align}
\alpha = \CovarianceOp{\IntegrandPSS,h}/\VarianceOp{h}
\end{align}

minimizes the variance of $\CVEstimator$.
This approach can be seen as a special case of our formulation where $\theta = \alpha$ and $\CVfunctionPSS(u, \theta) = \alpha h(u)$.

Let us analyze the error of $\CVEstimator$ when $\CVfunctionPSS(u, \theta)$ is given by least-squares regression.
Since the expected squared error of $\CVEstimator$ is equal to the variance of $\IntegrandPSS - \CVfunctionPSS$ (<Ref>), we have

\begin{align}
&\VarianceOp{\IntegrandPSS(u) - \CVfunctionPSS(u, \theta)}  =\ExpectedOp{\left( \IntegrandPSS(u) - \CVfunctionPSS(u, \theta) \right)^2} - \ExpectedOp{ \IntegrandPSS(u) - \CVfunctionPSS(u, \theta)}^2 \nonumber\\
&= \int_U \left( \IntegrandPSS(u) - \CVfunctionPSS(u, \theta) \right)^2 \diff u - (F - G)^2 = \RIntegral(\theta) - (F - G)^2.
\end{align}

where we used the fact that $\PDF(u) = 1$ in the primary sample space and substituted $\RIntegral(\theta)$.
The expected squared error of the estimator is thus bounded as

\begin{align}
\label{eq:upper_bound_residual}
\ExpectedOp{\left( \CVEstimator - \Integral \right)^2} = \frac{1}{N}\left(\RIntegral(\theta) - (F - G)^2\right) \leq \frac{1}{N} \RIntegral(\theta).
\end{align}

We can see that, by defining the control variate $\CVfunctionPSS(u)$ via least-squares regression, one can also minimize (the bound of) the expected squared error of the estimator.
In practice, we cannot directly minimize $\RIntegral(\theta)$, but we minimize its estimate $\REstimator(\theta)$, so this relationship does not exactly hold for regression with a finite number of samples as one can also see in Fig. <ref>.
Our experiments, however, show that minimizing $\REstimator(\theta)$ works well in practice when $N$ is not too small.

Note also that control variates directly minimize the variance by setting $\alpha$, while we minimize only its bound by least-squares regression.

While it is tempting to consider more general regression to directly minimize the variance in our method, it is the use of least-squares regression which allowed us to reveal clear connections among least-squares regression, control variates, and MC integration as we see in the followings.

§.§ Reduction to Monte Carlo integration

It is now trivial to show that our approach reduces to MC integration when regression happened to return a constant function $\CVfunctionPSS(u) = c$ as its best fit.
Substituting $\CVfunctionPSS(u) = c$ to <Ref>,

\begin{align}
\CVEstimator = c|U| + \frac{1}{N} \sum_{i=1}^{N} \left( \IntegrandPSS(u_i) - c \right) =  \frac{1}{N} \sum_{i=1}^{N} \IntegrandPSS(u_i) = \FEstimator.
\end{align}

It shows that, when $\CVfunctionPSS(u) = c$ is the solution to regression, our approach reduces to MC integration.

In this case, since we have $F = G$, the upper-bound in <Ref> becomes exact and

\begin{align}
\ExpectedOp{\left( \CVEstimator - \Integral \right)^2} = \ExpectedOp{\left( \FEstimator - \Integral \right)^2}  =  \RIntegral(c)
\end{align}

where $\RIntegral(c)$ is the residual for $\CVfunctionPSS(u) = c$.

Based on this reduction to Monte Carlo integration, we can show for the first time that the expected squared error of this approach is provably better than or equal to MC integration.
Deceivingly simple, but the important difference from prior work is that we include a constant function in least-squares regression by adding $c$ to $\CVfunctionPSS(u)$ and we minimize the residual $\RIntegral(\theta)$.
This allows us to have a "safe-guard" against performing worse than MC integration, when combined with least-squares regression.
Since $\RIntegral(\theta)$ is minimized for a model function including a constant,

\begin{align}
\label{eq:proof_mce_ours}
\RIntegral(\theta) \leq \RIntegral(c).
\end{align}

It thus follows that

\begin{gather}
\ExpectedOp{\left( \CVEstimator - \Integral \right)^2} \leq \frac{1}{N} \RIntegral(\theta) \leq \frac{1}{N} \RIntegral(c) = \ExpectedOp{\left( \FEstimator - \Integral \right)^2} \nonumber\\
\therefore \ExpectedOp{\left( \CVEstimator - \Integral \right)^2} \leq \ExpectedOp{\left( \FEstimator - \Integral \right)^2}.
\end{gather}

The two errors become equal only when $\RIntegral(\theta) = \RIntegral(c)$, that is, our estimator exactly reduces to MC integration when a constant function is the solution to regression.
This inequality also shows that, when a constant function is not the solution, it will always lead to a better or equivalent estimator $\CVEstimator$ than MC integration $\FEstimator$.
Our estimator is thus provably better than or equal to MC integration in terms of its expected squared errors.

In practice, since the solution to least-squares regression is not exact for any finite $N$, this theoretical property of reduction of errors is not guaranteed for any finite $N$.
Nonetheless, our numerical results demonstrate that reduction of error happens and stabilize well even for a practical range of $N$.

While regression is often used as an approach to find a control variate, none of the prior work shows how including a constant function in regression allows us to unify control variates and MC integration. For example, a very thorough discussion on control variates <cit.> immediately dismisses a constant control variate as redundant since it will not reduce variance at all.

Recognizing a problem of control variates that it could perform worse than MC integration, <cit.> proposed switching to conventional MC when the estimated variance of adaptive control variate is determined to be larger. Our approach does not need this explicit switching, as it is automatically done as a result of least-squares with theoretical guarantee.
Although our estimator looks similar to a general mixture estimator by <cit.>, their estimator does not revert back to MC integration.
The difference becomes apparent when we fit a mixture density including a constant function.
When regression returns a constant as the solution, the method of  also changes the sampling density to a constant.
On the other hand, our method reduces to MC integration while keeping the sampling density as the given mixture density. The two estimators are thus fundamentally different.

Some existing methods <cit.> use multiple different control variate functions to take advantage of their respective strengths. <cit.>, in particular, shows that the variance reduction depends on the optimal number of control variates used in the estimator.
Our work is orthogonal to this work.

§ POLYNOMIAL-BASED ESTIMATORS

While our formulation supports arbitrary integrable functions for $\CVfunctionPSS$, we elaborate a practical implementation with a polynomial function for $\CVfunctionPSS$ to demonstrate the strength of our method. In general, a polynomial function can be written as

\begin{align}
\label{eq:function_approximation}
\CVfunctionPSS(u) = \sum_{q=0}^{\numBasis} \Coefficient_q \Basis_q(u)\,
\end{align}

with monomials $\Basis$ defining a polynomial.
We can analytically integrate $\CVfunctionPSS$ as the weighted sum of the integrals of monomials as

\begin{align}
\label{eq:approximate_integrate}
\int_{U} \CVfunctionPSS(u) \diff u
=   \sum_{q=0}^{\numBasis} \Coefficient_q \int_{U}  \Basis_q(u) \diff u.
\end{align}

The advantage of polynomials is that it is simple and easily generalizable for different dimensions with different orders of polynomials.
Note that we do not claim that polynomials are optimal as control variates.
The integrand $\IntegrandPSS$ is, in general, far more complex than a polynomial, and it will not be well approximated by a polynomial in many cases.
Using polynomials in the conventional control variates might be seen as a bad approach because of this reason.

What is important here is that a polynomial is still more complex than a constant and includes a constant.
Our formulation thus provably guarantees improvement over MC integration in this case.
Therefore, even with a lower-order polynomial, it is expected that our method either outperforms or at least performs equally to MC integration.
We explain two different approaches for regression.

$f \leftarrow \{\}$; $u \leftarrow \{\}$; $\BasisInnerProductMatrix_0 \leftarrow 0^{M \times M}$; $b_0 \leftarrow 0^{M \times 1}$
$i \leftarrow 1$ $N$
$u_i \leftarrow$ random()
$f_i \leftarrow \IntegrandPSS(u_i)$
$(\BasisInnerProductMatrix_{i+1}, b_{i+1}) \leftarrow$ update_system($\BasisInnerProductMatrix_i$, $b_i$, $u_i$, $f_i$)

$\theta_{\mbox{min}} \leftarrow$ solve_linear_system($\BasisInnerProductMatrix_{N+1}$, $b_{N+1}$)
$G(\theta_{\mbox{min}}) + \frac{1}{N}(\sum_{i=1}^{N} f_i -  \CVfunctionPSS(u_i, \theta_{\mbox{min}}))$

Direct matrix estimator.

Equal-time comparisons between regression-based and conventional Monte Carlo integration. We demonstrate the effectiveness of our approach with polynomials of order 1 to 3. The rendering involves 3D integrations (Killeroos and Bathroom). relMSE values are scaled by $10^3$.

§.§ Direct matrix estimator

It is well known that regression with a polynomial can be formulated as a solution to a linear system.
The size of this linear system is proportional to the number of monomials $M$ and independent from the number of samples $N$, making it suitable for MC integration.

We can define a linear system by taking partial derivatives with respects to coefficients and setting all to zero

\begin{align}
\frac{\partial\REstimator(\theta)}{\partial c_q} = 0 \qquad \forall q \in \{0,1,\ldots,\numBasis\} \,.
\end{align}

Since $\REstimator(\theta)$ is the squared $\Lnorm{2}$ norm of the differences between $\IntegrandPSS(u_i)$ and $\CVfunctionPSS(u_i, \theta)$, the solution to this system is guaranteed to give a global minimum of $\REstimator(\theta)$.
In a matrix form, this linear system can be represented as

\begin{align}
\label{eq:least_squares_order_based}
%     \left\lVert
\begin{bmatrix}
\BasisInnerProduct{0}{0} &  \cdots & \BasisInnerProduct{0}{\numBasis}\\
\BasisInnerProduct{1}{0} &  \cdots & \BasisInnerProduct{1}{\numBasis}\\
\vdots \\
\BasisInnerProduct{\numBasis}{0} &  \cdots & \BasisInnerProduct{\numBasis}{\numBasis}\\
\end{bmatrix}
\begin{bmatrix}
\Coefficient_0 \\
\Coefficient_1 \\
\vdots \\
\Coefficient_\numBasis
\end{bmatrix}
\begin{bmatrix}
\Summation\IntegrandPSS(u_i)\Basis_0(u_i) \\
\Summation\IntegrandPSS(u_i)\Basis_1(u_i) \\
\vdots \\
\Summation\IntegrandPSS(u_i)\Basis_\numBasis(u_i)
\end{bmatrix}
\end{align}

where the function $\BasisInnerProduct{\cdot}{\cdot}$ represents

\begin{align}
\BasisInnerProduct{p}{q} = \Summation\Basis_p(u_i)\Basis_q(u_i)
\end{align}

which is precomputed for given bases $\Basis$.
Updating the full matrix requires $M^2 + M$ floating point multiplication and summation if we ignore the cost of evaluating the monomials.

Since we need to evaluate the resulting $\CVfunctionPSS(u)$ at each of $u_i$ after regression, we chose to store all the $N$ samples first, solve the linear system, and then evaluate $\CVfunctionPSS(u_i)$ afterward.
An iterative approach which does not require storing all the samples is possible under the matrix formulation <cit.>, but we leave it as future work and use the simpler alternative to store all the samples.

Note also that the size of the matrix stays the same regardless of the number of samples $N$.
We use complete orthogonal decomposition as a solver as it is fast for small matrices and does not require that the matrix is invertible (which may happen when $N$ is small).
Algorithm <ref> shows the details of this matrix-based approach.

§.§ Gradient descent estimator

Since our formulation does not assume a specific approach for regression, we can use an alternative approach based on gradient descent.
Rather than setting partial derivatives to zero to define a linear system, we use them as a gradient to iteratively update $\theta_k$ by starting from its initial guess $\theta_1$ and calculating

\begin{align}
\theta_{k+1} = \theta_{k} - \gamma \left(\cdots, \frac{\partial\REstimator(\theta_k)}{\partial c_{q,k}}, \cdots \right)
\end{align}

where $\gamma$ is a parameter which we set to $0.01$.
Using stochastic gradient descent, we can modify the above equation by observing that $u_i$ is a random sample used for Monte Carlo estimator $\REstimator(\theta)$ of $R(\theta)$ and use only one sample for $k=1, \cdots, N$ as

\begin{align}
\theta_{k+1} = \theta_{k} - \gamma \left(\cdots, \frac{\partial r(u_k, \theta_k) }{\partial c_{q,k}}. \cdots \right).
\end{align}

where we defined $r(u, \theta) = (\IntegrandPSS(u) - \CVfunctionPSS(u, \theta) )^2$ for simplicity.
This process can be repeated multiple times over the same set of samples to potentially reach a better solution for regression.

This gradient descent approach is quite general and can be used for more complex integrable functions than polynomials such as neural networks <cit.>.
In the case of polynomials, the computation of gradients is straightforward and only requires one evaluation of $\CVfunctionPSS(u)$ for each sample $u_i$ and the evaluations of the monomials.
We thus found that the gradient descent approach is generally less computationally costly than the matrix-based approach.
While more sophisticated algorithms may perform better, we have found that the above simple algorithm works well in our case.
Algorithm <ref> shows the details.

lr: learning rate, T: number of iterations
$f \leftarrow \{\}$; $u \leftarrow \{\}$
$i \leftarrow 1$ $N$
$u_i \leftarrow$ random()
$f_i \leftarrow \IntegrandPSS(u_i)$

$\theta_{\mbox{min}} \leftarrow 0^{M \times 1}$
$t \leftarrow 1$ $T$
$\theta_{\mbox{min}} \leftarrow$ gradient_descent($\theta_{\mbox{min}}$, lr, $u_i$, $f_i$)

$G(\theta_{\mbox{min}}) + \frac{1}{N}(\sum_{i=1}^{N} f_i - \CVfunctionPSS(u_i, \theta_{\mbox{min}}))$

Gradient descent estimator.

The gradient descent approach is general and updating the coefficient is less costly compared to solving a matrix.
The matrix approach can be used only when the residual is quadratic to the parameters, as in the case for polynomials.
In more general cases where the residual is not a quadratic function of the parameters, one can still use the gradient descent approach.
We advise to use the matrix for polynomials because it directly solves for the optimal solution.

At a low sample count, however, we found that the matrix approach can sometimes be numerically unstable as the matrix could be ill-conditioned.
We have observed this issue in particular for a higher-order polynomial in a high dimensional space, as the number of entries in the matrix will be large.
In this case, the gradient descent approach might be better suited.

Our regression is by no means perfect, and using more efficient regression techniques and implementation could lead to further performance improvement.

Convergence plot for equal-time comparisons in <ref> and <ref>. The metric is relMSE.
Our method converges faster than Monte Carlo estimation in all but has similar performance in the Bathroom scene.

§ EXPERIMENTS

Our implementation uses PBRT v3 <cit.> as a rendering system and Eigen <cit.> as a solver for linear systems.
All the results were generated on AMD EPYC 7702 64-Core or Intel Core i9-8950HK Processor using eight threads for each. For comparison, we compare our method against conventional MC integration and <cit.> using the relative mean squared error (relMSE):

\textstyle \nicefrac{1}{N} \sum_{i=1}^{N} \nicefrac{(I_i - R_i)^2}{(R_i^2 + 0.01)}, %\nonumber

where $I$ and $R$ represent the rendered and the reference image, respectively. The reference image $R$ for each scene is computed with 65,536 samples per pixel.

To handle a vector (i.e., RGB) integrand in rendering, we first estimate the RGB estimate $F_{\mathit{rgb}}$ by averaging RGB samples just as
F_{\mathit{rgb}} \approx \langle F_{\mathit{rgb}} \rangle = \nicefrac{1}{N} \sum_{i=1}^{N} \IntegrandPSS_{\mathit{rgb}}(u_i)
where $\IntegrandPSS_{\mathit{rgb}}(u_i)$ is the integrand that returns an RGB sample.
We then perform regression on the luminance value of the samples, resulting in the luminance valued $G$ as the analytical integral. The reconstruction of the RGB value based on our estimator is written as
\langle F^*_{\mathit{rgb}} \rangle = (\nicefrac{\langle F_{\mathit{rgb}} \rangle }{y(\langle F_{\mathit{rgb}} \rangle)}) \CVEstimator
where $y$ is the luminance function and $\CVEstimator$ estimate the luminance value of $F_{\mathit{rgb}}$ based on our method. In other words, we correct the noisy luminance value of $ F_{\mathit{rgb}}$ by our estimate.

§.§ Low-dimensional comparisons

We show in <Ref> and <Ref> an equal-time comparison between the conventional MC integration and our approach using polynomials of order 1 to 3 on three different scenes.

The Dining Room scene (<Ref>) is composed of a single triangle light source for direct illumination, making it a 2-dimensional integration problem. Our method shows significant variance reduction compared to MC integration as predicted.
Even in the penumbra region where the reduction of error is small, our method is no worse than MC integration.

The Killeroos scene (top row, <Ref>) is composed of a single spherical light and motion blur, which is a 3-dimensional integration problem. Despite the fact that motion blur introduces more discontinuities in the integrand—which makes it difficult to approximate by polynomials of any order—our method can still reduce the variance in most of the regions.
This result demonstrates that it is possible to achieve variance reduction with model functions that only roughly approximate the underlying integrand.

The Bathroom scene (bottom row, <Ref>) has multiple mesh light sources picked randomly using a random number, making it also a 3-dimensional integration problem. Moreover, each light in this scene is a mesh light with several triangles making discontinuities in the $2$D area light parameterization.
For this scene, the improvement in our method is limited but no worse than MC integration.

<Ref> shows the error convergence plots for the scenes described above.
As predicted by our theory, our method consistently outperforms MC integration once the solution to regression is almost converged.
At a fewer sample count than 10 spp, our method with higher order polynomials (order two and three) may perform worse than Monte Carlo integration due to the error in regression.
This is because the matrix tends to be ill-conditioned due to noise in the sampling.
Even though a linear model function (O1) is only slightly more complex than a constant model function (MC), our method with O1 significantly outperforms MC integration.
We also observe that the maximum possible reduction in variance is limited by the complexity of the model function. There is thus a trade-off between the maximum reduction of variance and its stability at a few sample count, when one chooses a model function. The order or number of basis also seems to impact the ovearhead.

Equal sample comparison of two MIS techniques (BRDF+light sampling) is performed using balance heuristics <cit.> and optimal weights <cit.>. The combination with our method uses an order-2 polynomial basis.
(a-d) Error maps and the relMSE values marked at the bottom shows improvement brought by our method.

§.§ Combination with multiple importance sampling

In <Ref>, we show that our method can be used in conjunction with MIS <cit.>. From different sampling strategies available during MIS, we simply perform the regression on weighted samples from each sampling strategy.

<Ref> shows the results using 64 samples per pixel of our method with two different importance sampling techniques: BRDF and light sampling. We compare the balance heuristic <cit.> and the optimal weights <cit.>.
Our method improves the quality of the result for both weighting schemes compared to Monte Carlo counterpart.
We also found that the combination of our method with optimal weight <cit.> keeps the strength of these two techniques, leading to the best variance reduction.
The error reduction brought by the optimal weights allows a closer regression to the integrand leading to a greater error reduction.
This result indicates that our method can take advantage of orthogonal variance reduction methods that simplify the shape of the integrand in primary sample space.

Convergence plots comparing our method using polynomial bases with a classical Monte Carlo estimator at different dimensionality. We observe that performance stays unchanged if the function complexity does not depend on the integration dimensionality (left column). Otherwise, our method experiences a slight decrease of improvement with the dimensionality but is still clearly advantageous compared to the Monte Carlo estimator (right column).

Equal sample comparison between classical Monte Carlo and our method using an 2 polynomial basis.
In this duck-cornell scene, only the upper part is directly illuminated by the light source.
The lower part has indirect-illumination. The error map (relMSE) shows clear improvements using our method. relMSE values marked at the bottom.

We study the impact of our method on the function complexity by increasing the path length: 1, 2 and 3 bounces between the camera and the light source. The dimensionality of the respective path length is 3D, 5D and 7D. At equal sample count, the improvements of our method are directly affected by the complexity of the integrand.

Equal-time comparison between our method and <cit.> where each pixel is considered independently.
Columns (a)-(c) show the control variate function for a pixel centered at each inset.
Columns (d)-(f) show the error images for each inset. Following , we ensure $30\%$ of the samples are used for building the control variate.
For both the dragon and the chopper-titan scene, our method show relatively less error ((d)-(f)).
is efficient in some cases. For dragon scene, column (e),  gives worse error than the classical MC.
The relMSE value of our method is lower than .
For better visualization, please see the html viewer in the supplemental.

[black,ultra thick] (3.1, 0.5) circle (0pt) node[anchor=north,rotate=0]
<cit.> ;
[black,ultra thick] (0.4, 3.5) circle (0pt) node[anchor=north,rotate=0]
Ours ;
[black,ultra thick] (-0.15, 0.90) circle (0pt) node[rotate=90, anchor=west]

[white,ultra thick] (3.1, 0.6) circle (0pt) node[anchor=north,rotate=0]
relMSE: $2.064\times 10^{-3}$ ;
[white,ultra thick] (1.4, 3.5) circle (0pt) node[anchor=north,rotate=0]
relMSE: $1.237\times 10^{-3}$ ;

[black,ultra thick] (-0.15, 0.90) circle (0pt) node[rotate=90, anchor=west]

[white,ultra thick] (3.5, 0.6) circle (0pt) node[anchor=north,rotate=0]
$6.174\times 10^{-3}$ ;
[white,ultra thick] (1.0, 3.5) circle (0pt) node[anchor=north,rotate=0]
$3.477\times 10^{-3}$ ;

[black,ultra thick] (-0.15, 0.90) circle (0pt) node[rotate=90, anchor=west]

[white,ultra thick] (3.5, 0.6) circle (0pt) node[anchor=north,rotate=0]
$2.412\times 10^{-3}$ ;
[white,ultra thick] (1.0, 3.5) circle (0pt) node[anchor=north,rotate=0]
$2.572\times 10^{-3}$ ;

Equal-sample comparison of our method using an $\Order{2}$ polynomial function basis with <cit.> method.

§.§ High-dimensional comparisons

Our method directly extends to higher dimensions.

We start the numerical analysis with simple analytic integrands. In <ref>, we show the convergence curves over 1D, 5D and 15D space for two integrands: (a) sum of sinusoids and (b) a multi-dimensional exponential function.
In both cases, our method shows significant improvements in variance reduction. Both 1 and 2 polynomials gives significant improvements. However, at low sample count (up to about 100 samples) in 15D, 2 polynomial performs poorly. This is because with increasing polynomial order, in higher dimensions, more samples are needed to better approximate the control variate.

We apply our method to a path tracing rendering.
<Ref> shows the results of our method for a duck-cornell scene rendered in path tracing with 3 bounces.
This corresponds to a 15D problem in our algorithm: a 3D random number for each next-event estimation and a 2D random number for the next-direction sampling at each bounce.
Improvements are visible in both directly- and indirectly-illumniated regions.
However, the improvements in the indirectly-lit regions are less due to the strong variations of the function. These variations come from the light making several bounces, which makes the control variate function less accurate than in the case of direct illumination.
One can of course consider having a more complex model function, but it will put an additional computational overhead on our method.

<Ref> shows the improvement due to our method on illumination coming from different number of bounces.
Our method provides the largest improvement on direct illumination (the top row), and the improvement generally decreases with more bounces (the middle and the bottom rows).
As we discussed above, it is not directly due to the increased number of dimensions for higher number of bounces.
We speculate that it is rather specific to the fact that the integrand in the primary sample space becomes more and more complex for higher number of bounces.
Similar observations have been made in prior work that attempt to approximate the integrand in the primary sample space <cit.>, so we believe that it is not specific to our regression approach.
We note that our method is still not performing worse than MC integration even for higher number of bounces.

§.§ Comparison to <cit.>

<cit.> proposed two adaptive control variate techniques using quadratures: one based on a quadrature per pixel and one with a global quadrature for the whole image. As we also treat each pixel independently, we have chosen to compare our method to the former method, which has been shown to be competitive.
Specifically, <cit.> construct an adaptive piecewise-polynomial control variate (per pixel) in the primary sample space by recursively splitting the space based on an error heuristic built upon nested quadrature rules and regions' volume.
We leave the possibility of extending our approach over a block of pixels as future work.

<Ref> provides an equal-time comparison of 's approach with our $\Order{2}$ polynomial basis method for a $2$D direct illumination integration in two scenes. The dragon scene consist of a simple object lit by a high-frequcency environment map. The chopper-titan scene contains glossy complex objects lit by a low-frequency environment map. For both scenes, we show the reference control variate function, the control variate functions constructed by  and our method (columns (a)-(c)). We also shows on columns (d)-(f) the error images of MC and these different technique.

For functions shown in column (a), 's method can produce very accurate function approximations, resulting in error reduction compared to MC method (d).
However, in columns (b, c), their adaptive control variate construction misses some critical function regions, resulting in a poor approximation. Moreover, stratifying the samples across regions with importance sampling increases this problem, resulting in higher error than MC (columns e,f).

In comparison, our method uses a crude approximation of control variate via a polynomial function but consistently outperforms the classical MC and is more robust compared to 's method.

<Ref> shows an equal-sample comparison between 's method and ours.
For the vw-van and the teapot scenes, our method significantly outperforms 's in difficult regions such as the bottom of the van and the reflections on the bottom side of the teapot. However, for the house scene, 's method performs better especially on the house's facades.
Despite such, it is noteworthy that the error reduction by our method is more uniform across the scenes.
's method, on the other hand, shows significant improvement upon classical MC in some regions, e.g., directly illuminated regions of the vw-van and the teapot, but performs worse than MC in some other regions.

w/o IS     w/ IS

Importance sampling (IS) performed in 's approach can increase the variance wrt classical MC due to bad IS strategy.

MC      no subdiv.     Ours $O2$ poly.

$3.33 \times 10^{-3}$     $2.69 \times 10^{-3}$     $2.057 \times 10^{-3}$

Without recursive space subdivision,  appears to be similar to our approach, but it is not. Our approach still performs better. All methods are rendered with 64 spp. relMSE values are at the bottom.

<Ref> demonstrates the impact of importance sampling in 's method.
We found that their importance sampling strategy can increase variance, which can be explained by closely investigating their adaptive subdivision scheme.
In the first row of <ref>, 's control variate function misses important integrand features. They used this approximation for importance sampling as well which may increase variance.

<Ref> shows another experiment with no adaptive subdivision in the method of  .
Without adaptive subdivision, their method falls back to a control variate estimator using a simple quadrature for a polynomial model function.
This setup looks deceivingly equivalent to ours, but they are different.

We compare our method with different basis functions at equal sample count. Please refer to the supplemental material for more scenes and bases.

method finds a polynomial via a quadrature rule with a fixed set of points within the domain for their control variates.
Consequently, their control variate cannot take advantage of any more samples than being required for the quadrature rule without splitting.
Our method, on the other hand, can take any number of samples to improve the least-squares fit without any splitting.
The least-squares regression further makes sure to reduce the variance bound as our formulation shows.
Both the error map and the relMSE values reported in <ref> confirm our observation.

In summary, due to the lack of robustness during the quadrature construction caused by their fixed sample locations and error heuristics,  's method does not consistently outperform the classical MC overall.
In the case of inaccurate subdivisions of the primary sample space, their method might result in inefficient importance sampling, leading to higher errors than MC.
By contrast, our method is provably never worse than the MC estimator despite that ours could result in rather coarse approximations of the function, limiting the error reduction in some regions.
As future work, 's adaptive model for $\CVfunctionPSS$ can potentially be useful in our formulation, but developing a practical least-squares regression algorithm for piecewise polynomials remains challenging.

§.§ Analysis

Choice of basis functions

We perform an experiment using different basis functions in the regression.
The results are shown in <ref>.
We compare the performance of order 2 polynomial, step functions, Gaussian mixture and Sine functions. All these basis have been chosen to contain the same number of parameters for a fair comparison.
The quality of the results obtained slightly varies according to the bases used and the scenes chosen.
Other results for additional bases and scenes is available in the supplemental materials.

Comparison of integration error ratio with Monte Carlo as a function of the number of lights in the scene. Our method is used here with polynomial basis of order 1, 2 and 3.
The drawing on the left schematizes the integration problem.

Regression using a direct matrix approach and stochastic gradient descent (SGD). At a low-sample count, SGD performs better. At a high-sample count, the direct matrix approach outperforms SGD. At high-order polynomial (order 7), the direct matrix approach has a very high overhead, making SGD an attractive alternative. Time is reported in seconds.

In general, while it is  challenging to define a function basis robust for all situations, our experiments show that polynomial basis works relatively well in many cases. As long as a constant function is included in least-squares regression, our method would be no worse than MC integration, no matter, what basis functions are chosen.

Impact of discontinuities on polynomial basis

Polynomial bases have difficulty in representing functions with many discontinuities.
In <ref>, we study the evolution of the integration error for a fixed sample budget of 1024 samples when integrating the contribution of multiple light sources over a diffuse plane (<ref>.a).
In this case, the number of discontinuities inside the $3$D integrand increases with the number of lights due to the random decision on the light selection (<ref>.b).
We study the impact of using 1, 2 and 3 polynomials with our method against the conventional Monte Carlo estimator.
The plot in <ref>(c) demonstrates that at a low number of lights, our method significantly improves over the conventional Monte Carlo estimation. However, this improvement decreases with the increase in the number of discontinuities. With a large number of lights, the error produced by our method becomes equivalent to that of the Monte Carlo estimator, confirming that the polynomial basis behaves equivalently to a constant model function, the conventional approach.

Computation time

We empirically found that the computation overhead of our method becomes significant for higher-order polynomial regression using the direct matrix approach.
For example, our method with low-order polynomials has negligible overhead over MC integration (<ref>) but a polynomial of order $5$ can cause $68\%$ overhead in terms of running time.
Similarly, higher-dimensional integration also requires significant overhead over MC integration.
However, we focus in demonstrating the versatility of our approach.
Improving the regression computation is out of the scope for this work.
Following above observation, we also investigate another approach for regression, such as gradient descent, which can be a better option for higher-order polynomials.
The gradient descent approach scales well thanks to its lower computational complexity ($O(M)$ vs $O(M^2)$) where $M$ is the number of monomials.

Matrix vs. Gradient descent solver

<Ref> shows the comparison between the gradient descent and the direct-matrix solver. The results are demonstrated with 32spp and 256spp on a bathroom scene taken from PBRT. We use three different polynomials of order 2, 5 and 7.
We compare the direct matrix approach based on complete orthogonal decomposition and stochastic gradient descent estimator (SGD) to solve the polynomial regression.
At a low sample count, SGD performs better than the direct matrix approach as the matrix tends to be ill-conditioned due to noise in sampling.
At a high sample count, the direct matrix approach outperforms SGD due to its ability to find a global minimum.
For a high-order polynomial (order 7),  the direct matrix approach becomes computationally expensive whereas SGD runs at at least $2\times$ faster speed while giving similar relMSE.

§ DISCUSSION

Incremental regression-based Monte Carlo integration with gradient descent. We render the Bathroom scene (256 spp) and compare with MC and non-incremental gradient descent estimator. While the incremental gradient descent has slightly higher error than its non-incremental counterpart, both achieves less error than MC integration. relMSE values are scaled by $10^{-2}$.

Incremental estimators

In our implementation, we assume that the $N$ sample estimator $\CVEstimator$ is constructed by first finding the solution to regression $\CVfunctionPSS$ using all the $N$ samples, and then evaluating $\CVEstimator$ with the same $N$ samples.
Therefore, when $N$ is changed to $N+1$, we need to solve regression again with $N+1$ samples first and then re-evaluate $\CVEstimator$ with the new $\CVfunctionPSS$ using $N+1$ samples.
In some applications, such as progressive rendering, being able to incrementally update both $\CVEstimator$ and $\CVfunctionPSS$ as we add a new sample is desirable.
Thanks to the use of control variates, our formulation allows such incremental update as follows.

Let us denote $\CVEstimator_N$ and $\CVfunctionPSS_N$ as the estimated value and the solution to regression with $N$ samples.
Using an online regression algorithm, one can obtain $\CVfunctionPSS_{N+1}$ based on $\CVfunctionPSS_N$ by adding a new sample $u_{N+1}$.
We define the incremental estimate $\CVEstimator_{N+1}$ as

\begin{align}
\resizebox{0.9\columnwidth}{!}{
$\CVEstimator_{N+1} = \frac{1}{N+1} \left(N \CVEstimator_{N} + \left( G_{N} + \IntegrandPSS(u_{N+1}) - \CVfunctionPSS_{N}(u_{N+1}) \right)\right)$
\end{align}

where $G_{N+1}$ is the integral of $\CVfunctionPSS_{N+1}$.
This estimator remains unbiased since

\begin{align}
& \ExpectedOp{\CVEstimator_{N+1}} \nonumber\\
& = \ExpectedOp{\frac{1}{N+1} \left(N \CVEstimator_N + \left( G_{N} + \IntegrandPSS(u_{N+1}) - \CVfunctionPSS_{N}(u_{N+1}) \right)\right)} \nonumber\\
&=\frac{1}{N+1} \left(N \Integral + \ExpectedOp{\left( G_{N} + \IntegrandPSS(u_{N+1}) - \CVfunctionPSS_{N}(u_{N+1}) \right)}\right) \nonumber\\
& = \frac{1}{N+1} \left(N \Integral +\Integral \right) = \Integral.
\end{align}

The expected squared error of incremental estimators is likely worse than that of non-incremental counterparts; the $N$ sample incremental estimator now combines the results of regression for all the $N-1, N-2, ..., 1$ samples which might be numerically unstable at the beginning.
The incremental estimators are, however, often computationally less expensive for the same number of samples than the non-incremental counterparts.
We leave the exact analysis of these differences as future work, but we show a preliminary result in Fig. <ref> where the incremental estimate is performed with the gradient descent estimator. The numerical error of the incremental estimate is slightly higher than the non-incremental version, but it still outperforms MC integration.

Differences from MCLS

Our regression-based Monte Carlo integration bears some similarities to Monte Carlo Least-Squares MCLS) <cit.>.

reformulated Monte Carlo integration as least-squares regression followed by the analytical integration of the resulting function.
MCLS also performs no worse than Monte Carlo integration as long as a constant function is included in regression.
There are two major differences between MCLS and our approach.
Firstly, MCLS considers a model function $g$ which integrates exactly to $F$ (i.e., $F = G$).

Such a family of model functions is far more restrictive than our model function which is allowed to integrate to any value.
Secondly,  claims that MCLS is derived under an entirely different framework than control variates, thus its connection to control variates is unclear.
In contrast, we show how control variates is closely related to the combination of Monte Carlo integration and least-squares regression.

When we use a model function $\CVfunctionPSS$ such that $\CVIntegral = \Integral$ (i.e., keeps the integral value), then our approach in fact reduces to MCLS.
MCLS is thus can be seen as a special case of our formulation.

Russian roulette
In our current implementation, we manually set all path lengths to be the same.
This is because, for regression it is necessary to fix the dimensionality of the search space.
Russian roulette, however may result in paths of different length/dimensionality.

One simple way to combine Russian roulette with our approach would be to perform regression up to a certain prescribed dimension and only invoke Russian roulette for paths that have dimensionality higher than the prescribed one.
The high-dimensional sample components can then be projected back into the lower-dimensional space.

§ LIMITATIONS AND FUTURE WORK

While our method often just works as a drop-in improvement over Monte Carlo integration in many cases, there are several limitations and future work.

Regression in our method incurs additional computation over MC integration.
Regression may or may not be worth the additional cost depending on the amount of error reduction in the end.
In general, regression becomes computationally more costly for a more complex model function, but a complex function typically leads to better error reduction.
There is a trade-off among the reduction of error, the complexity of a model function, and the computational overhead of the regression process.
Finding an optimal balance among those can potentially be achieved in an adaptive manner.
Our formulation does not specify how regression should be solved in practice, so there is a large room of design choices in its implementation.

Since our formulation is based on the primary sample space, applying our method to a higher number of bounces in light transport simulation is challenging.
The dimensionality of the primary sample space is directly proportional to the number of bounces, and regression in a high dimensional is not numerically stable.
It is being recognized that light transport simulation has a much lower dimensional structure in the original path space, as was utilized by path guiding methods <cit.>.
It would be interesting to study how our regression-based formulation would benefit from such a low dimensional structure of light transport.

Depending on how the regression is implemented, our method can have a very small bias.
When the same set of samples are used for both regression of the control variate and estimation of the difference, our estimator can be biased.
<cit.> and <cit.> have analyzed
the amount of bias introduced when estimating the $\alpha$ parameter and evaluating the difference estimator using the same set of samples in control variates.
In brief, the bias is negligible compared to the variance of the estimator and reduces at the order of $O(N^{-1})$ for $N$ samples as opposed to $O(N^{-1/2})$ of the error of Monte Carlo integration.
This makes our estimator consistent and the bias reduces faster than variance.
One can also always remove this bias by separating a set of samples into two halves, performing regression only on the other half, and taking the average of the two difference estimators, in a manner akin to cross-validation.
We, however, have not used this approach in the results in this paper since we have observed that bias is negligible.
Our analysis is not influenced by this bias since this bias happens only when the same set of samples are used in its implementation of regression, while our analysis does not assume any specific implementation of regression.

We thank the anonymous reviewers for their comments that helped shape the paper. We thank the following for scenes used in our experiments: Wing42 (Dinning-room), nacimus (Bathroom), Karl Li (PBRT-book), julioras3d (chopper-titan), Greyscalegorilla (vw-van), Benedikt Bitterli (Teapot), MtChimp2313 (House).
This project is supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC) RGPIN1507.

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

[Belcour, Xie, Hery, Meyer, Jarosz, and
NowrouzezahraiBelcour et al2018]
authorpersonLaurent Belcour, personGuofu
Xie, personChristophe Hery, personMark Meyer,
personWojciech Jarosz, and personDerek
Nowrouzezahrai. year2018.
Integrating clipped spherical harmonics
journalACM Transactions on Graphics (TOG)
volume37, number2 (year2018),

[Bitterli, Rousselle, Moon,
Iglesias-Guitián, Adler, Mitchell, Jarosz, and NovákBitterli
et al2016]
authorpersonBenedikt Bitterli, personFabrice
Rousselle, personBochang Moon, personJosé A.
Iglesias-Guitián, personDavid Adler, personKenny
Mitchell, personWojciech Jarosz, and personJan
Novák. year2016.
Nonlinearly Weighted First-Order Regression for
Denoising Monte Carlo Renderings.
volume35, number4 (dateJune
year2016), pages107–117.

[Clarberg and Akenine-MöllerClarberg and
authorpersonPetrik Clarberg and
personTomas Akenine-Möller. year2008.
Exploiting Visibility Correlation in Direct
volume27, number4 (year2008),

[Crespo, Jarabo, and MuñozCrespo
et al2021]
authorpersonMiguel Crespo, personAdrian
Jarabo, and personAdolfo Muñoz.
Primary-Space Adaptive Control Variates Using
Piecewise-Polynomial Approximations.
journalACM Trans. Graph. volume40,
number3, Article articleno25 (datejul
year2021), numpages15 pages.

authorpersonPhilip J Davis.
booktitleMethods of numerical integration.
publisherAcademic, addressNew York ; London.

[Fan, Chenney, Hu, Tsui, and LaiFan
et al2006]
authorpersonShaohua Fan, personStephen
Chenney, personBo Hu, personKam-Wah Tsui, and
personYu-Chi Lai. year2006.
Optimizing Control Variate Estimators for

[Glynn and SzechtmanGlynn and
authorpersonPeter W. Glynn and
personRoberto Szechtman. year2002.
Some New Perspectives on the Method of Control
Variates. In booktitleMonte Carlo and Quasi-Monte Carlo
Methods 2000, editorpersonKai-Tai Fang,
personHarald Niederreiter, and personFred J.
Hickernell (Eds.). publisherSpringer Berlin Heidelberg,
addressBerlin, Heidelberg, pages27–49.

[Guennebaud, Jacob, et alGuennebaud
et al2010]
authorpersonGaël Guennebaud,
personBenoît Jacob, et al
titleEigen v3.

[Guo, Bauszat, Bikker, and EisemannGuo
et al2018]
authorpersonJerry Guo, personPablo Bauszat,
personJacco Bikker, and personElmar Eisemann.
Primary sample space path guiding. In
booktitleEurographics Symposium on Rendering,
Vol. volume2018. The Eurographics Association,

authorpersonStefan Heinrich.
Multilevel Monte Carlo Methods. In
booktitleProceedings of the Third International Conference
on Large-Scale Scientific Computing-Revised Papers
(seriesLSSC '01). publisherSpringer-Verlag,
addressBerlin, Heidelberg, pages58–67.

[Hickernell, Lemieux, and OwenHickernell
et al2005]
authorpersonFred J. Hickernell,
personChristiane Lemieux, and personArt B. Owen.
Control Variates for Quasi-Monte Carlo.
journalStatist. Sci. volume20,
number1 (year2005), pages1 – 31.

[Hua, Gruson, Petitjean, Zwicker,
Nowrouzezahrai, Eisemann, and HachisukaHua et al2019]
authorpersonBinh-Son Hua, personAdrien
Gruson, personVictor Petitjean, personMatthias
Zwicker, personDerek Nowrouzezahrai, personElmar
Eisemann, and personToshiya Hachisuka.
A Survey on Gradient-Domain Rendering. In
booktitleComputer Graphics Forum,
Vol. volume38. Wiley Online Library, pages455–472.

[Kelemen, Szirmay-Kalos, Antal, and
CsonkaKelemen et al2002]
authorpersonCsaba Kelemen, personLászló
Szirmay-Kalos, personGyörgy Antal, and
personFerenc Csonka. year2002.
A Simple and Robust Mutation Strategy for the
Metropolis Light Transport Algorithm.
volume21, number3 (dateSept.
year2002), pages531–540.

authorpersonAlexander Keller.
Hierarchical Monte Carlo Image Synthesis.
journalMath. Comput. Simul. volume55,
number1–3 (dateFeb. year2001),

[Kondapaneni, Vévoda, Grittmann, Skrivan,
Slusallek, and KrivanekKondapaneni et al2019]
authorpersonIvo Kondapaneni, personPetr
Vévoda, personPascal Grittmann, personTomas
Skrivan, personPhilipp Slusallek, and
personJaroslav Krivanek. year2019.
Optimal Multiple Importance Sampling.
journalACM Transactions on Graphics (Proceedings of
SIGGRAPH 2019) volume38, number4
(dateJuly year2019), pages37:1–37:14.

[Kutz, Habel, Li, and NovákKutz
et al2017]
authorpersonPeter Kutz, personRalf Habel,
personYining Karl Li, and personJan Novák.
Spectral and decomposition tracking for rendering
heterogeneous volumes.
journalACM Transactions on Graphics (TOG)
volume36, number4 (year2017),

[Lafortune and WillemsLafortune and
authorpersonEric P. Lafortune and
personYves D. Willems. year1994.
The Ambient Term as a Variance Reducing Technique
for Monte Carlo Ray Tracing. pages163–171.

authorpersonW. W. Loh.
titleOn the Method of Control Variates.
Ph.D. Thesis. schoolStanford University.

[Moon, Carr, and YoonMoon
et al2014]
authorpersonBochang Moon, personNathan
Carr, and personSung-Eui Yoon.
Adaptive Rendering Based on Weighted Local
volume33, number5 (dateSept.
year2014), pages170:1–170:14.

[Moon, Iglesias-Guitian, Yoon, and
MitchellMoon et al2015]
authorpersonBochang Moon, personJose A.
Iglesias-Guitian, personSung-Eui Yoon, and
personKenny Mitchell. year2015.
Adaptive Rendering with Linear Predictions.
volume34, number4 (dateJuly
year2015), pages121:1–121:11.

[Müller, Gross, and NovákMüller
et al2017]
authorpersonThomas Müller, personMarkus
Gross, and personJan Novák.
Practical Path Guiding for Efficient
Light-Transport Simulation.
volume36, number4 (dateJune
year2017), pages91–100.

[Müller, Rousselle, Keller, and
NovákMüller et al2020]
authorpersonThomas Müller, personFabrice
Rousselle, personAlexander Keller, and personJan
Novák. year2020.
Neural Control Variates.
journalACM Trans. Graph. volume39,
number6, Article articleno243 (dateNov.
year2020), numpages19 pages.

authorpersonYuji Nakatsukasa.
titleApproximate and integrate: Variance reduction in
Monte Carlo integration via function approximation.

authorpersonHarald Niederreiter.
booktitleRandom Number Generation and Quasi-Monte
Carlo Methods.

[Owen and ZhouOwen and Zhou2000]
authorpersonArt Owen and personYi Zhou.
Safe and effective importance sampling.
journalJ. Amer. Statist. Assoc.
volume95, number449 (year2000),

authorpersonArt B. Owen.
booktitleMonte Carlo Theory, Methods and Examples.
publisherTo be published.

[Pajot, Barthe, and PaulinPajot
et al2014]
authorpersonAnthony Pajot, personLoic
Barthe, and personMathias Paulin.
Globally Adaptive Control Variate for Robust
Numerical Integration.
journalSIAM Journal on Scientific Computing
volume36, number4 (year2014),

[Pharr, Jakob, and HumphreysPharr
et al2016]
authorpersonMatt Pharr, personWenzel Jakob,
and personGreg Humphreys. year2016.
booktitlePhysically Based Rendering: From Theory to
Implementation (edition3 ed.).
addressSan Francisco, CA, USA.

[Rousselle, Jarosz, and NovákRousselle
et al2016]
authorpersonFabrice Rousselle,
personWojciech Jarosz, and personJan Novák.
Image-Space Control Variates for Rendering.
volume35, number6 (dateNov.
year2016), pages169:1–169:12.

[Rubinstein and MarcusRubinstein and
authorpersonReuven Y. Rubinstein and
personRuth Marcus. year1985.
Efficiency of Multivariate Control Variates in
Monte Carlo Simulation.

authorpersonKartic Subr.
Q-NET: A Network for Low-dimensional Integrals of
Neural Proxies. In booktitleComputer Graphics Forum,
Vol. volume40. Wiley Online Library, pages61–71.

authorpersonLloyd N. Trefethen.
booktitleApproximation Theory and Approximation
Practice (Other Titles in Applied Mathematics).
publisherSociety for Industrial and Applied Mathematics,

authorpersonEric Veach.
titleRobust Monte Carlo Methods for Light
Transport Simulation.
Ph.D. Thesis. schoolStanford University,
addressUnited States – California.

[Veach and GuibasVeach and Guibas1995]
authorpersonEric Veach and
personLeonidas J. Guibas. year1995.
Optimally Combining Sampling Techniques for Monte
Carlo Rendering, Vol. volume29. pages419–428.

[Vévoda, Kondapaneni, and
KřivánekVévoda et al2018]
authorpersonPetr Vévoda, personIvo
Kondapaneni, and personJaroslav Křivánek.
Bayesian online regression for adaptive direct
illumination sampling.
journalACM Transactions on Graphics (TOG)
volume37, number4 (year2018),

[Vorba, Karlík, Šik, Ritschel, and KřivánekVorba et al2014]
authorpersonJiří Vorba, personOndřej Karlík, personMartin Šik, personTobias
Ritschel, and personJaroslav Křivánek.
On-Line Learning of Parametric Mixture Models for
Light Transport Simulation.
volume33, number4 (dateAug.
year2014), pages101:1–101:11.

[Zheng and ZwickerZheng and Zwicker2019]
authorpersonQuan Zheng and personMatthias
Zwicker. year2019.
Learning to importance sample in primary sample
space. In booktitleComputer Graphics Forum,
Vol. volume38. Wiley Online Library, pages169–179.

[Zwicker, Jarosz, Lehtinen, Moon, Ramamoorthi,
Rousselle, Sen, Soler, and YoonZwicker et al2015]
authorpersonMatthias Zwicker, personWojciech
Jarosz, personJaakko Lehtinen, personBochang Moon,
personRavi Ramamoorthi, personFabrice Rousselle,
personPradeep Sen, personCyril Soler, and
personSung-Eui Yoon. year2015.
Recent Advances in Adaptive Sampling and
Reconstruction for Monte Carlo Rendering.
volume34, number2 (dateMay
year2015), pages667–681.