Ronak Molazem

Eigenvalue estimation of light transport using Monte Carlo algorithms

Keywords: multiple light sources

Constructing solutions for the discrete form of rendering equation has been extensively studied, due to the fact that digital computations are restricted to finite representations, i.e. matrices. The disadvantage of such representations is their computational cost for accurate estimations. On the other hand, the behavior of an operator acting on an infinite-dimensional space may not be the same as its matrix form in some states. In this study, we develop methods to estimate dominant eigenvalues of light transport based on the continuous form of light transport. We commence this investigation by formulating light transport as an integral operator defined on an infinite-dimensional Hilbert space. Applying all the tools offered by functional analysis formalism, we are able to discover properties such as continuity, closedness, compactness, and Hilbert-Schmidt characteristic. We also demonstrate the effect of the compactness property of light transport to compromise the accuracy of the finite and infinite representation of light transport. Furthermore, we study the resolvent of light transport and discuss obstacles on the way to building resolvent-based tools for the spectral study of light transport. The resolvent study is highly connected to the Fredholm theory. In this framework, operators are kernel operators. Fredholm derived an infinite series expansion for his famous determinant. We illustrate how it is possible to transform computations (related to Fredholm determinant) for light transport, when light transport is Hilbert-Schmidt, in terms of path tracing techniques. Finally, we design a few methods that leverage path tracing techniques normally used in the production of images, in order to numerically compute approximations of the dominant eigenvalues of the operator.

Estimation des valeurs propres du transport de la lumière à l'aide d'algorithmes de Monte Carlo

Ronak Molazem

Résumé

La construction de methodes discrètes d'approximation de l'équation de rendu a été largement étudiée. Celles-ci impliquent des calculs numériques utilisant des représentations de dimension finie (des matrices), dont la taille suit le besoin en precision du resultat. La convergence de ces methodes est intrinsequement liee a la compacite de l'operateur de transport lumineux. Dans cette these, nous développons des méthodes pour estimer les valeurs propres dominantes de l'operateur de transport de la lumière, basées sur la forme continue de l'operateur. Nous formulons le transport de la lumière comme un opérateur intégral défini sur un espace de Hilbert de dimension infinie. Grace aux outils offerts par l'analyse fonctionnelle, nous analysons tout d'abord des propriétés telles que la continuité, la fermeture, la compacité et la caractéristique de Hilbert-Schmidt de l'operateur. Nous démontrons également l'impact de la compacité du transport de la lumière sur la précision des approximations finies de l'equation du rendu. De plus, nous étudions le résolvent de l'operateur de transport et discutons des possibilites qu'il offre pour l'étude spectrale du transport de la lumière. L'étude du résolvent est fortement liée à la théorie de Fredholm, qui propose une extension en série infinie pour son célèbre déterminant. Nous illustrons comment il est possible de transformer les calculs (liés au déterminant de Fredholm) pour le transport de la lumière, lorsque l'operateur de transport de lumière est Hilbert-Schmidt, sous forme d'integrales d'ensembles de chemin lumineux. Enfin, nous proposons plusieurs méthodes qui tirent parti des techniques de Path Tracing, normalement utilisées dans la production d'images, afin de calculer numériquement des approximations des valeurs propres dominantes de l'opérateur de transport.

8.2 Validation of the method on the matrix . . . . . . . . . . . . .

List of Figures

The analytic continuation allows to extend a series expression

of a function beyond its initial convergence disk when carefully avoiding poles of the function. . . . . . . . . . . . . . . . . . . . . 6.2 Noise-free images of T 16 L and T 17 L of the "CornellBox.ifs" scene, computed using multiple initial paths. . . . . . . . . . . . . . . . 6.3 Noise-free images of T 9 L and T 10 L of the "ThreeSquaresNoSymmetry" scene, computed using multiple initial paths. . . . . . . . 6.4 Noise-free images of T 9 L and T 10 L of the "Contact" scene, computed using an initial path. . . . . . . . . . . . . . . . . . . . . . 6.5 Comparison of the method's performance on "ThreeSquaresNoSymmetry" by increasing the number of initial paths from 100 to 10000. 6.6 Performing the method on the "Contact" scene with a single initial path. .

List of Tables

5.1

Test scenes: For each scene the albedo is 0.8. . . . . . . . . . . . 5.2 Test scenes: For each scene the albedo is 0.8, except the Cornell box scene, which has albedos of 0.15, 0.295, 0.4 and 0.5. . . . . . 7.1 Computing c n coefficients and eigenvalues of the light transport operator of the "Contact" scene using both methods. . . . . . . . 7.2 Computing c n coefficients and eigenvalues of the light transport operator of the "TwoSquares" scene using both methods. . . . . . 7.3 Computing c n coefficients and eigenvalues of the light transport operator of the "Cubes0027" scene using both methods. . . . . . 7.4 Computing c n coefficients and eigenvalues of the light transport operator of the "Cubes0343" scene using both methods. . . . . . 7.5 Computing c n coefficients and eigenvalues of the light transport operator of the "TwoSquaresNoSymmetry" scene using both methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Computing c n coefficients and eigenvalues of the light transport operator of the "CornellBox.ifs.scaled" scene using both methods. 7.7 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Contact" scene using both methods. . . . . 7.8 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "TwoSquares" scene using both methods. . . 7.9 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Cubes0027" scene using both methods. . . . 7.10 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Cubes0343" scene using both methods. . . . 7.11 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "TwoSquaresNoSymmetry" scene using both methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.12 Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "CornellBox.ifs.scaled" scene using both methods. .

Global notations

S

Surfaces of the scene Ω

Half-sphere of outgoing directions (local frame) ℓ

The ray space S × Ω L 2

The space of square-integrable functions H

The space of light distributions B(H)

The space of linear bounded operators from H to itself r(•)

The spectral radius of an operator (Def. 5.1.10) < , > Albedo at x ∈ S ρ(x, ω o , ω i ) BRDF at x in directions ω o and ω i V (x, y)

Visibility function between points x and y κ(x, y)

Integral kernel of T b p : ℓ → S Function mapping (x, ω i ) to the point (and direction) seen from x in direction ω i , with p 2 (x, ω i ) = (x, ω i )

Figure 1: Notations used in this manuscript.

Chapter 1

Introduction

During the past few centuries, mathematical models simulating the physical behavior of light propagation and its interaction with surfaces are developed by physicians. With the aid of computers, it is possible to evaluate these models numerically, given the material properties and the description of the environment. The process of creating images using numerical methods is called image synthesis and has become an active research area in computer graphics since 1970. Early applications of shading models come from the pivotal effort of Phong [START_REF] Bui Tuong | Illumination for computer generated pictures[END_REF], Bouknight [START_REF] Bouknight | An algorithm for producing half-tone computer graphics presentations with shadows and movable light sources[END_REF], and Gouraud [START_REF] Gouraud | Computer display of curved surfaces[END_REF]. All these models were local i.e. illumination caused by other surfaces was not considered. It has been established that considering illumination from other surfaces, the global model achieves greater realism. This model essentially solves the equilibrium of light energy in a scene by determining the visible surfaces.

The central purpose of any global illumination algorithm is to provide a solution for the so-called rendering equation. For example, Kajiya [START_REF] James | The rendering equation[END_REF] illustrated that distribution ray tracing provides a solution to the rendering equation by the Monte Carlo algorithm, but it was Arvo [START_REF] Richard | Analytic methods for simulated light transport[END_REF] who developed the operator formalism for the rendering concept. Based on the same foundation, Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] built estimators for this equation using path tracing and the Markov chain concept, leading to constructing the Markov chain Monte Carlo estimator. For more restricted cases, considering purely diffuse materials, the radiosity method, derived by Gloral et al. [START_REF] Goral | Modeling the interaction of light between diffuse surfaces[END_REF], is a finite element method. For this special case, the type of solution is achieved by discretization, while Monte Carlo estimators build estimations without it. The cost of discretization increases dramatically for non-diffuse scenes, since not only patch elements but also directions are discretized. Therefore, building a light transport simulator without discretization requirement is beneficial from a practical point of view. All these approaches approximate the solution of the rendering equation. Evaluation of these approximations can be analyzed by studying the spectral properties of the integral operator, contained in the rendering equation, that transforms the incident light to the reflected one. This operator is called the light transport operator.

Functional analysis is a well-known tool to study operators. Applying the knowledge developed in this framework, it is possible to discover the true nature of this operator as an operator acting on an infinite-dimensional space (the space of light distributions). Observing the light transport operator from this window will help answering the questions challenging the solutions based on finite dimensional representations (such as the solution obtained by radiosity methods) associated to an operator which is not necessarily a matrix.

The eigenvalue problem of light transport

Path tracing algorithms provide the solution by generating paths between light sources and the surface points in the scene. These algorithms can be applied for any material model and geometry. In general, path tracing algorithms with acceleration structures (such as kd-trees, and octrees) build an image with highquality illumination effects faster. Despite the advent of GPUs to reduce the computation cost, the determination of the visibility between point surfaces still remains costly. Therefore, there are other families of approaches developed based on recorded data for interactive rendering, such as Precomputed Radiance Transfer (PRT) that records the transport matrix to handle complex illumination effects. PRT algorithms assume that a static scene is lit by infinitely far light sources i.e. environment maps. In the preprocessing step, the transfer functions are generated to transform the incident light into the output radiance. These transfer functions are usually represented in terms of some basis functions, such as Spherical Harmonics (SH), wavelet, spectral mesh basis, etc. Dimension. From the computational point of view, studying the concept of dimension is crucial to constructing any light transport simulator. Dimensionality might be also expressed as the number of necessary eigenvalues required to represent the solution in terms of the eigendecomposition of the operator (under the assumption that this decomposition exists). It has been shown that the local light transport operator is a low-dimensional operator.

Considering PRT methods, the number of basis elements, required to represent the solution, might be defined as the dimension of local light transport. For example, Mahajan et al. [START_REF] Mahajan | A theory of locally low dimensional light transport[END_REF] defined the dimension as the number of principal components required to approximate local light transport with a low dimensional subspace offered by PCA.

In application, studying the behavior of the scene when changing the lighting conditions is possible by the spectral study of the transport matrix, provided by discretization. Mahajan et al. [START_REF] Mahajan | A theory of locally low dimensional light transport[END_REF] developed the concept of locally low dimensional light transport by searching for the eigenvalues of the covariance matrix of the canonical cases. Considering all these approaches which attempt to study the light transport locally, there is a common requirement which is discretization. On the other hand, discretization is usually computationally expensive (even for Lambertian scenes).

Our contribution

In comparison to path tracing methods, the techniques based on discretization, such as radiosity methods, provide the solution in low rank approximations. However, the light transport operator is intrinsically an operator acting on an infinite-dimensional space, and may not be necessarily finite rank. Our main goal in this study is to develop the methods that compute eigenvalues without the need of discretization.

To achieve this goal, some necessary information about the spectrum of the operator, i.e. the spectral study, is indispensable. The spectral study of the global light transport provides the knowledge for understanding the behavior of this operator on a global scale since the eigenpair of the operator does not depend on specific parameters such as the location of light sources (or lighting condition) or the camera position. However, it may depend on the geometry of the scene and the material model.1

Spectral study of the light transport operator

Based on the continuous representation of the light transport operator, we apply functional analysis to obtain the necessary information on the spectrum of light transport, which will guide us to design our methods. In the following, we address some important facts about the light transport operator that are proved using functional analysis.

Eigenvalues of the light transport operator. When implementing the algorithms to extract eigenvalues, knowing that eigenvalues are real or complex, helps to simplify the complexity. For example, in the methods based on polynomials (Chapter 8), choosing an appropriate algorithm to extract the roots is highly dependent on such information. Therefore, similar to the proof of Baranoski et al. [START_REF] Gladimir | Eigen-Analysis for Radiosity Systems[END_REF] 2 , we prove that the light transport operator in Lambertian scenes has real eigenvalues, Theorem 4.A.3. Compactness of the light transport operator. Facing an operator acting over an infinite-dimensional space, the first question rises to the head is that to what extent does the operator behave like a matrix? From a mathematical point of view, this question is answered by checking the compactness property of the operator. Using this property, one may establish whether an estimation of the solution, provided by an algorithm, based on the discrete representation, converges to the actual solution. Considering light transport operator, we prove that the operator is not compact given an arbitrary geometry, Chapter 4. However, once the edges removed from the geometry, the operator becomes compact.

In other words, the lack of compactness explains the failure of bounding the approximation error of radiosity methods for arbitrary lighting conditions of a scene.

Resolvent study. The resolvent formalism of an operator provides a framework to apply complex analysis to obtain information about the spectrum of an operator on Banach spaces3 . The integral of the resolvent defines a projection onto the related eigenspace. We study the resolvent of the light transport operator and discuss the obstacles that prevent designing any resolvent-based algorithm to extract eigenvalues.

Methods

To design the methods for obtaining eigenvalues, we seek to avoid discretization. Therefore, it is possible to generalize the methods to non-Lambertian scenes. However, for each method, it is required to check that the necessary (and theoretical) conditions hold for that specific method, which is not an easy task. Therefore, to simplify the task, we design and validate our methods for Lambertian scenes.

Searching for eigenvalues of an operator is a well-known problem. This task becomes more challenging when the operator acts on an infinite-dimensional space, i.e. it is not a matrix. Therefore, inspired by mathematical approaches that support the generalization of the eigenvalue problem to the operators over infinite-dimensional spaces, we build algorithms that estimate the dominant eigenvalues using Monte Carlo approaches. Fredholm theory. Fredholm theory studies Fredholm integral equations. The existence and the uniqueness of the solution of these integral operators are equivalent to obtaining a non-zero values for the so-called Fredholm determinant. Depending on the characteristics of the kernel and the integration domain, several approaches are offered by mathematicians. The formalism of the Fredholm determinant is also generalized to some specific class of compact operators, such as Hilbert-Schmidt and trace class operators. Considering the light transport operator for Lambertian scenes, we are able to compute the regularized determinant, which is the extension of the Fredholm determinant to the class of Hilbert-Schmidt operators, for a specific geometry (i.e., a geometry without edges).

Regression method. Inspired by regression methods, we estimate the characteristic polynomial of the light transport operator. Because of the sensitivity of extracting roots of polynomials to the coefficients, we repeat the process and obtain the dominant eigenvalues from histograms.

Power iteration method. Similar to the matrix, under some specific assumptions on the spectrum of the operator, it is possible to reach the first dominant eigenvalue by acting the operator on the output vector repeatedly.

Thesis organization

Chapter 3 reviews the basics of the transport theory, and includes the mathematical formulation of the operators included in the rendering equations. Based on the operator formulation of light transport, introduced by Arvo [START_REF] Richard | Analytic methods for simulated light transport[END_REF], in Chapter 4, we first review the functional properties of this operator, which were also studied by Arvo [START_REF] Richard | Analytic methods for simulated light transport[END_REF] and Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF]. In addition to these primary properties, we discuss the compactness property of this operator. In Chapter 5, we discuss the resolvent of the light transport operator and demonstrate how the method is dependent on extra information about the eigenvalues and the set of eigenvectors. In Chapter 6, we perform the power iteration method by Metropolis approach to extract the dominant eigenvalue of the light transport operator. Chapter 7 represents specific coefficients that all participate to calculate the regularized determinant of the light transport operator. In Chapter 8, we build histograms of the roots of a polynomial, which is constructed from the data provided by path tracing.

Chapter 2

Related work

The low-dimensionality of the light transport operator has been widely studied (e.g. Mahajan et al. [START_REF] Mahajan | A theory of locally low dimensional light transport[END_REF] analyzed the dependency of a local light transport on the block and cluster size) and used in many global illumination algorithms. The class of Precomputed Radiance Transfer (PRT) methods is a typical example of the application of this property. However, there are only a few studies on the spectrum of the light transport operator. Baranoski et al. [START_REF] Gladimir | Eigen-Analysis for Radiosity Systems[END_REF] illustrated that the eigenvalues of a finite rank approximation of I -T are real and positive for Lambertian scenes. Ashdown et al. [START_REF] Ashdown | Eigenvector radiosity[END_REF] interprets the eigenvectors of a radiosity matrix as the eigenvectors of a conductance matrix. Building on the continuous representation of light transport (operator formalism), which was introduced by Arvo [Arv95b], and similar to the functional study of this operator by Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF], we start reviewing the very basic functional properties of the operator. Next, we go further to analyze the compactness property of this operator. In Section 2.1, we address a few references used for the functional study of the light transport operator. Section 2.2 lists some applications of low-dimensionality of the local light transport. In Section 2.3, we review a few applications based on sparseness property of the light transport operator. Existence of this property permits exploring the structure of the transport matrix using a small number of rows and columns.

Spectrum of operators on infinite-dimensional vector spaces

Spectral theory refers to a wide range of studies regarding the eigenpairs of an operator which may not be necessarily a matrix. Therefore, many studies have been trying to find situations in which the spectral properties of the operator are known. Most of these theorems are developed based on essential properties, such as boundedness (otherwise, at least to be a closed operator), compactness and self-adjoint (or at least a normal operator).

Unfortunately, the light transport operator does not have these nice properties (e.g. as a kernel operator, its kernel is not continuous or even bounded, it is neither compact nor normal). These facts about the light transport operator limit the previous work studies on operators we can exploit. However, we still benefit from the boundedness and the (local) low-dimensionality property of this operator.

Regarding the functional properties of the light transport operator, we rely upon the knowledge provided by Cheverry et al. [CR19] as a handbook of spectral theory to build possible foundations of our methods. Dealing with a non self-adjoint operator, we also refer to the book by Gohberg et al. [START_REF] Gohberg | Introduction to the theory of linear nonselfadjoint operators[END_REF] to study the Schmidt representation of non self-adjoint operators. The spectral study of an operator is connected to the resolvent of that operator, see Kato et al.'s book [START_REF] Kato | Perturbation theory for linear operators[END_REF] for the resolvent theory of a finite rank operator and Dimov et al. [START_REF] Dimov | Parallel resolvent monte carlo algorithms for linear algebra problems[END_REF] for the application of the resolvent to the eigenvalue problem. Resolvent theory is highly related to Fredholm theory1 . Zemyan et al. [START_REF] Stephen M Zemyan | The classical theory of integral equations: a concise treatment[END_REF] provided several algorithms to compute the Fredholm determinant of Fredholm integral equations over R 2 . However, Bornemann et al. [START_REF] Bornemann | On the numerical evaluation of fredholm determinants[END_REF] addressed the numerical methods of computing infinite Fredholm determinants (i.e. extension from matrices to compact operators).

Dimensionality reduction

Mahajan et al. [START_REF] Mahajan | A theory of locally low dimensional light transport[END_REF] provided the very first paper to analyze the dimensionality of light transport versus local patch/block sizes. They demonstrated that light transport is locally low-dimensional concluded from eigenvalues of the covariance matrix. Inspired by this fact, several approaches contained in the category of the PRT methods, have been developed.

Precomputed radiance transfer (PRT).

It is referred to a class of algorithms supporting complex and expensive illumination effects by precomputing the transport matrix. However, these early methods assumed that a static scene is lit by infinitely far away light sources, i.e. environment maps. Ramamoorthi et al. [START_REF] Ramamoorthi | An efficient representation for irradiance environment maps[END_REF] illustrated that the unshadowed irradiance can be represented with low-order Spherical Harmonics (SH). As in the original paper of the PRT methods, Sloan et al. [START_REF] Sloan | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF] introduced a real-time approach for shading diffuse and glossy objects lit by environment maps. Focusing on lowfrequency illumination, they used spherical harmonics to represent both incident radiance and transfer functions. They showed that a few element vectors (9-25) suffice to support lighting effects of diffuse objects. Next, Ng et al. [START_REF] Ng | All-frequency shadows using non-linear wavelet lighting approximation[END_REF] extended this precomputation idea to all-frequency environments by applying an orthonormal basis (Haar basis) over a cube-map domain.

Clustered principal component analysis (CPCA) for rendering.

PRT methods record the transfer matrix over many surface points which forms a high-dimensional signal over a set of points. Sloan et al. [START_REF] Sloan | Clustered principal components for precomputed radiance transfer[END_REF] presented a method to compress the storage by approximating this surface signal with an affine subspace, provided by principal component analysis (PCA). Liu et al. [START_REF] Liu | All-frequency precomputed radiance transfer for glossy objects[END_REF] applied the BRDF factorization to CPCA method to support allfrequency shadowing effects of glossy materials that yields a great reduction in the number of rows of transfer matrices. These approaches propose a way for dimensionality reduction over recorded transfer functions, but Nowrouzezahrai et al. [START_REF] Nowrouzezahrai | Eigentransport for efficient and accurate allfrequency relighting[END_REF] introduced a new approach, that is building a new basis directly from the samples of transfer functions. This basis consists of eigenvectors of the transfer function samples and these basis elements are applicable for allfrequency diffuse relighting issues.

