<figure><img src="image_1.png" /></figure>

HAL
open science

# Spectral Theory of Light Transport Operators

Cyril Soler, Kartic Subr

► To cite this version:

Cyril Soler, Kartic Subr. Spectral Theory of Light Transport Operators. ACM Transactions on Graphics, 2025, 45 (2), 10.1145/3774756. hal-05407513

HAL Id: hal-05407513

https://inria.hal.science/hal-05407513v1

Submitted on 10 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spectral Theory of Light Transport Operators

CYRIL SOLER, INRIA - Grenoble University, France

KARTIC SUBR, University of Edinburgh, UK

Light Transport Operators (LTOs) represent a fundamental concept in computer graphics, modeling single bounces of light within a virtual environment as linear operators on infinite dimensional spaces. While the LTOs play a crucial role in rendering, prior studies have primarily focused on spectral analyses of the light field rather than the operators themselves. This article presents a rigorous investigation into the spectral properties of the LTOs. Due to their non-compact nature, traditional spectral analysis techniques face challenges in this setting. However, many practical rendering methods effectively employ compact approximations, suggesting that non-compactness is not an absolute barrier. We show the relevance of such approximations and establish various path integral formulations of their spectrum. These findings enhance the theoretical understanding of light transport and offer new perspectives for improving rendering efficiency and accuracy.

CCS Concepts: • **Computing methodologies** → **Rendering**; • **Mathematics of computing** → **Functional analysis**; **Integral equations**.

Additional Key Words and Phrases: Light Transport Operator, Compactness, Fredholm Equations

**ACM Reference Format:**

Cyril Soler and Kartic Subr. 2026. Spectral Theory of Light Transport Operators. *ACM Trans. Graph.* 45, 2, Article 112.2092 (April 2026), 20 pages. https://doi.org/10.1145/3774756

## 1 Introduction

Light Transport Operators (LTOs) are a mathematical construct in computer graphics, that formalizes light propagation within a virtual environment as linear operators acting on a given light distribution. These operators encapsulate one bounce of the input light distribution capturing its intricate interactions with geometric structure and material reflectance properties in the environment. During rendering, the light energy at each pixel of a virtual camera is estimated by evaluating the sum of an infinite series of increasing powers (multiple bounces) of the operators applied to the distribution of emitted light.

Despite the foundational role of transport operators in rendering, most prior studies have concentrated on spectral analyses of the light field [Durand et al. 2005; Mahajan et al. 2007] rather than the operator that governs its evolution. Consequently, the spectral properties of the LTOs remain underexplored. In this work, we use the term *spectral* strictly to

refer to harmonic analysis of the LTOs and not to wavelength-dependent phenomena. Furthermore, our analysis is restricted to monochromatic light (in the absence of fluorescence, light transport operators are defined independently across wavelengths). In this article, we define a *scene* to mean a combination of geometry and their material properties—light sources are excluded since they are considered part of the input space on which the LTOs operate (See Figure 1).

The LTOs inherently act in infinite-dimensional function spaces, posing substantial challenges for both theoretical investigation and computational implementation. Traditional approaches often involve discretization to make the problem tractable. Recent findings by Soler et al. [2022] reveal that LTOs are generally non-compact, leading to non-uniform convergence with respect to the number of elements of the discretization. Non-compactness complicates the application of classical spectral analysis results from finite-dimensional settings. Nevertheless, practical techniques such as Galerkin approximations [Baranoski et al. 1997] and precomputed radiance transport [Sloan et al. 2002] circumvent these issues by approximating the LTO with compact operators, suggesting that non-compactness is not an absolute limitation. This observation motivates deeper inquiry into the relationships between finite-dimensional approximations and the infinite-dimensional operators.

Our theoretical framework aims to bridge these gaps, connecting the eigenspectrum of the LTOs to measures of specific light path sets within a scene. By establishing these connections, we contribute to both the theoretical understanding of light transport and the development of practical methodologies for spectral analysis. We envision that our findings may have future direct implications for improving the efficiency and accuracy of rendering algorithms. In this article, we present a detailed analysis of the spectral properties of the infinite-dimensional LTOs, with the following key contributions:

• A novel and rigorous mathematical characterization of the eigenspectrum of the LTOs;

• an investigation into how compact approximations of light transport operators can be used to approximate their spectral properties; and

• identification of connections between the eigenvalues of LTO's and path space integration, setting foundations for fresh perspectives on using eigenanalysis in rendering.

This article presents an analysis of the light transport operator rather than of the light distribution. The operator is defined specifically by structure within a scene and so the specific choice of emitters within the scene is largely irrelevant. Occasionally, we exemplify our discussion with specific

Authors' Contact Information: Cyril Soler, INRIA - Grenoble University, Montbonnot, France, cyril.soler@inria.fr; Kartic Subr, University of Edinburgh, Edinburgh, UK, k.subr@ed.ac.uk.

<figure><img src="image_3.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License.

© 2025 Copyright held by the owner/author(s). ACM 1557-7368/2026/04-ART112.2092 https://doi.org/10.1145/3774756
scenes that we chose to illustrate (only showing geometry), rather than render, to avoid misleading readers.

We separate the treatment of related work (Section 2) and the mathematical background needed (Section 3) to explain our work. Before performing a spectral analysis of operators that are non-compact in general, we explain how they permit compact approximations which lend themselves to analysis (Section 4). The core results surrounding spectral analysis of LTO's are presented as theorems (Section 5). We further explore how these spectra relate to different aspects of light transport (Section 6). Next, we briefly describe practical methods for computing the eigenvalues of the LTOs, providing insights into the computational challenges and opportunities in spectral analysis (Section 7). We conclude with a discussion of the implications of our findings and outline future research directions (Section 8).

## 2 Previous work

This article pertains to the spectral analysis of linear operators over infinite dimensional spaces. Although the operators that we analyze have a long history in computer graphics, there has been little theoretical analysis of spectral properties. Yet, several practical approaches have exploited their low-dimensional approximations.

**Light transport operators and their spectrum.** We analyze operators that are minorily adapted from standard formulations [Arvo 1996; Veach 1997]. These formulations have provided insight to various applications such as inverting transport equations to compensate for global illumination [Ng et al. 2012], unification of forward and inverse rendering [Bai et al. 2010], analysis of compactness [Soler et al. 2022], and so on. Their spectrum however has received little attention.

Baranoski [1997] examined the spectrum of the radiosity matrix in the discrete Lambertian case. Ashdown [2001] reformulated the radiosity transport matrix as an electrical network and studied its eigenvectors via analysis of the conductance matrix of the resulting network. Machida [2014] showed that, in an infinite and isotropic medium, the singular eigenfunctions of the radiative transport equations can be independently solved along each dimension.

**Spectral analysis of linear operators.** Although this is a vast area of study within mathematics, light transport operators exhibit challenging characteristics–such as having asymmetric unbounded and discontinuous kernels– restricting the applicability of standard results from the theory of linear operators. Despite this, we draw from a range of theoretical tools and concepts including trace-class operators [Brislawn 1988], non self-adjoint operators [Gohberg and Krein 1978], perturbation theory for linear operators [Kato 1995] and Fredholm equations [Zemyan 2012]. The connections between the eigenvalues of an operator and its discretized counterpart can be found in the work of Boffi [2010], Sun [2016], and Chatelin [2011]. We provide specific citations to the relevant sections of these works in our derivations.

**Dimensionality of light transport matrices.** Discrete light transport operators have been observed empirically to be conducive to dimensionality reduction. Precomputed Radiance Transfer methods [Nowrouzezahrai et al. 2007; Sloan et al. 2002, 2005; Wang et al. 2007] and modular radiance transfer [Loos et al. 2011] utilized this property to compress transport matrices. Wang [2009] exploited this property to reconstruct the matrix from a limited set of images. The solution space of light distributions obtained via low-dimensional approximations of the transport matrix is limited by the low rank of the matrix. Some methods attempt to adjust spatial sampling to optimize for this [Belcour et al. 2022; Huang and Ramamoorthi 2010] or use low rank approximations to parts of the transport matrix [Garg et al. 2006]. Dimensionality reduction has also been applied to simplify computation in the ‘many-light’ setting with several point lights [Hašan et al. 2007; Ou and Pellacini 2011]. Mahajan et al. [2007] presented a theory of locally low dimensionnal light transport as the number of principal components required to represent the output of the transport matrix under varying conditions. This is analogous to the problem of estimating the number of eigenvalues of the light transport matrix necessary for local approximations. Lessig [2010] relates the problem of finding the effective dimension of local bandlimited and radially symmetric light transport to the spatio-spectral concentration problem. While all the above works operate on discretizations of the operator as matrices, the non-compactness of these operators [Soler et al. 2022] precludes finite dimensional uniform approximations whose spectra are not representative of the original operator's spectrum. In this article, we analyze the spectrum of the light transport operator without discretization.

## 3 Background

### 3.1 Definitions: Linear operators and their properties

We study linear operators over Hilbert spaces. We denote spaces, operators and space elements using upper case caligraphic, upper case bold, and roman letters respectively.

**Norm and convergence of linear operators.** Let $\mathcal{H}$ be a Hilbert space. The norm of a linear operator $\mathbf{A}: \mathcal{H} \to \mathcal{H}$ is defined as $\|\mathbf{A}\| = \sup_{\|f\|=1} \|\mathbf{A}f\|$. The remaining definitions in this section apply to operators for which the supremum exists, or bounded linear operators.

A sequence of linear operators $\{\mathbf{A}_n\}$ is said to be uniformly converging to an operator $\mathbf{A}$ when

$$ \lim_{n \to \infty} \| \mathbf{A} - \mathbf{A}_n \| = 0. \qquad (1) $$

A sequence of linear operators $\{\mathbf{A}_n\}$ is said to be strongly converging to an operator $\mathbf{A}$ when

$$ \forall f \in \mathcal{H} \quad \lim_{n \to \infty} \| \mathbf{A}f - \mathbf{A}_n f \| = 0. \qquad (2) $$

Uniform convergence implies strong convergence [Chatelin 1981].
<figure><img src="image_4.png" /><figcaption>Fig. 1. An illustrated overview of the central research question tackled in this article. (a) The radiance distribution in a scene $L$ (magenta) at equilibrium satisfies the rendering equation $L = E + TL$, where $T$ is a linear operator that transports radiance by one bounce and $L_e$ is the emitted radiance distribution (yellow). (b) The operator $T$ is fully specified by the geometry of the scene along with the reflectance distributions associated with all surfaces in the scene. It is independent of the choice of emitted distribution $L_e$. (c) If the scene (and hence $L$) were discretized then $T$ could be represented as a matrix and its spectrum (e.g. eigenvalues and eigenvectors) could be analyzed empirically. In this article, we present mathematical results about the spectrum of the infinite-dimensional operator $T$ without discretization.</figcaption></figure>

**Symmetrizable operators.** The adjoint of operator $\mathbf{A}: \mathcal{H} \to \mathcal{H}$ is the unique operator $\mathbf{A}^*$ for which

$$ \forall (f,g) \in \mathcal{H}^2 \quad \langle \mathbf{A}f, g \rangle = \langle f, \mathbf{A}^*g \rangle. $$

If $\mathbf{A}^* = \mathbf{A}$, then $\mathbf{A}$ is self-adjoint and enjoys special properties (See Lax [Lax 2002], chap.28). If there exists a non-negative self-adjoint operator $\mathbf{S}$ such that $\mathbf{SA}$ is self-adjoint, then $\mathbf{A}$ is known as *symmetrizable*. In this situation, $\mathbf{A}$ is self-adjoint for the dot product

$$ \langle f, g \rangle_S = \langle f, \mathbf{S}g \rangle. \quad (3) $$

Symmetrizable operators in separable Hilbert spaces share many properties with self-adjoint operators.

The definition of $\langle \cdot, \cdot \rangle_S$ in Equation 3 requires $\mathbf{S}$ to be bounded. Alternatively, if $\mathbf{A}$ can be symmetrized by a positive self-adjoint operator $\mathbf{S}'$ with a bounded inverse on $\mathcal{H}' \subset \mathcal{H}$ (See Silbertein [1962] part II, Th. 9.1) then $\mathbf{A}$ is a self-adjoint operator on $\mathcal{H}'$ with dot product $\langle f, g \rangle_{S'} = \langle \sqrt{\mathbf{S}'}f, \sqrt{\mathbf{S}'}g \rangle$, so it still shares the properties of symmetrizable operators in this subspace.

**Compact operators.** A compact operator maps any bounded sequence in the domain into a sequence containing a converging subsequence. The family of compact operators itself forms a closed set, meaning that the limit of a uniformly converging sequence of compact operators is always a compact operator ([Gohberg and Krein 1978] p67). Practically speaking, compact operators behave like “infinite dimensional matrices”. Light transport operators are generally not compact [Soler et al. 2022].

**Trace-class and Hilbert–Schmidt operators.** An orthogonal basis $\{e_i\}_{i \ge 0}$ of $\mathcal{H}$ that generalizes to infinite-dimensional spaces is known as a Schauder orthogonal basis: every element in $\mathcal{H}$ can be expressed as a–possibly infinite–linear combination of $e_i$. Such a basis can be used to define the trace and the Hilbert–Schmidt norm

$$ \mathrm{Tr}(\mathbf{A}) = \sum_{i=0}^{\infty} \langle \mathbf{A}e_i, e_i \rangle \quad \text{and} \quad \|\mathbf{A}\|_{\mathrm{HS}} = \sum_{i=0}^{\infty} \|\mathbf{A}e_i\|^2 \quad (4) $$

respectively. In finite-dimensional spaces, these correspond to the trace of the corresponding matrix and its Frobenius

norm. Trace-class operators are those for which Tr($\mathbf{A}$) converges (See Gohberg [1978] §8)<fnref n="1" /> and Hilbert–Schmidt operators are operators for which the Hilbert–Schmidt norm converges. All trace-class operators are Hilbert–Schmidt and all Hilbert–Schmidt operators are compact (See Simon [Barry 2000], sec.5). The product of two Hilbert–Schmidt operators is trace-class (following Gohberg [1978] Cor. 4.1 and §4.3).

**Integral operator.** If $\mathcal{H}$ is a space of integrable functions over a domain $S$, an operator $\mathbf{A}: \mathcal{H} \to \mathcal{H}$ is an integral operator if

$$ \forall f \in \mathcal{H} \quad (\mathbf{A}f)(\mathbf{x}) = \int_S \kappa(\mathbf{x}, \mathbf{y})f(\mathbf{y})d\mathbf{y}, $$

where $\kappa: S \times S \to \mathbb{C}$, is known as the integration kernel. A necessary and sufficient condition for an integral operator to be Hilbert–Schmidt is that $\int_S \int_S \kappa(\mathbf{x}, \mathbf{y})^2 d\mathbf{x}d\mathbf{y} < \infty$ (Gohberg [2012] Chp.IX), which is automatically verified when $\kappa$ is bounded over its domain. The trace of trace-class integral operators can be computed (Barry [2000] Th. 3.9)

$$ \mathrm{Tr}(\mathbf{A}) = \int_S \kappa(\mathbf{x}, \mathbf{x})d\mathbf{x}. \quad (5) $$

Convergence of this integral does not, by itself, imply that $\mathbf{A}$ is trace-class.<fnref n="2" />

## 3.2 Spectra of operators

In finite-dimensional spaces, linear operators are matrices. The spectrum of a $n \times n$ matrix $\mathbf{M}$ is the set of exactly $n$ eigenvalues $\lambda_i$ (accounting for multiplicity) such that $\mathbf{M} - \lambda \mathbf{I}$ is singular, and corresponds to the roots of the characteristic polynomial of $\mathbf{M}$. In infinite-dimensional spaces these definitions are slightly more complex:

**Spectrum, resolvent set and resolvent operator.** The spectrum of an operator $\mathbf{A}$ which we denote $\sigma(\mathbf{A})$ is the set of complex scalars $\lambda$ for which $\mathbf{A} - \lambda \mathbf{I}$ is bijective and therefore cannot be inverted. The resolvent set $\rho(\mathbf{A}) = \mathbb{C} \setminus \sigma(\mathbf{A})$ is