Spectral mesh basis in PRT.

Wang et al. [START_REF] Wang | Precomputed radiance transfer for real-time indirect lighting using a spectral mesh basis[END_REF] presented a new method for real-time rendering of multibounce indirect lighting. The key idea of their approach is using a new basis, spectral mesh basis, to approximate the surface signal, which is defined over the mesh. Namely, these basis elements are actually the eigenvectors of a matrix, mesh Laplacian, and provided by singular value decomposition (SVD) of the matrix. To reduce the SVD cost, the mesh is divided into submeshes and the basis is exploited for each submesh.

Modular radiance transfer (MRT).

Extending PRT algorithms to local lighting [START_REF] Wang Kristensen | Precomputed local radiance transfer for real-time lighting design[END_REF] increases the amount of recorded data. Based on direct-to-indirect transfer ([WZH07] and [START_REF] Hašan | Direct-to-indirect transfer for cinematic relighting[END_REF]), Loos et al. [LAM + 11] introduced a novel PRT method, in which the scene dependency is eliminated and the basis is more appropriate for direct and indirect lighting tasks. In details, the solution of the transport operator is approximated using SVD decomposition for the direct-to-indirect map over simple shapes, called shape dictionaries. Then, the full transport solution is estimated by combining precomputed values of shape dictionaries at run-time.

Matrix reconstruction from sparse sampling

It has been observed that the transport matrix is sparse [START_REF] Ng | All-frequency shadows using non-linear wavelet lighting approximation[END_REF]. Namely, it is possible to recover the structure of this matrix using a small number of rows and columns.

Many-light problems.

Matrix row-column sampling algorithms search for a solution to the problem of shading sampled surface points lit by a number of light sources. Hašan et al.

[HPB07] offered a GPU-based method for many-lights problem. The authors developed this algorithm based on the key observation that the transport matrix is close to low rank. This fact is concluded from the decreasing order of singular values of a large matrix (containing all samples interactions with an individual light in columns and shaded samples by all lights in rows), and also based on the fact that the illumination of the scene can be represented linearly in terms of other light sources.

Nÿstrom kernel method.

Wang et al. [WDT + 09] applied Nÿstrom method to reconstruct the transport matrix from a small number of images from a static scene. This method, which works well with complex illuminations and occlusions, can be used for imagebased relighting applications. This approach is originally inspired by Nÿstrom method [START_REF] Williams | Using the nyström method to speed up kernel machines[END_REF] and used for reconstructing a low-rank matrix from a small number of rows and columns. Applying this method to light transport matrix enables full reconstruction of the matrix from the known rows and columns.

Light-Slice approach.

Ou et al. [START_REF] Ou | Lightslice: matrix slice sampling for the many-lights problem[END_REF] introduced the Light-Slice algorithm to explore both local and global structure of the transport matrix. The light-slice algorithm solves the many-lights problem by sampling matrix slices2 . This method works better than similar algorithms, such as light-cuts and mrcs, for large and complex environments.

Sparse sampling.

Huang et al. [START_REF] Fu | Sparsely precomputing the light transport matrix for real-time rendering[END_REF] showed that the cost of precomputation step can be reduced by adaptive and sparse sampling. In their method, a subset of vertices, called dense vertices, are selected for which more angular dimensions are sampled. Then, for the remaining vertices, called sparse vertices, a few angular samples are used.

Chapter 3

The physics of light transport

The aim of rendering is to produce images with a correct distribution of light in all point locations in the environment. In practice, it means to simulating the propagation of light into the environment. Dedicated algorithms require a deep insight into the physical properties of light, materials, and the interactions between these two. In particular, we are able to describe phenomena such as scattering, reflection, refraction, and absorption, which are physical descriptions of light interactions with the environment. A rendering algorithm forms an image given a full description of objects, light sources, material properties, the shape of objects, and the camera. Therefore, it is sufficient to characterize the related physical quantities. These quantities reformulate the behavior of light and the environment that leads to designing a model. In this chapter, we review key radiometric quantities that build the foundation of rendering algorithms in computer graphics.

A touch of optics

Light is electromagnetic radiation. Electromagnetic radiation can be represented as a wave, which then could be characterized by a wavelength. The wavelength of radiation varies from radio waves up to gamma-ray. There are different sub-fields of optics such as ray optics, wave optics, and photon optics. Ray optics is the most useful area of optics in computer graphics. Ray optics (or geometrical optics) describes the macroscopic properties of light. This is helpful to get insights into the basic laws, but not sufficient by itself. Wave optics (or physical optics) studies the key concepts of optics for developing technologies. Finally, it is photon optics (or quantum optics) that describes the interaction of light with objects at atomic scales.

Radiometry

In this section, we address the most important quantities, very well-known in computer graphics. Radiometry deals with measuring electromagnetic energy. Considering the energy of light, it is possible to measure the amount of energy for an arbitrary wavelength. In comparison to radiometry, photometry measures the visual sensation of the light spectrum viewed by the eye. The human eye is sensitive to a specific range of waves, which is between 330 nm and 770 nm. It is worth mentioning that different light sources with different wavelengths appear to have different colors. It was Pierre Bouguer who established photometry in 1760. He asked observers to compare different light sources. He concluded that the human eye is not capable of recognizing the amount of brightness while observing different light sources, however, it could detect whether two sources are equally bright [START_REF] Edgar | Pierre Bouguer's Optical Treatise on The Gradation of Light[END_REF]. It is now known that different spectra have different brightness to the human eye, and that is because of the differences in the responses received by our photo-receptors from different wavelengths. In computer graphics and image synthesis, radiometric quantities are used. These quantities can be computed from radiometric measurements. Radiant power or flux. The radiant power, also called flux, indicates the amount of energy flowing through a surface per unit of time. This quantity is usually denoted by ϕ and is expressed in watts (joules/sec). Radiance. Radiance is the most important quantity in image synthesis. Considering global illumination problems, this is the quantity characterizing the appearance of objects in the scene. It is defined as power per unit projected area per unit solid angle (watts/(steradian.m 2 )), denoted by L(x, ω). That is,

L(x, ω) = d 2 ϕ cos θdS x dω , (3.1)
see Figure 3.3. Note that radiance is a five-dimensional quantity varying with position, x, and direction ω.

Remark 3.2.1 Note that in equation (3.1) we used the same notation ω for the direction and the solid angle, which are not exactly the same. But, since the direction is actually contained in a differential solid angle around ω, without loss of generality, we keep using the same notation in this dissertation.

Here, we briefly recall two properties of radiance:

1. Radiance is constant along a ray, under the assumption of no loss due to absorption or scattering. One may conclude this property from the conservation law of energy.

2. The response of a sensor contained in a display device is proportional to the radiance received by the sensor.

Irradiance. Irradiance is defined as the incident radiant power per unit surface area (watts/m2 ) E = dϕ dS .

One may formulate

dϕ = [ Ωx L(x, ω) cos θdω]dS x (3.2)
where cos θdω is often called projected solid angle1 , in which θ indicates the angle between the differential solid angle and the normal at surface point x, Figure 3.3. Therefore, using equation (3.2)

E(x) = Ωx L(x, ω) cos θdω. (3.3)
Radiosity. This quantity is similar to irradiance, E(x), with the difference that radiosity refers to the amount of energy per unit area that leaves a surface. Namely,

B(x) = Ωx L(x, ω) cos θdω. (3.4)
Radiosity is also known as radiant exitance. 2 We limit our focus on the previous radiometric quantities.

Distribution functions

Reflection refers to the event in which light leaves a surface after some interaction with the surface. A distribution function is defined from S × Ω to R + that returns the amount of reflected light at a given surface point in some direction. In Section 3.5, we will briefly review the functional properties of this space.

Remark 3.2.2 In this dissertation, we consider both incident and outgoing light distributions in the upper hemisphere for notation simplification.

The Bidirectional Reflection Distribution Function (BRDF)

The BRDF characterizes the interaction of light with materials. Consider the differential irradiance incident through dω i . Then, the differential radiance reflected from a local outgoing direction

ω o is proportional to dE(x, ω i ), Figure 3.1. Namely, dL(x, ω o ) ∝ dE(x, ω i ) (3.5)
As a conclusion of equation (3.5), increasing incident light energy per unit area increases the reflected light. It is also worth mentioning that it is possible to increase the incident radiance by increasing the solid angle subtended by a source. The BRDF ρ(x, ω o , ω i ) is defined as the constant of proportionality of equation (3.5):

ρ(x, ω o , ω i ) = L(x, ω o ) L(x, ω i ) cos θdω i . (3.6)
The BRDF is a four-dimensional function, one may denote it by ρ(θ i , ϕ i ; θ o , ϕ o ), where θ i (θ o ) and ϕ i (ϕ o ) are polar angle and azimuthal angle respectively. From physical aspect, the BRDF is actually a mathematical description of light interaction with different materials. This characterization is needed to compute the appearance of materials correctly.

In the general case, due to the subsurface scattering, light can leave a surface at a point that differs from the incident point. The function dealing with this phenomenon is called the Bidirectional Surface Scattering Reflectance Distribution Function (BSSRDF) [NRH + 77]. In this dissertation, we do not consider the scattering phenomenon. We just focus on the BRDF3 .

Some properties of the BRDF.

1. Range. The BRDF is always positive.

2. Reciprocity. The value of the BRDF function does not change with reversing directions:

ρ(x, ω o , ω i ) = ρ(x, ω i , ω o ),
this property is also called Helmholtz reciprocity.

3. Anisotropy. The BRDF is anisotropic in general. Namely, the value of this function varies by rotating a given surface, without changing ω o and ω i . This phenomenon can easily be observed on brushed metals or some fabrics in which the distribution of normals at the microscale is not uniform.

4. Energy conservation. The law of conservation ensures that the total amount of reflected light in all directions is less than or equal to the total amount of light power incident on the surface. The necessary and sufficient condition to fulfill the conservation law is

∀ω i Ωx ρ(x, ω o , ω i ) cos θdω o ≤ 1, see Dutre et al.'s book [DBB18].
To design models that characterize the BRDF accurately, it is important to ensure that energy conservation and Helmholtz reciprocity properties are satisfied. More specifically, bidirectional global algorithms compute images by tracing paths from both light sources and cameras. So, using both types of paths is feasible under the assumption of reversibility of paths, i.e. Helmholtz property.

BRDF examples

Depending on the material, the appearance of objects changes. One may characterize the appearance of a surface either as diffuse or glossy. In this section, we introduce some common BRDF examples. Diffuse (Lambertian) reflection. For a diffuse surface, the BRDF value is constant in any outgoing direction. This is because of the uniform reflection over the entire hemisphere. More precisely,

ρ(x, ω o , ω i ) = ρ d π , (3.7)
where ρ d is the fraction of the total incident energy reflected from a surface , called albedo and belongs to [0, 1]. Factor π in the denominator compensates for the integral of the cosine over outgoing directions. Specular reflection. An ideal specular surface reflects the incident light in a specific direction, which makes an equal angle with the incident direction about the normal in a plane containing both the incident direction and the normal.

It is worth mentioning here that there is no such material, ideal specular, in the real world. However, some materials may behave very similar to it, like a combination of an ideal diffuse reflection with an ideal specular that becomes the glossy reflection. Glossy reflection. In practice, one may consider a BRDF model as a combination of three components: ideal specular reflection, ideal diffuse (Lambertian) reflection, and glossy reflection.

As was pointed out before, the perfect diffuse and ideal specular reflection models are just the simplest mathematical descriptions of surface behavior. Therefore, a realistic BRDF contains a component that results in more complex light interactions. This component is variously called: glossy reflection, specular, narrow, and/or directional diffuse.

In computer graphics, this component, say glossy reflection, is studied by microfacet theory. In this theory, a surface is considered as a union of little reflective facets, reflecting light perfectly. This model consists of several terms that characterize the distribution of microfacets, accounting for self-shadowing and Fresnel reflection4 . The analytical models of reflection include Phong's model (Phong et al. [START_REF] Bui Tuong | Illumination for computer generated pictures[END_REF]) and physically based models, such as Cook-Torrance (Cook et al. [START_REF] Robert | A reflectance model for computer graphics[END_REF] and He et al. [START_REF] Xiao D He | A comprehensive physical model for light reflection[END_REF]).

The rendering equation

Reflectance equation. It is known that reflectance behaves linearly, meaning that incident lights from different directions independently contribute to the resulting outgoing radiance in a given direction. Therefore, the total reflected light from a surface point can be computed as an integral over all possible incident directions:

L(x, ω o ) = Ωx ρ(x, ω o , ω i )L(x, ω i ) cos θdω i , (3.8)
this equation is called the reflectance equation.

Direct illumination

Having the incident light distribution and the BRDF, it is possible to compute the reflected distribution using the reflectance equation. The most straightforward case is considering no occlusion, and supposing that all light arriving at a surface comes from a distant light source. We may then put one more assumption: all light leaving the sources will arrive at the surface. Thus, we expect no shadow. In this case, instead of considering a whole environment, surfaces are considered individually. This model is often called local or direct illumination. Suppose that the source is a point light. It can be shown [START_REF] Michael F Cohen | Radiosity and realistic image synthesis[END_REF] that the irradiance from a single point light source x s that arrives at point x is

E(x) = Φ cos θ 4πr 2 xxs .
Now, let ω s be the direction that connects x to x s , that is ω s := x s -x. Then the radiance reaching x is

L(x, ω i ) Φ 4πr 2 xxs δ(cos θ i -cos θ s )δ(cos Φ i -cos Φ s ).
(3.9)

Therefore, substituting equation (3.9) in (3.8) gives

L(x, ω o ) = Φ 4π cos θ s ρ(x, ω o , ω s ). (3.10)
It is also possible to model direct illumination with linear and area light sources.

To read more on the subject, check [START_REF] Nishita | Continuous tone representation of three-dimensional objects taking account of shadows and interreflection[END_REF] and [START_REF] Amanatides | Ray tracing with cones[END_REF].

Indirect illumination

Indirect/Global illumination is slightly more difficult than direct illumination since in this case the entire environment is involved. To add shadows, a ray tracer is required. In a ray tracer, it is possible to test whether a point on a surface is visible from another point or not. This information is essential to form shadows. In an indirect illumination model, the radiance of a surface depends on the incident light that comes from all possible surfaces. That is the reason why the indirect illumination model is able to capture inter-reflections of the environment.

Using the property that radiance is constant along a ray, it is possible to express the incident radiance at surface point x from another point y,

L(x, ω i ) = L(y, ω o )V (x, y), (3.11)
where V (x, y) is called the visibility function, returning 1 when there is no object between x and y otherwise 0, see Figure 3.2. Considering the dependency characterized by equation 3.11, one may reformulate the equation in terms of an integral over surface points. This task can be done by a change of variable from dω i to dS y , that originates from the connection of the solid angle sustained by point y to its projected area

dω i = cos θ ′ r 2 xy dS y .
(3.12)

Finally, substituting equation (3.12) into the reflectance equation leads to

L(x, ω o ) = S ρ(x, ω o , ω i )L(y, ω o ) cos θ ′ r 2 xy V (x,

The radiosity equation

The power emitted from a diffuse emitter5 can be computed as

Φ = S Ωx L(x, ω o ) cos θdω o dS x = L S Ωx cos θdω o dS x = LπA, (3.15)
where A is the total area of this emitter. Since radiosity is defined as the exitant power per unit surface area, so

B(x) = πL(x). (3.16)
Therefore, the radiosity equation can be derived from the rendering equation.

The BRDF does not depend on the incident and exitant directions, so it can be taken out of the integral,

L(x) = ρ d (x) π S G(x, y)V (x, y)L(y)dS y + L e (x).
(3.17)

Multiplying equation (3.17) by π yields,

B(x) = ρ d (x) S G(x, y)V (x, y) B(y) π dS y + B e (x). (3.18)
Discretization of radiosity equation.

The radiosity methods belong to the class of Galerkin finite element methods. They approximate the solution of the rendering equation by a function that belongs to a finite-dimensional vector space. The simplest of these methods uses the space of piecewise constant functions over a tessellation of the scene.

In the classic radiosity equation, the goal is to compute the average radiosity, let's say B i on each patch6 s i . The self-emitted radiosity B e i and ρ i , the reflectivity of patch s i , are the inputs of the algorithm. The algorithm returns B i associated to each patch. Mathematically speaking, a system of equations formed by inputs is constructed, in which the solution is the pixel values. The average radiosity emitted from patch s i is

B i = 1 A i si Ωx L(x) cos θdω o dS x . (3.19)
Since the environment is diffuse, L(x) can be taken out of the inner integral,

B i = 1 A i si L(x)( Ωx cos θdω o )dS x = 1 A i si L(x)πdS x = 1 A i si B(x)dS x .
(3.20)

The radiosity equation is an integral equation. To solve such equations, the problem is usually reduced to an approximation problem. That is to say, the radiosity across an element (patch) is approximated by a set of basis functions, also called shape functions [START_REF] Cecil Zienkiewicz | The finite element method[END_REF].

In the following, we will derive the so-called classical radiosity system of equations starting with the radiosity integral equation.

Integrating both sides of equation (3.18) yields

1 A i si B(x)dS x = 1 A i si S ρ d (x)G(x, y)V (x, y) B(y) π dS y dS x + 1 A i si B e (x)dS x , (3.21) which is equivalent to 1 A i si B(x)dS x = 1 A i si j sj ρ d (x)G(x, y)V (x, y) B(y) π dS y dS x + 1 A i si B e (x)dS x .
(3.22) Now, under the assumption that the radiosity across each patch is constant, which is very unlikely to happen in the real world, equation (3.22) becomes

B ′ i = j B ′ j 1 A i si sj ρ d (x) G(x, y)V (x, y) π dS y dS x + B e i , (3.23)
where B ′ i is a constant value of radiosity across patch s i . Now, if we suppose that ρ d (x) is also constant over each patch, denoted by ρ i for patch i, then equation (3.23) becomes

B ′ i = ρ i j F ij B ′ j + B e i .
(3.24)

In equation (3.24), F ij coefficients are called form factors:

F ij = 1 A i si sj G(x, y)V (x, y) π dS y dS x . (3.25)
To derive an equation closer to the real case, using

B i + ϵ B (x) in equation (3.22)
turns equation (3.23) to the following

B i = ρ i j F ij B j + B e i + M, (3.26)
where M denotes the error term of the approximation.

To read more on the classic radiosity method, refer to Gloral et al. [START_REF] Goral | Modeling the interaction of light between diffuse surfaces[END_REF], Cohen et al. [START_REF] Michael | The hemi-cube: A radiosity solution for complex environments[END_REF] and Nishita et al. [START_REF] Nishita | Continuous tone representation of three-dimensional objects taking account of shadows and interreflection[END_REF].

Radiosity methods have developed more and more ever since. New classes of radiosity algorithms are based on sampling strategies. For example, one branch of such algorithm is called stochastic relaxation. These methods solve linear systems of equations using iterative algorithms.

The solution of a linear system may be achieved using Monte Carlo algorithms, see Forsythe et al. [START_REF] George | Matrix inversion by a monte carlo method[END_REF]. This new category of methods is called discrete random walk radiosity.

It is worth mentioning that there exists another class of methods using Monte Carlo approaches, in which the integral of the rendering equation is solved directly. This procedure does not rely on the assumption of piecewise constant reflectivity and works with every kind of reflectance model. One may summarize the classic radiosity method in the following steps:

1. Tessellating the scene geometry, i.e. discretizing the scene into patches.

Computing form factors using equation (3.25).

3. Solving the radiosity system of equations, equation (3.24).

Turning the resulting solution into display colors by tone mapping.

Despite the fact that radiosity is an object-based method, there are major problems such as discretization and form factor computation. In the discretization step, the size of the patches is important. It should be small enough to have an acceptable estimate of the illumination of the surfaces. On the other hand, the number of patches should be taken into consideration. The number of patches should not exceed the storage capacity. For this problem, hierarchical refinement and clustering methods are proposed. The other challenge is computation of form factors. Calculation of form factors is an expensive task since the cost grows quadratically with the number of patches.

Considering the error rising in equation (3.26), one may consider the radiosity method as an unreliable approach. In Chapter 4, we show the connection of this error bound with the compactness property of light transport. Namely, when the operator is compact, this error can be made arbitrarily small independent of the solution. However, we will see that this operator is not compact in general.

Monte Carlo Methods

Monte Carlo methods are techniques using sampling strategies to estimate integrals. These methods were first developed to simulate neutron transport by researchers like Stanislaw Ulam (1946), John von Neumann (1947), and Nicholas Metropolis (1949), while developing nuclear weapons. However, random sampling had been used earlier in an experiment conducted by Comte de Buffon (1777), to compute a specific probability. Later, it was Laplace who expressed that this repetitive experiment could be used to estimate the value of π.

The Monte Carlo approach is used to compute the value of an integration, which is provided by the expected value of a random variable. One major advantage of Monte Carlo methods is the simplicity that can be applied to a large variety of stochastic problems, such as transport problems. However, one may also mention a big disadvantage of Monte Carlo techniques as their convergence rate 1 √ N with N number of samples. This problem demands a search for optimization approaches, or variance reduction techniques. Despite all efforts until today, Monte Carlo techniques are still slow. However, for some specific problems, such as computing high-dimensional integrals or integrals with non smooth integrand, these techniques are the only practicable methods. For more details on Monte Carlo methods, refer to Kalos and Whitlock [START_REF] Malvin | Monte carlo methods[END_REF], Hammersley and Handscomb [START_REF] Hammersley | Monte carlo methods[END_REF], Spanier and Gelbard [START_REF] Spanier | Monte Carlo principles and neutron transport problems[END_REF]. For quasi-Monte Carlo methods, see Niederreiter [START_REF] Niederreiter | Random number generation and quasi-Monte Carlo methods[END_REF].

Monte Carlo integration

Let f : Ω -→ R be a real-valued function defined on a measure space. Consider the following integral,

I = Ω f (x)dµ, (3.27)
where µ is the corresponding measure.

The goal is to estimate I. To reach this goal using the Monte Carlo method, it is required to draw samples in Ω according to a probability distribution function p(x). Then, the following expression is an estimator of I using N samples

< I >:= 1 N N i=1 f (x i ) p(x i ) . (3.28)
Using linearity property of expectation, it is easy to show that

E(< I >) = I.
The variance of this estimator is

σ 2 = 1 N Ω ( f (x) p(x) -I) 2 p(x)dµ. (3.29)
Equation (3.29) indicates a dependency on N . That is, variance decreases as N increases. It actually indicates that the error of the estimation, which is proportional to variance, is decreasing by using four times more samples.

Bias

The expected value of the estimation error, denoted by B[< I >], is defined as the difference between the expectation of < I > and the actual value, i.e.

B[< I >]

:= E[< I >] -I.
This quantity is called Bias7 . If B[< I >] = 0, then the estimator is called unbiased, otherwise biased. This notation is very important in classifying Monte Carlo approaches. Using an unbiased estimator, we are more or less able to control the error. Despite the essential accuracy offered by an unbiased algorithm, these methods are expensive in practice.

Monte Carlo steps

A Monte Carlo algorithm includes the following steps:

1. Drawing samples from a probability distribution.

2. Evaluating the function at sample points.

3. Computing the estimator using resulting values of the previous step and equation (3.28).

Operator formulation

Ray space.

The four-dimensional ray space ℓ = S × Ω is the Cartesian product of the set of surface points and the unit outgoing directions, defined in the upper hemisphere centered at each surface point.

Space of light distributions.

Definition 3.5.1 (Hilbert space) A vector space which is a complete metric space with respect to the metric induced by the inner product in the vector space.

Let H be the space of linear real-valued light distribution functions defined from ℓ to R + . One may define an inner product on H as follows

< L 1 , L 2 >= S Ωx L 1 (x, ω)L 2 (x, ω) cos θdS x dω, (3.30)
where θ indicates the angle between the normal at point x and the outgoing direction ω, see Figure 3.3. Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] showed that the induced norm by this inner product is complete. We denote this norm by ∥ • ∥ H .

Measure on H.

Definition 3.5.2 (Measure space) A Measure space is a triple (S, S, µ), where S is a set, S is a σ-algebra8 of the subsets of S, and µ : S → [0, ∞] is a measure. Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] defined a measure on the ray space, called throughput measure and denoted by µ. For any given ray

(x, ω) ∈ ℓ dµ(x, ω) = dS ⊥ x dω,
where dS ⊥ x is the projected area measure, see Figure 3.3.

The local reflectance operator.

Definition 3.5.3 We define the local reflectance operator K x : L(Ω x ) -→ L(Ω x ) as the following integral operator

(K x L i )(ω o ) = Ωx ρ(x, ω o , ω i )L i (x, ω i ) cos θdω i , (3.31)
acting on an incident distribution, arriving at point x, and resulting another distribution,

K x L i = L(ω o ) at point x.
Note that this operator is defined locally, which means we are interested to interpret the interaction of light at fixed point x.

The global reflectance operator.

Definition 3.5.4 Similar to Definition 3.5.3, this operator transforms an incident distribution to another one. However, this operator, denoted by K, is not restricted to a fixed point and defined from H to itself. Therefore, these two operators are connected with the following equation

K x L(ω) = KL(x, ω).
The reparametrization operator.

Before defining the reparametrization operator, it is required to define the raycasting function. Similar to Veach's definition [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF], we define p :

ℓ -→ S such that p(x, ω i ) = x + t int ω i ,
where

t int = min{t | x + t ω i ∈ S, t ∈ (0, ∞)}.
In practice, this function indicates the closest visible point to x in a global direction corresponding to local direction ω i . Definition 3.5.5 Arvo [START_REF] Richard | Analytic methods for simulated light transport[END_REF] defined an operator G to represent an incoming radiance at a point as transported radiance reflected from another (the nearest point), i.e.

(GL)(x, ω i ) = L(y = p(x, ω i ), ω o ) = L i (x, ω i ), (3.32)
where ω i and ω o are the results of localizing W at point x and -W at point y in which W = yx, see Figure 3.4.

Note that this definition implies that the output distribution after action of G will be the same distribution, but at a different ray. That is why we call it the reparametrization operator.

The light transport operator.

Definition 3.5.6 The light transport operator T : H → H is defined as composition of K and G i.e. T = KG, that is The radiant exitance transport operator. Definition 3.5.7 Let B denotes the radiant exitance field which is defined by equation (3.4). For this specific case, the purely diffuse material, the integral of equation

∀(x, ω o ) ∈ ℓ : (T L)(x, ω o ) = Ωx ρ(x, ω o , ω i )L i (x, ω i ) cos θdω i . ( 3
(3.35) simplifies to T b : B -→ B (T b B)(x) = S κ(x, y)B(y)dS y (3.36) with κ(x, y) = ρ d (x) cos θ cos θ ′ πr 2 xy V (x, y). (3.37)
We call this operator the radiant exitance transport operator.

Part I

On the functional study of the light transport operator

Chapter 4

Functional properties of the light transport operator

In this chapter, first, we review the basic functional properties of the light transport operator in Sections 4.1, 4.2, and 4.3. For the completeness of the context, some of the known proofs and extra information on the light transport operator in a very specific situation are delivered in the appendix of the chapter. In addition to these basic properties, we also discuss the compactness property of the light transport operator, Section 4.4. We finish the chapter by discussing the impact of this property on bounding the approximation error of algorithms, developed based on the discretized transport operator (matrices), and the difficulty of computing the eigenvalues of this operator in the absence of this property.

Boundedness

Definition 4.1.1 Suppose (X, ∥.∥ X ) and (Y, ∥.∥ Y ) are normed spaces and T :

(X, ∥.∥ X ) → (Y, ∥.∥ Y ) is linear. T is said to be bounded if ∃M ∀x ∈ X ∥T (x)∥ Y ≤ M.
Notation. We denote the space of bounded linear operators defined from a Hilbert space H to itself by B(H).

Definition 4.1.2 (Operator norm) Let (X, ∥.∥ X ), (Y, ∥.∥ Y ) be two normed spaces, and T : (X, ∥.∥ X ) -→ (Y, ∥.∥ Y ) be a linear map. Then, the operator norm is defined ∥T ∥ = sup

∥L∥ X =1 ∥T L∥ Y . (4.1)
The light transport operator is bounded. Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] proved that ∥K∥ < 1. It can be shown that the reparametrization operator is also bounded, i.e. ∥G∥ ≤ 1 (see Section 4.A.1). Therefore, The light transport operator is also bounded

∥T ∥ = ∥K∥∥G∥ < 1. (4.2)

Continuity

Definition 4.2.1 (Continuous operator) A linear map T : (X, ∥.∥ X ) → (Y, ∥.∥ Y ) from one Hilbert space to another is continuous if for any x ∈ X and given ϵ > 0, there exists δ > 0 s.t.

∥x -x ′ ∥ X < δ ∥T x -T x ′ ∥ Y < ϵ.
Using the boundedness property, the following proposition [START_REF] Rudin | Functional analysis, mcgrawhill. Inc[END_REF] implies that the light transport operator is also continuous.

Proposition 4.2.1 Suppose (X, ∥.∥ X ) and (Y, ∥.∥ Y ) are normed spaces. For a linear map T : (X, ∥.∥ X ) → (Y, ∥.∥ Y ), the following statements are equivalent

1. T is bounded; 2.
T is continuous at 0;

3. T is continuous.
The light transport operator is continuous. Since the light transport operator is linear and bounded, the previous proposition implies its continuity.

Normality

Before bringing the definition of a normal operator, it is required to recall the definition of the adjoint operator.

Definition 4.3.1 (Adjoint operator) Let T ∈ B(H).

There exists a unique operator,

T * ∈ B(H) such that ∀x, y ∈ H < T x, y > H =< x, T * y > H . Definition 4.3.2 (Self-adjoint operator) The operator T ∈ B(H) is self- adjoint if T = T * .
Adjoint of G and K operators. Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] proved that operator G is self-adjoint, and K is also self-adjoint under the assumption of symmetric BSDF. In the definition of K and K x , we use the BRDF instead of the BSDF, therefore, based on Remark 3.2.2 and the Helmoltz reciprocity property of the BRDF, we conclude the same for both of these operators, Section 4.A.

T T * = T * T.
The light transport operator is not normal. By a counterexample, we show that the light transport operator is not normal.

Proof 4.3.1 Consider the following reflection model (perfect specular reflection) for a scene that only contains two surfaces with different reflectance properties (Figure 4.1):

ρ(x, ω o , ω i ) = δ(ω i -ω o )F r (ω o ) cos θ ,
where F r indicates the fraction of light reflected. For an arbitrary L, we show

T T * L ̸ = T * T L.
We know T T * = K 2 and T * T = GK 2 G. Define L ′′ := K 2 L, so

L ′′ (x, ω o ) = Ωx ρ(x, ω o , ω i )KL(x, ω i ) cos θdω i = Ωx ρ(x, ω o , ω i )( Ωx ρ(x, ω i , ω ′′ )L(x, ω ′′ ) cos θ ′′ dω ′′ ) cos θdω i = F 2 r (ω o )L(x, ω o ). (4.3)
On the other hand,

T * T L(x, ω o ) = GK 2 GL(x, ω o ) = GK 2 L(y, ω y o ) = G(F 2 r (ω y o )L(y, ω y o )) = F 2 r (ω y o )GL(y, ω y o ) = F 2 r (ω y o )L(x, ω o ).
(4.4)

Since F r has different values on the surfaces, the values of equations (4.3) and (4.4) are not the same.

Compactness

Definition 4.4.1 (Compact operator) Let (X, ∥.∥ X ), (Y, ∥.∥ Y ) be two normed spaces. A linear map T : (X, ∥.∥ X ) → (Y, ∥.∥ Y ) is called compact if for each bounded sequence (u n ) ∈ X N , (T u n ) has a convergent subsequence.
The following theorem expresses an important property of compact operators. We refer to this theorem frequently in this manuscript. Proof 4.4.1 The integral operator defined in equation 3.33 is partial, because it is only performed over directions instead of space and directions. Therefore, this operator is not compact [START_REF] As Kalitvin | On the theory of partial integral operators[END_REF].

Theorem 4.4.2 concludes that the operator is not compact in general. However, by restricting the material to purely diffuse, we can still address the condition under which the radiant exitance transport operator is compact, Theorem 4.4.4. The following theorem provides the sufficient condition of a kernel integral operator over L 2 functions to be compact.

Theorem 4.4.3 (Hilbert-Schmidt operator) Let X and Y be two measure spaces with µ and σ as their corresponding measures respectively. Consider κ(x, y) be a measurable function on X × Y with the following property

X×Y |κ(x, y)| 2 dµ(x)dσ(y) < ∞ then (T L)(x) = Y κ(x, y) L(y)dσ(y), is a compact map from L 2 (Y, dσ) to L 2 (X, dµ) [CR19]
. Such an operator is called Hilbert-Schmidt operator. Proof 4.4.2 The bounded property of the kernel implies that it is also squareintegrable:

S 2 |κ(x, y)| 2 dS x dS y ≤ M 2 S 2 dS x dS y = M 2 A 2 ,
where A indicates the area of S. Theorem 4.4.3 concludes the compactness.

In addition to Theorem 4.4.4, providing a sufficient condition for the compactness of radiant exitance operator, we also illustrate that an unbounded kernel implies non-compactness of the radiant exitance operator. Proof 4.4.3 An arbitrary geometry is a composition of joint surfaces, sharing a common edge, and disjoint surfaces. The kernel of T b becomes unbounded for arbitrary close points with different normals i.e. near a common edge. Therefore, a geometry with at least one edge is considered as a geometry with an unbounded kernel. In what follows, without loss of generality, we proceed with the proof for a simple scene with two orthogonally abutting half-planes. Let x 0 be the common origin of both planes, lied on the common edge, and S 0 be a region of the unit area on plane P 1 (See Figure 4.2). We construct a sequence of regions in this plane by scaling S 0 by a factor of 2 n , and let s n indicate the area of S n . Next, we build a sequence of radiant exitance functions, B n , as constant functions over S n regions, i.e.

B n (y) = 2 n ∀y ∈ S n Each B n has a unit L 2 norm by construction ∥B n ∥ 2 2 = Sn B 2 n (y)dS y = 4 n s n = 1.
We will show that sequence (T b B n ) does not have any converging subsequence. Choose point x in P 2 . Then

T b B n (x) = ρ 2 n F n (x),
where F n (x) and ρ denote the form factor from x to S n and the albedo of plane P 1 respectively. Since F n (x) = F n-1 (2x), we have

T b B n (x) = ρ 2 n F n (x) = ρ 2 n F n-1 (2x) = 2T b B n-1 (2x) = 2 n T b B 0 (2 n x).
We define d p as the sequence of distances between T b B n and

T b B n+p i.e. d p = ∥T b B n -T b B n+p ∥ 2 .
By definition:

d 2 p = S |T b B n (x) -T b B n+p (x)| 2 dS x = S |2 n T b B 0 (2 n x) -2 n T b B p (2 n x)| 2 dS x .
(4.5)

Let y = 2 n x. Expressing the integral of equation (4.5) in terms of y gives

d 2 p = S |T b B 0 (y) -T b B p (y)|dS y , (4.6)
where dS y = 4 n dS x was used to derive the equation.

To compute the limit of (d p ) sequence, it is required to evaluate the integrand of equation (4.6) when p tends to infinity. We know that

T b B p (y) = 2 p T b B 0 (2 p y) = ρ 2 p F 0 (2 p y).
Using a large value of p for 2 p y, makes this point far away from x 0 . Therefore, the form factor F 0 (2 p y) is equal to the multiplication of s 0 and the point-topoint form factor from 2 p y to x 0 . Orthogonality of P 1 and P 2 yields the cosine between the direction x 0 -2 p y and the normal at y tends to 0 and the cosine between 2 p yx 0 and the normal at x 0 tends to 1, therefore

lim p→∞ d p = ∥T b B 0 ∥ 2 . (4.7)
Let (T b B q(n) ) be a convergent subsequence of (T b B n ), in which q : N -→ N is a strictly increasing map. Convergence implies that this sequence is also Cauchy, i.e. for an arbitrary ϵ, there exists δ ϵ s.t.

∀n 1 , n 2 > δ ϵ ∥T b B q(n1) -T b B q ( n2) ∥ 2 < ϵ (4.8) Let n > δ ϵ then ∥T b B q(n) -T b B q(n+p) ∥ 2 = d q(n+p)-q(n) < ϵ.
However, equation (4.7) implies that the the limit of d p is a non negative value, on the other hand because q is strictly increasing, therefore d q(n+p)-q(n) can not be arbitrary small, which contradicts equation (4.8).

We provide an experiment to illustrate the result of Theorem 4.4.5 practically. In this experiment, the sequence of S n are disk emitters. The differences of the transported radiance to P 2 plane are given in a tabulation, see

Discussion

Bounding the approximation error and compactness. A compact operator can be interpreted as the limit of a sequence of square matrices, say (T n ), Theorem 4.A.1, which is obviously independent of the choice of the initial vector. Mathematically speaking, this kind of convergence is characterized in terms of the operator norm in the space of the bounded operators with the following definition.

Definition 4.5.1 (Operator-norm convergence) Let (T n ) be a sequence of bounded operators, defined from (X, ∥.∥ X ) to itself. This sequence converges to T in the operator norm, if

lim n→∞ ∥T n -T ∥ = 0. (4.9)
The operator norm convergence is stronger than the pointwise convergence defined below.

Definition 4.5.2 (Pointwise convergence)

∀L ∈ X lim n→∞ ∥(T n -T )L∥ X = 0. (4.10)
Pointwise convergence involved in approximation methods that use discretization, such as radiosity (Goral et al. [START_REF] Goral | Modeling the interaction of light between diffuse surfaces[END_REF]) and matrix row-column sampling ([HPB07]), explains the reason why such algorithms have struggled with bounding the approximation error given an arbitrary distribution. Definition 4.5.2 indicates that (T n ) converges to T , but the approximation error ∥(T n -T )L∥ X depends on L.

Designing methods in the absence of compactness. Eigen-solvers are developed based on matrix representation, which implicitly implies the compactness property. To be able to use these algorithms for extracting the eigenvalues of an operator acting over an infinite-dimensional space, the compactness property is a sufficient condition to guarantee the decreasing order of eigenvalues with finite algebraic multiplicity (Theorem 4.4.1). The absence of this property for the light transport operator (even for Lambertian scenes) makes the eigenvalue extraction task more and more challenging.

4.A Appendix

4.A.1 G is bounded.
Let the scene be closed. So, G is defined everywhere. By the inner product defined in H, Definition 3.30, we have

∥GL∥ 2 H = S Ωx [(GL)(x, ω i )] 2 dω i dS ⊥ x = S Ωy L 2 i (y, ω o )dω o dS ⊥ y = ∥L i ∥ 2 H . (4.11) Therefore, ∥G∥ = sup ∥L∥ H =1 ∥GL∥ H = 1. (4.12) 4.A.2 G is self-adjoint. < GL i , L ′ i > H = S Ωx (GL)(x, ω i )L ′ i (x, ω i )dω i dS ⊥ x = S Ωy L i (y, ω o )GL ′ i (y, ω o )dω o dS ⊥ y =< L i , GL ′ i > H , (4.13)
where For the given two incident distributions L i , L ′ i we have