<footnote n="1">¹Another definition of the trace involving the polar decomposition of $\mathbf{A}$ is sometimes found.</footnote>

<footnote n="2">²These are generally proved when $\kappa$ is continuous and can be extended to domains where $\kappa$ is continuous almost everywhere.</footnote>
the complement of the spectrum, for which the resolvent operator $R(\mathbf{A}; \lambda) = (\lambda\mathbf{I} - \mathbf{A})^{-1}$ exists. When $\lambda \in \rho(\mathbf{A})$, the resolvent operator is a bounded bijection over $\mathcal{H}$ (See Jeribi [2021] Sec 1.2.2).

**Point spectrum, continuous spectrum and residual spectrum.**

The non-bijectivity of $\mathbf{A} - \lambda\mathbf{I}$ for $\lambda \in \sigma(\mathbf{A})$ reflects different situations: it might be non injective, or injective but not surjective in which case the range of the operator may be dense or not in $\mathcal{H}$. This defines the following partition of the spectrum:

$\sigma_p$: the point-spectrum is the set of $\lambda \in \mathbb{C}$ such that $\mathbf{A} - \lambda\mathbf{I}$ is not injective;

$\sigma_r$: the residual spectrum is the set of $\lambda \in \mathbb{C}$ such that $\mathbf{A} - \lambda\mathbf{I}$ is injective, but does not have a dense range;

$\sigma_c$: the continuous spectrum is the set of $\lambda \in \mathbb{C}$ such that $\mathbf{A} - \lambda\mathbf{I}$ is injective with a dense range, but not bounded below.

The eigenvalues of $\mathbf{A}$ are the complex numbers $\lambda$ such that there exists $\Lambda \in \mathcal{H}$ for which $\mathbf{A}\Lambda = \lambda\mathbf{A}$, which corresponds to the elements of $\sigma_p(\mathbf{A})$. As such, we're interested in this article in the point-spectrum of light transport operators.

**Intuition about operator spectra.** Although the point-spectrum of an operator is a generalization of the eigenvalues of a matrix to infinite-dimensional spaces, the general properties of the point-spectrum of operators cannot be extrapolated from those of matrices [Mortad 2022]. Some classical examples of operators and their point-spectrum include the Fourier transform, which point spectrum is the finite set $\{1, i, -1, -i\}$, and the unidirectional right-shift operator $\mathbf{S}_r$ (resp. left-shift operator $\mathbf{S}_l$) over square-summable sequences in $\mathbb{C}$, for which $\sigma_p(\mathbf{S}_r) = \emptyset$ whereas $\sigma_p(\mathbf{S}_l)$ is the open unit disc [D'Aniello and Maiuriello 2022]. In general, $\sigma_p$ may be either empty, finite, countable infinite, or even uncountable. Even “good looking” integral operators over $L^2[0, 1]$ may have any compact subset of $[0, 1]$ as a point-spectrum [Kalisch 1972]. Like matrices however, the span of eigenfunctions (including those corresponding to a zero eigenvalue) is not guaranteed to be the entire space. Except when standard theorems apply, proofs for every particular operator need a careful derivation.

One advantage of studying compact operators is that their spectrum boils down to a countable point-spectrum, with zero as the only possible accumulation point (Simon [2000] Th. 1.1). As such, they behave like “infinite-dimensional matrices”. This also explains why they are the closure of the set of finite dimensional operators, in operator norm. In this article, we analyse the spectral properties of light transport operators.

**Infinite-dimensional analogue of SVD.** For bounded operators over infinite dimensional Hilbert spaces, the notion of singular value decomposition (SVD) is replaced by a continuous analogue [Crane and Gockenbach 2020], due to the fact that such operators may not have a countable set of singular values. For compact operators however, the singular values are countable (See Gohberg et al [Gohberg and Krein 1978]

Chap.2, Sec.2) and the SVD remains a discrete sum known as the *Schmidt expansion*. For a compact operator $\mathbf{A}$ it takes the form of

$$ \mathbf{A}f = \sum_{k=0}^{\infty} s_k \langle v_k, f \rangle u_k, \quad (6) $$

where $\{u_k\}_{k \ge 0}$ and $\{v_k\}_{k \ge 0}$ are complete orthogonal sequences in the operator value space standing for the left and right singular vectors respectively, and the $s_k$ are the singular values (also named s-numbers) of $\mathbf{A}$.

When an operator is self-adjoint in a specific Hilbert space (self-adjointness depends on the dot product), its s-numbers are equal to its eigenvalues up to the sign, which in turn do not depend on the dot-product.

### 3.3 Definitions: Light transport operators

**Assumptions and spaces considered.** We limit our formulation to monochromatic light transport with surface reflection only. The space of radiance functions is $\mathcal{H} = L^2(S \times \Omega)$, where $S$ is the set of surfaces (with bounded area) in a scene and $\Omega$ is the upper hemisphere of directions. $\mathcal{H}$ is a separable Hilbert space with the inner product [Soler et al. 2022]

$$ \langle l_1, l_2 \rangle_{\mathcal{H}} = \int_S \int_{\Omega} l_1(\mathbf{x}, \omega) \bar{l}_2(\mathbf{x}, \omega) \cos \theta' d\omega' d\mathbf{x}. \quad (7) $$

The cosine in Equation 7 imbues valuable properties to some of the operators defined below. The notation in this inner product is depicted in Figure 3.

**Reflectance operator K.** The reflectance equation [Kajiya 1986] expresses radiance leaving a point $\mathbf{x}$ in direction $\omega$ as:

$$ L(\mathbf{x}, \omega) = \int_{\Omega} \rho(\mathbf{x}, \omega, \omega') L_i(\mathbf{x}, \omega') \cos \theta' d\omega', \quad (8) $$

where $L(\mathbf{x}, \omega)$ is the exitant radiance at $\mathbf{x}$ in direction $\omega$, $L_i$ is incident radiance and $\rho$ is the **bidirectional reflectance distribution function (BRDF)** at $\mathbf{x}$. Figure 3 summarizes the geometric notations. For consistency we consider “mirror” reflections to be using a very sharp true function as opposed to a distribution.

Equation 8 defines an integral operator $\mathbf{K}: \mathcal{H} \to \mathcal{H}$ which maps an incident radiance field into an exitant field: $L = \mathbf{K}L_i$. Since the integral is only over $\Omega$, the reflectance operator $\mathbf{K}$ is a *partial integral operator*. We may express $\mathbf{K}$ as an operator tensor product of $\mathbf{K}_\mathbf{x}: L^2(\Omega) \to L^2(\Omega)$ at $\mathbf{x}$ and the identity $\mathbf{I}: L^2(S) \to L^2(S)$ where for all $\mathbf{x} \in S$, $\mathbf{K}_\mathbf{x}(L_i(\mathbf{x}, .)) = (\mathbf{K}L_i)(\mathbf{x}, .)$. As will be seen later on, $\mathbf{K}_\mathbf{x}$ is a well behaved integral operator over $\mathcal{O} = L^2(\Omega)$.

**Reparameterization operator G.** When radiance is constant along directions in free space (absence of volumetric scattering), the incident radiance field $L_i(\mathbf{x}, \omega')$ is equal to the exitant radiance at some location $\mathbf{y} = p_S(\mathbf{x}, \omega')$ along the direction expressed within that point's own local frame as $\omega'' = p_\Omega(\mathbf{x}, \omega')$. This reparameterization of the radiance field is captured via an operator $\mathbf{G}: \mathcal{H} \to \mathcal{H}$ where

$$ L_i(\mathbf{x}, \omega') = (\mathbf{G}L)(\mathbf{x}, \omega') = L(p_S(\mathbf{x}, \omega'), p_\Omega(\mathbf{x}, \omega')). $$
<table><caption>Fig. 2. Notation used in this article.</caption><tr><td>C</td><td>complex plane</td></tr><tr><td>S</td><td>surfaces in the scene</td></tr><tr><td>Ω</td><td>hemisphere of outgoing directions (local frame)</td></tr><tr><td>L</td><td>4D domain S × Ω</td></tr><tr><td>H</td><td>space of radiance distributions L<sub>2</sub>(L)</td></tr><tr><td>B</td><td>space of spatial distributions L<sub>2</sub>(S)</td></tr><tr><td>O</td><td>space of directional distributions L<sub>2</sub>(Ω)</td></tr><tr><td>⟨., .⟩<sub>H</sub></td><td>dot product weighted by cosine (Eq. 7)</td></tr><tr><td>||<sub>H</sub></td><td>norm induced by ⟨., .⟩<sub>H</sub></td></tr><tr><td>K : H → H</td><td>global reflectance operator</td></tr><tr><td>K<sub>x</sub> : O → O</td><td>local reflectance operator at x</td></tr><tr><td>G : H → H</td><td>re-parameterization operator</td></tr><tr><td>T : H → H</td><td>light transport operator (T = KG)</td></tr><tr><td>T<sub>b</sub> : B → B</td><td>radiant exitance transport operator</td></tr><tr><td>ρ(x, ω, ω')</td><td>BRDF at x in directions (ω, ω') ∈ Ω<sup>2</sup></td></tr><tr><td>ρ̄(x)</td><td>albedo at x (Eq. 11)</td></tr><tr><td>v(x, y)</td><td>binary visibility function between points x and y</td></tr><tr><td>κ<sub>b</sub>(x, y)</td><td>integration kernel of T<sub>b</sub></td></tr><tr><td>κ(x, ω, y, ω')</td><td>integration kernel of T</td></tr><tr><td>p<sub>S</sub>/Ω</td><td>functions mapping (x, ω) to the point and direction seen from x in direction ω (Fig. 3)</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Notation for Equations 7,8,9 and 12. We deviate from convention in the literature [Veach 1997] by always considering outgoing directions in the local coordinate system. Consequently $ω'' \neq -ω'$ in general.</figcaption></figure>

**Transport operator T.** The light transport operator T combines the reparameterization and reflectance operators as follows:

$$ (\mathbf{T}\mathbf{L})(\mathbf{x}, \omega) = \int_{\Omega} \rho(\mathbf{x}, \omega, \omega')(\mathbf{G}\mathbf{L})(\mathbf{x}, \omega') \cos \theta' \, d\omega'. \quad (9) $$

Thus $\mathbf{T}: \mathcal{H} \to \mathcal{H}$ transports a radiance field through one-bounce. Since we consider outgoing directions as the upper hemisphere in the local coordinate system at each point, our definitions of $\mathbf{K}$ and $\mathbf{G}$ differ slightly from those in the literature [Arvo 1996; Veach 1997]. However the operator form of the rendering equation remains as

$$ L = E + \mathbf{T}\mathbf{L} = E + \mathbf{K}\mathbf{G}\mathbf{L}, \quad (10) $$

where $E \in \mathcal{H}$ is the emitted radiance field and $L \in \mathcal{H}$ is the unknown radiance.

Lambertian transport operator $T_b$. When all materials in a scene are Lambertian, the BRDF $\rho$ is reduced to a direction-less albedo $\bar{\rho}(\mathbf{x})$:

$$ \bar{\rho}(\mathbf{x}) = \int_{\Omega} \rho(\mathbf{x}, \omega, \omega') \cos \theta \, d\omega = \pi \rho(\mathbf{x}, \dots, \dots). \quad (11) $$

Also, reflected radiance is independent of directions and radiant exitance $B(\mathbf{x}) = \int_{\Omega} L(\mathbf{x}, \omega) \cos \theta \, d\omega$, is considered instead. Rewriting Equation 9 in this space defines the radiant exitance transport operator $T_b: L^2(S) \to L^2(S)$ [Sillion and Puech 1994] as

$$ (\mathbf{T}_b B)(\mathbf{x}) = \int_S v(\mathbf{x}, \mathbf{y}) \frac{\bar{\rho}(\mathbf{x})}{\pi} \frac{\cos \theta' \cos \theta''}{\|\mathbf{x} - \mathbf{y}\|^2} B(\mathbf{y}) \, d\mathbf{y}, \quad (12) $$

where the integral is over surfaces rather than angle and $v(\mathbf{x}, \mathbf{y})$ is sometimes referred to as the “visibility function”, which restricts the integral to points $\mathbf{y}$ that are visible from $\mathbf{x}$.

Kernel expressions of $\mathbf{T}$ and $T_b$. We define a geometric term

$$ g(\mathbf{x}, \mathbf{y}) = v(\mathbf{x}, \mathbf{y}) \frac{\cos \theta' \cos \theta''}{\|\mathbf{x} - \mathbf{y}\|^2}. $$

From Equation 12 we write the kernel expression of $T_b$ as:

$$ \mathbf{T}_b B = \int_S \kappa_b(\mathbf{x}, \mathbf{y}) B(\mathbf{y}) d\mathbf{y} \quad \text{with} \quad \kappa_b(\mathbf{x}, \mathbf{y}) = \frac{\bar{\rho}(\mathbf{x})}{\pi} g(\mathbf{x}, \mathbf{y}). $$

To express $\mathbf{T}$ as a proper kernel-integral operator, the integral needs to be defined over the entire domain $S \times \Omega$. This requires a Dirac distribution to make the integral over $S$ “select” pairs of elements of $\mathcal{L}$ that can actually exchange light:

$$ (\mathbf{T}\mathbf{L})(\mathbf{x}, \omega) = \int_{S \times \Omega} L(\mathbf{y}, \omega'') \kappa(\mathbf{x}, \omega, \mathbf{y}, \omega'') \mathbf{d}\mathbf{y} \, d\omega'', \quad (13) $$

where the operator “kernel” is the distribution

$$ \kappa(\mathbf{x}, \omega, \mathbf{y}, \omega') = g(\mathbf{x}, \mathbf{y}) \rho(\mathbf{x}, \omega, p_{\Omega}(\mathbf{y}, \omega'')) \delta_{\mathbf{x}}(p_{S}(\mathbf{y}, \omega'')). $$

## 3.4 Review: Properties of light transport operators

**General properties.** In the absence of refractive materials, operators $\mathbf{K}$ and $\mathbf{G}$ are both self-adjoint due to the cosine in the dot product of Eq. 7 [Veach 1997] (See Appendix A for a proof using our own notations). From conservation of energy, we know that $\|\mathbf{K}\|_2 < 1$ (See Arvo [1996] Sec. 6.3.3) and $\|\mathbf{G}\|_2 \le 1$, where the equality occurs when the scene is “closed”. There, Arvo also proves

$$ \forall p \ge 1 \quad \|\mathbf{T}\|_p < 1 \text{ and } \|\mathbf{T}_b\|_p < 1. \quad (14) $$

**Non-compact operators.** Neither $\mathbf{K}$ nor $\mathbf{T}$ are compact as they are both partial integral operators [Kalitvin and Zabrejko 1991; Soler et al. 2022]. The product of a bounded operator and a compact operator is always compact (see Kato [1995] p158), but $\mathbf{T} = \mathbf{K}\mathbf{G}$ is known to be non-compact. Therefore $\mathbf{G}$ cannot be compact.
**Conditionally-compact operators.** The operator $\mathbf{K}_x$ is a Hilbert-Schmidt (hence compact) operator, due to its bounded kernel (the sliced BRDF in Eq. 8). Thus, the non-compactness of $\mathbf{K}$ is due to the “larger” space that it acts on ($\mathcal{H}$ instead of $O$).

$\mathbf{T}_b$ is compact only in scenes where $\kappa_b$ is everywhere bounded, which imposes conditions of not having abutting edges (or corners) or contact points between surfaces. At an abutting edge, it is indeed possible to construct a bounded sequence of light distributions with increasing spatial frequency content whose image under $\mathbf{T}_b$ does not contain any converging subsequence, proving that $\mathbf{T}_b$ is non compact. Fortunately, $\mathbf{T}_b$ coincides with a compact operator everywhere away from such edges [Soler et al. 2022].

## 4 Compact approximations of operators

In this section we explain that although neither $\mathbf{T}$ nor $\mathbf{T}_b$ is compact in general, they admit-not necessarily finite-dimensional-strongly converging compact approximations in the sense of Eq. 2. We also draw connections between these compact approximations and popular “fixes” adopted in the practical computation of light transport in order to avoid arbitrary large values of the kernel. For example, to avoid speckles (white spots) in path-traced images and images rendered with virtual point lights [Keller 1997], multiple authors proposed to cap the geometric factor [Kollig and Keller 2006], or use a pre-integrated point-to-virtual-surface element geometry factor [Hašan et al. 2009]. Another example is vertex merging [Georgiev et al. 2012] as used in Galerkin approximations within non-Lambertian scenes, where couples of surfaces and directional elements are paired despite their position or directions not matching exactly.

Strong convergence in this context implies that the radiance function resulting from the action of a non-compact transport operator can be obtained as the limit of the value resulting from an approximating sequence of compact operators. Since compact operators form a closed set in the operator norm, these approximations will not be uniformly converging (operator-norm convergence in the sense of Eq. 1). In practice a sequence of operators $\mathbf{T}_n$ that is strongly converging but not uniformly converging to the operator $\mathbf{T}$, allows to find for any $\epsilon > 0$ a sequence of light distributions $\{L_n\}$ for which $\|\mathbf{T}_n L_n - \mathbf{T} L_n\| > \epsilon$ for any $n \ge 0$. We will see that this difference materializes as compact approximations filtering out high frequencies where the original operator would actually preserve them.

**Compact approximations of $\mathbf{T}_b$ and $\mathbf{T}$.** As recalled in Section 3.4, the radiant exitance operator $\mathbf{T}_b$ is not compact since its kernel is not bounded at abutting edges. One way to overcome this is by parameterizing the kernel via a small threshold $\epsilon > 0$ so that

$$ \kappa_{b,\epsilon}(\mathbf{x}, \mathbf{y}) = \frac{\bar{\rho}(\mathbf{x})}{\pi} f_{\epsilon}(\mathbf{x}, \mathbf{y}), \quad (15) $$

where $f_{\epsilon}(\mathbf{x}, \mathbf{y})$ is the capped point-to-point geometry factor<fnref n="3" />

$$ f_{\epsilon}(\mathbf{x}, \mathbf{y}) = v(\mathbf{x}, \mathbf{y}) \min \left( \frac{1}{\epsilon}, \frac{\cos \theta' \cos \theta''}{\|\mathbf{x} - \mathbf{y}\|^2} \right). \quad (16) $$

Because $\kappa_{b,\epsilon}$ is bounded, its use as the integration kernel (instead of $\kappa_b$) leads to a Hilbert-Schmidt operator $\mathbf{T}_{b,\epsilon}$, which is therefore compact.

In the non-Lambertian case, there are two causes for the operator $\mathbf{T}$ not being compact: abutting edges and the partial integration which causes the $\mathbf{1}_x$ function to appear in Equation 13. One possible way to overcome the latter while ensuring strong convergence is to replace $\mathbf{1}_x$ in $\kappa$ by a Gaussian with parameterized sharpness such as

$$ \kappa_{\epsilon}(\mathbf{x}, \omega, \mathbf{y}, \omega'') = g_{\epsilon}(\omega \cdot p_{\Omega}(\mathbf{y}, \omega'')) \rho(\mathbf{x}, \omega, p_{\Omega}(\mathbf{y}, \omega'')) f_{\epsilon}(\mathbf{x}, \mathbf{y}), $$

where $g_{\epsilon}$ is a normalized Gaussian of variance $\epsilon$ centered at 1. In Appendix B we prove the strong convergence of $\mathbf{T}_{\epsilon}$ (resp. $\mathbf{T}_{b,\epsilon}$) to $\mathbf{T}$ (resp. $\mathbf{T}_b$) in the $L_1$ norm for bounded light distributions. While we cannot prove strong convergence for the $L_2$ norm, strong convergence in $L_1$ for bounded light distributions is sufficient in all arguments that require strong continuity in this article. Note in particular that eigenfunctions of $\mathbf{T}$ (resp. $\mathbf{T}_b$) are necessarily bounded because of energy conservation.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Effect of Equation 15 on the first bounce of indirect illumination in the Cornell box. The flattening of curves along geometrically continuous regions is a practical consequence of the energy reduction in high frequencies by compact approximations.</figcaption></figure>

Finite rank Galerkin approximations of $\mathbf{T}_b$ and $\mathbf{T}$ [Sillion et al. 1991], are alternative examples of compact approximations of these operators with strong convergence when their basis functions form a complete Schauder basis for $\mathcal{H}$ [Chatelin 2011]. This explains why finite element methods can yield a converging approximation to solve equation 10 for a specific emission function $E$ and also explains why (due to non-uniform convergence) adaptive meshing is required for controlled error across different choices of $E$.

<footnote n="3">³This arbitrary choice of a capping method leads a simpler convergence proof.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 5. A toy scene with a compact (no abutting edges in the scene) radiant exitance operator that attenuates high frequencies. It is impossible to design an emission function E such that the transported radiance $L = T_b L$ contains a discontinuity. The estimate of the solution E “compensates” for the filtering by introducing arbitrarily high spatial frequencies with large intensities. Due to the strong convergence of the Galerkin approximation, this experiment provides a valid insight into the behavior of the non-discretized operator.</figcaption></figure>

**Compact approximations favor low frequencies.** Although there is no general definition of spatial or angular “frequency” of a function, we adopt common use of this terminology as the non-zero coefficients of functions projected onto a basis with local support in the Fourier domain. For example, spatial or directional wavelets [Peers et al. 2009], Fourier bases [Durand et al. 2005] or spherical harmonics [Soler et al. 2015]. If $\mathbf{A}: \mathcal{X} \to \mathcal{X}$ is a compact operator over a separable Hilbert space, and $\{\varphi_n\}_n$ is any orthonormal Schauder basis of $L^2(\mathcal{X})$, we prove in appendix D that

$$ \lim_{n \to \infty} \| \mathbf{A} \varphi_n \| = 0. \quad (17) $$

In the specific situation of a basis with local frequency support, Eq. 17 directly proves that $\mathbf{A}$ attenuates higher frequencies.<fnref n="4" /> This result applies to compact approximations $T_{b,\epsilon}$ and $T_\epsilon$. This effect is illustrated on Figure 4.

Equation 17 does not imply that compact operators are low-pass filters. Rather, that the action of the operator results in stronger attenuation of higher frequencies. Thus, Lambertian light transport tends be “low frequency” away from edges. For example, consider assigning an emission distribution E to a planar emitter some distance from a parallel planar receiver. If we wish the receiver to contain a step edge in the spatial distribution of reflected radiance L (e.g. as seen in a hard shadow), solving for E such that $L = T_b E$ is not possible. An attempt to approximate E would contain arbitrarily high spatial frequencies at large energy to compensate for Equation 17. This is illustrated in Figure 5. This theory underpins the reason that numerical approximations of light transport inherently attenuate high spatial and angular frequencies.

Since compact approximations have eigenvalues that converge to 0, when the corresponding eigenfunctions form a complete system, they necessarily increase in spatial (or angular) frequency with the index of the corresponding eigenvalue. We prove this in Appendix C.

**Convergence of eigenelements.** While the spectrum of a parameter dependent approximation is generally “continuous” w.r.t. to norm-convergence (see Conway [1991] p431), less can be expected when only strong convergence is satisfied. Generally speaking, strong convergence of a holomorphic family of parameter-dependent compact operators (such as $T_\epsilon$ and $T_{b,\epsilon}$) may display very unexpected behaviors. For example, the approximation could have an infinite number of eigenvalues even though the operator (when $\epsilon = 0$) might have no eigenvalues at all, or an uncountable point-spectrum, etc (See Kato [1995], VII.4, p371).

A useful result however is that when a parameter dependent approximation $\mathbf{A}_\epsilon$ of an operator $\mathbf{A}$ strongly converges to $\mathbf{A}$, then the isolated eigenvalues of $\mathbf{A}$ are the limit of eigenvalues of $\mathbf{A}_\epsilon$ [Chatelin 1981, 2011], neglecting multiplicity. This theory simplifies even further when the operator and its approximations are both self-adjoint and densely defined, in which case the convergence properties cover both eigenvalues and eigenprojections (See Chatelin [1981] definition 2.1 and Section 2.3). Consequently, both eigenvalues and eigenprojections of $T_{b,\epsilon}$ converge toward those of $T_b$ (We show in Sec. 5.1 that $T_b$ is self-adjoint in a particular Hilbert space).

Contrary to $T_b$, having no formal proof that the point-spectrum of $\mathbf{T}$ is countable, the previous conditions do not apply. We are currently limited to observing the convergence of eigenvalues in a number of typical geometric configurations, which tends to suggest that the point-spectrum of $\mathbf{T}$ is countable (See Figure 6). Numerical experiments indeed show that Galerkin approximations of $\mathbf{T}$ are stable (In the sense of Chatelin [Chatelin 1981]). A formal proof of such stable convergence (left for future work) would allow to apply Proposition 2.2 of F.Chatelin’s paper [Chatelin 1981] in order to show that eigenvalues of these approximations actually converge to isolated elements of the point-spectrum of $\mathbf{T}$.

The strong convergence of approximations $T_\epsilon$ and $T_{b,\epsilon}$ also mean that partial eigenfunction expansions of these operators strongly converge, which validates the use of eigenfunctions and eigenvalues of these compact approximations to represent the behavior of the operators over finite dimensional sets of light distributions.

**Usage for eigenfunction expansions.** Expanding elements of $\mathcal{H}$ with the eigenfunctions of an operator requires the family of eigenfunctions of that operator to be dense in $\mathcal{H}$ (accounting for eigenvalue 0). This is the case for $T_b$ (See Section 5),

<footnote n="4">⁴Note that a sequence of functions with decreasing frequency would not be complete, so the same argument cannot be used to show that compact operators filter low frequencies.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Convergence and stability of the largest eigenvalues of a Galerkin approximation of the non-Lambertian transport operator $\mathbf{T}$ in a toy scene made of two squares with the MERL gold-paint material with an increasing number of elements. Left: as opposed to the Lambertian case, the eigenvalues of $\mathbf{T}$ are not necessarily real, because of the inherent non-symmetry of $\mathbf{T}$. Zooming on the red squared region of the complex plane (center) where we show four eigenvalues computed with an increasing number of surface elements, suggests that they display consistent convergence. Since pairs of surface elements are used to define directional finite elements, the effective matrix size is the square of the indicated numbers (e.g. max size here is $300^2 \times 300^2$). Note that eigenvalues come in this scene by pairs of opposite signs because of the symmetry of the scene.</figcaption></figure>

but this is yet to be proven for $\mathbf{T}$. Regardless, the compactness of $\mathbf{T}_\epsilon$ is not sufficient to ensure completeness of its eigenfunctions.<fnref n="5" /> The eigenfunctions $\Lambda_i$ of $\mathbf{T}_b$ (resp. $\mathbf{T}_{b,\epsilon}$) are

$$ \forall L \in \mathcal{B} \quad L = \sum_{i=0}^{\infty} \langle \Lambda_i, L \rangle_S \Lambda_i. \quad (18) $$

For $\mathbf{T}_\epsilon$, we can at least say that a similar equation holds for elements selected inside the span of its eigenfunctions.

In summary we have justified in this section that the eigen-decomposition of the compact approximations of light transport can be used to estimate the countable part of the operators' point spectrum, and that the largest eigenvalues correspond to the energy exchanges of lowest frequency.

## 5 Point-Spectrum of transport operators

In this section, we analyze the point-spectra of operators $\mathbf{T}$, $\mathbf{K}$, $\mathbf{G}$ and $\mathbf{T}_b$. For each operator, we define the elements that are contained within its point-spectrum as a theorem accompanied by a proof. Figure 9 provides a table summarizing our results.

### 5.1 Spectrum of the radiant exitance transport operator

Although $\mathbf{T}_b$ is the infinite-dimensional analog of the radiosity matrix, its non-compactness in the general case precludes direct extensions of the properties of the latter. However, a common “trick” to make the radiosity matrix symmetric positive definite [Baranoski et al. 1997] is useful toward proving that $\mathbf{T}_b$ is symmetrizable.

Assuming that the albedo $\bar{\rho}$ is never null, let $\mathbf{S}: \mathcal{B} \to \mathcal{B}$ be the trivial self-adjoint linear operator defined as

$$ (\mathbf{S}\mathbf{B})(\mathbf{x}) = \frac{1}{\bar{\rho}(\mathbf{x})} \mathbf{B}(\mathbf{x}). \quad (19) $$

This operator is also positive since

$$ \forall B \in \mathcal{B} \quad \langle \mathbf{S}B, B \rangle = \langle \frac{1}{\bar{\rho}}B, B \rangle = \int_S (\rho(\mathbf{x}))^{-1/2} |B|^2 d\mathbf{x}. $$

It is therefore possible to define a new dot-product

$$ \forall f, g \in L^2(S) \quad \langle f, g \rangle_S = \langle f, Sg \rangle, $$

and we have

$$ \langle \mathbf{T}_b f, g \rangle_S = \langle \mathbf{T}_b f, \mathbf{S}g \rangle = \langle f, \mathbf{S}\mathbf{T}_b g \rangle = \langle f, \mathbf{T}_b g \rangle_S. $$

The operator $\mathbf{S}\mathbf{T}_b$ is trivially self-adjoint over $\mathcal{B}$, and $\mathbf{S}\mathbf{T}_b = \mathbf{T}_b\mathcal{S}$. In other words, $\mathbf{T}_b$ is symmetrizable on $\mathcal{B}$, which also means that it is self-adjoint on the Hilbert space $\mathcal{B}_S = (L^2(S), \langle \cdot, \cdot \rangle_S)$. We obtain the following result:

**THEOREM 5.1.1:** *The point-spectrum of $\mathbf{T}_b$ has the following properties:*

(1) $\sigma(\mathbf{T}_b) = \sigma_p(\mathbf{T}_b)$;

(2) $\sigma_p(\mathbf{T}_b) \subset (-1, 1)$;

(3) $\sigma_p(\mathbf{T}_b)$ is countable.

(4) eigenfunctions of $\mathbf{T}_b$ for different eigenvalues are orthogonal w.r.t to the dot product $\langle \cdot, \cdot \rangle_S$;

(5) eigenfunctions of $\mathbf{T}_b$ span $\mathcal{B}$.

(6) the first eigenfunction of $\mathbf{T}_b$ is positive everywhere

**PROOF:** Property (2) is the consequence of the point-spectrum of self-adjoint operators being real (See Lax [Lax 2002], chap.28) combined with Equation 14. Other properties are direct consequences of the spectral theorem: A self-adjoint operator $\mathbf{A}$ on Hilbert space has no residual spectrum, proving (1), and its eigenfunctions for different eigenvalues are orthogonal (4) (See for instance Reed [1972], Th. VI.8). Because $\mathcal{B}$ is separable these eigenfunctions can therefore only be countable (3)– in a separable space every orthogonal sequence is countable. From the multiplication form of the spectral theorem ([Reed and Simon 1972], Th. VII.3), we deduce that the eigenfunctions of $\mathbf{T}_b$ form a Schauder basis of $\mathcal{B}$ (5). The last property

<footnote n="5">⁵A notable counterexample in finite dimensions is that of defective matrices.</footnote>
is a consequence of $\mathbf{T}_b$ being a "positive operator" (in the sense of Lax [Lax 2002], p253). □

Self-adjointness alone does not imply that eigenfunctions of $\mathbf{T}_b$ span $\mathcal{B}$. The operator's point-spectrum needs to be countable which happens in this case when the space is separable, or more generally when the operator is compact. It is possible that 0 belongs to $\sigma_p(\mathbf{T}_b)$, for instance when the geometry includes some finite surface that is not visible from any other surface. If 0 is in the point spectrum, a complete system of eigenfunctions of $\mathbf{T}_b$ includes its corresponding eigenfunctions. Finally, a complete system of eigenfunctions with no null eigenvalue does not make the operator invertible, since the span of eigenfunctions is only dense in $\mathcal{B}$ and not equal to $\mathcal{B}$.

**Distribution of the eigenvalues.** When $\mathbf{T}_b$ is compact, in addition to being countable, $\sigma_p(\mathbf{T}_b)$ can only have 0 as an accumulation point. When it is not compact, eigenvalues are still countable (since the operator is self-adjoint and the space is separable), and will therefore necessarily accumulate at extra non zero values within $(-1, 1)$. These eigenvalues correspond to eigenfunctions for which most the energy is concentrated next to an abutting edge. This explains why compact approximations of $\mathbf{T}_b$ are practical for image synthesis: while they remove non negligible eigenvalues, the corresponding energy exchanges are mostly local.

When $\mathbf{T}_b$ is compact, it is also Hilbert-Schmidt as are its compact approximations $\mathbf{T}_{b,\epsilon}$. Since it is self-adjoint in $\mathcal{B}_S$, we can apply Eq. 4 to an orthogonal basis of eigenfunctions in $\mathcal{B}_S$ to get $\sum_i \lambda_i^2 < \infty$. This brings the following upper bound on the asymptotic behavior of its eigenvalues: $\lambda_n = o(1/\sqrt{n})$. Figure 7 demonstrates this using a simple scene with no abutting edges. Determining the actual asymptotic behavior of the spectrum is non-trivial and depends on the specific operator (see for instance Jeribi [2021] Sec 11.4).

<figure><img src="image_10.png" /><figcaption>Fig. 7. Toy example of a scene where $\mathbf{T}_b$ is compact showing the upper bound on the asymptotic decrease of eigenvalues.</figcaption></figure>

Finally, while finite multiplicity of eigenvalues is not a limitation in this article, it is fairly straightforward to construct

scenes where the eigenvalues of $\mathbf{T}_b$ have multiplicity greater than 1: given a closed scene $S_1$ (no light leaks), we build a new scene by duplicating the geometry of $S_1$ multiple times, to obtain as many closed and separate environments. Every eigenpair $(\lambda, \Lambda)$ of $\mathbf{T}_b$ in $S_1$ will give birth to as many eigenpairs in the new scene where each closed subpart will have $\Lambda$ while other subparts are 0. Having an infinite number of scaled copies which total surface area stays finite would allow to construct a scene were eigenvalues have infinite multiplicities, making the operator non-compact even if the original scene had no abutting edges. This however is outside of physical constraints.

**Singular value decomposition of $\mathbf{T}_b$.** Since $\mathbf{T}_b$ is self-adjoint in $\mathcal{B}_S$, and given that its eigenvalues $\{\lambda_k\}_{k \ge 0}$ are real and countable, and its eigenfunctions $\{\Lambda_k\}_{k \ge 0}$ form a complete orthogonal basis of $\mathcal{B}_S$ (See Section 5.1), $\mathbf{T}_b$ has the following "Schmidt expansion" in $\mathcal{B}_S^6$:

$$ \mathbf{T}_b L = \sum_{k=0}^{\infty} |\lambda_k| \langle L, \Lambda_k \rangle_S \underbrace{\operatorname{sign}(\lambda_k) \Lambda_k}_{u_k}. $$

The sign of $\lambda_k$ has been transferred to the left singular functions $u_k$ in order to keep the singular values non-negative. This expansion is however not uniformly converging since the eigenvalues of $\mathbf{T}_b$ may have non zero accumulation points when the operator is not compact. In $\mathcal{B}$ however, there is no such trivial expression for a Schmidt expansion of $\mathbf{T}_b$, and numerical computations bring singular values that are different from its eigenvalues, mostly because the eigenfunctions of $\mathbf{T}_b$ are not orthogonal in this space.

## 5.2 Spectrum of the reflectance operators

The local and global scattering operators $\mathbf{K}_x$ and $\mathbf{K}$ have different spectra. While they share eigenvalues, the eigenvalues of the latter have infinite dimensionality.

**THEOREM 5.2.1:** *The spectrum of $\mathbf{K}_x$ is equal to its pointspectrum, and has the following characteristics:*

(1) $\sigma_p(\mathbf{K}_x) \subset (-1, 1)$

(2) the eigenfunctions of $\mathbf{K}_x$ span $\mathcal{O}$

**PROOF:** This operator is Hilbert-Schmidt (therefore compact) and self-adjoint. Its eigenvalues are consequently real and equal to its singular values, and its eigenfunctions—including eigenvalue 0, if present—span the half-space of directions $\mathcal{O}$. Furthermore the eigenvalues of $\mathbf{K}_x$ lie in the open interval $(-1, 1)$ since energy conservation imposes $\|\mathbf{K}_x\| < 1$. Since $\mathbf{K}_x$ is compact, the only possible accumulation point for its eigenvalues is zero, which may or not be part of the pointspectrum even if the rank of $\mathbf{K}_x$ is infinite. □

Figure 8 shows the largest eigenvalues of $\mathbf{K}_x$ for some elements of the MERL database. Since the integral kernel of $\mathbf{K}_x$ is the BRDF, we sometimes refer to its eigenvalues as the "eigenvalues of the BRDF".

<footnote n="6">⁶Note that the term *schmidt expansion* is not perfectly appropriate here since the convergence is not uniform.</footnote>
**THEOREM 5.2.2:** *The point-spectrum of K has the following characteristics:*

**case 1:** if there is no measurable area where the BRDFs all share the same eigenvalue then $\sigma_p(\mathbf{K}) = \emptyset$;

**case 2:** else if the scene can be partitioned into a finite number of measurable areas with respective BRDFs $\rho_i$, then $\sigma_p(\mathbf{K}) = \bigcup_{i \ge 0} \sigma_p(\rho_i)$. Its eigenvalues have infinite dimensionality and its eigenfunctions span $\mathcal{H}$. This case is typical of common light transport scenes;

**case 3:** else, the point spectrum of $\mathbf{K}$ will be limited to the material eigenvalues that are shared across measurable subsets of $S$, each having infinite dimensionality. In this case the eigenfunctions of $\mathbf{K}$ do not span $\mathcal{H}$.

**PROOF:** Because $\mathbf{K}$ is a tensor product between $\mathbf{K}_x$ and an identity operator over $L^2(S)$, the existence of an eigenpair $(k, \phi)$ of $\mathbf{K}$ also implies that $\mathbf{K}_x\phi(\mathbf{x},.) = k\phi(\mathbf{x},.)$ wherever $\phi(\mathbf{x},.) \ne 0$, which means that $(k, \phi(\mathbf{x},.))$ is also an eigenpair of the BRDF at $\mathbf{x}$. And, in order to be a non-zero function in $\mathcal{H} = L^2(S \times \Omega)$, the property $\phi(\mathbf{x},.) \ne 0$ should remain valid over a measurable subset $A \subset S$. On the contrary, when the BRDFs across $A \subset S$ all share the same eigenvalue $k$ with eigenfunction $\rho_x$, any function of the form $(\mathbf{x}, \omega) \mapsto \rho_x(\omega)f(\mathbf{x})$ over $A$ and $0$ elsewhere is an eigenfunction of $\mathbf{K}$, and from Th. 5.2.1, $L^2(A \times \Omega)$ is spanned by these eigenfunctions. □

<figure><img src="image_11.png" /><figcaption>Fig. 8. Top: Point-spectrum of some reflectance operators from the MERL database. Bottom: eigenfunctions 0, 6 and 20 of the gold-paint material, computed by projecting the BRDF onto a basis of spherical harmonics up to degree 30 (green means negative). As theory predicts (See Lax [2002], chap.28), only the first eigenfunction is entirely positive.</figcaption></figure>

## 5.3 Spectrum of the reparametrization operator

Let $\mathcal{L}_1 \subset \mathcal{L}$ denote the set of points $(\mathbf{x}, \omega) \in \mathcal{L}$ for which $p_x(\mathbf{x}, \omega)$ is defined, $\mathcal{H}_1 = L^2(\mathcal{L}_1)$ and $\mathbf{P}_1: \mathcal{H} \to \mathcal{H}_1$ the orthogonal projection on $\mathcal{H}_1$. We also call $\mathbf{G}_1: \mathcal{H}_1 \to \mathcal{H}_1$ the restriction of $\mathbf{G}$ to $\mathcal{H}_1$. Note that because $p_x$ is reversible wherever it is defined, using $\mathcal{H}_1$ for the output space of $\mathbf{G}_1$ makes sense. With these settings, we have

$$ \mathbf{G} = \mathbf{G}_1\mathbf{P}_1 \text{ and } \mathbf{G}_1^2 = \mathbf{I}. \qquad (20) $$

Similarly to $\mathbf{G}$ operator $\mathbf{G}_1$ is self-adjoint, and $\|\mathbf{G}_1\| = 1$ which makes it an isometry on $\mathcal{H}_1$. Its eigenvalues are therefore $\pm 1$. Furthermore, to any eigenpair $(\lambda, V_1)$ of $\mathbf{G}_1$ corresponds an eigenpair $(\lambda, V)$ of $\mathbf{G}$ with $V(\mathbf{x}, \omega) = 0$ for any $(\mathbf{x}, \omega) \notin \mathcal{L}_1$.

Finally, any function $V \in \mathcal{H}$ such that $P_1V = 0$ satisfies $GV = 0$. Therefore 0 is also in the point-spectrum of $\mathbf{G}$. No other eigenvalue can exist, as it would automatically be an eigenvalue of $\mathbf{G}_1$. In closed scenes, $P_1 = \mathbf{I}$, $\mathcal{L}_1 = \mathcal{L}$, and $\mathbf{G}_1 = \mathbf{G}$, and the part of the point-spectrum corresponding to the nul eigenvalue vanishes.

**THEOREM 5.3.1:** *Spectrum of the reparametrization operator*

(1) $\sigma_p(\mathbf{G}) = \{-1, 1\}$ in closed scenes;

(2) $\sigma_p(\mathbf{G}) = \{-1, 0, 1\}$ otherwise;

(3) the eigenfunctions of $\mathbf{G}$ form a complete family of $\mathcal{H}$.

**PROOF:** In order to prove (3), let $\{\varphi_n\}_{n \ge 0}$ (resp. $\{\psi_n\}_{n \ge 0}$) be an orthogonal Schauder basis of $\mathcal{H}_1$ (resp. $\mathcal{H} \setminus \mathcal{H}_1$).<fnref n="7" /> Assuming that both functions can be extended by 0 in $\mathcal{L}$, we define the sequence:

$$ G_{3n} = \frac{1}{2}(\varphi_n + \mathbf{G}\varphi_n) $$

$$ G_{3n+1} = \frac{1}{2}(\varphi_n - \mathbf{G}\varphi_n). $$

$$ G_{3n+2} = \psi_n. $$

Since $\mathbf{G}^2 = \mathbf{I}$ we obtain

$$ \mathbf{G}\mathbf{G}_{3n} = G_{3n} \quad \mathbf{G}\mathbf{G}_{3n+1} = -G_{3n+1} \quad \mathbf{G}\mathbf{G}_{3n+2} = 0. $$

After noting that $\varphi_n = G_{3n} + G_{3n+1}$ we see that $\{G_n\}_{n \ge 0}$ spans $\mathcal{H}$, which makes it a complete—yet not orthogonal inside each eigenspace—family of eigenfunctions of $\mathbf{G}$. A fully orthogonal family of eigenfunction can further be obtained by Gram–Schmidt orthogonalization in each eigenspace. □

## 5.4 Spectrum of the radiance transport operator

Unlike the point-spectrum of $\mathbf{T}_b$ (Section 5.1), much less can be determined about the point-spectrum of $\mathbf{T}$ because of the combination of the following properties:

(1) $\mathbf{T}$ is not compact;

(2) $\mathbf{T}$ is not normal;

(3) $\mathbf{T}$ cannot be symmetrized.

<footnote n="7">⁷These spaces are separable, and any such combination of spatial and directional wavelets would work.</footnote>
The non-compactness of $\mathbf{T}$ has been discussed in Section 3.4. To demonstrate non-normality, since $\mathbf{K}$ and $\mathbf{G}$ are self-adjoint, we have

$$ 
\begin{aligned}
\mathbf{T}\mathbf{T}^* &= \mathbf{KGGK} = \mathbf{K}^2, \\
\mathbf{T}^*\mathbf{T} &= \mathbf{GKKG} = \mathbf{GK}^2\mathbf{G}.
\end{aligned}
$$

That is, reflecting twice is not equivalent to transporting, reflecting twice and transporting back. Note that neither is $\mathbf{T}$ quasi-normal nor hyponormal, which would yield beneficial properties [Conway 1991]. (3) is a side effect of the Dirac in Equation 13, which precludes multiplication of $\mathbf{T}$ by a trivial operator (such as $\mathbf{S}$ in Eq. 19) to make it self-adjoint.

**THEOREM 5.4.1:** *The spectrum of $\mathbf{T}$ has the following properties:*

(1) $\sigma_p(\mathbf{T})$ is not empty; the largest eigenvalue is positive with a positive eigenfunction;

(2) $\sigma_p(\mathbf{T})$ lies in the disk of radius $\|\mathbf{K}\|$;

(3) the eigenfunctions of $\mathbf{T}$ for non-zero eigenvalues do not necessarily span $\mathcal{H}$.

**PROOF:** Since $\mathbf{T}$ is a bounded operator it has a non-empty spectrum (See Kato [1995] Ch3,§6.2), but this does not guarantee the same property for its point-spectrum. But since its kernel is positive (assuming non-zero reflectance), the Krein-Rutman theorem [Phát and Dieu 1994] applied to the cone of light distributions of constant sign ensures that when the scene is closed we have $\sigma_p(\mathbf{T}) \neq \emptyset$ and its maximal eigenvalue is a positive real value of unit multiplicity, associated to an all-positive eigenfunction. When the scene is not closed, restricting $\mathbf{T}$ to $\mathcal{H}_1$ preserves the conditions to apply this theorem, hence proving (1). Note that the same applies to $\mathbf{T}_b$ which coincides with $\mathbf{T}$ when materials are Lambertian. (2) is straightforward: let $(V, \lambda)$ be an eigenpair of $\mathbf{T}$ (which means $\|V\| = 1$), then $\|\mathbf{T}V\| = |\lambda|$ and $\|\mathbf{T}V\| \le \|\mathbf{K}\|\|\mathbf{G}\|$. Using Equation 14 we get $|\lambda| \le \|\mathbf{K}\|$.

Since $\mathbf{T} = \mathbf{K}\mathbf{G}$, any eigenfunction $V$ of $\mathbf{T}$ is in the range of $\mathbf{K}$, so is the span of the eigenfunctions, which proves (3). The extreme case is in Lambertian scenes, where eigenfunctions for non-zero eigenvalues are directionally constant, and therefore cannot represent any non-directionally constant function in $\mathcal{H}$.

The following might be useful in practical contexts where one set of eigenfunctions might be computed easily, either for transporting light or importance:

**THEOREM 5.4.2:** *The spectrum of $\mathbf{T}$ and $\mathbf{T}^*$ are connected via the operators $\mathbf{G}$ and $\mathbf{K}$:*

(1) The eigenfunctions of $\mathbf{T}$ and $\mathbf{T}^*$ are mutually orthogonal for different eigenvalues;

(2) $\mathbf{G}$ turns the former into the later, while $\mathbf{K}$ turns the later into the former;

(3) $\sigma_p(\mathbf{T}) = \sigma_p(\mathbf{T}^*)$.

**PROOF:** Orthogonality is a classical result. In order to prove (2) and (3), let $(\lambda, L)$ be an eigenpair of $\mathbf{T}^*$. We have $\mathbf{G}\mathbf{K}\mathbf{L} = \lambda\mathbf{L}$. Left multiplying by $\mathbf{G}$ leads to $\mathbf{K}\mathbf{L} = \lambda\mathbf{G}\mathbf{L}$. Reporting the former equation in the later, we get $\mathbf{K}(\mathbf{G}\mathbf{K}\mathbf{L}) =$

$\lambda\mathbf{G}(\mathbf{G}\mathbf{K}\mathbf{L})$, which simplifies to $\mathbf{T}(\mathbf{K}\mathbf{L}) = \lambda(\mathbf{K}\mathbf{L})$. Conversely, if $(\lambda, V)$ is an eigenpair of $\mathbf{T}$, we have $\mathbf{G}\mathbf{K}\mathbf{L} = \lambda\mathbf{L}$. Left-multiplying by $\mathbf{G}$ directly leads to $\mathbf{T}^*(\mathbf{G}\mathbf{L}) = \lambda(\mathbf{G}\mathbf{L})$. □

Note that the above properties do not imply that eigenvalues are real, and indeed, numerical experiments tend to show that even in the simplest cases they are not (See Figure 6). Finally, we lack a formal proof that $\sigma_p(\mathbf{T})$ is countable. It seems to be a reasonable conjecture given the stable convergence of eigenvalues of its Galerkin approximations [Chatelin 2011] (See Figure 6). We leave such a proof for future work.

## 5.5 Schmidt expansion of radiance transport operator

Interestingly, $\mathbf{T}$ admits a simple explicit (yet non-uniformly converging) “Schmidt expansion” in scenes with a finite number of non-spatially-varying materials $m_i$<fnref n="8" /> which we briefly recall here for completeness [Soler et al. 2022]. This result will be used to derive the spectral properties of $\mathbf{T}$. We call $\rho_{ij}$ the j-th eigenvalue of material $m_i$ and its corresponding eigenfunction $r_{ij} \in \mathcal{O}$. We define $\{\phi_k^i\}_{k \ge 0}$ to be an orthogonal basis of the subset of $S$ where the material is $m_i$ (wavelets are such an example). Finally we name $g: j \mapsto (n_1, n_2)$ an arbitrary bijection between $\mathbb{N}$ and $\mathbb{N}^2$, and use $j = g(n_1, n_2)$. Thus, we have

$$ \mathbf{T}\mathbf{L} = \sum_{i=1}^{n} \sum_{j>0} |\rho_{ij}| \langle L, \varphi_j^i \rangle \psi_j^i, \quad (21) $$

where for all $(n_1, n_2) \in \mathbb{N}^2$,

$$ \psi_j^i(\mathbf{x}, \omega) = \phi_{n_1}^i(\mathbf{x}) r_{n_2}^i(\omega) \quad (22) $$

$$ \varphi_j^i(\mathbf{x}, \omega) = \text{sign}(\rho_{ij}) \mathbf{G} \psi_j^i(\mathbf{x}, \omega). \quad (23) $$

The above equation takes the form of Equation 6 after re-ordering the elements in the sum.

## 6 Light transport formulations of the spectrum

We now derive connections between the eigenvalues of compact approximations of the light transport operators and light paths and their integrals. In some cases, the connection is direct such as the structure of the operators’ output space being governed by its eigenelements (Section 6.1). In Section 6.2 we discover a relationship between eigenvalues of the operator and linear combinations of multiple bounces of light. In Section 6.3 we explain how the measure of single light paths of fixed length, as a pure geometric quantity (with no specific light distribution involved), is also closely related to the eigenvalues. Finally in 6.4 we show that eigenfunctions can be expressed as an integral over solutions of a parameter-dependent generalized light transport problem.

All results derived in this section remain valid of any compact approximation of $\mathbf{T}$ and $\mathbf{T}_b$, including such as those described in Section 4. We use the generic notation $\tilde{\mathbf{T}}$ to designate them. Note that while some of the following results

<footnote n="8">⁸This is the most general situation, assuming that surfaces with spatially-varying-materials can further be separated in sub-regions with constant materials</footnote>
<table><caption>Fig. 9. Summary of spectral properties of the different operators involved in light transport (See definitions in Section 3.1). Completeness accounts for eigenfunctions corresponding to eigenvalue 0.</caption><thead><tr><th>operator</th><th>space</th><th>normality</th><th>class</th><th>point-spectrum</th><th>completeness</th></tr></thead><tbody><tr><td>T<sub>b</sub></td><td>B</td><td>symmetrizable</td><td>bounded</td><td>countable</td><td>yes</td></tr><tr><td>T<sub>b,ε</sub></td><td>B</td><td>symmetrizable</td><td>Hilbert-Schmidt</td><td>countable → 0</td><td>yes</td></tr><tr><td>K<sub>x</sub></td><td>O</td><td>self-adjoint</td><td>Hilbert-Schmidt</td><td>countable → 0</td><td>yes</td></tr><tr><td>K</td><td>H</td><td>self-adjoint</td><td>bounded</td><td>countable (∞-dim)</td><td>yes</td></tr><tr><td>G</td><td>H</td><td>self-adjoint</td><td>unitary</td><td>{0, 1, -1}</td><td>yes</td></tr><tr><td>T</td><td>H</td><td>none</td><td>bounded</td><td>countable (conjectured)</td><td>conjectured</td></tr><tr><td>T<sub>ε</sub></td><td>H</td><td>none</td><td>Hilbert-Schmidt</td><td>countable → 0</td><td>conjectured</td></tr></tbody></table>

would be straightforward for matrices, careful proofs are still required for these infinite dimensional operators.

## 6.1 Structure of the operator's output space

Projecting $\tilde{T}$ onto a finite $N$-dimensional subspace of its output space, results in a matrix whose eigenvalues estimate the largest elements of $\sigma_p(\tilde{T})$. This formulation is inspired from the Monte-Carlo approximations of large matrices proposed by Halko [2011].

Let $\mathcal{H}_N$ be the $N$-dimensional subspace of $\mathcal{H}$ spanned by the $N$ first elements of a complete sequence of functions of finite support in the frequency domain (wavelets, spherical harmonics, etc). Consider $p$ light distributions $L_1, \dots, L_p$ in this space sampled from a Gaussian distribution. For some $n$, let $q_1, \dots, q_n$ be an orthogonal basis of the $n$-dimensional span $Q$ of $\{\tilde{T}L_1, \tilde{T}L_2, \dots, \tilde{T}L_p\}$ (most of the time $n = p$). We also define

$$Q: \mathbb{R}^n \to Q \\ (a_1, \dots, a_n) \mapsto a_1 q_1 + \dots + a_n q_n.$$

Using this notation, the adjoint $Q^*: Q \to \mathbb{R}^n$ of $Q$ computes the $a_i$ coordinates for any $L \in Q$. Both $Q$ and $Q^*$ can be plunged into $\mathcal{H}$ so that $QQ^*$ is the orthogonal projection on $Q$, whereas $Q^*Q$ is the identity over $\mathbb{R}^n$. Given this, the following theorem holds.

**THEOREM 6.1.1:** For any $\epsilon > 0$ there exists $(n, N) \in \mathbb{N}^2$ such that the $n \times n$ matrix $M = Q^* \tilde{T} Q$ satisfies

$$(Mv = \lambda v) \Rightarrow \| \tilde{T} Q v - \lambda Q v \| \le \epsilon,$$

with probability $1-10^{-n}$. That is, the eigenvalues of $M$ statistically approximate some eigenvalues of $\tilde{T}$ using sufficiently many input light distributions $L_i$.

**PROOF:** Let $\epsilon > 0$ and $T_N = P_N \tilde{T}$ where $P_N$ is the orthogonal projection on $\mathcal{H}_N$. Since $\tilde{T}$ is compact, there exists $N$ such that $\|T_N - \tilde{T}\| \le \epsilon/3$. We therefore have

$$\| \tilde{T} - QQ^* \tilde{T} \| \le \frac{\epsilon}{3} + \| T_N - QQ^* T_N \| + \| QQ^* (T_N - \tilde{T}) \|.$$

The 3rd term above verifies $\|QQ^*(T_N - \tilde{T})\| \le \|QQ^*\|\|T_N - \tilde{T}\| \le \epsilon/3$ since, as an orthogonal projection $\|QQ^*\| = 1$. The second term is majorated using the framework of Halko [2011] (Sections 4.3 and 4.4) which we rephrase here in our context:

**PROPOSITION 6.1.1 (HALKO, 2011).** Let $A$ be a $N \times P$ matrix, and $\omega_0, \omega_1, \dots$ a sequence of random Gaussian vectors, and $Q$ the matrix built as above. As soon as $n$ consecutive vectors verify $\|(I - QQ^*)A\omega_k\| \le \epsilon/(10\sqrt{2/\pi})$, then with probability $1 - 10^{-n}$ we have $\|(I - QQ^*)A\| \le \epsilon$.

We apply this result to the operator $T_N P_N$, that is a finite rank operator over the finite dimensional space $\mathcal{H}_n$ which therefore can be assimilated to a matrix, to determine the existence of $n$ such that $\|T_N - QQ^*T_N\| \le \epsilon/3$. In summary, there exists $N$ and $n$ such that

$$\| \tilde{T} - QQ^* \tilde{T} \| \le \epsilon. \quad (24)$$

Let $(\lambda, v)$ be an eigenpair of $M$. We have

$$\| \tilde{T} Q v - \lambda Q v \| \le \| \tilde{T} Q v - QQ^* \tilde{T} Q v \| + \| \underbrace{Q Q^* \tilde{T} Q v}_{QMv} - \lambda Q v \|.$$

The first term on the right hand side is less than $\epsilon$ because of Eq. 24 and since $Mv = \lambda v$, the second term is equal to 0, which yields $\| \tilde{T} Q v - \lambda Q v \| < \epsilon$. $\square$

This theorem should be understood as follows: $\mathcal{H}_N$ being the space where the $L_i$ are sampled, choosing $n$ large enough will reduce the approximation error $\|\tilde{T}Qv - \lambda Qv\|$ down to the threshold imposed by $N$. However, using a very large $N$ will require $n$ to be large as well in order to reduce $\|T_N - QQ^*T_N\|$. Halko [2011] conjectures that $n$ should be slightly larger than $N$ in order to satisfy Eq. 24. In other words, controlling the error requires sampling a small number of low frequency light distributions, or a larger number of higher-frequency light distributions.

As defined in the theorem, the matrix $M = Q^* \tilde{T} Q$ resembles a product of matrices, but it is not. Operator $Q$ may be viewed as an "infinite matrix" with $n$ columns each of which is one of the light distributions $q_i$. However, we may estimate $M$ by sampling these light distributions at a finite (yet sufficiently large) number of points and directions.

## 6.2 Linear combinations of powers of $\tilde{T}$

Here we derive that the minimization of a linear combination of powers of $T$ results in a polynomial whose roots estimate the largest elements of the point-spectrum of $T$. While such a connection seems obvious for a matrix (choosing $1 + P(x)$ to be proportional to the characteristic polynomial of that
matrix) the following theorem shows that it holds for our compact approximations.

**THEOREM 6.2.1:** For any light distribution $L \in \mathcal{H}$ with a finite support on the sequence of eigenfunctions of $\tilde{\mathbf{T}}$ there exists a finite set of coefficients $\alpha_1, ..., \alpha_n$ such that

$$L + \sum_{k=1}^{n} \alpha_k \tilde{\mathbf{T}}^k L = 0. \quad (25)$$

Furthermore, the roots of polynomial $1 + \sum_k \alpha_k x^k$ are the corresponding eigenvalues of $\tilde{\mathbf{T}}$.

**PROOF:** If $\{\Lambda_i\}$ is the sequence of eigenfunctions of $\tilde{\mathbf{T}}$ (which we do not assume to be complete nor orthogonal), and $L$ has a finite support on this basis, then there exist $n+1$ coefficients $\beta_i \neq 0$ (without loss of generality we assume no gap in the related $\Lambda_i$) such that

$$L = \sum_{p=0}^{n} \beta_p \Lambda_p. \quad (26)$$

Let $\alpha_0, ..., \alpha_n \in \mathbb{C}$. Applying $\tilde{\mathbf{T}}^k$ to equation 26, and exchanging sums (both are finite) we immediately get

$$\sum_{k=0}^{n} \alpha_k \tilde{\mathbf{T}}^k L = \sum_{p=0}^{n} \beta_p P(\lambda_p) \Lambda_p. \quad (27)$$

where $P$ is the polynomial $\sum_{k=0}^n \alpha_k X^k$, hence defining operator $P(\tilde{\mathbf{T}})$ (See Conway [1990] Sec.4.10).<fnref n="9" />

Because of equation 27 the necessary and sufficient condition for $P(\tilde{\mathbf{T}})L = 0$ is that $P(\lambda_p) = 0$ for each $p$. Such a polynomial being defined up to a factor, we choose $\alpha_0 = 1$. Polynomial $\prod_{i=0}^n (1 - \frac{1}{\lambda_i} X)$ is in this case the unique solution, whose expansion provides the $\alpha_i$ coefficients for Equation 25. □

In practice, this theorem proves useful to estimate eigenvalues because of the increasing oscillatory nature of eigenfunctions of $\tilde{\mathbf{T}}$, which means that any finite-frequency light distribution will have a limited support on the sequence of eigenfunctions, or at least very small values of $\beta_n$ beyond some $N$. Consequently even when the coefficients do not exactly vanish, optimizing

$$(\alpha_1, ..., \alpha_n) = \operatorname{argmin} \|L + \sum_{k=0}^{n} \alpha_k \tilde{\mathbf{T}}^k L\|^2 \quad (28)$$

produces a polynomial whose roots will approach the largest eigenvalues. We illustrate this in Figure 10.

Theorem 6.2.1 is related to theorem 6.1.1 in the following way: let's choose $\{L_k\}$ (in theorem 6.1.1) to belong to the span $\mathcal{M}$ of the first eigenfunctions of $\tilde{\mathbf{T}}$, and call $C$ the characteristic polynomial of matrix $\mathbf{M}$. Because $\mathcal{M}$ is stable by $\tilde{\mathbf{T}}$ we have $C(\mathbf{M}) = \mathbf{Q}^* C(\tilde{\mathbf{T}}) \mathbf{Q}$. Since $C(\mathbf{M}) = 0$ (Cailey-Hamilton theorem), we also have $\mathbf{Q}^* C(\tilde{\mathbf{T}}) \mathbf{Q} = 0$.

Consequently, for any $L \in \mathcal{M}$ in theorem 6.2.1 the corresponding polynomial $P$ is $C$ (or a factor of $C$ depending on the support of $L$ over the $\{\Lambda_i\}$).

## 6.3 The measure of circular light paths

In this subsection we reveal an unexpected connection between the measure of circular light paths in a scene and the eigenvalues of $\tilde{\mathbf{T}}$.

**THEOREM 6.3.1:** The eigenvalues of $\tilde{\mathbf{T}}$ are the reciprocals of the poles of the entire function $d(\lambda) = \sum_{n=0}^{\infty} (-1)^n a_n \lambda^n$, where

$$a_n = -\frac{1}{n} \sum_{k=2}^{n} (-1)^k p_k a_{n-k}, \quad (29)$$

and $p_n$ is the measure of circular light paths of length $n$ in the scene, with $p_0 = a_0 = 1$ and $p_1 = a_1 = 0$.

**PROOF:** In Section 4 we saw that the compact approximations $\mathbf{T}_\epsilon$ and $\mathbf{T}_{b,\epsilon}$ of $\mathbf{T}$ and $\mathbf{T}_b$ also happen to be Hilbert-Schmidt. Assuming $\tilde{\mathbf{T}}$ is one of these, we can express the Hilbert-Carleman determinant (See Gohberg et al. [2012] p176) of $\tilde{\mathbf{T}}$ as

$$\operatorname{Det}(\mathbf{I} + \lambda \tilde{\mathbf{T}}) = \sum_{n=0}^{\infty} a_n \lambda^n \quad \text{with} \quad (30)$$

$$a_n = \frac{1}{n!} \begin{vmatrix} 0 & n-1 & 0 & \dots & 0 & 0 \\ \operatorname{Tr}(\tilde{\mathbf{T}}^2) & 0 & n-2 & 0 & \dots & 0 \\ \operatorname{Tr}(\tilde{\mathbf{T}}^3) & \operatorname{Tr}(\tilde{\mathbf{T}}^2) & 0 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \operatorname{Tr}(\tilde{\mathbf{T}}^n) & \operatorname{Tr}(\tilde{\mathbf{T}}^{n-1}) & \operatorname{Tr}(\tilde{\mathbf{T}}^{n-2}) & \dots & \operatorname{Tr}(\tilde{\mathbf{T}}^2) & 0 \end{vmatrix}$$

with $a_0 = 1$ and $a_1 = 0$. Since $\tilde{\mathbf{T}}$ is Hilbert-Schmidt, all operators $\tilde{\mathbf{T}}^k$ are also trace-class for $k \ge 2$, hereby justifying the use of the trace in Equation 30. Just like the more classical Fredholm determinant (Gohberg et al. [2012] Chp.VI), the Hilbert-Carleman determinant is an entire function over $\mathbb{C}$ which zeros are the $\lambda = -\frac{1}{\lambda_i}$, where $\lambda_i$ are the eigenvalues of $\tilde{\mathbf{T}}$.

For $k \ge 2$, denoting by $\kappa_k$ the kernel of the trace-class operator $\tilde{\mathbf{T}}^k$, thanks to Equation 5 we have [Brislawn 1988]

$$\begin{align} \operatorname{Tr}(\tilde{\mathbf{T}}^k) &= \int_S \kappa_k(\mathbf{x}, \mathbf{x}) d\mathbf{x} \nonumber \\ &= \int_S \int_{S^{k-1}} \kappa(\mathbf{x}, \mathbf{x}_1) \kappa_{k-1}(\mathbf{x}_1, \mathbf{x}) d\mathbf{x}_1 d\mathbf{x} \nonumber \\ &\dots \nonumber \\ &= \int_{S^k} \kappa(\mathbf{x}_1, \mathbf{x}_2) \kappa(\mathbf{x}_2, \mathbf{x}_3) \dots \kappa(\mathbf{x}_k, \mathbf{x}_1) d\mathbf{x}_1 \dots d\mathbf{x}_k \tag{31} \end{align}$$

In the non lambertian case the $\mathbf{x}_i$ in these equations should be understood as couples of points and directions. Because of the arbitrarily sharp Gaussian $g_\epsilon$ in the kernel expression of $\tilde{\mathbf{T}}$, only the ones that actually constitute a correct light path

<footnote n="9">⁹Note that self-adjointess is not required here</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Illustration of Th.6.2.1, where we show that the largest eigenvalues of $\tilde{T}$ can be approximated as roots of a polynomial which coefficients $\alpha_k$ zero the first bounces of light $\tilde{T}^k L$.</figcaption></figure>

are non zero. In other words for both Lambertian and non-Lambertian scenes the above equations prove that $\mathrm{Tr}(\tilde{T}^k)$ corresponds to the integral of all circular light paths of length $k$ in the scene. Henceforth, we use the notation $p_k = \mathrm{Tr}(\tilde{T}^k)$. Rewriting the determinant in Equation 30 using Laplace expansion along the first column leads to the recursive expression of Equation 29 (See Appendix E for further details). □

The simplest method to compute $p_k$ is Monte-Carlo integration, via a random starting point and direction followed by importance sampling of the path up to length $k$. The last vertex is then connected to the first.<fnref n="10" /> Doing so, the probability density of sampling the circular path $\pi_k = (\mathbf{x}_1, ..., \mathbf{x}_k)$ is

$$P(\pi_k) = \frac{1}{S} \kappa(\mathbf{x}_1, \mathbf{x}_2) \dots \kappa(\mathbf{x}_{k-1}, \mathbf{x}_k),$$

and we have

$$p_k = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^{N} \frac{1}{P(\pi_k)} \kappa(\mathbf{x}_1, \mathbf{x}_2) \kappa(\mathbf{x}_2, \mathbf{x}_3) \dots \kappa(\mathbf{x}_k, \mathbf{x}_1) \\ = \lim_{N \to \infty} \frac{S}{N} \sum_{i=1}^{N} \kappa(\mathbf{x}_k, \mathbf{x}_1).$$

Even though the calculation of $p_n$ is stable, estimating eigenvalues via a finite number of these path integrals converges slowly. The roots of polynomials are generally highly sensitive to noise in their coefficients, especially when these roots are evenly spaced [Wilkinson 1959]. This is why in a practical calculation scenario truncating $\mathrm{Det}(\mathbf{I} + \lambda \tilde{T})$ to its first terms causes instability when determining its roots. Figure 11 shows practical values of $p_n$ and $\lambda_n$ in a simple scene where the values vanish rapidly with increasing $n$.

## 6.4 Parameter-Dependent light transport solutions

The previous subsections derived connections between path integrals and eigenvalues. We now explain how it is also possible to express eigenfunctions as a sum of solutions to a parameter-dependent light transport problem in the complex domain.

<figure><img src="image_13.png" /><figcaption>Fig. 11. The measure of n-bounce circular light paths $p_n$ shown for n up to 16 on an example scene (inset). The scene is composed asymmetrically of 3 Lambertian squares (albedo 0.5) so that $\kappa$ is bounded, making $T_b$ compact.</figcaption></figure>

**THEOREM 6.4.1:** Any eigenfunction $\Lambda_i$ of $\tilde{T}$ corresponding to an eigenvalue $\lambda_i$ with unit geometric multiplicity can be expressed as a path integral

$$\Lambda_i = -\frac{1}{2i\pi} \oint_{\Gamma} L(\alpha) d\alpha, \qquad (32)$$

where $L(\alpha)$ is the solution to the “generalized light transport problem”

$$\alpha L = E + \tilde{T} L, \qquad (33)$$

and $\Gamma \subset \mathbb{C}$ is any simple closed curve enclosing the sole eigenvalue $\lambda_i$, and $E$ any light distribution such that $\langle \Lambda_i, E \rangle \neq 0$.

**PROOF:** The explanation to this intriguing theorem relies on the fact that because $L(\alpha) = (\alpha - \tilde{T})^{-1}E = R(\tilde{T}; \alpha)E$, Eq. 32 corresponds to the Riesz integral of the resolvent, that is known to be a projection onto the union of eigenspaces of the operator which corresponding eigenvalues are enclosed by $\Gamma$ (See for instance Gohberg [1978] ch.I). The integrated “light” distribution in Equation 32 will therefore be equal

<footnote n="10">¹⁰Note that implementation efficiency is not an issue at this point, and we acknowledge that better sampling strategies that reuse the same path for multiple $p_i$ are available</footnote>
to an (unnormalized) eigenfunction of $\tilde{\mathbf{T}}$, or 0 depending on whether $\langle E, \Lambda_i \rangle \neq 0$ or not. $\square$

Equation 33 can be viewed as a generalization of the traditional light transport equation 10 by introducing a parameter $\alpha$. While the existence of a solution to this equation is guaranteed by the Fredholm alternative theorem [Zemyan 2012] (since $\alpha \in \Gamma \subset \rho(\tilde{\mathbf{T}})$), its solution need not be positive nor even real. Thus, it would not represent a physically meaningful distribution of light. Nevertheless, this is helpful to effectively compute an eigenfunction of $\tilde{\mathbf{T}}$.

A practical, yet non computationally efficient way to compute $\Lambda_i$ is to apply a Galerkin approximation to $\tilde{\mathbf{T}}$ and solve for $(\alpha - \tilde{\mathbf{T}})L(\alpha) = E$ using a biconjugate gradient. Computing this solution using the traditional Monte-Carlo light transport paradigm is however not feasible, since the series obtained by expanding $(\alpha - \tilde{\mathbf{T}})^{-1}$ only converges for $|\alpha| > \|\tilde{\mathbf{T}}\|$. This precludes any curve $\Gamma$ that would only encircle a strict subset of the point-spectrum of $\tilde{\mathbf{T}}$.

# 7 Practical computation of the spectrum

In this section, we suggest a couple of practical methods to compute the eigenvalues of operators $\mathbf{T}_b$ and $\mathbf{T}$.

## 7.1 Lambertian operator $\mathbf{T}_b$: Galerkin is best

Following the work of Chatelin [1981], since the point-spectrum of $\mathbf{T}_b$ is countable, a converging approximation of the eigenvalues of $\mathbf{T}_b$ can be obtained by computing the eigenvalues of an $n$-elements Galerkin approximation with matrix $\mathbf{A}_n$ of that operator. The scene is partitioned into $n$ surfaces of area $S_i$ and the light transport operator is projected onto a finite element basis. For a basis of piecewise constant functions, the matrix entry $a_n^{i,j}$ of the projected transport operator is [Sil- lion and Puech 1994]

$$a_n^{i,j} = \frac{1}{S_i} \int_{S_i} \int_{S_j} \kappa_b(\mathbf{x}_i, \mathbf{x}_j) \, d\mathbf{x}_i \, d\mathbf{x}_j \quad (34)$$

where $S_i$ the area of surface element $i$ and $\kappa_b$ the kernel of $\mathbf{T}_b$ (See Section 3.3). The eigenvalues of the matrix can then be calculated using a standard eigensolver.

## 7.2 Non-Lambertian operator $\mathbf{T}$: Monte-Carlo approach

In non-Lambertian scenes, a Galerkin approach would require discretization of surfaces and directions so that Equation 34 integrates over pairs of surface element-direction pairs $\ell_i = S_k \times \Omega_l$

$$a_n^{ij} = \frac{1}{\mu(\ell_i)} \int_{\ell_i} \int_{\ell_j} \kappa(\mathbf{x}, \omega, \mathbf{y}, \omega'') \mathrm{d}\mathbf{x}\mathrm{d}\omega\mathrm{d}\mathbf{y}\mathrm{d}\omega'' \quad (35)$$

where $\mu$ is the measure associated with $\ell$. Unfortunately, this approach is computationally impractical. Compared to the Lambertian case where the eigensolver would be applied on a matrix that is about 1000 × 1000, the non-Lambertian case for comparable accuracy would be $10^6 \times 10^6$ or larger. The matrix in Eq. 35 is very sparse, and even more when reflectance

functions are highly glossy. Computing $a_{ij}$ is however expensive because the correlation between the geometric factor in $\kappa$ and the BRDF cannot be neglected while maintaining physical consistency.

Since Equation 35 is an integral, $a_n^{i,j}$ is the expected value of an appropriately scaled random variable:

$$a_n^{ij} = \mu(\ell_i) \mathbb{E}(\kappa(\mathbf{x}, \omega, \mathbf{y}, \omega'')).$$

Consequently we form the following random matrix:

$$\mathbf{A}_n = \frac{2\pi S}{n} \left\| \kappa(\mathbf{x}_i, \omega_i, \mathbf{y}_j, \omega_j'') \right\|_{ij},$$

where points and directions are randomly sampled. In practice, we suggest to select $n$ random points $\mathbf{x}_i$, and outgoing directions for each point are implicitly defined by the connection to all other points at the cost of a bias. This also ties sparsity to the glossiness of reflectance functions only, since in this formulation only the directional integration is implicitly considered. For compact approximations $\mathbf{T}_e$ of $\mathbf{T}$ (see Section 4) the eigenvalues of $\mathbf{A}_n$ form an estimator of the eigenvalues of $\mathbf{T}_e$ [Chatelin 1981]. The estimator is biased since the eigenvalues of $\mathbf{A}_n$ do not linearly depend on its coefficients, but the bias converges to 0 when $n$ tends to infinity. Note that because of the lack of a proof that $\sigma_p(\mathbf{T})$ is actually countable and $\mathbf{T}_e$ strongly converges to $\mathbf{T}$, we can only ensure that this method computes converging approximations of the isolated elements of $\sigma_p(\mathbf{T})$ (See Section 4).

We use an eigensolver that only requires computing the product of vectors by $\mathbf{A}_n$ or its transpose [Lehoucq et al. 1998], which is performed efficiently as a triple loop over the $\mathbf{x}_i$. Figure 12 shows an example of using this technique on a non-Lambertian scene. The method presented in this section shares similarity with random submatrix sampling as proposed by Frieze [2004] to compute singular values of large matrices.

# 8 Discussion

## 8.1 Discussion of results

*Challenges specific to light transport operators.* Despite the large array of mathematical tools to analyse the spectra of linear operators, light transport operators $\mathbf{T}$ and $\mathbf{T}_b$ present unique challenges. Although their non-compactness can be mitigated by approximations, these approximations inherit some of the bad properties of the operators: non normality (or any sub-normality class), discontinuous integration kernel, non-compact resolvent. These severely limit the set of applicable methods. Although they are not positive operators or symmetric, the fact that they are integral operators with a strictly positive kernel allows some of the analyses presented.

*Relationship between T and its adjoint.* A curious result from Theorem 5.4.1 is that the adjoint $\mathbf{T}^*$ has the same point spectrum (eigenvalues) as $\mathbf{T}$. The adjoint operator has a practical significance in light transport, as it represents the transport of importance or tracking light backwards from the eye to the scene. Our result confirms that the two operators
<figure><img src="image_14.png" /><figcaption>Fig. 12. Left: validation of the Monte-Carlo estimation of eigenvalues of Section 7.2 on the two-squares scene of Figure 6 (Lambertian version with albedo 0.8). When increasing the number of points n, eigenvalue histograms show sharper peaks that exactly align with the reference eigenvalues. Middle: Monte-Carlo estimation of eigenvalues in the CornellBox scene. Right: 3D histogram in the complex plane of Monte-carlo estimation of eigenvalues in the non-Lambertian scene of Figure 6.</figcaption></figure>

share similar spectral properties and therefore similar computational challenges with respect to convergence. Although $\mathbf{T}$ is not self-adjoint, $\mathbf{T} = \mathbf{KG}$ where $\mathbf{K}$ and $\mathbf{G}$ are self-adjoint operators. As a result, $\mathbf{T}$ and $\mathbf{T}^*$ have mutually orthogonal eigenfunctions and operators $\mathbf{K}$ and $\mathbf{G}$ can be used to transform one into the other and vice versa. This could potentially have practical implications in precomputed light transport where one of them is diagonalized.

**Relationship between $\mathbf{T}$ and $\mathbf{T}_b$ for diffuse scenes.** In Lambertian scenes, any eigenfunction $\Lambda_b$ of $\mathbf{T}_b$ can be turned into an eigenfunction of $\mathbf{T}$ as $\Lambda(\mathbf{x}, \omega) = \Lambda_b(\mathbf{x})$. Conversely, any eigenfunction of $\mathbf{T}$ with a non-zero eigenvalue will have a constant directional component (by means of $\mathbf{T}\Lambda = \lambda\Lambda$) and is similarly tied to an eigenfunction of $\mathbf{T}_b$. Additionally, 0 is an eigenvalue of $\mathbf{T}$ with infinite dimensionality, corresponding to all directional functions of mean value zero.

Generally $\mathbf{T}$ is not compact as a partial integral operator, it is therefore counterintuitive that it may have the pointspectrum of a compact operator plus eigenvalue 0 ($\mathbf{T}_b$ is compact when the scene has no abutting edges). In Lambertian scenes however $\mathbf{T}$ can be restricted to the space of directionally constant functions for which a single value represent the whole function, which in effect makes $\mathbf{T}$ a non partial integral operator. Note that in this case, the eigenspaces of $\mathbf{T}$ (including eigenvalue 0) are complete and its eigenvalues are countable.

**Specific case of $\rho_{ij} \ge 0$.** An interesting observation is that because of Equation 21, when all the $\rho_{ij}$ are positive, $\mathbf{T}$ can be left-symmetrized by the operator

$$ \mathbf{R}: \mathbf{T}(\mathcal{H}) \rightarrow \mathbf{T}(\mathcal{H}) \\ \mathbf{RTL} = \sum_{i=1}^{n} \sum_{j>0} \langle L, \varphi_j^i \rangle \psi_j^i, \qquad (36) $$

which in this case makes $\mathbf{RT}$ self-adjoint (the proof involves showing that $\langle \varphi_j^i, \psi_l^k \rangle = \langle \varphi_l^k, \psi_j^i \rangle$ which only happens when $\varphi_j^i$ does not include the negative sign of some of the $\rho_{ij}$). Since $\mathbf{R}$ has eigenvalues $|\rho_{ij}|^{-1}$, it is bounded below by 1 (which means it has a bounded inverse) and we can in this case apply Th. 9.1 of Silberstein [1962] to deduce that $\mathbf{T}$ is

symmetrizable over $\mathbf{T}(\mathcal{H})$. Consequently its eigenvalues are real, and its eigenfunctions are orthogonal for the dot product $\langle f, g \rangle_R = \langle \sqrt{\mathbf{R}}f, \sqrt{\mathbf{R}}g \rangle$ [Silberstein 1962]. Furthermore orthogonality implies countability because $\mathbf{T}(\mathcal{H}) \subset \mathcal{H}$ is separable. Since $\mathbf{G}$ is invertible on $\mathcal{H}$, we know that $\mathbf{K}(\mathcal{H}) = \mathbf{T}(\mathcal{H})$ and $\ker(\mathbf{T}) = \mathbf{G}^{-1}(\ker(\mathbf{K}))$. Since $\mathbf{K}$ is self-adjoint, its eigenfunctions (including those in $\ker(\mathbf{K})$ span $\mathbf{H}$, and therefore eigenfunctions of $\mathbf{T}$ span $\mathcal{H}$, proving a strong result for the spectrum of $\mathbf{T}$.

This happens at least in Lambertian scenes, where the only eigenvalues of $\mathbf{K}$ correspond to the albedo of all materials in the scene. In classical scenes, it seems that BRDFs have both positive and negative eigenvalues (See Figure 8). Although the effect of flipping the sign of negative eigenvalues of BRDFs is unclear, this would preserve the frequency content of material reflectances [Durand et al. 2005]. Perhaps this case would be useful in the context of analysis of light fields via a symmetrizable operator.

**Estimating eigenfunctions from $\alpha_k$.** The result from Theorem 6.2.1 means that if $L$ belongs to the span of a finite set of eigenvectors of a compact approximation of $\mathbf{T}$ or $\mathbf{T}_b$ the roots of polynomial $P$ given by the theorem gives the corresponding eigenvalues. Given one such $\lambda_i$, one can also expand $Q(x) = P(x)/(x - \lambda_i)$ and obtain $Q(\mathbf{T})L = \beta_i\Lambda_i$. While a random $L$ will generally not be in the span of a finite set of eigenfunctions of $\mathbf{T}$, the above expression only approximate an eigenfunction, assuming that eigenvalues decrease fast enough for the contribution of $\{\Lambda_i\}_{i>n}$ to be small enough.

**Total response of the generalized transport system.** The measures of circular light paths $p_k$ defined in Section 6.3 can be used as a measure of the total response of the generalized transport equation $\alpha L = E + \mathbf{T}L$ to border conditions ($E$ in this expression). This quantity is generally defined as

$$ d = \int_{\mathcal{L}} r_{\alpha}(\mathbf{x}, \mathbf{x}) \, d\mathbf{x}, $$

where $r_\alpha$ is the integral kernel [Brislawn 1988] of the resolvent operator $\mathbf{R}(\mathbf{T}; \alpha)$. Recalling the calculation leading to Eq. 31, the resolvent of $\tilde{\mathbf{T}}$ is clearly not trace class due to the terms
αI + ˜T. Dropping the first two terms still allows to define a regularized total response ˜d, which turns into

$$ \tilde{d} = \alpha^{-1} p_2 + \alpha^{-2} p_3 + \alpha^{-3} p_4 + \dots \quad (37) $$

This gives a formal proof of the intuitive following fact: for compact approximations of T, the larger the measure of circular light paths, the more energy will be in the solution of the generalized equation.

## 8.2 Open questions

**Eigenfunctions of T.** The spectrum of T diagonalizes the operator and so estimation of the eigenvalues and eigenfunctions of T would simplify the calculation of higher powers (multiple-bounce transport), inversion (relighting, inverse rendering), and so on. In this article, we have analyzed the spectra (eigenvalues and eigenfunctions) mathematically and estimated eigenvalues numerically. However, stable methods for the estimation of the eigenfunctions of T remain elusive.

**Quantifying scene complexity using $p_n$.** The complexity of a scene with respect to computing light transport solutions is difficult to quantify. There are a variety of factors that make computation challenging, including number of light sources, reflectances in the scene, visibility, geometric elements, and so on. We conjecture that the measure of circular light paths in a scene may be useful for quantitative assessment of the interplay between these factors. As explained in Section 6.3, $p_n$ is a single number per length of circular paths $n$ obtained as the trace of $T^n$. Its distribution over $n$ encodes information about the relative contributions of different path lengths to the light transport solution. Thus, scenes where $p_n$ has a heavy or long tail may be more challenging to solve accurately than scenes where $p_n$ has a sharp peak for lower values of $n$.

**Is the point-spectrum of T countable?** Despite our efforts, we were unable to derive a formal proof for this challenging question. Since $\mathcal{H}$ is a separable Hilbert space, completeness of the eigenfunctions of T would prove countability of its point-spectrum. However T does not meet the conditions for Gohberg theorems on completeness of eigenfunctions ([Gohberg and Krein 1978] Chp.V), nor is it dissipative, quasinormal or hyponormal [Gohberg and Krein 1978]. Exploiting the characteristics of the convergence of $T_\epsilon$ toward T through asymptotic perturbation theory may be an interesting direction to consider (See Chatelin [2011], and Kato [1995] Chp.8).

**New class of light transport algorithms?** Solving the generalized light transport equation, $\lambda L = E + TL$ for $|\lambda| < \|T\|$, with path tracing could be an important step toward practical calculation of the eigenfunctions of T, using the Riesz integral in Equation 32. As explained in Section 6.4, a generic formulation of its solution as an integral over light paths remains elusive. Such a formulation may reveal fresh perspective of the rendering equation and potentially lead to a new class of light transport algorithms.

## 9 Conclusion and future work

In this article we have analyzed the spectral characteristics of light transport operators $K, K_x, G, T_b$ and $T$. We demonstrated that, although neither $T_b$ nor $T$ are compact, these two operators can still be strongly approximated by compact operators. Further, we highlight that these approximations actually correspond to the ones used in practice while implementing path tracing. We have also revealed various connections between the point spectrum of these operators and path integrals in a scene, combining circular light paths, powers of the operator and projection on the operators' output space. We identify important questions that remain unanswered and leave them as future work in this area.

In conclusion, we believe that there is much to be discovered at the confluence of stochastic approximations of large matrices [Dimov et al. 2015; Frieze et al. 2004; Halko et al. 2011; Kobayashi et al. 2001] and light transport estimation.

### Acknowledgments

Kartic Subr was supported by a Royal Society University Research Fellowship.

### A Self-Adjointness of G and K

Although $K$ and $G$ are long known to be self-adjoint [Veach 1997], we derive this using our notations. For two points $x$ and $y$ that are visible to each other (via local directions $\omega'$ and $\omega''$ as in Figure 2),

$$ \langle GL_1, L_2 \rangle = \int_S \int_{\Omega} (GL_1)(x, \omega') \overline{L_2}(x, \omega') \cos \theta' \, dx \, d\omega'. $$

Now we change variables from $(x, \omega')$ into $(y, \omega'')$. We have $d\omega' = dy \cos \theta''/r^2$ and $d\omega'' = dx \cos \theta'/r^2$, from which we get $dy \cos \theta''d\omega'' = dx \cos \theta'd\omega'$. Using the fact that $(GL_1)(x, \omega') = L_1(y, \omega'')$ and $\overline{L_2}(x, \omega') = (GL_2)(y, \omega'') = \overline{GL_2}(y, \omega'')$, the above equation turns into

$$ \langle GL_1, L_2 \rangle = \int_S \int_{\Omega} L_1(y, \omega'') \overline{GL_2}(y, \omega'') \cos \theta'' \, dy \, d\omega'' \\ = \langle L_1, GL_2 \rangle. $$

Similarly,

$$ \langle KL_1, L_2 \rangle = \int_S \int_{\Omega} (KL_1)(x, \omega) \overline{L_2}(x, \omega) \cos \theta \, dx \, d\omega \\ = \int_S \int_{\Omega} \left( \int_{\Omega} \rho(x, \omega, \omega') L_1(x, \omega') \cos \theta' \, d\omega' \right) \overline{L_2}(x, \omega) \cos \theta \, dx \, d\omega $$

Since $\rho(x, \omega, \omega') = \rho(x, \omega', \omega)$ due to Helmholtz reciprocity and $KL_2 = \overline{KL_2}$, the two definite integrals over $\Omega$ can be exchanged to yield

$$ \langle KL_1, L_2 \rangle = \langle L_1, KL_2 \rangle. $$

Both operators are self-adjoint thanks to the cosine in the dot product. □
## B Strong convergence of $T_\epsilon$ and $T_{b,\epsilon}$ in $L_1$ norm

Let $L$ be a bounded light distribution. For any $\epsilon > 0$ we denote by $S_\epsilon(\mathbf{x})$ the set of points $\mathbf{y}$ for which $\kappa_b(\mathbf{x}, \mathbf{y}) > \frac{1}{\epsilon}$, and $M_\epsilon$ the set

$$M_\epsilon = \left\{ \mathbf{x} \in S \;\exists \mathbf{y} \in S \; \kappa(\mathbf{x}, \mathbf{y}) \ge \frac{1}{\epsilon} \right\}.$$

Clearly for all $\mathbf{x} \in S$ we have $S_\epsilon(\mathbf{x}) \subset M_\epsilon$. Using successively $S_\epsilon \subset S$, Hölder's inequality and the non-negativity of $\kappa_b$ we get

$$ \begin{aligned} |((T_b - T_{b,\epsilon})L)(\mathbf{x})| &\le \int_{S_\epsilon(\mathbf{x})} \left| (\kappa_b(\mathbf{x}, \mathbf{y}) - \frac{1}{\epsilon})L(\mathbf{y}) \right| \mathrm{d}\mathbf{y} \\ &\le \|L\|_{\infty} \int_{S_\epsilon(\mathbf{x})} \left| \kappa_b(\mathbf{x}, \mathbf{y}) - \frac{1}{\epsilon} \right| \mathrm{d}\mathbf{y} \\ &\le \|L\|_{\infty} \int_{S_\epsilon(\mathbf{x})} \kappa_b(\mathbf{x}, \mathbf{y}) \mathrm{d}\mathbf{y}. \\ &\le \|L\|_{\infty} \int_{M_\epsilon} \kappa_b(\mathbf{x}, \mathbf{y}) \mathrm{d}\mathbf{y}. \end{aligned} \quad (38) $$

Therefore

$$ \begin{aligned} \|(T_b - T_{b,\epsilon})L\|_1 &\le \|L\|_{\infty} \int_S \int_{M_\epsilon} \kappa_b(\mathbf{x}, \mathbf{y}) \mathrm{d}\mathbf{y} \mathrm{d}\mathbf{x} \\ &\le \|L\|_{\infty} \int_{M_\epsilon} \underbrace{\int_S \kappa_b(\mathbf{x}, \mathbf{y}) \mathrm{d}\mathbf{x}}_{\le 1} \mathrm{d}\mathbf{y} \\ &\le \|L\|_{\infty} \mu(M_\epsilon) \end{aligned} $$

In classical physically meaningful scenes, where edges and corners constitute a non measurable set of points, the measure $\mu(M_\epsilon)$ verifies $\lim_{\epsilon \to 0} \mu(M_\epsilon) = 0$, therefore

$$\lim_{\epsilon \to 0} \|(T_b - T_{b,\epsilon})L\|_1 = 0.$$

The strong convergence of $T_\epsilon$ to $T$ is immediate after noticing that for any directional function $\omega \mapsto h(\omega)$ we have whenever $h$ is continuous at $\omega$,

$$\lim_{\epsilon \to 0} \int_{\Omega} h(\omega')g_{\epsilon}(\omega \cdot \omega') \mathrm{d}\omega' = h(\omega). \quad (39)$$

Once again, the "proper" geometry of scenes ensuring that light distributions are only non continuous at a set of points which measure is 0, and applying the preceding proof for the remaining term $f_\epsilon(\mathbf{x}, \mathbf{y})$, we deduce the strong convergence of $T_\epsilon$ to $T$ almost everywhere. $\square$

## C Eigenfunctions increase in frequency

In this proof $T$ denotes a compact approximation of one of the transport operators. Let $\{\varphi_n\}_{n \ge 0}$ be an orthogonal basis of the space with non decreasing frequency content (e.g. directional/spatial wavelets). Let $\{\Lambda_k\}_{k \ge 0}$ be the eigenfunctions of $T$ assumed to form a complete system.

Assume, first, that eigenfunctions are orthogonal. Then, for all $n$

$$\varphi_n = \sum_{k=0}^{\infty} \alpha_{nk} \Lambda_k \quad \text{with} \quad \alpha_{nk} = \langle \varphi_n, \Lambda_k \rangle.$$

The convergence of that sum implies that

$$\forall n \in \mathbb{N} \quad \lim_{k \to \infty} \alpha_{nk}^2 = 0.$$

Selecting the first $N+1$ values of $n$ allows us to apply the limit to the supremum of the $\alpha_{nk}$:

$$\forall N \in \mathbb{N} \quad \forall n \in [0, N] \quad \limsup_{k \to \infty} \sup_{n \le N} \alpha_{nk}^2 = 0,$$

which exactly states that when $k$ increases, the coefficients of $\Lambda_k$ over $\{\varphi_n\}$ uniformly converge to 0, over every interval $[0, N]$. Since $\sum_n \alpha_{nk}^2 = 1$, that means most of the energy in $\Lambda_k$ lies beyond $N$.

Now if eigenfunctions are not orthogonal, we can still apply Gram-Schmidt orthogonalization to $\{\Lambda_k\}$ (which remains valid for countably infinite sequences) to obtain an orthogonal sequence $\{\Lambda_k'\}_{k \ge 0}$ such that

$$\varphi_n = \sum_{k=0}^{\infty} \alpha'_{nk} \Lambda'_k \quad \text{with} \quad \alpha'_{nk} = \langle \varphi_n, \Lambda'_k \rangle.$$

Therefore, $\Lambda'_k$ increases in frequency with $k$. But since $\Lambda'_k$ is built using Gram-Schmidt orthogonalization, we know that for all $k$, $\Lambda_k$ is a linear combination of $\Lambda'_0, \dots, \Lambda'_k$, and thus the increase in frequency applies to $\Lambda_k$. $\square$

## D Proof of Equation 17

For any $N \in \mathbb{N}$ we call $\mathbf{P}_N : \mathcal{H} \to \mathcal{H}$ the projection on the finite dimensional space spanned by $\{\varphi_n\}_{n<N}$, and we define the finite rank operator $\mathbf{A}_n$ to be

$$\mathbf{A}_n = \mathbf{A} \mathbf{P}_n.$$

Since $\mathbf{A}$ is compact and the sequence $\{\varphi_n\}_{n>0}$ is complete, we know that $\mathbf{A}_n$ converges to $\mathbf{A}$ in the operator norm [Chatelin 1981]:

$$\forall \epsilon > 0 \quad \exists p \in \mathbb{N} \quad \forall n \ge p \quad \|\mathbf{A}_n - \mathbf{A}\| < \epsilon.$$

That means

$$\forall \epsilon > 0 \quad \exists p \in \mathbb{N} \quad \forall n \ge p \quad \forall L \in \mathcal{A} \quad \|(\mathbf{A}_n - \mathbf{A})L\| < \epsilon \|L\|.$$

We now restrict the above equation to $L = \varphi_{n+1}$ for every possible $n > p$. We have in this case $(\mathbf{A} - \mathbf{A}_n)\varphi_{n+1} = \mathbf{A}\varphi_{n+1}$ and $\|\varphi_{n+1}\| = 1$, which gives after shifting $n+1$ to $n$:

$$\forall \epsilon > 0 \quad \exists p \in \mathbb{N} \quad \forall n > p \quad \|\mathbf{A}\varphi_n\| < \epsilon. \quad \square$$

## E Proof of Equation 29

Let $\mathbf{A}_n$ be the matrix in Eq. 30 and let $c_n = \det(\mathbf{A}_n)$. In order to prove Eq. 29 one needs to notice that when expanding this determinant along the first column, the $k$-th minor $m_k$ is

$$m_k = \left| \frac{\mathbf{D}_{n,k}}{\mathbf{B}_{n,k}} \frac{\mathbf{D}_{n,k}}{\mathbf{A}_{n-k}} \right|,$$

where $\mathbf{D}_{n,k}$ is a $(k-1) \times (k-1)$ lower-triangular matrix with $n-1, n-2, \dots, n-k+1$ on its diagonal, making the bottom-left submatrix $\mathbf{B}_{n,k}$ irrelevant in this calculation. Following this, we get

$$c_n = \sum_{k=2}^{n} (-1)^{k+1} p_k m_k = \sum_{k=2}^{n} (-1)^{k+1} \frac{(n-1)!}{(n-k)!} p_k c_{n-k}$$
Substituting $n! a_n = c_n$ directly leads to Eq. 29.

## References

James Arvo. 1996. *Analytic Methods for Simulated Light Transport*. Ph. D. Dissertation. Yale University.

Ian Ashdown. 2001. *Eigenvector radiosity*. Ph. D. Dissertation. University of British Columbia. https://doi.org/10.14288/1.0051491

Jiamin Bai, Manmohan Chandraker, Tian-Tsong Ng, and Ravi Ramamoorthi. 2010. A Dual Theory of Inverse and Forward Light Transport. In *Computer Vision - ECCV 2010 (Lecture Notes in Computer Science)*, Kostas Daniilidis, Petros Maragos, and Nikos Paragios (Eds.). Springer, Berlin, Heidelberg, 294–307. https://doi.org/10.1007/978-3-642-15552-9_22

Gladimir V. G. Baranoski, Randall Bramley, and Jon G. Rokne. 1997. Eigen-Analysis for Radiosity Systems. In *Proceedings of the Sixth International Conference on Computational Graphics and Visualization Techniques*. Vilamoura, Algarve, Portugal.

Simon Barry. 2000. *Trace Ideals and Their Applications: Second Edition* (second edition ed.). Mathematical Surveys and Monographs, Vol. 120. American Mathematical Society. https://bookstore.ams.org/surv-120-s

Laurent Belcour, Thomas Deliot, Wilhem Barbier, and Cyril Soler. 2022. A Data-Driven Paradigm for Precomputed Radiance Transfer. Proc. *ACM Comput. Graph. Interact. Tech.* 5, 3 (July 2022), 26:1–26:15. https://doi.org/10.1145/3543864

Daniele Boffi. 2010. Finite element approximation of eigenvalue problems. *Acta Numerica* 19 (May 2010), 1–120. https://doi.org/10.1017/S0962492910000012 Publisher: Cambridge University Press.

Chris Brislawn. 1988. *Kernels of Trace Class Operators*. Proc. Amer. Math. Soc. 104, 4 (1988), 1181–1190. https://doi.org/10.2307/2047610 Publisher: American Mathematical Society.

F. Chatelin. 1981. The Spectral Approximation of Linear Operators with Applications to the Computation of Eigenelements of Differen- tial and Integral Operators. *SIAM Rev.* 23, 4 (Oct. 1981), 495–522. https://doi.org/10.1137/1023099

Françoise Chatelin. 2011. *Spectral Approximation of Linear Operators*. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611970678

John Conway. 1991. *The Theory of Subnormal Operators*. Mathematical Surveys and Monographs, Vol. 36. American Mathematical Society, Providence, Rhode Island. https://doi.org/10.1090/surv/036

J. B. Conway. 1990. *A Course in Functional Analysis* (2 ed.). Springer-Verlag, New York. https://www.springer.com/gp/book/9780387972459

Daniel K. Crane and Mark S. Gockenbach. 2020. The Singular Value Expansion for Arbitrary Bounded Linear Operators. *Mathematics* 8, 8 (Aug. 2020), 1346. https://doi.org/10.3390/math8081346 Number: 8 Publisher: Multidisciplinary Digital Publishing Institute.

Ivan Dimov, Sylvain Maire, and Jean Michel Sellier. 2015. A new Walk on Equations Monte Carlo method for solving systems of linear algebraic equations. *Applied Mathematical Modelling* 39, 15 (Aug. 2015), 4494–4510. https://doi.org/10.1016/j.apm.2014.12.018

Frédéric Durand, Nicolas Holzschuch, Cyril Soler, Eric Chan, and François X. Sillion. 2005. A Frequency Analysis of Light Transport. In *ACM SIGGRAPH 2005 Papers (SIGGRAPH '05)*. ACM, New York, NY, USA, 1115–1126. https://doi.org/10.1145/1186822.1073320

Emma D'Aniello and Martina Maiuriello. 2022. On the spectrum of weighted shifts. *Revista de la Real Academia de Ciencias Exactas, Físicas y Naturales. Serie A. Matemáticas* 117, 1 (Oct. 2022), 4. https://doi.org/10.1007/s13398-022-01328-z

Alan Frieze, Ravi Kannan, and Santosh Vempala. 2004. Fast monte-carlo algorithms for finding low-rank approximations. *J. ACM* 51, 6 (Nov. 2004), 1025–1041. https://doi.org/10.1145/1039488.1039494

Gaurav Garg, Eino-Ville Talvala, Marc Levoy, and Hendrik P. Lensch. 2006. Symmetric Photography: Exploiting Data-sparseness in Reflectance Fields. In *Proceedings of the 17th Eurographics Conference on Rendering Techniques (EGSR '06)*. Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 251–262. https://doi.org/10.2312/EGWR/EGSR06/251-262

Iliyan Georgiev, Jaroslav Krivánek, Tomáš Davidovič, and Philipp Slusallek. 2012. Light transport simulation with vertex connection and merging. *ACM Trans. Graph.* 31, 6 (Nov. 2012), 192:1–192:10. https://doi.org/10.1145/2366145.2366211

Israel Gohberg, Seymour Goldberg, and Nahum Krupnik. 2012. *Traces and Determinants of Linear Operators*. Birkhäuser. Google-Books-ID: 6DMDCAAAQBAJ.

Israel Gohberg and Mark Grigor evich Kreïn. 1978. *Introduction to the Theory of Linear Nonselfadjoint Operators*. American Mathematical Soc. Google-Books-ID: GOJFyTscyrkC.

N. Halko, P. G. Martinsson, and J. A. Tropp. 2011. Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. *SIAM Rev.* 53, 2 (2011), 217–288. https://doi.org/10.1137/090771806

Miloš Hašan, Jaroslav Krivánek, Bruce Walter, and Kavita Bala. 2009. Virtual spherical lights for many-light rendering of glossy scenes. *ACM Transactions on Graphics* 28, 5 (Dec. 2009), 1–6. https://doi.org/10.1145/1618452.1618489

Miloš Hašan, Fabio Pellacini, and Kavita Bala. 2007. Matrix Row-column Sampling for the Many-light Problem. In *ACM SIGGRAPH 2007 Papers (SIGGRAPH '07, Vol. 26)*. ACM, New York, NY, USA, 26.1–26.10. https://doi.org/10.1145/1275808.1276410

Fu-Chung Huang and Ravi Ramamoorthi. 2010. Sparsely Precomputing the Light Transport Matrix for Real-Time Rendering. *Comput. Graph. Forum* 29, 4 (2010), 1335–1345. https://doi.org/10.1111/j.1467-8659.2010.01729.x

Aref Jeribi. 2021. *Perturbation Theory for Linear Operators: Density and Bases with Applications*. Springer, Singapore. https://doi.org/10.1007/978-981-16-2528-2

James T. Kajiya. 1986. The rendering equation. *ACM SIGGRAPH Computer Graphics* 20, 4 (Aug. 1986), 143–150. https://doi.org/10.1145/15886.15902

Gerhard K. Kalisch. 1972. On Operators on Separable Banach Spaces with Arbitrary Prescribed Point Spectrum. *Proc. Amer. Math. Soc.* 34, 1 (1972), 207–208. https://doi.org/10.2307/2037930 Publisher: American Mathematical Society.

A. Kalitvin and P.P. Zabrejko. 1991. On the Theory of Partial Integral Operators. *Journal of Integral Equations and Applications* 3, 3 (Sept. 1991), 351–382. https://doi.org/10.1216/jiea/1181075630

Tosio Kato. 1995. *Perturbation Theory for Linear Operators*. Springer Science & Business Media. Google-Books-ID: 8j12kN_D3BwC.

Alexander Keller. 1997. Instant radiosity. In *Proceedings of the 24th annual conference on Computer graphics and interactive techniques (SIGGRAPH '97)*. ACM Press/Addison-Wesley Publishing Co., USA, 49–56. https://doi.org/10.1145/258734.258769

Mei Kobayashi, Georges Dupret, Oliver King, and Hikaru Samukawa. 2001. Estimation of Singular Values of Very Large Matrices using Random Sampling. *Computer & Mathematics with Applications* 42, 10-11 (2001), 1331–1352.

Thomas Kollig and Alexander Keller. 2006. Illumination in the Presence of Weak Singularities. In *Monte Carlo and Quasi-Monte Carlo Methods 2004*, Harald Niederreiter and Denis Talay (Eds.). Springer, Berlin, Heidelberg, 245–257. https://doi.org/10.1007/3-540-31186-6_15

Peter D. Lax. 2002. *Functional Analysis*. Wiley-Interscience. Open Library ID: OL7618836M.

Richard B. Lehoucq, Danny C. Sorensen, and Chao Yang. 1998. *ARPACK Users' Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods*. SIAM. Google-Books-ID: iMUEa23N_CQC.

Christian Lessig and Eugene Fiume. 2010. On the Effective Dimension of Light Transport. In *Proceedings of the 21st Eurographics Conference on Rendering (EGSR'10)*. Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 1399–1403. https://doi.org/10.1111/j.1467-8659.2010.01736.x

Bradford J. Loos, Lakulish Antani, Kenny Mitchell, Derek Nowrouzезahrai, Wojciech Jarosz, and Peter-Pike Sloan. 2011. Modular Radiance Transfer. In *Proceedings of the 2011 SIGGRAPH Asia Conference (SA '11)*. ACM, New York, NY, USA, 178:1–178:10. https://doi.org/10.1145/2024156.2024212

Manabu Machida. 2014. Singular eigenfunctions for the three-dimensional radiative transport equation. *JOSA A* 31, 1 (Jan. 2014), 67–74. https://doi.org/10.1364/JOSA.A.31.000067

Dhruv Mahajan, Ira Kemelmacher Shlizerman, Ravi Ramamoorthi, and Peter Belhumeur. 2007. A Theory of Locally Low Dimensional Light Transport. In *ACM SIGGRAPH 2007 Papers (SIGGRAPH '07, Vol. 26-3)*. ACM, New York, NY, USA, 62.1–62.10. https://doi.org/10.1145/1275808.1276454

Mohammed Hichem Mortad. 2022. *Counterexamples in Operator Theory*. Birkhauser. https://link.springer.com/book/10.1007/978-1-4506-0173-6
3-030-97814-3

Tian-Tsong Ng, Ramanpreet Singh Pahwa, Jiamin Bai, Kar-Han Tan, and Ravi Ramamoorthi. 2012. From the Rendering Equation to Stratified Light Transport Inversion. *Int. J. Comput. Vision* 96, 2 (Jan. 2012), 235–251. https://doi.org/10.1007/s11263-011-0467-6

Derek Nowrouzezahrai, Patricio Simari, Evangelos Kalogerakis, and Eugene Fiume. 2007. Eigentransport for Efficient and Accurate All-frequency Relighting. In *Proceedings of the 5th International Conference on Computer Graphics and Interactive Techniques in Australia and Southeast Asia (GRAPHITE '07)*. ACM, New York, NY, USA, 163–169. https://doi.org/10.1145/1321261.1321290

Jiawei Ou and Fabio Pellacini. 2011. LightSlice: Matrix Slice Sampling for the Many-lights Problem. In *Proceedings of the 2011 SIGGRAPH Asia Conference (SA '11)*. ACM, New York, NY, USA, 179:1–179:8. https://doi.org/10.1145/2024156.2024213

Pieter Peers, Dhruv K. Mahajan, Bruce Lamond, Abhijeet Ghosh, Wojciech Matusik, Ravi Ramamoorthi, and Paul Debevec. 2009. Compressive light transport sensing. *ACM Trans. Graph*, 28, 1 (Feb. 2009), 3:1–3:18. https://doi.org/10.1145/1477926.1477929

Vũ Ngoc Phát and Trinh Cong Dieu. 1994. On the Krein-Rutman theorem and its applications to controllability. In *Proceedings of the American Mathematical Society*, Vol. 120. American Mathematical Society, 495–500. https://doi.org/10.1090/S0002-9939-1994-1182706-X ISSN: 0002-9939, 1088-6826 Issue: 2 Journal Abbreviation: Proc. Amer. Math. Soc.

MICHAEL Reed and BARRY Simon. 1972. *Methods of Modern Mathematical Physics*. Vol. I Functional Analysis. Academic Press. https://doi.org/10.1016/B978-0-12-585001-8.50001-5

J. P. O. Silberstein. 1962. Symmetrisable operators. *Journal of the Australian Mathematical Society* 2, 4 (Dec. 1962), 381–402. https://doi.org/10.1017/S14467887000027427

François Sillion and Claude Puech. 1994. *Radiosity and Global Illumination*. Morgan Kaufmann.

François X. Sillion, James R. Arvo, Stephen H. Westin, and Donald P. Greenberg. 1991. A global illumination solution for general reflectance distributions. In *Proceedings of the 18th annual conference on Computer graphics and interactive techniques (SIGGRAPH '91)*. Association for Computing Machinery, New York, NY, USA, 187–196. https://doi.org/10.1145/122718.122739

Peter-Pike Sloan, Jan Kautz, and John Snyder. 2002. Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments. *ACM Transactions on Graphics* 21, 3 (July 2002), 527–536. https://doi.org/10.1145/566654.566612

Peter-Pike Sloan, Ben Luna, and John Snyder. 2005. Local, deformable precomputed radiance transfer. In *ACM SIGGRAPH 2005 Papers (SIGGRAPH '05)*. Association for Computing Machinery, New York, NY, USA, 1216–1224. https://doi.org/10.1145/1186822.1073335

Cyril Soler, Mahdi M. Bagher, and Derek Nowrouzezahrai. 2015. Efficient and Accurate Spherical Kernel Integrals Using Isotropic Decomposition. *ACM Trans. Graph*. 34, 5 (Nov. 2015), 161:1–161:14. https://doi.org/10.1145/2797136

Cyril Soler, Ronak Molazem, and Kartic Subr. 2022. A Theoretical Analysis of Compactness of the Light Transport Operator. In *ACM SIGGRAPH 2022 Conference Proceedings (SIGGRAPH '22)*. Association for Computing Machinery, New York, NY, USA, 1–9. https://doi.org/10.1145/3528233.3530725

Jiguang Sun and Aihui Zhou. 2016. *Finite Element Methods for Eigenvalue Problems*. CRC Press. https://doi.org/10.1201/9781315372419 Journal Abbreviation: Finite Element Methods for Eigenvalue Problems Publication Title: Finite Element Methods for Eigenvalue Problems.

Eric Veach. 1997. *Robust Monte-Carlo Methods for Light Transport Simulation*. Ph. D. Dissertation. Stanford University.

Jiaping Wang, Yue Dong, Xin Tong, Zhouchen Lin, and Baining Guo. 2009. Kernel Nystrom method for light transport. In *ACM SIGGRAPH 2009 papers (SIGGRAPH '09)*. Association for Computing Machinery, New York, NY, USA, 1–10. https://doi.org/10.1145/1576246.1531335

Rui Wang, Jiajun Zhu, and Greg Humphreys. 2007. Precomputed Radiance Transfer for Real-time Indirect Lighting Using a Spectral Mesh Basis. In *Proceedings of the 18th Eurographics Conference on Rendering Techniques (EGSR'07)*. Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 13–21. https://doi.org/10.2312/EGWR/EGSR07/013-021

J. H. Wilkinson. 1959. The evaluation of the zeros of ill-conditioned polynomials. Part I. *Numer. Math.* 1, 1 (Dec. 1959), 150–166. https://doi.org/10.1007/BF01386381

Stephen M. Zemyan. 2012. *The Classical Theory of Integral Equations: A Concise Treatment* (1 ed.). Birkhäuser Boston.