L ′ i (x, ω i ) = GL ′ i (y,
< L ′ i , K x L i > H = ℓ L ′ i (x, ω o )K x L i (ω o )dω o dS ⊥ x = ℓ L ′ i (x, ω o ) Ωx ρ(x, ω o , ω i )L i (x, ω i ) cos θdω i dω o dS ⊥ x = S Ωx Ωx L ′ i (x, ω o )ρ(x, ω o , ω i )L i (x, ω i ) cos θdω i dω o dS ⊥ x = S Ωx L i (x, ω i ) cos θ Ωx L ′ i (x, ω o )ρ(x, ω o , ω i )dω o dω i dS ⊥ x =< K x L ′ i , L i > (4.14)
To obtain equation 4.14, the Fubini's theorem and the Helmoltz property were used. By definition, K x L(ω) = KL(x, ω). Therefore, K is also self-adjoint.

4.A.4 K x is compact.

The following theorem provides a necessary and sufficient condition for compactness of a linear operator [CR19]. Applying this theorem we conclude compactness of K x .

Theorem 4.A.1 Let (L k ) be a basis for a Hilbert space H, and T be a linear map from H to itself. Define

ϵ n = sup L∈span(L k ),∥L∥ H =1 ∥T L∥ H , T n = n k=0 < ., L k > T L k , (4.15) then • ϵ n = ∥T -T n ∥,
• T is compact if and only if lim n→∞ ϵ n = 0.

Theorem 4.A.2 The local reflectance operator is compact.

Proof 4.A.1 Recall the spherical harmonics, (y m p ) -p≤m≤p , 0≤p≤∞ , constructing an orthonormal basis for the space of functions defined over sphere. Therefore, for L be a radiance distribution at fixed point x

L(ω) = ∞ p=0 m=p m=-p y m p (ω)c m p , ω ∈ Ω x , (4.16)
where (c m p ) m,p are the coefficients of projecting L onto each basis element, (y m p ) m,p . One may reformulate equation (4.16) in terms of the spherical basis parametrized with one index set i ∈ {0, • • • , n 2 }, where i = p(p + 1) + m. Next, approximating the infinite series of equation (4.16) with n 2 coefficients, i.e.

L n (ω) = n 2 -1 i=0 y i (ω) c i .
(4.17)

So the previous equation is the n-th order approximation of L. Let

K n := n 2 -1 i=0 < ., y i > K x y i , since (y i ) n 2
i=0 is an orthonormal basis, we have

∥K x -K n ∥ → 0,
in terms of operator norm. Therefore, applying the previous proposition implies the compactness of K.

4.A.5 The eigenvalues of the radiant exitance operator

Considering Lambertian scenes, Baranoski et al. [START_REF] Gladimir | Eigen-Analysis for Radiosity Systems[END_REF] proved that the eigenvalues of I -T , which they called it the radiosity matrix, are real and positive. Similar to this proof, we provide the proof in terms of continuous (operator) representation of T .

Theorem 4.A.3 The eigenvalues of the radiant exitance operator are real.

Proof 4.A.2 We first prove that in Lambertian scenes, K is a positive operator. Let L i be a directional distribution of light at a point:

< KL i , L i >= ℓ cos θ ′ L i (x, ω o ) Ωx cos θρ(x, ω o , ω i )L i (x, ω i )dω i dω o dS x
When the material is Lambertian, we have ρ(x, ω o , ω i ) =1 π ρ d (x) ≥ 0. Therefore, in this case

< KL i , L i >= 1 π S ρ d (x) Ωx cos θL i (x, ω)dω 2 dS x > 0.
Let (λ, Λ) be an eigenpair of T . It follows that (λ, Λ) is an eigenpair of GK:

GKΛ = λΛ
Since K is a positive operator, there exists an operator C such that K = C * C 1 . Thus we have GC * CΛ = λΛ.

Acting C on the both sides of the previous equation leads to

(CGC * )CΛ = λCΛ
Since G is Hermitian, CGC * is Hermitian too, and λ is one of its eigenvalues.

It is therefore necessarily a real number, and so is λ.

4.A.6 Trace class property of the radiant exitance operator

Continuity of κ(x, y) at x = y.

Theorem 4.A.4 The kernel function κ is bounded on any curved/flat surface (more specifically the surface needs to be of class G 2 ). It is however continuous in only one very specific situation: flat surfaces and portions of spheres. In the general case no value can be attributed to κ(x, x).

Proof 4.A.3 In order to find the value in the limit, we consider the case of a point x on a surface that of class G 2 at this point and compute a series expansion of κ(x, y) at a nearby point y. Without loss of generality, we use 2D coordinates to represent both points using an order 2 differentiable parametrization of the surface in the neighborhood of x, along the tangent plane at x to the surface of the scene. We call h the distance between x and the projection of y on this tangent plane. Figure 4.4 shows a cross-section of the surface containing both points (their normals are in the same plane, since we consider the limit as x and y are arbitrarily close to each other). With these notations, the coordinates at x are (0, 0), and the coordinates at y are (h, α 1 h + α2 2 h 2 + α3 6 h h + ...). Because the normal is vertical at x, we have α 1 = 0. The value α 2 = 0 corresponds to a locally flat surface (or at least locally flat surfaces fall in this case). Coefficient α 2 is actually equal to the surface linear curvature at x. When α 2 ̸ = 0, the surface is locally equivalent to a spherical cap, because it has locally constant curvature. We have therefore θ = θ ′ and

cos θ cos θ ′ = ( 1 2 α 2 h 2 + O(h 3 )) 2 h 2 + ( 1 2 α 2 h 2 + O(h 3 )) 2 . So κ(x, y) = ρ(x) cos θ cos θ ′ ∥x -y∥ 2 = ρ(x) (cos θ) 2 ∥x -y∥ 2 = ρ(x) ( 1 2 α 2 h 2 + O(h 3 )) 2 h 2 + ( 1 2 α 2 h 2 + O(h 3 )) 2 2 .
Assuming that α 2 ̸ = 0, the limit of the kernel can therefore be studied when h tends to 0. We have

lim h→0 κ(x, y) = ρ(x) ( 1 2 α 2 h 2 + O(h 3 )) 2 h 2 + ( 1 2 α 2 h 2 + O(h 3 )) 2 2 = ρ(x) 1 4 α 2 2 .
Now if the surface around x is differentiable, the curvature α 2 of a cross section will depend on the direction v of that cross-section, involving the shape operator ∇n(x) of the surface at x. Writing h = [h cos β, h sin β], see Figure 4.4, we get

lim h→0 κ(x, x + h) = ρ(x) 4 h ⊺ ∇n(x)h ∥h∥ 2 = ρ(x) 4 (k 1 cos 2 β + k 2 sin 2 β) 2 ,
where k 1 and k 2 are two principal curvatures. This expression obviously depends on the direction θ at which y approaches x, which proves that even if κ is bounded, it has no limit for x = y when the two principal curvatures are not equal.

In summary, we have proved that:

• On differentiable surfaces the kernel of the operator is bounded, but generally not continuous;

• It is however continuous whenever the two principal curvatures of the surface are equal.

Specific case of the interior of a sphere. When

k 1 (x) = k 2 (x) = k(x)
, κ becomes a continuous function at y = x which makes T b a trace-class operator in this case. In particular the trace of T b can be computed by applying the theorem for trace class operators with continuous kernels:

Tr(T b ) = S κ(x, x)dS x = S ρ(x) 4 (k(x)) 2 dS x .
Inside a sphere of radius r, the two principal curvatures are constant and equal to 1/r. Applying the above formula we get

Tr(T b ) = S 1 4r 2 ρ(x)dS x = ρ d , (4.18)
where ρ d is the average albedo inside the sphere.

Chapter 5

A resolvent study of light transport

In this chapter, the main goal is to study the resolvent of the light transport operator and explain how it relates to the eigenspaces of this operator. This connection is expressed by an integral of the resolvent and the projections into the eigenspaces of the isolated eigenvalues. In Section 5.1, we start with the resolvent of the matrix, then we update the required definitions according to the case of operators over infinite-dimensional spaces. Based on this background, we conclude the derivation regarding the resolvent of the light transport operator still holds for its isolated eigenvalues. In Section 5.2, we represent the resolvent of the light transport operator as a series in terms of the powers of the operator. Speaking of power series, we also discuss the convergence areas in terms of spectral radius. Finally, Section 5.4 elaborates the difficulties and obstacles of applying the resolvent of the light transport operator to extract eigenpairs.

Tools of the trade

In this section, we review the preliminary definitions and theorems required to proceed with this chapter. For operators acting on finite-dimensional vector spaces (matrices), the books by Brown et al. Definition 5.1.3 (Isolated singular point) A singularity z 0 is called isolated if there exists at least a disk D(z 0 , ϵ) such that f is analytic on D(z 0 , ϵ) \ {z 0 }.

Theorem 5.1.1 (Cauchy's integral formula) Let f be an analytic complex function defined over a subset of the complex plane, denoted by U . Let γ be the positively oriented (oriented counterclockwise) boundary of U , then for any

a ∈ U f (a) = 1 2πi γ f (z) z -a dz.
(5.1)

The proof of the theorem uses the Cauchy integral theorem. Since f is infinitely differentiable,

f n (a) = n! 2πi f (z) (z -a) n+1 dz, (5.2)
this formula is often referred as Cauchy's differentiation formula.

Theorem 5.1.2 (Laurent series) Let f be a function which is analytic over an annular domain r < |z -z 0 | < R, and Γ be a positively oriented closed curve around z 0 . Then, f (z) has the following representation

f (z) = ∞ n=0 a n (z -z 0 ) n + ∞ n=1 b n (z -z 0 ) n , (5.3)
where

a n = 1 2πi Γ f (z) dz (z -z 0 ) n+1 ,
(5.4)

and b n = 1 2πi Γ f (z) dz (z -z 0 ) -n+1 .
(5.5) Definition 5.1.4 (Residue) Under the same assumptions of Theorem 5.1.2, where z 0 is an isolated singular point of f , the residue of function

f at z = z 0 is referred to b 1 coefficient in equation (5.3), i.e. Γ f (z) dz = 2πib 1
, we denote this coefficient with Res z=z0 f (z).

Theorem 5.1.3 (Cauchy's residue theorem) Let Γ be a positively oriented closed curve the interior of which contains a finite number of singular points

{z 1 , • • • , z k } of f , then Γ f (z) dz = 2πi i=k i=1
Res z=zi f (z).

Theorem 5.1.4 (Cauchy-Hadamard theorem) Consider the analytic expression of f , f (z) = ∞ n=0 a n (z -a) n , then the radius of convergence, r, of this series is defined as

1 lim sup n→∞ (|a n | 1 n ) .

Analytic continuation for a complex function

The aim of this section is to recall the Analytic Continuation technique. Let f (z) = ∞ n=0 a n (z -z 0 ) n be the Taylor series of f (z) around z 0 . This series is absolutely convergent in a disk D(z 0 , r 0 ). In addition to have the analytic expression of f (z) valid on this disk, it is also possible to extend this analytic form in a larger area. This is exactly the task of the analytic continuation process.

In the following of this section, we recall this process using an example then, apply this tool for a specific function, the resolvent function.

Example 5.1.1 Consider function f (z) =1 z . This function has the following expansion series

f (z) = ∞ n=0 (-1) n (z -1) n ,
which is absolutely convergent in the disc D(1, 1). The aim is to express f (z) around other point, say a, in D(-1, 1). That is

f (z) = ∞ n=0 a n (z -a) n .
Since this is supposed to be the Taylor series of f (z) around z = a, we know

a n = f n (a) n! .
Applying Cauchy's differentiation formula for a circle, centered at a with a radius less than 1 -|a -1|, we obtain

a n = (-1) n a -n-1 1 . So f (z) is represented as follows f (z) = ∞ n=0 a n (z -a) n = ∞ n=0 (-1) n a -n-1 (z -1) n = a -1 ∞ n=0 a n (1 - z a
) n , (5.6)

using Cauchy-Hadamard theorem, the radius of convergence of the series (5.6)

Figure 5.1: The analytic continuation allows to extend a series expression of a function beyond its initial convergence disk when carefully avoiding poles of the function.

is |a|, see Figure 5.1. Therefore f (z) might be represented in a larger area than D, that is D ′ (a, |a|).

Choosing the points carefully, avoiding poles of the function, this procedure can be repeated, which leads to extending the analytic form of f (z) on the complex plane, see Figure 5.1.

Resolvent of an operator over a finite-dimensional vector space

This section is dedicated to explain the concept of the resolvent of an operator defined on a finite vector space, X. The main goal is to obtain a decomposition of the resolvent in terms of infinite series, Laurent series. This series is expressed in terms of an eigenvalue of the operator, since the singularities of the resolvent are exactly the eigenvalues of the operator. At the end of this section, we bring a proposition stating a line integral of the resolvent in terms of the sum over projections of the resolvent into the related eigenspace.

Definition 5.1.5 (Resolvent of an operator) Let T : X -→ X be a linear map. The operator-valued function R T (z)

R T (z) := (T -z) -1 (5.7)
is called the resolvent of T . This operator is defined when z is not an eigenvalue.

Proposition 5.1.1 (Resolvent formula) Let z 1 , z 2 to not be eigenvalues of

T . Then (z 1 -z 2 )R T (z 1 )R T (z 2 ) = R T (z 1 ) -R T (z 2 ) (5.8)
This equality holds for the more general case, such as T be a bounded operator(or unbounded but closed) on a Banach space. For the proof of the proposition refer to [CR19].

Applying equality (5.8), one may write

R T (z) = [1 -(z -z 0 )R T (z 0 )] -1 R T (z 0 ) = ∞ n=0 (z -z 0 ) n R n+1 T (z 0 ).
(5.9)

This infinite series is absolutely convergent inside D(z 0 , r) by Theorem 5.1.4, where r = ( lim

n-→∞ ∥R n T (z 0 )∥ 1 n ) -1 .
(5.10)

Note that the value of r in equation above is also related to the spectral norm of R T , see Definition 5.1.6. This connection is concluded from the Gelfand spectral radius theorem [START_REF] Nevanlinna | Spectrum, resolvent and power boundedness[END_REF]. This theorem expresses the spectral norm of the matrix in terms of the limit of

∥T n ∥ 1 n as n → ∞: ∥T ∥ spr = lim n→∞ ∥T n ∥ 1 n .
Therefore, r = ∥R T (z 0 )∥ -1 spr .

(5.11) Definition 5.1.6 (Spectral radius of a matrix) Let T be a k × k matrix with the set of eigenvalues {λ 1 , • • • , λ k }. The spectral radius norm is defined

∥ • ∥ spr := max 1≤i≤k |λ i |.
Note that since R T (z) is analytic on D(z 0 , ∥R T (z 0 )∥ -1 spr ) with Taylor series expansion on the right-hand-side of (5.9), so

d n R T (z) dz n = n!R n+1 T (z 0 ) n = 1, 2, • • • (5.12)
Neumann series for the resolvent. To expand R T (z) in terms of T ,

R T (z) = -(z -T ) -1 = -z -1 (I -z -1 T ) = - ∞ n=0 T n z n+1 , (5.13) this series is convergent if |z| > ∥T ∥ spr .
One may reformulate equation (5.3) by

f (z) = ∞ -∞ (z -z 0 ) n A n , (5.14)
where

A n = 1 2πi Γ z -n-1 f (z)dz.
Let {λ 1 , λ 2 , • • • , λ m } be the set of eigenvalues of T . Note that the eigenvalues of T are singularities of R T (z). Therefore, it is possible to express the resolvent with the Laurent series around a specific eigenvalue. Without loss of generality, suppose λ i = 0 for some i. So, one may formulate R T (z) as follows

R T (z) = ∞ -∞ z n A n ,
where

A n = 1 2πi Γ z -n-1 R T (z)dz, (5.15)
where Γ is a positive-oriented circle that only contains z = 0. Let Γ ′ contain Γ, it can be shown [Kat13]

A n A m = α n + α m -1 2πi Γ z -n-m-2 R T (z)dz = (α n + α m -1)A n+m+1 , (5.16)
where

α n = 1 n ≥ 0,
and

α n = 0 n < 0.
To obtain equation (5.16), one may start from

A n A m = ( 1 2πi ) 2 Γ ′ Γ z -n-1 z ′-m-1 R T (z)R T (z ′ )dz dz ′ = ( 1 2πi ) 2 Γ ′ Γ z -n-1 z ′-m-1 (z ′ -z) -1 [R T (z ′ ) -R T (z)]dz dz ′ ,
(5.17)

where the RHS of the second line is obtained using resolvent formula (Equation (5.8)). Then, the following equations are also applied

1 2πi Γ z -n-1 (z ′ -z) -1 dz = α n z ′-n-1 1 2πi Γ ′ z ′-m-1 (z ′ -z) -1 dz = (1 -α m )z -m-1 .
(5.18)

Equations listed in (5.18) are derived using Cauchy's integral formula. Therefore, using equation (5.16)

A 2 -1 = -A -1 , which concludes that -A -1 is a projection 2 . Let D := -A -2 , so A -n = -D n-1 , n ≥ 2 2 A
linear map from a vector space to itself is called projection if P 2 = P . and

A n = S n+1 , n ≥ 0
, where S := A 0 . Finally, the Laurent series of R T (z) around a given eigenvalue λ i is

R T (z) = -(z -λ i ) -1 P i - ∞ n=1 (z -λ i ) -n-1 D n i + ∞ n=0 (z -λ i ) n S n+1 i (5.19)
, where

P i = Γi R T (z)dz
in which Γ i is a circle containing just λ i . Note that P i indicates -A -1 , where the subscript i indicates the index of the eigenvalue, considered to from the Laurent series around it. Using the Laurent series of the resolvent, it is possible to form a decomposition of X using P i as follows

X = M i ⊕ M ′ i , (5.20)
where

M i = P i X , M ′ i = (I -P i )X.
Applying the Laurent expansion of the resolvent around λ i and Cauchy's residue theorem, it is possible to project a vector, belonging to X, into the eigenspace associated to λ i , we denote it by Λ i . The following theorem expresses this fact.

Lemma 5.1.1 Let Γ i be a closed curve with positive direction only containing

λ i . Then Γi R T (z)dz = -2πiP i (5.21)
Proof 5.1.1 From Cauchy's residue theorem, we have

Γ R T (z)dz = 2πiRes z=λi R T (z) = 2πiRes z=λi (-(z -λ i ) -1 P i - ∞ n=1 (z -λ i ) -n-1 D n i + ∞ n=0 (z -λ i ) n S n+1 i ) = -2πi(P i ).
(5.22) Applying Cauchy's residue theorem, it is possible to derive an equation in terms of the resolvent and projections onto the eigenspaces. The following proposition illustrates this fact.

λ i < ., Λ i > Λ i .
(5.24)

Based on the fact that eigenvalues of (T -z) -1 are related to the eigenvalues of T , i.e. (λ i -z) -1 are R T (z)'s eigenvalues, and on the other hand, eigenvectors of T and R T (z) are the same, applying (5.24), one may build a decomposition for R T (z) as follows

R T (z) = i=k i=1 1 λ i -z < ., Λ i > Λ i .
(5.25)

Consider the same assumptions as Proposition 5.

1.2. Residue of R T (z) at z = λ i is Res z=λi R T (z) = lim z→λi [(z -λ i )R T (z)] = -< ., Λ i > Λ i .
(5.26)

On the other hand, we know

Γi R T (z)dz = 2πi Res z=λi R T (z) = -2πi < ., Λ i > Λ i = -2πiP i .
(5.27)

So, integrating the resolvent on Γ confirms equality in (5.23).

Resolvent of an operator over an infinite-dimensional vector space

Before starting a discussion about the resolvent of light transport, we need to clarify which definitions are generalized to the case when X is infinitedimensional vector space, such as a Hilbert space. This attention is required, since light transport is defined from H to itself, and this space is infinitedimensional.

Definition 5.1.7 (The resolvent set) The resolvent set is the set of z ∈ C such that the operator T -z : H -→ H is a bijective map3 . We denote this set with ρ(T )4 .

Definition 5.1.8 (spectrum) The spectrum of T is the set sp(T ) := C \ ρ(T ).

It is worth mentioning that since H is infinite-dimensional then, the spectrum of T does not contain only eigenvalues. In other words, each point of the spectrum is not necessarily an eigenvalue. Therefore, it is required to generalize some definitions.

Definition 5.1.9 (Eigenvalue) An eigenvalue λ is a number such that the kernel of T -λ is not empty. The set of eigenvalues of T is called point spectrum, and denoted by sp p (T ).

Definition 5.1.10 (Spectral radius of a general operator) The definition of the spectral radius of a general operator, we denote it by r(T ), is defined similar to the matrix case, however in this case, the maximum absolute value is taken over the entire spectrum, i.e. Applying the resolvent formula for light transport enables to represent its resolvent with Taylor expansion, as equation (5.9). Note that this expansion is valid for both norms ∥T ∥ and r(T ), however these two norms do not necessarily coincide, unlike the matrix case in which these two norms coincide when the matrix is normal [CR19]. In general, the spectral norm of a bounded operator is less than the operator norm using Gelfand's formula, which is also valid for an operator defined on H.

Proposition 5.1.4 (Gelfand's Formula) For a bounded operator R

r(R) = inf n∈N * ∥R n ∥ 1 n ≤ ∥R∥.
(5.30)

Resolvent of light transport

Since the light transport operator is bounded, then all the mathematical properties reviewed in Section 5.1.3 remain valid for it. Therefore,the resolvent of light transport also verifies equation (5.12), and the Taylor representation in (5.9) is valid on D(z 0 , r(R T (z 0 )) -1 ). To be able to extend the representation of the resolvent on a wider area of the complex plane, one may think of analytic continuation. To make analytic continuation for the resolvent of light transport practically feasible, it is required to use an equation containing powers of T . Practically speaking, to build any algorithm related to the resolvent of the light transport, images of powers of T are only available materials. The key to bring powers of T into the play is the Neumann series of resolvent. However, to be able to use the Neumann series of the resolvent z must be out of D(0, ∥T ∥). In Section 5.2.1 we will elaborate this idea. This is also possible to build the Laurent series of the light transport operator around an isolated eigenvalue. In general the eigenvalues of this operator may not be isolated, however, Proposition 5.1.2 is still valid for set of isolated eigenvalues of the operator [START_REF] Chatelin | Spectral approximation of linear operators[END_REF].

Analytic continuation for the resolvent of light transport

When T is acting on a finite-dimensional space then, the only primary information needed to do analytic continuation for R T , is the location of eigenvalues. However, when the operator defined on an infinite-dimensional space, it becomes even more complicated. Since, the singularities are the whole spectrum, not only eigenvalues.

To be able to use the Neumann series of the resolvent in its Taylor representation, z must be outside of D(0, ∥T ∥). Therefore, to start the analytic continuation, let z 0 ∈ C \ D(0, ∥T ∥). Equation (5.9) provides the analytic expression of R T around z 0 , and plugging the Neumann series of R T in this equation gives

R T (z) = ∞ n=0 (z -z 0 ) n R n+1 T (z 0 ) = ∞ n=0 (z -z 0 ) n [- ∞ m=0 T m z m+1 0 ] n+1 = ∞ n=0 (z -z 0 ) n c n , (5.31)
where

c n = (-1) n+1 [ ∞ m=0 T m z m+1 0 ] n+1 (5.32)
The convergence radius of the series is

r 0 = 1 r(R T (z 0 )) = 1 sup λ∈sp(T ) |λ -z 0 | -1 = 1 dist(z 0 , sp(T )) -1 = dist(z 0 , sp(T ))
(5.33)

where

dist(z 0 , sp(T )) = inf λ∈sp(T )
|z 0 -λ|.

Equation (5.33) illustrates the dependency of the convergence disk on the nearest point of the spectrum. Therefore, before choosing a new point for the next step of the analytic continuation, it is essential to know the nearest spectrum element. Unfortunately, in general, no primary information on the spectrum of light transport is available.

Dealing with operators over infinite-dimensional spaces, the compactness property usually helps to simplify problems. For this particular problem, under the assumption of compactness, computation of the convergence radius of the circle in any step of analytic continuation simplifies to searching the nearest eigenvalue. From our study, we know the state in which the exitance radiant operator is compact. So we finish this section for this special case, and under the assumption of knowing the location of the first (largest) eigenvalue, λ 1 . Under compactness property the convergence radius of the first circle is

r 0 = |z 0 -λ 1 |.
(5.34)

Then, for the next step of analytic continuation, we chose a ∈ D(z 0 , r 0 ). To compute the new coefficients of the new representation of the resolvent, we need to evaluate n-th order derivatives of the resolvent at z = a. From Taylor series representation around point z = a, we have

a n = 1 n! d n R T (a) dz n = R T (a) n+1 = [- ∞ m=0 T m a m+1 ] n+1 = (-1) n+1 [ ∞ m=0
T m a m+1 ] n+1 .

(5.35) Therefore, R T (z) = ∞ n=0 a n (z -a) n is the new representation of the resolvent in the new disk D(a, r 1 ) where

r 1 = |a -λ 1 |.

Application of the resolvent to extract eigenvalue

As it is pointed out in the previous section, equation (5.23) holds under the assumption of having isolated eigenvalues. Under this only assumption, the information it brings is the projection coefficients of L ∈ H onto the related eigenspaces, and it can be theoretically used to compute images of the eigenfunctions. However, considering some extra assumptions on T may bring us more information. Namely, if we suppose that the set of eigenfunctions of T is complete and the first eigenvalue is dominant, we are able to extract an estimate of the dominant eigenvalue. In theory, this is possible using both equations (5.24) and (5.26) together. In other words

∀L ∈ H : λ 1 ≈ -2πi T L Γ1 R T (z)Ldz .
(5.36)

To have decreasing eigenvalues in magnitude, the compactness property of the operator is sufficient (See Theorem 4.4.1). So, we assume compactness as a primary assumption on the operator. And for completeness, one may address the following conditions:

• When T is self-adjoint [Rud73].
• when T is normal [START_REF] Arch | Linear operator theory in engineering and science[END_REF].

In practice, to use the resolvent of light transport in this way, it is required to check whether the listed assumptions are provided or not. From this study, we know that the light transport is not compact in general. Therefore, in general, there is no information on the set eigenvalues to even be sure that they are discrete or isolated. The operator is not neither self-adjoint nor normal. Therefore, in the absence of completeness of light transport's eigenfunctions, decomposing an arbitrary light distribution using eigenfunctions is hopeless. Except, if the distribution belongs to the space spanned by eigenfunctions.

Discussion

In spite of the theoretical beauty of the resolvent study of an operator, it is impossible to apply the resolvent of light transport for any eigenpair-related objective. This is basically because of the nature of the resolvent concept and light transport. The resolvent concept is highly dependent on the singularities of the operator. On the other hand, the light transport operator does not have any nice properties, leading to extract information on the spectrum of T , such as compactness or normality. The only information, extracted from our study on light transport, is that the eigenvalues of the radiant exitance operator are real. When the operator is not compact, it is very challenging to obtain information about eigenvalues. Since the light transport operator is not compact in general, there is no mathematical proof to claim that eigenvalues are isolated, or there is an decreasing order on the set of eigenvalues 5 . About the analytic continuation of the resolvent, even under the assumption of compactness, it it essential to have information of the location of eigenvalues, which is not available for an arbitrary scene. Under the assumption of having this information too, the only way to perform the analytic continuation practically, is to use powers of T . On the other hand, there is a limitation of using Neumann series of T in the expression of the resolvent namely, z ∈ C\D(0, ∥T ∥), and computing the coefficients of each new representation of the resolvent suffers from numerical instability.

It is worth mentioning that, the Neumann series of the resolvent diverges inside D(0, ∥T ∥), but the resolvent obeys the Laurent series expression around each eigenvalue. However, in this case there is no practical way to do the job. That is also an obstacle to investigate the idea described in Section 5.3. Chapter 6

Scene

Stable computation of the largest eigenvalue

In this chapter, we estimate the dominant eigenvalue using a Monte Carlo formulation of the power iteration method. In this method, a sequence of vectors that converges to the eigenvector associated with the dominant eigenvalue is generated. To apply this method to the light transport operator, it is required to generate the images of the powers of light transport given an initial light distribution L. Since producing an image of long light paths, say for a large k, would be noisy with path tracing, we use Metropolis light transport (MLT) for T k L, with an additional bounce of path tracing for T k+1 L. Using this strategy, both images of T k L and T k+1 L will have the same scaling factor. Applying these two images, the dominant eigenvalue is obtained by computing the ratio of the two images.

6.1 Method's strategy

Power iteration for a matrix

The power iteration1 method is a simple eigenvalue algorithm. In the field of numerical analysis, eigenvalue analysis is referred to as the method in which the aim is computing the eigenvalues of a given matrix. This method starts with a random vector, say b 0 , and the goal is to achieve the eigenvector associated with the largest eigenvalue in magnitude. Namely, a sequence of vectors, (b k ), is generated, then under some sufficient conditions, the sequence converges to the eigenvector. The power iteration method is stated as follows:

1. Build a random vector to start the iteration process.

2. Generate vectors using the following formula

b k+1 = Ab k ∥Ab k ∥ , k ≥ 1. (6.1)
If the largest eigenvalue is strictly greater than others in magnitude, and b 0 has at least one non-zero component in the direction of the associated eigenvector, then (b n ) converges to that eigenvector [START_REF] Van Loan | Matrix computations (johns hopkins studies in mathematical sciences)[END_REF]. And the sequence

b t k Ab k b t k b k
converges to the dominant eigenvalue.

Power iteration for light transport

The light transport operator is an operator over an infinite-dimensional space, therefore it does not act necessarily as a matrix. To apply the power iteration method for this operator, it is required to adapt the arguments regarding the convergence issue of this method. We consider the light transport operator when it is compact, therefore all the eigenvalues are isolated. Define (a n ) as the sequence of the ratio of two consecutive powers, that is

a n := ∥T n+1 L∥2 ∥T n L∥2
. Since this sequence depends on L, its convergence also depends on the representation of this vector. If the set of eigenvectors is complete, then each light distribution can be expressed in terms of eigenfunctions, (Λ i ), provided by eigendecomposition:

T L = ∞ i=1 λ i < L, Λ i > Λ i . (6.2)
Using the equation above, we have

T k L = ∞ i=1 λ k i < L, Λ i > Λ i . (6.3) Therefore, lim n→∞ a n = |λ 1 |. (6.4)
Let I(T n+1 L) and I(T n L) indicate the correspondent images of T n+1 L and T n L respectively. Since the image is the result of a linear operation, projection, the sequence ∥I(T n+1 L)∥2 ∥I(T n L)∥2 also converges to λ 1 :

lim n→∞ ∥I(T n+1 L)∥ 2 ∥I(T n L)∥ 2 = |λ 1 |. (6.5)

Computing required images

To be able to obtain |λ 1 | from equation (6.5), it is required to have the images of T k L and T k+1 L for a large value of k. Since computing these images using path tracing methods consumes a large amount of time, and on the other hand the rendered image for a large k is noisy, the Metropolis light transport (MLT) method sounds a better option. However, in comparison to path tracing, the energy of the rendered image, the norm of the image vector is exact up to an unknown factor. Namely, if the exact energy of an image be E, then each run of MLT produces an image with a new unknown scaled energy, say αE. For our purpose, paying attention to this fact is very crucial. Therefore, we make a modification on the MLT algorithm, and also add a new step of path tracing to render the required images with the same α.

In Section 6.2.1, we review the original MLT algorithm, and in Section 6.2.2, we describe the way to compute the required images with the same scaling factor for the given Lambertian scene.

A brief reminder of Metropolis light transport (MLT) method

It was first Veach [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF] who applied the Metropolis algorithm, developed in computational physics, for transport problems. For the case of light transport, the aim is to compute an integral, D f (x)dµ x 2 , which is the energy transported by a set of light paths x ∈ D to the image plane, i.e. the value of this integral determines the illumination on the screen. The algorithm produces a Markov chain of light paths that follows a probability density proportional to f . It does so by mutating light paths and accepting them according to an acceptance probability that is computed as the ratio of the paths contributions normalized by their respective sampling probabilities. The advantage of the Metropolis method is that the rendered image is unbiased and it is usable for any scene with any material model and geometry. The Metropolis algorithm may be started with either a single light path or multiple ones. Generating an image with a single light path assumes that the ergodicity 3 of the set of mutation strategies used, which is not necessarily the case, depending on the implementation. Contrariwise, if multiple light paths are used, they need to be sampled using the same probability density, in order to keep the balance between the different related components (due to the mutations) of the total set of light paths. One solution to do this is to generate initial light paths proportional to f . Generating multiple light paths ensures building different families of light paths. Considering the case with multiple initial paths, the method is separately performed for each initial path in parallel, but the result of mutations on each initial path is recorded in a common image. Let X 0 be an initial path generated by a sampling method. Recall that MLT does not depend on X 0 , but it is important to generate the samples proportional to f . For example, using any importance sampling strategy provides us with this goal. Once X 0 is built, the mutation step is repeated long enough to reach the equilibrium distribution4 . Peskun [START_REF] Peter H Peskun | Optimum monte-carlo sampling using markov chains[END_REF] illustrated that to obtain this distribution faster, the best strategy is to maximize a(x → y), the acceptance probability of y given x, and a(y → x). It is sufficient to compute a(x → y) as follows a(x → y) = min{1, f (y)P (y → x) f (x)P (x → y) }, (6.6)

where P (• → •) is the tentative transition function5 , and can be chosen arbitrary.

Producing the required images with the same scaling factor

To compute the image of T k L using MLT first, it is required to generate at least one initial path, X 0 , with k number of bounces. This initial path starts from the given static light source, denoted by x s , and makes k number of intersections in the scene. Note that, the initial path is generated such that it is not connected to the camera at this stage, i.e. X 0 is of the form

X 0 = (x s , x 1 , • • • , x k ).
The reason behind this modification is that if the algorithm only selects the paths that can be connected to the camera after k bounces, then the result will be biased. Since the paths that can not be connected to the camera after k intersections can still contribute to the second image after making one more intersection in the scene. Connecting the last vertex of light paths to the camera is done in the step of projecting the contribution of paths into the screen plane, or in the stage of recording pixel values.

Once the initial path is generated, then the mutation process starts operating. Speaking of the mutation, we perform the point perturbation strategy on the given path, which preserves the length of the path [START_REF] Veach | Robust Monte Carlo methods for light transport simulation[END_REF]. Using this strategy, we keep the calculation of the acceptance probability simple. However, to ensure covering all necessary families of paths (ergodicity), more than one initial path is required. We explain the strategy for an initial path, and then generalize it for multiple seeds.

Generating images using a single seed path. We proceed the following three steps:

• Step 1. Producing a set of paths with length k, set D, by mutating one initial path or more.

•

Step 2. Adding one more vertex to each path belonging to D using path tracing.

• Step 3. Recording pixel values corresponding to each related path of length k and k + 1. After this step, we will have two separate images of T k L and T k+1 L, both proportional to the unbiased image, with a common unknown factor.

In the remainder of this section, we elaborate each of these steps in detail.

Step 1. To start the MLT algorithm, either with a single seed or multiple seeds, it is required to generate at least one initial path according to a probability density function proportional to f . In terms of practice, we sample X 0 using an importance sampling technique. First, we sample a direction at x s , according to the source emittance distribution, to trace a new point, x 1 . For the next k -1 points, at each current intersection point x t , we sample a direction according to the BRDF, to get a new surface point, x t+1 .

Once the initial path is available, then we mutate the path, using point perturbation, many times. To review this mutation technique briefly, consider an index p ∈ {1, • • • , k}. The result of the mutation on the current path X j will be a new path which has the same vertices, except the point with index p, x p , we denote the index of this new point by p ′ . That is, if X j+1 indicates the new path, then

X j = (x s , • • • , x p , • • • , x k ) , X j+1 = (x s , • • • , x p ′ , • • • , x k )
, under the assumption that p ̸ = k. If p = k, then

X j = (x s , • • • , x p ) , X j+1 = (x s , • • • , x p ′ )
Note that, x p ′ is chosen such that V (x p ′ , x p+1 ) = 1, where V (•, •) is the visibility function. Therefore, one may summarize this approach with the following steps:

1. Choosing an index p ∈ {1, • • • , k}.

2. Sampling a direction from x p-1 to trace x p ′ .

3. Checking the visibility condition on x p ′ to decide whether accept or deny this point.

4. If x p ′ is accepted then, subpath x p ′ x p+1 contributes in the calculations with the following factor

ρ(x p-1 → x p ′ → x p+1 )ρ(x p ′ → x p+1 → x p+2 )G(x p ′ , x p+1 ). (6.7)
Recall that, in our implementation, the paths are not connected to the camera in the first place. Therefore, we do not involve any coefficient related the last vertex of the path while computing the acceptance probability of the following cases: p = k and p + 1 = k. We will take care of this coefficient in the projection step. So, to compute the acceptance probability, we separate the following three cases:

• p = k: a(X j -→ X j+1 ) = 1 (6.8) • p + 1 = k: a(X j -→ X j+1 ) = ρ(x k-2 → x p ′ → x k )G(x p ′ , x k ) ρ(x k-2 → x p → x k )G(x p , x k ) (6.9) • p + 1 < k: a(X j -→ X j+1 ) = ρ(x p-1 → x p ′ → x p+1 )ρ(x p ′ → x p ′ +1 → x p+2 )G(x p ′ , x p ′ +1 ) ρ(x p-1 → x p → x p+1 )ρ(x p → x p+1 → x p+2 )G(x p , x p+1
) .

(6.10)

To explain the terms, appearing in equation (6.9) and (6.10), we look back at the mutation strategy used. The acceptance probability is computed in terms of factors that form equation (6.7).

Step 2 and 3. In this part, we add the contribution of each related path to the corresponding image, considering some factors coming from first step of this method.

In the remainder of this section, we discuss factors that are required to be multiplied by the contribution of the paths of length k and k + 1.

Image of T k L. Let {X 0 , • • • , X d }
be the sequence of paths generated by repeating mutation operation on paths, starting from X 0 as an initial path. Let X j be the current light path, and X j+1 be the new path as a result of one mutation operation. To reduce the variance, the new path is always accepted. Therefore, the contribution of both old and new paths are scaled with factors 1 -pr mu and pr mu respectively, where

pr mu := min{1, |a(X j → X j+1 )|}. (6.11)
Given X, let (u X , v X ) denote the corresponding pixel indices in the image plane, I, which are determined by the last vertex of X. Therefore,

I(u Xj+1 , v Xj+1 ) = pr mu c x k xe C(X j+1 ) (6.12)
, and

I(u Xj , v Xj ) = (1 -pr mu ) c x k xe C(X j ). (6.13)
In both equations (6.12) and (6.13), C(X) = f (X) p(x=X) indicates a function that returns the contribution value of a path, which is sampled according to a probability density function p(x). c x k xe represents the factor, corresponding to the fact of connecting x k to the camera, computed as follows

c x k xe = ρ(x k-1 → x k → x e )G(x k , x e )V (x k , x e ),
where x e denotes the camera 6 . So what remains is to compute C(X j ) and C(X j+1 ).

In our implementation, each vertex of X 0 is sampled using importance sampling. Namely, vertex x t+1 is traced from x t according to the BRDF at this point. Therefore, the contribution of X 0 is

C(X 0 ) = Π t=2 ρ(x t-1 → x t → x t+1 ) Π t=2 ρ(x t-1 → x t → x t+1 ) = 1. (6.14)
Equation (6.14) is true since at each surface point, light reflects with the probability equal to the BRDF value. Similarly, C(X j ) = C(X j+1 ) = 1. Image of T k+1 L. Constructing the second image requires paths of length k +1.

Since in this method we use the same paths that already contributed to the first image, by sampling a new vertex x k+1 . Similar to other vertices along the given path X j , we sample a direction at vertex x k to trace x k+1 . We denote this new path with Y j , Figure .

Similar to equations (6.12) and (6.13), the following values are added to the second image

I(u Yj+1 , v Yj+1 ) = pr mu c x k+1 xe C(Y j+1 ), (6.15)
, and

I(u Yj , v Yj ) = (1 -pr mu ) c x k+1 xe C(Y j ) (6.16)
, where

c x k+1 xe = ρ(x k → x k+1 → x e )G(x k+1 , x e )V (x k+1 , x e ).
To compute C(Y j ) (C(Y j+1 )) from C(X j ) (C(X j+1 )), it is supposed that light is reflected at x k . To consider this fact, it is required to scale these values with this probability. In the specific case of Lambertian scenes, this value is equal to the albedo, so

C(Y j ) = C(X j )ρ d (x k ) = ρ d (x k ), (6.17) similarly C(Y j+1 ) = C(X j+1 )ρ d (x k ) = ρ d (x k ). (6.18)
Generating images using multiple initial paths. Image of T k L. In this case, we repeat the three steps, expressed earlier for each initial path, for each initial path individually.

To ensure ergodicity, we need to sample all initial paths, {X (i) 0 } i , using a common probability distribution. This distribution needs to be proportional to the image on the screen. To do this, we sample points uniformly in the scene and then form paths in the way that any two consecutive vertices are not on the same surface i.e. subpath x t x t+1 belongs to

X (i) 0 if V (x t , x t+1 ) = 1 , < n xt , n xt+1 >< 0.
Then, each initial path is mutated in parallel and contributes to the same image using equations (6.13), (6.12), (6.16) and (6.15).

The only difference is the value of C(X (i) 0 ). Since each initial path is built explicitly, so C(X (i) 0 ) is the multiplication of point-to-point form factors, i.e.

C(X

(i) 0 ) = Π t=k-1 t=1 ρ(x t-1 → x t → x t+1 )G(x t , x t+1 ). (6.19) Since each X (i)
j+1 is obtained by mutating X (i) j , and variable pr mu has already accounted for the modification applied on

X (i) j , so C(X (i) j ) = C(X (i) j+1 ) = C(X (i) 0 ). (6.20)
The values of C(Y

(i) j ) (C(Y (i) j+1
)) are computed using its connection with C(X

(i) j ) (C(X (i) j+1 )): C(Y (i) j ) = C(X (i) j )ρ d (x k ).
(6.21)

Results

To obtain the dominant eigenvalue, we experimented the power iteration method on our test scenes. For each scene with a static light source, the images of pairs of powers of T have been computed. Then, the graph of the ratio of images Figure 6.3: Noise-free images of T 9 L and T 10 L of the "ThreeSquaresNoSymmetry" scene, computed using multiple initial paths.

T k L and T k+1 L, for k = 2, • • • ,

Discussion

In this section, we discuss two subjects related to the theory and implementation parts respectively.

Power iteration method

To be able to analyze the power iteration algorithm for the light transport operator, the compactness property is required to limit the spectrum set to the set of eigenvalues. We also make an extra assumption that the set of eigenvectors span H. So, every light distribution L can be written in terms of eigenvectors i.e.

L = ∞ i=1 c i Λ i .
Figure 6.4: Noise-free images of T 9 L and T 10 L of the "Contact" scene, computed using an initial path.

Figure 6.5: Comparison of the method's performance on "ThreeSquaresNoSymmetry" by increasing the number of initial paths from 100 to 10000.

b k+1 : = T b k = T k L = T k ( ∞ i=1 c i Λ i ) = ∞ i=1 c i λ k i Λ i = λ k i [c 1 Λ 1 + ∞ i̸ =1 c i ( λ i λ 1 )Λ i ]. (6.22) Equation (6.22) indicates that that sequence (b k ) converges to Λ 1 if | λi λ1 | < 1. Therefore, if λ 1 = -λ 2 , then (b k ) does not converge to Λ 1 . Depending on index k, (b k ) flip-
flops in the direction of two vectors Λ 1 +Λ 2 and Λ 1 -Λ 2 . This reason may explain the method's behavior performed on the "TwoSquaresNoSymmetry" scene. However, this event does not happen for "TwoSquares". This scene is symmetric and if the light source lies in the middle of two surfaces, it may cause orthogonality of L to one of the eigenvectors Λ 1 or Λ 2 . Therefore, the corresponding coefficient (c 0 or c 2 ) is zero, which avoids the method's failure. Remark 6.4.1 Note that the compactness and the completeness properties are sufficient conditions to explain how the method works. Therefore, the method may still work for the scenes without these properties, such as the "Cornell-Figure 6.6: Performing the method on the "Contact" scene with a single initial path.

Box.ifs" scene.

Sampling strategies for long light paths

The problem of choosing a sampling strategy for long light paths is the following:

1. Sampling a path explicitly, sampling vertices uniformly then connecting them, is costly because of visibility calculation, in particular, in scenes where the couples of points have a high probability of not seeing each other.

2. In the case following importance sampling, it is hard to get a full light path of length k in open scenes, because the geometry is easily be missed (e.g. the "TwoSquares" scene has this problem).

It is generally better to try both strategies, assuming that their respective probability densities are correctly set. This the reason that we have used both strategies, using the importance sampling technique for a single initial path, and building multiple initial paths explicitly. As it is expected, using more that one initial light path ensures ergodicity. Considering the "ThreeSquaresNoSymmetry" scene, the curves of the eigenvalue estimation in Figure 6.5 illustrates more accurate estimation versus increasing the number of initial light paths. However, for a very specific scene among our test scenes, "Contact" scene, the values computed by equation (6.19) become very small, so it will increase the numerical instability, Figure 6.1a. Since this scene is very simple, using even one initial path still estimates the dominant eigenvalue, Figure 6.6.

Chapter 7

Calculation of eigenvalues via regularized determinant

It is known that the Fredholm theory is a useful tool in several fields such as operator theory, mathematical physics and quantum field theory. Therefore, the numerical methods to evaluate the so-called Fredholm determinant is very critical. But choosing a numerical method to fulfill the task is highly dependent on the properties of the kernel, such as continuity, singularity, and the domain of the kernel function. For an operator defined on a Hilbert space, based on the assumption of Hilbert-Schmidt property of the operator, the eigenvalues of the operator are the reciprocals of a determinant, called the regularized determinant, which is not exactly the Fredholm determinant but relates to it. Similar to the Fredholm determinant, this determinant is also expressed in terms of some specific coefficients, denoted by c n and p n in this chapter. To evaluate p n coefficients for the light transport operator, we derive an interesting equation in terms of path integrals. Therefore, it is possible to compute the determinant using the Monte Carlo algorithms. Despite of exploring all the connections between these coefficients and eigenvalues of the light transport operator, our implementation suffers from a numerical obstacle. In Section 7.1, we recall the original definition of the Fredholm determinant, and then we address briefly other similar constructions of the determinant for the specific case of the Hilbert-Schmidt operator. In Section 7.2, we explain how to compute the regularized determinant of light transport. In Section 7.3, c n and p n coefficients with the eigenvalues of the light transport operator evaluated for our test scenes are provided. We finish the chapter by a discussion on the numerical results, Section 7.4.

Fredholm determinant and its generalizations

It was Fredholm himself who introduced the determinant defined by (7.2) for the first time, and then the subject started to grow and it has been generalized

to the operators over infinite-dimensional spaces. In addition to the great interest of the theoretical point of view, the numerical evaluation of the Fredholm determinant has a great significance. Regardless of the theoretical development, there is no general framework for numerical methods to compute the Fredholm determinant since implementations depend on the kernel's property. For example, Nystrom-type method seems fit to the case when the kernel is smooth, or some other ways are suitable when the operator is self-adjoint, for example Gaussian unitary ensemble (GUE) problem, see [START_REF] Gaudin | Sur la loi limite de l'espacement des valeurs propres d'une matrice aléatoire[END_REF] and [START_REF] Madan | Random matrices[END_REF].

When the operator is trace class, there is a projection-based method. This method requires an orthogonal basis which projects the Hilbert space to some finite subspaces, check Bornemann et al.'s paper [START_REF] Bornemann | On the numerical evaluation of fredholm determinants[END_REF]. In Section 7.1.1, we start with the Fredholm integral equation of the second kind under the assumption that the unknown solution function is an integrable function on (a, b). In Section 7.1.2, we list equivalent infinite-dimensional determinants that are applicable when the operator is trace class. In Section 7.1.3, we review the regularized determinant, which is used in the place of Fredholm operator when the operator is Hilbert-Schmidt.

The original definition

Consider the following equation

ϕ(x) = f (x) + µ b a κ (x, y)ϕ(y)dy , x ∈ (a, b) (7.1)
, which is called a Fredholm equation of the second kind when the kernel is bounded and continuous. Fredholm showed that equation (7.1) is uniquely solvable if and only if d(µ) with the following expression

d(µ) = det(I+µT ) = ∞ n=0 µ n n! (a,b) n det      κ(t 1 , t 1 ) κ(t 1 , t 2 ) • • • κ(t 1 , t n ) κ(t 2 , t 1 ) κ(t 2 , t 2 ) • • • κ(t 2 , t n ) . . . . . . . . . . . . κ(t n , t 1 ) κ(t n , t 2 ) • • • κ(t n , t n )      dt 1 • • • dt n (7.
2) be a non-zero. d(µ) is known as Fredholm determinant. Searching for the solution of equation (7.1), ϕ(x), this unknown function is assumed to be Riemann integrable. One may also make extra assumptions on f (x) and κ such as continuity on the interval (a, b) and on the square (a, b) × (a, b) respectively. Define c n coefficients as

c n = (a,b) n det      κ(t 1 , t 1 ) κ(t 1 , t 2 ) • • • κ(t 1 , t n ) κ(t 2 , t 1 ) κ(t 2 , t 2 ) • • • κ(t 2 , t n ) . . . . . . . . . . . . κ(t n , t 1 ) κ(t n , t 2 ) • • • κ(t n , t n )      dt 1 • • • dt n , (7.3)
which reformulates equation (7.2) into

d(µ) = ∞ n=0 µ n n! c n .
Let B 0 (x, y) := κ(x, y), and c 0 = 1. For n ≥ 1 define B n (x, y) as a recursive linear combination of iterated kernel

B n (x, y) = c n κ(x, y) -n b a κ(x, y)B n-1 (y, z)dy, (7.4)
be a family of functions defined in terms of c n . With this setting, it can be shown [START_REF] Stephen M Zemyan | The classical theory of integral equations: a concise treatment[END_REF] c n = b a B n-1 (x, x)dx. (7.5)

Fredholm determinant of a trace class operator

Let H be a Hilbert space and T : H -→ H be a trace class operator. One may find several equivalent reformulations of

d(µ) = det(I + µT ),
for µ ∈ C(̸ = 0). We address some of these constructions as follows:

• Gohberg et al. defined d(µ) as an infinite product [START_REF] Gohberg | Introduction to the theory of linear nonselfadjoint operators[END_REF], namely

d(µ) = N (T ) n=1 (I + µλ n ), (7.6)
where

(λ n ) N (T )
n=1 , counting multiplicity, is either finite or infinite set of eigenvalues. The infinite product on the right side of equation (7.6) is locally uniformly convergent. He called this determinant as characteristic determinant, and showed that it is bounded [START_REF] Gohberg | Introduction to the theory of linear nonselfadjoint operators[END_REF] i.e.

|d(µ)| ≤ ∞ j=1 (1 + |µ|s j (T )) ≤ exp |µ|∥T ∥ J 1 .
• Consider (T n ) as a sequence of finite-rank operators that converges to T , in the trace norm, ∥ • ∥ J11 . Then, the sequence of (det(I + µT n | ran(Tn) )) converges locally uniform to d(µ) [START_REF] Gohberg | Classes of linear operators[END_REF]. Note that the definition of d(µ) as this limit is independent of the choice of (T n ).

• Dunford et Schwartz d(µ) in terms of the trace of powers of T , tr(T n ),

det(I + µT ) = exp( ∞ n=1 (-1) n+1 n µ n tr(T n )), (7.7)
which converges in {µ ∈ | |µ| < 1 |λ1| }. This definition comes from the theorem called Plemelj-Smithies formulas [START_REF] Gohberg | Traces and determinants of linear operators[END_REF]. Note that using analytic continuation, d(µ) can extended for all µ ∈ C.

The regularized determinant for a Hilbert-Schmidt operator

Let H be a Hilbert space and T : H -→ H be a Hilbert-Schmidt operator. In comparison to the trace class property,

N (T ) n=1 λ 2 n is finite but not necessarily N (T ) n=1 λ n .
Therefore, the definition of d(µ) by equation (7.6) is not meaningful, since it may not converge. In replace, one may define d 2 (µ) as follows

d 2 (µ) := det 2 (I + µT ) := N (T ) n=1 (1 + µλ n ) exp -µλn .
This definition introduced by Hilbert (1904) and Carleman (1921), and it is also known as Hilbert-Carleman determinant. Originally, det 2 (I + T ) is defined as det[(I + T ) exp(-T )] [GGK96]. For this case, equation (7.7) turns to

det 2 (I + µT ) = exp( ∞ n=2 (-1) n+1 n µ n tr(T n )), 2 .
For a specific case, when H is L 2 (a, b), there exists an expression for det 2 (I +µT ) similar to the one for d(µ) in equation (7.2). However in this case κ(x, x) is set to zero while computing the determinant, that is

∞ n=0 µ n n! (a,b) n det      0 κ(t 1 , t 2 ) • • • κ(t 1 , t n ) κ(t 2 , t 1 ) 0 • • • κ(t 2 , t n ) . . . . . . . . . . . . κ(t n , t 1 ) κ(t n , t 2 ) • • • 0      dt 1 • • • dt n . (7.8)
Hilbert proved that the expression above converges to det 2 (I + µT ) [START_REF] Hilbert | Basic principles of a general theory of linear integral equations[END_REF].

Calculation of regularized determinant of light transport

To be able to reach eigenvalues of light transport using equations developed in Fredholm theory, it is essential to choose an appropriate equation based on the property of the operator. Considering the light transport operator, acting on a Hilbert space, we can only use equation (7.8) in a very specific situation, when the operator is Hilbert-Schmidt. This property holds for the Lambertian light transport operator when we remove the edges from the scene.

Computing c n coefficients of the regularized determinant

To evaluate d 2 (µ) using equation (7.8) numerically, it is required to evaluate the interior integral, which is c n coefficients, defined by equation (7.3). To calculate these coefficients, we only consider the Lambertian materials. Recall the kernel expression of the light transport operator, equation (3.35). Since the material model is Lambertian, then

∀x i , x j ∈ S : κ(x i , x j ) = ρ d (x i ) cos θ cos θ ′ πr 2 x i x j V (x i , x j ). (7.9)
To compute c n coefficient, given n:

1. n number of surface points are sampled in the scene, 2. using equation (7.9), the matrix containing values of the kernel evaluated at a pair of surface points, selected in the previous step, is constructed, 3. the determinant of the matrix, which is constructed in the previous step, is computed.

The integral of cyclic light paths

In this section we define a quantity, p n , which brings a very interesting interpretation for general scenes. Under the assumption of Hilbert-Schmidt property for the light transport operator, we derive a new equation, indicating the relationship between c n and p n coefficients.

We define Π n (x) : S -→ S as a path integral given by the following formula

∀x ∈ S : Π n (x) := S n-1 κ(x, x 2 ) • • • κ(x n-1 , x n )κ(x n , x)dx 2 • • • dx n (7.10) in which (x i ) n i=2 ⊂ S.
In Lambertian scenes, this quantity can be interpreted as the integral of the amount of energy returning to a surface point x, after propagation along all paths starting from vertex x and bouncing n -1 times in the scene. Then define

p n := S Π n (x)dx (7.11)
, which evaluates the density of all cyclic light paths3 with length n in the scene.

The relationship between c n and p n coefficients. Starting from equation (7.4), we have

B n (x, x) = -n S κ(x, s)B n-1 (s, x)ds (7.12)
, in which we put κ(x, x) = 0. One may use equation (7.4) to expand B n-1 (s, x), so

B n (x, x) = -nc n-1 S κ(x, s)κ(s, x)ds + n(n + 1) S S κ(x, s)κ(s, t)B n-2 (t, x)dtds = -nc n-1 Π 2 (x) + n(n + 1) S S κ(x, s)κ(s, t)B n-2 (t, x)dtds, (7.13)
if we continue applying equation (7.4), and using B 0 (x, x) = κ(x, x), then

B n (x, x) = n k=1 (-1) k Π k+1 (x) n! (n -k)! c n-k . (7.14)
Now plugging the previous expression of B n (x, x) in equation (7.5) yields

c n+1 = x∈S ( n k=1 (-1) k Π k+1 (x) n! (n -k)! c n-k )dx = n k=1 (-1) k p k+1 c n-k n! (n -k)! .
(7.15) Therefore, equation (7.15) is another way to compute c n coefficients from the integrals of cyclic light paths.

Results

In this section, we provide the values of c n , computed using equation (7.3) in which the diagonal of the matrix is set to zero, and p n coefficients, using equation (7.11). The eigenvalues of the light transport operator are also calculated solving for the roots of equation (7.2). c n coefficients. A number of surface points are sampled in the scene. To construct the matrix in equation (7.8), it is required to evaluate the kernel at pairs of sampled points. To do this task, we construct the matrix coefficients in two ways:

• Method 1. Using equation (7.9) directly.

• Method 2. To decrease the numerical instability caused by the infinitely large values of the kernel evaluations for pairs of points that are located arbitrarily close, the kernel values, which are point-to-point form factors, can be approximated by an analytic formula, which is the point-to-disk form factors [START_REF] Mhn Naraghi | Radiation view factors from differential plane sources to disks-a general formulation[END_REF]. Using this approximation, the kernel values are bounded. p n coefficients. To compute p n coefficients using the Monte Carlo algorithm, set of N points are sampled uniformly in the scene for the starting points of the cyclic paths. By doing the change of variable from point surfaces to the projected solid angles in equation (7.10), and growing the path according to the cosine-weighted BRDF distribution, the estimation of the p n values simplifies to

p n = A N N i=1 Π n (x i 1 ) P (Π n (x i 1 )) = A N N i=1 κ(x n , x i 1 )
(7.16)

Discussion

Considering the fact that the light transport operator is acting on infinitedimensional space, the generalization of the Fredholm determinant to the operators over infinite-dimensional spaces is required. On the other hand, the Fredholm determinant is generalized for the operators with specific properties, such as compactness, or even stronger assumptions on the operator, such as being Hilbert-Schmidt or trace class. Therefore, the background theory to ensure the convergence of the calculations to the real values is limited to these specific classes of operators. We know that the light transport operator is not even compact in general (Theorem 4.4.5). So, we restrict the discussion to a very specific scene, the scene without edges, since the operator is Hilbert-Schmidt in this situation (Theorem 4.4.4), and the regularized determinant is allowed to be used for the operator. Among our test scenes, the scenes in which the operator is Hilbert-Schmidt are the "TwoSquares" scene, the "Cubes0027" scene, and the "Cubes0343" scene. Although the theory ensures extracting eigenvalues from this determinant, only the estimated eigenvalues of the "TwoSquares" scene are close to the reference values. The failure of the method for the two other scenes may originate from calculating lots of visibility calculations between the cubes conjugated with the fact that computing the roots of a polynomial is highly sensitive to the polynomial coefficients [START_REF] Forman S Acton | Numerical methods that work[END_REF].

Chapter 8

Eigenvalue extraction based on a polynomial approach using path tracing

Under the assumption of the compactness property, we know that the eigenvalues of the light transport operator are isolated, and have a decreasing order( tend to zero). On the other hand, the spectrum of the operator is contained in the unit disk. Benefiting from these properties, it is possible to consider a finite number of these eigenvalues and construct a polynomial of a finite degree and compute its roots by numerical algorithms. We develop a method to approximate the eigenvalues of the light transport operator by computing the roots of the (estimated) characteristic polynomial, depending on the initial distribution. To compensate for this dependency, we form histograms of the roots and obtain better estimates of the eigenvalues. In Section 8.1, we explain the theoretical background of the method. In Section 8.2, we perform the method on a matrix with the given properties such as, size, rank, explicit form of the matrix, and eigenvalues. In Section 8.3, we perform the method for light transport both in screen-space and object-space. However, since the method is originally designed based on histograms, we focused on the performance of the object-space run. Therefore, only a static scene is used for the former and different test scenes for the latter. We finish the chapter by a discussion on the designed method, Section 8.4.

Method's strategy

Suppose that T be the light transport operator when it is compact. Therefore, the eigenvalues of this operator are isolated. Let {λ 1 , • • • , λ m } be the set of first m dominant eigenvalues of T and T m : ⊕Λ i -→ ⊕Λ i indicates a finite rank approximation of the light transport operator of rank m. Let p m (x) = m k=0 α k x k denote the characteristic polynomial of T m . From Cayley-Hamilton theorem1 , we have p m (T ) = 0. Let (p i ) 1≤i≤n be a set of points sampled in the scene, and J be a matrix in which the columns of J are filled by the evaluations of T k L at these points, i.e. t ik := (T k+1 L)(p i ):

J =    v 1 t 11 • • • t 1m . . . . . . . . . v n t n1 • • • t nm    , (8.1) in which v = (T L)(p i ) 1≤i≤n . Practically, it is possible to estimate α = (1, α 1 , • • • , α m ) t
2 by solving the following system

   v 1 t 11 • • • t 1m . . . . . . . . . v n t n1 • • • t nm       1 . . . α m    = 0 (8.2)
, which the solution can be approximated by the following optimization problem

α * := argmin α ∥Jα∥ 2 .
The optimization problem above is equivalent to the following problem

α * := argmin α ∥J ′ α -v∥ 2 (8.3)
, where J ′ contains the same columns as J, except the first one. One may solve the previous problem by differentiating ϕ(α) := ∥J ′ α -v∥ 2 i.e.

ϕ(α) = (J ′ α -v) t (J ′ α -v) = α t J ′t J ′ α -2α t J ′t v + ∥v∥ 2 dϕ dα = 2J ′t J ′ α -2J ′t v, therefore α * = (J ′t J ′ ) -1 (J ′t v
) is an approximation of the coefficients of p m (x). So, depending on v, some of the eigenvalues of T m can be obtained by extracting the roots of p

* m (x) = k=m k=0 α * k x k .

Validation of the method on the matrix

To evaluate the method's performance for Lambertian scenes, we examine the method on two matrices A and B. We build these matrices such that their eigenvalues are real, have a magnitude less than 1, and follow a decreasing order.

Let A be a full-rank matrix of dimension 50 with the following eigenvalues It is clear that this sequence has the properties expressed earlier.

λ i = 0.8e -i , i = 0, 1, • • • , m -1.
Extracting the roots of a polynomial is highly sensitive to its coefficients [START_REF] Forman S Acton | Numerical methods that work[END_REF], therefore, we do not expect to approximate the eigenvalues by using only one individual v. But, one may repeat the process several times to plot a histogram of the computed roots. Therefore, the picks in the histogram are better approximations of the related eigenvalues. We construct 10 4 samples of J ′ , by computing A k v for k = 0, . . . , 25, while randomly changing v. We denote these samples with (J ′ j ) 1≤j≤M . For each J ′ j , the α * vector is computed then, the roots of polynomial p 25 (x) are calculated.

Remark 8.2.1 There is a remarkable difference between this matrix simulation and light transport, rather than just the dimension. It concerns the representation of light transport. In our method, we do not use any matrix representation of the light transport operator to compute the columns of J ′ . Mathematically speaking, we keep using the operator (continuous) representation of light transport. This means that a path tracer computes each power of T for a given initial light distribution L using a Monte Carlo algorithm. In this way, we avoid discretizing scenes3 . However, in this simulation, we produce A k v vectors by matrix-vector multiplication using the explicit form of the matrix directly.

In Figure 8.1a, the histogram of the roots of p 25 (x) for matrix A is given. Reading eigenvalues with index i ≥ 5 from the histogram in Figure 8.1a is not possible. Therefore, we look at the data from which the histogram is plotted. We observe that the estimation of eigenvalues becomes less accurate after i = 5, see Table 8 smaller when index i grows. Moreover, the magnitude of each eigenvalue is also decreasing.

We also examine the method on another full-rank matrix, B. We set the eigenvalues of this matrix to be ( 1 2 i ). We repeat the process on this matrix, where m = 20. However, the two peaks corresponding to λ 5 and λ 6 are less sharp in the histogram of Figure 8.1b, but the data of this experiment illustrates that the method still could catch these values, see Table 8 Looking at these values and the number of extracted eigenvalues by the method may imply the dependency of the process on the root's properties.

.2. Let d i = λ i -λ i+1 . In

The method's performance on light transport

This method can be performed both on screen-space and object-space. For the former, a static light source and camera are used to produce a sequence of images using a path tracer. Each image of this sequence is computed by generating a large number of light paths, starting from the light source and making k number of intersections in the scene before reaching the camera. The disadvantage of using a static light source is that p m will depend on the choice of L, so it will make the estimation of this polynomial biased by the choice of the initial distribution. To illustrate this fact, in Section 8.3.1, we perform the method on a single scene ("CornellBox" scene).

Considering the latter case, in the previous section, it was shown how the method estimates the suitable eigenvalues, the ones which are not too close, by constructing a histogram of the computed roots from the samples of J ′ , (J ′ j ) 1≤j≤M , where M is a large positive value. Implementing this method for the light transport operator in the object-space contains two important key parameters: the number of rows of J ′ j , n, and M . Unlike the case of the matrix, in which the rank of the operator is finite, no information is available on the rank of the global light transport operator. On the other hand, applying the large values of n and M increases the cost. So, the question that arises here is What are the minimum possible values of these two parameters to achieve good approximations? Answering this question requires insight into the connection between the scene properties such as the material and the geometry, with the rank of the global light transport operator, moreover, the level of the importance between n and M in comparison to each other. Despite the interest in answering such a question, we do not discuss this issue here, and the question remains open. In Section 8.3.2, the results of applying this method on the object-space is disposed. In this case, we do not compute the full images but only for a limited collection of points (p i ) 1≤i≤n . For this case, the point source, denoted by x s , is changed randomly to build new samples of J ′ j , and the histogram of the roots of p m (x) is plotted for the test scenes. We experiment with the "CornellBox" scene in which the location of the light source is fixed in the ceiling of the scene, Figure 8.4. To render each image of T k L, given k, lots of light paths of the form x = (x s , x 1 , • • • , x k , x e ) contribute to the final image. To generate each path, first, x s is sampled on the area light, then the path is grown by the Russian roulette technique. Adding a value to each corresponding pixel, is the task of projection step. That is, the reflected radiance from the last vertex of the suitable4 path is multiplied by a factor, which we call it sensor factor. Therefore, this operation is linear, and feasible on the image-space. The images of T k L that fill matrix J ′ are given in Figure 8.35 . The roots of the polynomial obtained from J ′ and the reference eigenvalues are given in the following table.

Results with a single light source in image space

Results with multiple light sources

For each row of the matrix J ′ , one point, say p i , is sampled in the scene, then the powers of T L are computed at this point. Computing (T k L)(p i ) means to collect the energy, transported along the paths of length k, at this point. For a given k, since the last vertex of all paths must be p i , it is required to connect the point x k to p i . In practice, it means to multiply the contribution of the subpath (x s , • • • , x k ) by the following factor

c x k p i = ρ(x k-1 → x k → p i )G(x k , p i )V (x k , p i ). (8.4)
In our implementation, we grow the path applying Russian roulette, that is growing the path stops until it is absorbed with an absorption probability. Let this absorption probability be 1 -ρ(x). Therefore, light reflects with the probability which is exactly equal to ρ(x). Starting from x s , to trace x 1 , a direction is sampled according to the cosineweighted BRDF distribution. For the next intersection points, at each current intersection point x t , we sample a direction according to the BRDF. Therefore, the contribution of the generated path with h number of bounces is

C(x s , x 1 , • • • , x h ) = Π h-1 t=1 ρ(x t-1 → x t → x t+1 ) Π h-1 t=1 ρ(x t-1 → x t → x t+1 ) . (8.5)
Since our method requires evaluating the contribution of the paths with k bounces in the scene, so the values, that are provided by equation (8.5), must be multiplied by the related factors. For each generated path by the Russian roulette technique, we construct the paths with different lengths 2 to k. Then, we evaluate the contribution of each path accordingly. Therefore, generating several light paths, and using the contribution measurements of all paths with different lengths, will lead to obtain J ′ j matrix, needed to compute α * . So what remains is to compute the contribution of paths with different lengths correctly. Let (x s , • • • , x h ) be a path generated by the Russian roulette method. If the length of this path h is less than k, then it adds no energy to the estimate of (T k L)(p i ). But if k ≤ h then, the contribution of the path is non-zero. We separate the two cases k < h and k = h respectively:

• when k < h, means light is reflected at x k , but we need to cut the path at this point, and then connect it to p i , see Figure 8.5a. So • if k = h, means the path is absorbed at x k . It only requires to connect it to p i . Therefore,

C(x s , • • • , x k , p i ) = c x k p i ρ(x k ) (8.6) (
C(x s , • • • , x k , p i ) = c x k p i 1 -ρ(x k ) (8.7) In equation (8.6), C(x s , • • • , x k , p i ) is divided by ρ(x k )
to compensate for the reflection at x k , however, in equation (8.7), C(x s , • • • , x k , p i ) is divided by 1 -ρ(x k ) to compensate for the absorption. The remainder of this section provides the pseudocode algorithm of the method. The method is implemented in the framework of a function, see Algorithm 2. There is a key function, we called it "pointSourcePathTrace", which computes the rows of the matrix J ′ in each iteration. The pseudocode of this function is also committed in Algorithm 1 separately. There are other functions applied in these two functions, such as "growPath", "computeAlphaValues" and "zroots". We briefly explain these functions.

• "growPath(x s )". This function generates a path, starting from a given surface point x s , by the Russian roulette technique.

• "computeAlphaValues(J ′ )". This function solves equation (J ′t J ′ )α = J ′t v for α using conjugate gradient method.

• "zroots(α * )". This function computes the roots of a polynomial of a finite degree with Laguerre's method6 .

In Table 8

Discussion

There are several methods to compute the roots of polynomials, such as Newton's method, Laguerre's method, Lehmer method, Root-squaring methods and many others. Regardless of such a simple task, there is no general method to extract all roots of a polynomial. This fact is because of the dependency of this procedure on many factors, especially root properties, like being real or complex, isolated or not, the algebraic multiplicity, distance of the roots with respect to each other, etc. Therefore, choosing an approach depends on the primary information (in the case of existence) on such factors and the purpose of the experiment. Computing the roots in our method is no exception among others. These factors influence the results produced by this method. The other disadvantage of extracting roots by machines is that, a very small change in the coefficients of a polynomial may change the expected roots notably. Wilkinson et al. [START_REF] James H Wilkinson | The evaluation of the zeros of ill-conditioned polynomials. part i[END_REF] showed that making a little change on a coefficient of a polynomial can even lead to get false complex roots in the place of reference real ones.

It is worth noting that, the problems do not only limit to the change of coefficients, but also on the roots itself. That is, most of methods do not obtain the roots at once, some roots are computed earlier. For example, if the polynomial contains complex roots, then the real roots are computed first. In such methods, if one root is not calculated accurately enough, then it may also interrupt extracting the rest of roots to search for. Wilkinson also revealed that the failure originates from the regular distance of the roots with respect to each other. The worst-case scenario is when no primary information is available, which is generally the case. In this situation, it is required to have a system equipped with the decision making ability, to choose the method depending on the problem. Designing such a system is far more complicated than the original issue, which is computing the roots.

Considering the light transport operator, we know that the eigenvalues are real in Lambertian scenes (See Theorem 4.A.3). However, there is no information available on the multiplicity of the eigenvalue. Such information is precious, since extracting multiple roots cause troubles for numerical algorithms, such as Newton's method7 . For a special geometry of the scene, in which the edges are removed, the light transport operator is compact. Therefore, the non-zero eigenvalues are isolated with finite algebraic multiplicity, but no detection of the multiplicity. Among our test scenes, the light transport acting on "TwoSquares-NoSymmetry", "TwoSquares", "Cubes0027" and "Cubes0343" is compact. In our implementation, the task of extracting roots is done by Laguerr's algorithm. This algorithm is able to obtain real and complex roots from complex coefficients. The advantage of this algorithm is that the method always converges to some root of the polynomial independent of the initial choice. This algorithm is also able to extract the roots with non-trivial multiplicity. However, this method is developed under some assumptions namely, the first root be isolated and the other roots to have some distance with the first root. The isolation assumption of all roots is not included in these assumptions, but very close space of roots increase the error of estimating roots using this algorithm. This originates from the stop condition of the iteration loop. Let q(x) be a polynomial of degree m. Then this algorithm estimates a root, we denote it by r, by choosing trials rk+1 = rk -a. Variable a is defined as the following fraction n qn ] 2 -q ′′ n qn . The iteration loop stops when a is small. Therefore, even considering Lambertian scenes, when the operator is not compact, root properties stay covert, so using Laguerr's method becomes a blind action. Considering the lack of the primary information on the eigenvalues of light transport, the limited theoretical background behind Laguerr's algorithm, and the dependency of our method on different factors, especially the position of all eigenvalues, it is not straightforward to conclude that which parameters exactly cause the behavior of histograms. Nevertheless, we finish the chapter by addressing a few observations from the histograms.

• Similar to the matrix case, discussed in Section 8.2, Laguerr's method is not accurate when the roots of the polynomial are closely spaced. This factor is also one of the items causes failure of detecting such eigenvalues of the "CornellBox.ifs.scaled" scene and the "Contact" scene.

• The accuracy of the method does not depend only on the magnitude of the first value. For example, the first eigenvalue of the "CornellBox.ifs.scaled" scene and the "Cubes0343" scene are similar however, the peaks of the "CornelBox.ifs.scaled" scene's histogram coincide the reference ones much better than the ones of the "Cubes0343" scene.

Chapter 9

Conclusion

In the simulation of light transport, it is essential to reduce the computation cost, which requires information on the dimension of the global light transport operator. The dimensionality might be defined as the dominant eigenvalues of the light transport operator. Even though local light transport has already been studied, the dimensionality of global light transport has not been studied extensively, up to our knowledge. In this project, we studied the light transport operator as an operator acting on an infinite-dimensional space, and proposed numerical methods, using path tracing algorithms, to search for eigenvalues of the operator in Lambertian scenes. Theory. To design such algorithms, spectral information of the light transport operator is required. This necessity originates from the fact that the set of points in which operator T -z does not have an inverse is not necessarily limited to the set of eigenvalues. The literature on the spectral study of bounded operators offers the compactness property to simplify this issue. This property ensures that the operator behaves similarly to a matrix, and the spectrum of the operator becomes limited to a decreasing set of eigenvalues with an only accumulation point at 0. In formulating the theories developed in this manuscript, this property showed up either as a necessary or a sufficient condition. Therefore, following Arvo's work [START_REF] Arvo | The role of functional analysis in global illumination[END_REF], while studying the functional properties of this operator (Chapter 4), we studied the compactness property. We proved that for diffuse surfaces, abutting edges of the geometry prevent the operator to be compact.

Methods. Although lack of this property made the task of finding eigenvalues more challenging, we implemented our methods based on the continuous representation of the operator.

Using a Monte Carlo formulation of the power iteration method, in Chapter 6, we estimated the dominant eigenvalue using images of the powers of the operator.

In Chapter 5, we discussed the resolvent of the light transport operator by passing the theory from finite-dimensional vector spaces to infinite-dimensional Hilbert spaces. We showed how the resolvent theory over infinite-dimensional 102 spaces connects the integral of the resolvent to the eigenspaces of corresponding isolated eigenvalues. If the set of eigenfunctions of the light transport operator is complete, we can project an arbitrary distribution into the eigenspaces of the isolated eigenvalue by a line integral of the resolvent (at least in theory). Therefore, proving this property (completeness of the set of eigenfunctions) is beneficial. We also described the resolvent of light transport in terms of the power series of this operator. Speaking of power series, to extend the representation of the resolvent to wider regions of the complex plane, we also discussed the analytic continuation of the resolvent. However, performing analytic continuation on the resolvent of light transport demands spectral information from the operator. Namely, to be able to apply the powers of the operator, it is necessary to know the location of the largest spectrum element in magnitude, belonging to the unit disk. In Lambertian scenes, when the operator is compact, searching for this element is reduced to the set of eigenvalues on the real line (since eigenvalues are real in this case). It could be possible to use estimates of the dominant eigenvalues, obtained in Chapter 6. But calculating coefficients that represent the resolvent in terms of power series suffer from numerical instability, due to sequential calculations of the coefficients from one region (circle) to another in the complex plane.

The generalization of the Fredholm determinant to a compact operator can be applied to estimate the eigenvalues of the operator by computing the poles of this determinant. However, the extension of this theory has developed for specific classes of compact operators, namely, Hilbert-Schmidt and trace class. We know that the light transport operator is not compact in general. However, in Lambertian scenes, it becomes Hilbert-Schmidt by removing the edges of the geometry. Therefore, the Fredholm theory provides a formula to compute the coefficients participating in the Fredholm determinant using kernel evaluations at pairs of points. Among our test scenes, there are three in which the operator is compact ("TwoSquares", "Cubes027", and "Cubes0343"), but the method only works for the "TwoSquares" scene. This failure may originate from intrinsic properties of the two other scenes which require multiple visibility calculations while computing kernel evaluations. Among the operators related to the rendering equation, K x is compact, therefore it is possible to solve the Fredholm equation of this integral to obtain its eigenvalues. On the other hand, the eigenpair of this operator is related to singular value decomposition (SVD) of the light transport operator1 . Therefore, extracting eigenvalues of the local reflectance operator is helpful to represent the SVD expansion of T .

In chapter 8, we estimated the characteristic polynomial of a finite rank operator (which is an approximation of the light transport operator) using a matrix filled by the estimates of powers of the operator at measurement points. The background theory of this method is developed under the assumption of the compactness property, which is a sufficient condition to have isolated eigenvalues with a decreasing order. However, the method still works for the "CornellBox" scene, in which the operator is not compact. Our intuition is that even though we do not know about the rank of global light transport and the operator is not compact, but this operator has a countable set of isolated and dominant eigenvalues allowing computing the eigenvalues of this scene. The performance of this method might be improved using the Metropolis algorithm. However, in this way, one needs to scale each column with a factor. So, it is required to compute the factor for each consecutive column of the matrix.

In developing the supportive theory of the methods experimented on in our test scenes, the compactness property helped us to simplify the problems. However, without this sufficient condition, some of our methods still work. This observation suggests the existence of special spectral properties which are always fulfilled with the compactness property. These properties might be the countability and the isolation of eigenvalues of the light transport operator.

1

  Notations used in this manuscript. . . . . . . . . . . . . . . . . . 3.1 Bidirectional Reflection Distribution Function (BRDF). . . . . . 3.2 Expressing the incident radiance at point x in terms of radiance at point y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The projected area measure is defined by dS ⊥ x = cos θ dS x . . . . . 3.4 Reparametrization operator: representing the incoming radiance at point x in terms of the outgoing radiance from the nearest point. 4.1 Notations used in Proof 4.3.1. . . . . . . . . . . . . . . . . . . . . 4.2 T b is not compact on a scene with abutting edges. . . . . . . . . 4.3 Illustration of Theorem 4.4.5: the radiant exitance operator T b is not compact. Left: rendered scene with diffuse emitter S n . Center: the radiant exitance computed on the receiver plane for n = 7. Right: table contains the pairwise L 2 distances between T b B n and T b B n+p . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Continuity of the kernel at x = y. . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Illustration of the method's performance on two pre-constructed matrices with decreasing eigenvalues. . . . . . . . . . . . . . . . . 8.2 The distance of consecutive eigenvalues of two pre-constructed matrices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Images of T k L of the "CornellBox" scene rendered by a path tracer. 8.4 An example of a light path that contributes to the image of T 3 L. 8.5 Building paths of length k from the generated paths by the Russian roulette technique. . . . . . . . . . . . . . . . . . . . . . . . 8.6 Performing the method on test scenes. . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Comparison of the estimated eigenvalues of matrix A from the histogram with the reference values. . . . . . . . . . . . . . . . . 8.2 Comparison of the estimated eigenvalues of matrix B from the histogram with the reference values. . . . . . . . . . . . . . . . . 8.3 Biased approximation of the eigenvalues of the "CornellBox" scene using the images of T k L ( k = 1, • • • , 5). . . . . . . . . . . . . . . 8.4 The results of the method in object-space. . . . . . . . . . . . . .

H

  Dot product weighted by cosine (Eq. (3.30))∥ • ∥ H Norm induced by ⟨ , ⟩ H L : ℓ → R Light distribution function. L ∈ H K : H → H Global reflectance operator K x : L(Ω x ) → L(Ω x )Local reflectance operator at x G(x, y) Geometric term between points x and y G : H → H reparametrization operator T : H → H Light transport operator (T = KG) T b : B → B Radiant exitance transport operator ρ d (x)

Figure 3

 3 Figure 3.1: Bidirectional Reflection Distribution Function (BRDF).

Figure 3 . 2 :

 32 Figure 3.2: Expressing the incident radiance at point x in terms of radiance at point y.

Figure 3 . 3 :

 33 Figure 3.3: The projected area measure is defined by dS ⊥ x = cos θ dS x .

Figure 3 . 4 :

 34 Figure 3.4: Reparametrization operator: representing the incoming radiance at point x in terms of the outgoing radiance from the nearest point.

  .33) Equation (3.33) can be reformulated in terms of point surfaces instead of solid angle i.e. (T L)(x, ω o ) = S κ(x, y)L(y, ω o )dS y (3.34) with κ(x, y) = ρ(x, ω o , ω i ) cos θ cos θ ′ r 2 xy V (x, y), (3.35) where cos θ ′ r 2 xy is the Jacobian of change of variables from dω i to dS y , and V (•, •) is the visibility function.

  3.Adjoint of the light transport operator (importance transport). The light transport operator is reformulated as a composition of two other operators, i.e. T = KG. By the property, related to the adjoint of the composition of operators in Hilbert spaces[START_REF] Rudin | Functional analysis, mcgrawhill. Inc[END_REF], we have T * = G * K * = GK, which is concluded from the self-adjointness property of K and G. Definition 4.3.3 (Normal operator) The operator T ∈ B(H) is normal if

Figure 4

 4 Figure 4.1: Notations used in Proof 4.3.1.

  Theorem 4.4.4 If the kernel of T b is bounded over S 2 , then the operator is compact.

  Theorem 4.4.5 If κ is not bounded, then T b is not compact.

Figure 4

 4 Figure 4.2: T b is not compact on a scene with abutting edges.

  Figure 4.3.

Figure 4 . 3 :

 43 Figure 4.3: Illustration of Theorem 4.4.5: the radiant exitance operator T b is not compact. Left: rendered scene with diffuse emitter S n . Center: the radiant exitance computed on the receiver plane for n = 7. Right: table contains the pairwise L 2 distances between T b B n and T b B n+p .

  ω o ) was used to re-express the integral of the first line of equation (4.13) in terms of point y, ω i and ω o are the same as in Definition 3.5.5. 4.A.3 K x and K are self-adjoint.

Figure 4 . 4 :

 44 Figure 4.4: Continuity of the kernel at x = y.

  [START_REF] Brown | Complex variables and applications[END_REF] and Kato et al.[START_REF] Kato | Perturbation theory for linear operators[END_REF], and for the resolvent of an operator acting on infinite-dimensional vector spaces the book byCheverry et al. [CR19] are used. Definition 5.1.1 (Open disk) The open disk D(z, r) is referred to the region of the complex plane which is defined as follows D(z 0 , r) = {z ∈ C | |z -z 0 | < r}. Definition 5.1.2 (Singular point) Let f be a complex function defined over D ⊂ C. Point z 0 ∈ D is called singular, if f can not be expressed analytically at z 0 .

  Let T be bounded. Then z ∈ ρ(T ) if and only if T -z is bijective. Therefore, (T -z) -1 is bounded and well-defined for z ∈ ρ(T ). Moreover, sp(T ) ⊂ D(0, ∥T ∥), (5.28) [CR19]. Based on Proposition 5.1.3, the resolvent of a bounded operator over H is defined for z ∈ ρ(T ). Similar to the finite-dimensional case, because the operator (T -z) is invertible, it has the infinite series expression R T (z) = -|z| > ∥T ∥. Remark 5.1.2 The resolvent formula holds for T : H -→ H, for z 1 , z 2 ∈ ρ(T ), [CR19].

  Figure 6.2: Noise-free images of T 16 L and T 17 L of the "CornellBox.ifs" scene, computed using multiple initial paths.

  The histogram of the eigenvalues of a pre-constructed matrix with the reference values (0.8e -i ). The histogram of the eigenvalues of a pre-constructed matrix with the reference values ( 1 2 i ).

Figure 8 . 1 :

 81 Figure 8.1: Illustration of the method's performance on two pre-constructed matrices with decreasing eigenvalues.

  Figure 8.2, d i values for matrices A and B are given.

  Figure 8.3: Images of T k L of the "CornellBox" scene rendered by a path tracer.

Figure 8

 8 Figure 8.4: An example of a light path that contributes to the image of T 3 L.

  a) Cutting the path at x k and connecting it to p i . (b) Connecting the point x k to p i .

Figure 8 . 5 :

 85 Figure 8.5: Building paths of length k from the generated paths by the Russian roulette technique.

for k = 2 ,

 2 Figure 8.6: Performing the method on test scenes.

  (a) "CornellBox.ifs.scaled" Scene (b) "TwoSquaresNoSymmetry" Scene (c) "Contact" Scene (d) "Cubes027" Scene (e) "Cubes0343" Scene (f) "TwoSquares" Scene

  Proposition 5.1.2 Let Γ be a closed curve with positive direction containing isolated eigenvalues λ 1 , • • • , λ k inside. It can be shown Since it is assumed that eigenvalues are isolated, therefore, it is possible to form closed curves those contain each eigenvalue separately, and every eigenvalue lies at the center of the associated circle. Using 5.1.1 for each eigenvalue separately, concludes equation (5.23). It is worth mentioning that the derivation of equation (5.23) does not depend on an extra assumption on T except, having isolated eigenvalues. Considering stronger assumptions on a matrix, such as normality or being diagonalizable allows concluding equation (5.23) simpler, i.e. without using Laurent series.

		i=k
	R T (z)dz = -2πi	P i .	(5.23)
	Γ	i=1
	Proof 5.1.2 Remark 5.1.1 Example 5.1.2 Let T be diagonalizable. Therefore
	i=k
	T =
	i=1

Table 5

 5

		"TwoSquares" "TwoSquaresNoSymmetry"	"Contact"
	κ	Continuous	Continuous		Not continuous
		Bounded	Bounded		Unbounded
	Scene
		"Cubes0027"	"Cubes0343"	"CornellBox"
	κ	Not Continuous Not Continuous Not continuous
		Bounded	Bounded		Unbounded

.1: Test scenes: For each scene the albedo is 0.8.

Table 5 .

 5

2: Test scenes: For each scene the albedo is 0.8, except the Cornell box scene, which has albedos of 0.15, 0.295, 0.4 and 0.5.

Table 7 .

 7 1: Computing c n coefficients and eigenvalues of the light transport operator of the "Contact" scene using both methods. Table7.2: Computing c n coefficients and eigenvalues of the light transport operator of the "TwoSquares" scene using both methods.

		cn		Eigenvalues		Reference Eigenvalues
		Method 1	Method 2	Method 1	Method 2
	n = 0	1.000000000000	1.000000000000	1.2636945069	2.1786475943	0.2979420180
	n = 1	0.000000000000	0.000000000000	1.2636945069	2.1786475943	0.2979950679
	n = 2	-3.325305689225	-9.654203629623	0.2983084146	0.2190373153	0.2873354311
	n = 3	0.000000000000	-0.000000000000	0.2983084146	0.2190373153	0.2873464849
	n = 4	2.717504977530	9.228603131626	0.2983084146	0.1421630612	0.2793885915
	n = 5	0.000000000000	-0.000000000000	0.2983084146	0.1421630589	0.2793980238
	n = 6	-9.536002701504	-6.758301880979	0.0759941015	0.1421630601	0.2727184421
	n = 7	-0.000000000001	-0.000000000000	0.0759941015	0.1421630597	0.2727202517
	n = 8	2.944601263605	3.750401401166	0.0000000000	0.0000000000	0.2663676942
		cn		Eigenvalues		Reference Eigenvalues
		Method 1	Method 2	Method 1	Method 2
	n = 0	1.000000000000	1.000000000000	0.1621222881	0.1621751199	0.1620320938
	n = 1	0.000000000000	0.000000000000	0.1621222881	0.1621751199	0.1620320938
	n = 2	-0.057301200393	-0.057339884526	0.0334121462	0.0335610978	0.0333833732
	n = 3	-0.000000000000	-0.000000000000	0.0334121460	0.0331730098	0.0333816687
	n = 4	0.001530773040	0.001532971797	0.0334121461	0.0335610912	0.0333816687
	n = 5	0.000000000000	-0.000000000000	0.0334121461	0.0331730163	0.0333833732
	n = 6	-0.000029844986	-0.000029730419	0.0092171512	0.0088483792	0.0089627591
	n = 7	-0.000000000000	-0.000000000000	0.0092171512	0.0088483792	0.0089627591
	n = 8	0.000000203247	0.000000202278	0.0053963737	0.0057280355	0.0055606790

Table 7 .

 7 3: Computing c n coefficients and eigenvalues of the light transport operator of the "Cubes0027" scene using both methods. Table7.4: Computing c n coefficients and eigenvalues of the light transport operator of the "Cubes0343" scene using both methods. Table7.5: Computing c n coefficients and eigenvalues of the light transport operator of the "TwoSquaresNoSymmetry" scene using both methods. Table7.6: Computing c n coefficients and eigenvalues of the light transport operator of the "CornellBox.ifs.scaled" scene using both methods. A indicates the scene area, and P denotes the probability density function to sample cyclic light paths. Similar to the c n coefficients, in Method 1, we use directly kernel values at pairs of points to compute p n values, in Method 2, the estimation of these values are used in the calculation. Note that due to the relationship between c n and p n coefficients by equation (7.15), it is also possible to compute c n from p n . For all our test scenes, the column titled as "c n from p n " contains these values. Table7.7: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Contact" scene using both methods.Table7.8: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "TwoSquares" scene using both methods. Table7.9: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Cubes0027" scene using both methods.

		cn		Eigenvalues		Reference Eigenvalues
		Method 1	Method 2	Method 1	Method 2
	n = 0	1.000000000000	1.000000000000	0.3511334245	0.2775410936	0.1565540933
	n = 1	0.000000000000	0.000000000000	0.3511334244	0.2459133874	0.0994781132
	n = 2	-0.276343122240	-0.278730786957	0.3475862150	0.2545465924	0.0953148799
	n = 3	0.007329412940	0.007540001169	0.3475862150	0.2545465904	0.0926948370
	n = 4	0.228419718834	0.222793371663	0.2820318993	0.2334714561	0.0915003017
	n = 5	-0.017332102881	-0.018117951668	0.2820318993	0.2334714541	0.0895453232
	n = 6	-0.263808779790	-0.272502148583	0.2826110841	0.1578578081	0.0879342769
	n = 7	0.048396818600	0.052066549738	0.2826110846	0.1578578084	0.0737271983
	n = 8	3.825454796271	0.334991712935	0.0611937985	0.1878533838	0.0722159153
		cn		Eigenvalues		Reference Eigenvalues
		Method 1	Method 2	Method 1	Method 2
	n = 0	1.000000000000	1.000000000000	1.0363432047	1.2415905679	0.3505501843
	n = 1	0.000000000000	0.000000000000	1.0363431855	1.0166310190	0.2476734734
	n = 2	-5.205916661381	-5.182023579929	0.9019082987	1.0166310054	0.2469152473
	n = 3	0.295064978695	0.306798943340	0.9019082987	1.0912010202	0.2465260685
	n = 4	80.258544188411	78.925546806727	0.8764974121	1.0912010215	0.1972150667
	n = 5	-15.308959730771	-12.272793542900	0.8764974080	1.0264440391	0.1959398298
	n = 6	-1634.094540949442	-2216.335405989682	0.9897568270	1.0264440279	0.1943451506
	n = 7	834.139205163719	1857.323790146492	0.9897568238	0.7134296635	0.1768256458
	n = 8	26463.081001032875	45867.097969776034	0.0652208256	0.0540458117	0.1752493375
		cn		Eigenvalues		Reference Eigenvalues
		Method 1	Method 2	Method 1	Method 2
	n = 0	1.000000000000	1.000000000000	0.1244746359	0.1245973535	0.1244895625
	n = 1	0.000000000000	0.000000000000	0.1244746359	0.1245973535	0.1244895625
	n = 2	-0.033250376741	-0.033300067914	0.0271942234	0.0270582658	0.0271874040
	n = 3	-0.000000000000	-0.000000000000	0.0271942234	0.0270582658	0.0271874040
	n = 4	0.000428066864	0.000426715049	0.0183179597	0.0184898859	0.0184309427
	n = 5	-0.000000000000	0.000000000000	0.0183179597	0.0184898859	0.0184309427
	n = 6	-0.000003452045	-0.000003435042	0.0074953708	0.0049126138	0.0050972532
	n = 7	-0.000000000000	-0.000000000000	0.0074953708	0.0049126137	0.0050972532
	n = 8	0.000000008709	0.000000008619	0.0000000000	0.0046015566	0.0044293806

Table 7 .

 7 10: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "Cubes0343" scene using both methods.

		cn from pn		pn		Eigenvalues
		Method 1	Method 2	Method 1	Method 2	Method 1	Method 2
	n =	1.000000000000	1.000000000000	1.000000000000	1.000000000000	1.1000928389	1.1631356201
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	1.0756078855	1.1631355857
	n =	-5.210606089970	-5.210790443404	5.210606089970	5.210790443404	1.0545023393	1.1447176044
	n =	0.302767329542	0.302594756490	0.151383664771	0.151297378245	1.0545023393	1.1447175722
	n =	81.067671352881	81.073486929993	0.063929353603	0.063920700870	1.0351373416	1.0188967352
	n =	-15.535566211644	-15.526468744297	0.010018612496	0.010046246214	1.0351373416	1.0188967346
	n =	-2091.539631230130	-2091.770349300102	0.003203685718	0.003215610492	0.9023724343	1.0052980202
	n =	833.393159904196	832.895413020379	0.000890292848	0.000895596109	0.9023724240	1.0052980005
	n =	75143.352967791128	75154.797648109001	0.000290081501	0.000293446770	0.9102283847	0.0753217116

Table 7 .

 7 11: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "TwoSquares-NoSymmetry" scene using both methods.

		cn from pn		pn		Eigenvalues
		Method 1	Method 2	Method 1	Method 2	Method 1	Method 2
	n =	1.000000000000	1.000000000000	1.000000000000	1.000000000000	0.1198089324	0.1127072930
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	0.1198089321	0.1127072926
	n =	-0.033270522719	-0.033271666352	0.033270522719	0.033271666352	0.1198089324	0.1127072965
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	0.1198089323	0.1127072925
	n =	0.000429321568	0.000430682084	0.000481910246	0.000481721544	0.0978065694	0.1051197717
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	0.0978065692	0.1051197713
	n =	-0.000002882675	-0.000011735242	0.000007445569	0.000007514570	0.0978065693	0.1053575006
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	0.0978065694	0.1053574955
	n =	0.000036624733	-0.000018131278	0.000000108105	0.000000120506	0.0931735515	0.1053574955

Table 7 .

 7 12: Computing p n , c n coefficients, and eigenvalues of the light transport operator of the "CornellBox.ifs.scaled" scene using both methods.

		cn from pn		pn		Eigenvalues
		Method 1	Method 2	Method 1	Method 2	Method 1	Method 2
	n =	1.000000000000	1.000000000000	1.000000000000	1.000000000000	2.1713602868	2.1660484073
	n =	0.000000000000	0.000000000000	0.000000000000	0.000000000000	2.1684989334	2.1631873046
	n =	-20.582239380717	-20.481654891761	20.582239380717	20.416548917761	2.0847585063	2.0796579898
	n =	0.178850361615	0.177936933362	0.089425180807	0.088968466681	2.0847585459	2.0796580294
	n =	1268.830874347934	1256.457537196904	0.342476571244	0.339504019859	2.0821364073	2.0770359614
	n =	-36.589350752259	-36.226214537492	0.009252450356	0.009092253285	2.0821363673	2.0770359214
	n =	-130154.572883096058	-128255.009634325063	0.007707471209	0.007451536397	1.8081862240	1.8037606255
	n =	7847.178572910076	7731.666955233658	0.000803148431	0.000787517193	1.8081861617	1.8037605634
	n =	18661115.733795121312	18298792.824101842940	0.000319090656	0.000295352848	1.8062892201	1.8018632427

  .1. Note that the distances of consecutive eigenvalues of A become Figure 8.2: The distance of consecutive eigenvalues of two pre-constructed matrices.

	A		B
	d 1	0.5056	d 1	0.5
	d 2	0.1860	d 2	0.25
	d 3	0.0684	d 3	0.125
	d 4	0.0251	d 4	0.0625
	d 5	0.0092	d 5	0.03125
	d 6	0.0034	d 6	0.0156
	(a) di values of matrix A.	(b) di values of matrix B.

Table 8 .

 8 1: Comparison of the estimated eigenvalues of matrix A from the histogram with the reference values.

				A
	Estimated eigenvalue	0.8	0.294		0.108	0.038	0.012
	Reference	0.8	0.294304	0.108268	0.0398297 0.0146525
	eigenvalue
	Table 8.2: Comparison of the estimated eigenvalues of matrix B from the his-
	togram with the reference values.
				B
	Estimated	0.998	0.498	0.248	0.124	0.06	0.03	0.014
	eigenvalue
	Reference	1	0.5	0.25	0.125	0.0625 0.03125 0.015625
	eigenvalue

Table 8 .

 8 3: Biased approximation of the eigenvalues of the "CornellBox" scene using the images ofT k L ( k = 1, • • • , 5).

	Roots

  .4, the value of M , n, and N used for each scene, and the correspondent histogram are addressed.Algorithm 1 pointSourcePathTrace functionx ← growPath(x s ) cumulRef ← 1 for each vertex x t ∈ {x 1 , • • • , x h } do if V (x t , p i ) ̸ = 0 then compute ρ(x t ) W ← p i -x t localize W at x t localize -W at p i compute ρ(x t-1 → x t → p i ) compute ρ(x t → p i → -) * ρ(x t-1 → x t → p i ) * G(x t , p i ) end if cumulRef ← cumulRef * ρ(x t → p i → -) end forAlgorithm 2 Polynomial method with path tracing for each x s do generate random measurement points (p i ) 1≤i≤n generate random directions for the measurement pointsJ ′ ← zero matrix for g = 1, • • • , N do x ← growPath(x s )for each p i do sample a direction for p i J ′ i• ← pointSourcePathTrace (x s , p i ) add J ′ i• to the i-th row of J ′

	J ′ i(t+1) ← cumulRef
	if ∥cumulRef ∥ ̸ = 0 then
	J ′ i(t+1) ←	J ′ i(t+1) ∥cumulRef ∥
	end if
	if k = h then J ′ J ′ i(t+1) i(t+1) ← 1-ρ(xt)
	else J ′ i(t+1) ←	J ′ i(t+1) ρ(xt)
	end if J ′ i(t+1) ←	J ′ i(t+1) 2
	J ′ i(t+1) ← J ′ i(t+1)

  Table 8.4: The results of the method in object-space.

		Results on Different Scenes
	Scene	"CornellBox.ifs.scaled"	"TwoSquares"	"TwoSquaresNoSymme-
				try"
	M	10 5	10 5	10 5
	n	10 3	200	10 3
	N	10 5	10 5	10 5
	Figure	8.6a	8.6f	8.6b
		Results on Different Scenes
	Scene	"Contact"	"Cubes0027"	"Cubes0343"
	M	10 5	10 5	10 5
	n	10 3	500	10 3
	N	10 5	10 5	10 5
	Figure	8.6c	8.6d	8.6e

It should be emphasized here that the eigenpair of the global light transport is different from the eigenpair of the local light transport, which is obtained from matrices provided by discretization.

However, the radiosity matrix in this paper is I -T .

A Banach space is a complete normed vector space.

This field includes theoretical and practical algorithms to calculate Fredholm determinant.

These slices are similar (clustered) image pixels, determined based on their geometry.

Interpreted as the projected differential area of Ωx onto the base of it.

In photometry it is called luminosity.

Note that our study does not include transparent surfaces. In that case, one may use the Bidirectional Scattering Distribution Function (BSDF).

Fresnel equations are used to calculate the reflectivity and transmissivity of light when it reaches the interface between two surfaces. This reflection is called Fresnel reflection.

A diffuse emitter emits equal light in all directions and points.

These patches are often triangles or convex quadrilaterals.

Bias does not appear as noise in the image. It is often detectable in the form of blurring or incorrect surface shading.

A σ-algebra on a set S is a collection of S's subsets which is closed under countable unions, countable intersections and complement.

See Treves' book [Tre16] for a proof.

Since this example is a very known case to explain analytic continuation, to avoid repetition, the details of computing an coefficients are skipped, for more information on the subject check Ahlfors et al.'s book[START_REF] Ahlfors | Complex analysis mcgraw-hill[END_REF].

A map is called bijective if it is both injective and surjective.

Note that ρ(T ) is not the same as ρ d (x) which is the albedo at surface point x.

However, eigenvalues of light transport matrices, provided by descritization, are always decreasing.

This algorithm is also known as Von Mises iteration.

x indicates the light path, and µ x is the corresponding measure defined on the space of light paths.

Considering two arbitrary light paths, ergodicity means that one light path can be reached by performing a finite set of mutations on the another one.

When a Markov chain reaches an equilibrium distribution, or a stationary distribution then, it will maintain its behavior for future times.

This function returns the probability density of y = X j given x = X j-1 .

Note that x k in both equations (6.12) and (6.13) represents the last vertex of X j+1 and X j respectively. However, these vertices may not be the same point.

This norm is defined as the sum of singular values of the operator.

Note that this series starts with index 2, since the trace of T n for n ≥ 2 is finite i.e. all powers of T n for n ≥ 2 are trace class operators.

By cyclic paths, we mean the light paths that have the same endpoints. For example, (x, x 2 , • • • , xn, x) is a cyclic light path. Since the sequence of points starting from x returns to x, as the final vertex.

This theorem states that evaluating the characteristic polynomial of a matrix in the matrix itself gives the zero matrix.

t indicates the transpose.

However, note that to provide the reference eigenvalues for our test scenes, we could not skip discretization.

By the word suitable we mean the path whose the last vertex of it is visible from camera position.

It is worth mentioning that rendering full images even with a static light source using path tracing is very costly.

Laguerre's method is a fairly simple algorithm to compute the roots of a polynomial. Using this algorithm, it is almost guaranteed to converge to some root of the polynomial[START_REF] Forman S Acton | Numerical methods that work[END_REF].

Newton's method is an appropriate method to extract real isolated roots. However, this method is very likely to miss multiple roots.

Since T is not compact, the existence of the SVD for light transport is not guaranteed. But, we proved, under a limited set of configurations, that this expansion exists.

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisor Cyril Soler for the continuous support of my Ph.D study, for his patience, motivation, and immense knowledge. His guidance helped me in all the time of research and writing of this thesis. I could not have imagined having a better advisor and mentor. Besides, I would like to thank my thesis committee: Mr Bruno Lévy, Derek Nowrouzezahrai, Valérie Perrier, Mathias Paulin, and Christian Lessig, for their insightful comments and encouragement. I would also like to thank my family for their continuous support during the time.

Part II

Computing eigenvalues using

Monte Carlo algorithms

Information on test scenes

We implemented and tested all methods described in the manuscript on various Lambertian test scenes. These scenes are extremely simple in order to compare the results of each method with the reference values.

Reference matrix calculation uses an eigensolver. Each scene is first tessellated. Then, surface-to-surface form-factor matrix is computed using Monte Carlo integration. Finally, the reference eigenvalues are computed using Arpack++.

Below is a small description for each scene.

• "TwoSquares": This is the most basic scene possible: only two squares facing each other. The scene only contains two polygons, the eigenvalues come by pairs of opposite signs. On this scene, the light transport kernel is continuous over its entire domain, see Table 5.1.

• "TwoSquaresNoSymmetry": Same scene than above, but with a geometric setup that prevents any symmetry. Here again eigenvalues come by pair of opposite signs. The kernel is here again continuous, see Table 5.1.

• "Contact": Elementary scene with a discontinuous light transport kernel, see Table 5.1.

• "Cubes0027": 27 cubes stacked regularly, but with a random orient. The light transport kernel is continuous, see Table 5.2.

• "Cubes0343": 343 cubes stacked regularly (7×7×7), with a random orient. The light transport kernel is continuous, see Table 5.2.

• "CornellBox": A classical Cornell Box that has been re-meshed in order to avoid contact surfaces. Contact surfaces indeed results in the kernel to be null on finite areas. Albedo is 0.295 (resp. 0.15) on the right face (resp. left face), and 0.4 everywhere else, see Table 5.2. 63