Spectral Laplace Transform of Signals on Arbitrary Domains

The Laplace transform is a central mathematical tool for analysing 1D/2D signals and for the solution to PDEs; however, its definition and computation on arbitrary data is still an open research problem. We introduce the Laplace transform on arbitrary domains and focus on the spectral Laplace transform, which is defined by applying a 1D filter to the Laplacian spectrum of the input domain. The spectral Laplace transform satisfies standard properties of the 1D and 2D transforms, such as dilation, translation, scaling, derivation, localisation, and relations with the Fourier transform. The spectral Laplace transform is enough general to be applied to signals defined on different discrete data, such as graphs, 3D surface meshes, and point sets. Working in the spectral domain and applying polynomial and rational polynomial approximations, we achieve a stable computation of the spectral Laplace Transform. As main applications, we discuss the computation of the spectral Laplace Transform of functions defined on arbitrary domains (e.g., 2D and 3D surfaces, graphs), the solution of the heat diffusion equation, and graph signal processing.


Introduction
Due to the increasing availability of data, which is supported by ongoing technological advances in the acquisition, storage, and processing, several transformations (e.g., the Fourier transform, the Laplace transform, and the Fuzzy transform) have been proposed to solve problems that spread from signal analysis to the solution of partial differential equations, from the analysis to the approximation of signals. From a general perspective, a transform is typically defined as a linear operator between functional spaces, and its discretisation reduces to a matrix-vector multiplication. Main examples include the definition of the Fourier and Laplace transforms as integral operators induced by a complex and a real exponential kernel, respectively. Transformations are crucial to address different signal processing and analysis B Giuseppe Patané patane@ge.imati.cnr.it 1 CNR-IMATI, Genova, Italy 2 RAISE Ecosystem, Genova, Italy tasks, such as multi-rate signal processing, multi-scale representations and frames, signal modelling and adaptive filtering.
The Laplace transform [79] (LT, for short) has been extensively applied to signal analysis and the solution of PDEs on 2D domains, such as boundary value problems in two variables (Sect. 2). The definition and computation of the LT on arbitrary domains are challenging for several reasons. On the one hand, the solution to PDEs on curved volumes and surfaces is still an open research problem, due to the complexity of partitioning the input domain and constructing proper functional spaces. On the other hand, the definition of the LT on arbitrary domains must generalise the properties of the LT from the 1D case to the nD case, also taking into account the geometry and topology of the input domain.
Instead of applying patch operators for signals defined on arbitrary data, which cannot be extended to graphs by combining the vertex-frequency analysis framework [85] with classical spectral analysis, we work in the spectral domain. Filtering the Laplacian spectrum, we introduce the spectral LT (Sect. 3) and show that the main properties of the 1D and 2D LTs, such as dilation, translation, scaling, derivation, and localisation, still apply to a signal defined on an arbitrary domain. Defining the LT of spectral wavelets and filters provides a "compact" representation of the input signal, which allows us to analyse any signal (not necessarily defined in 1D or 2D) and achieve a stable computation of the spectral LT through polynomial representations and geometric series.
The spectral LT is particularly useful for the solution of the heat equation at small scales, where it is generally difficult to achieve a high approximation accuracy and a numerically stable computation, due to the discretisation of the temporal variable (e.g., θ -method [1,70]) and heuristic criteria for the selection of the number of Laplacian eigenpairs with respect to time and target accuracy. A smaller scale generally requires a higher number (i) of eigenfunctions with larger frequencies for the spectral approximation and (ii) of iterations for the θ -method, which introduces oscillations in the approximate solution. In contrast, the LT removes the time variable and transforms the input PDE in a Helmholtz-type boundary value problem, which is split into the solution to a homogeneous equation and the computation of a particular solution. Furthermore, differentiation and integration are converted to function multiplication and division through the LT.
For the efficient and numerically stable computation of the spectral LT (Sect. 4), we show a low approximation error between the input filter and its polynomial or rational polynomial approximations, which corresponds to a low error between the corresponding spectral LTs. As the first approach, we approximate the input filter with a polynomial, showing that the spectral LT can be iteratively computed. The second approach approximates the 1D LT with a rational polynomial and maps this approximation on the corresponding spectral LTs. Both methods allow us to compute the spectral LT without evaluating the Laplacian spectrum, which is computationally expensive and numerically unstable in the case of large data. As the main advantage, the spectral LT is intrinsic by construction, is independent of the discretisation of the input domain, and involves only a small set of Laplacian eigenvalues and eigenfunctions. Indeed, the proposed approach supports the analysis of any type of signal and data representation, e.g., meshes, point clouds, and graphs.
In the experimental analysis (Sect. 5), we evaluate the spectral LT of functions defined on 3D domains and graphs, also discussing their main properties, computation, and applications to the solution of the heat diffusion equation on arbitrary domains (e.g., 3D surfaces), graph signal processing, and 3D shape correspondence.

Related Work
We briefly recall previous work on the 1D LT (Sect. 2.1) and its link with the Fourier transform (FT, for short), the Laplacian spectrum (Sect. 2.2) and the spectral filters that will be applied to the definition of the spectral LT. Finally, we review LT-based and meshless solvers of the heat equation (Sect. 2.3).

1D Laplace and Fourier Transforms
The Laplace operator of a 1D function f : R + → R is defined as the integral operator The LT is linear and bijective (Lerch's theorem) in the space of continuous functions. In our discussion, we focus on the case of a real variable s. Since differentiation and integration are converted to function multiplication and division through the LT, this transform is commonly applied to solve time-depending PDEs, such as the heat diffusion equation. Then, the solution is mapped back to the initial space through the inverse LT [25,79]. Recalling that the Fourier transform (FT, for short) of a square-integrable function f : 2 , the FT decomposes signals into frequencies and maps signals localised in the time domain to signals spread out across the frequency domain, and vice versa. The continuous FT is equivalent to the LT, where the integration domain is R instead of [0, +∞) and with imaginary argument s := iω, i.e., On the one hand, the FT represents a signal as a sum of modes of vibration. On the other hand, the LT converts time-depending functions into functions of complex angular values. Applying the variable change ω := −iω in Eq. (1), we get that (L f )(ω) = (F f )(−iω); indeed, computing the LT is equivalent to evaluating the FT.

Spectral Signal Processing
Laplacian Spectral Operators For the definition of spectral LT, we introduce the Laplace-Beltrami operator and its spectrum, which provides a generalisation of the Fourier basis to arbitrary domains [76,77]. Let M be a compact and connected domain and L 2 (M) the space of square-integrable functions on M equipped with the scalar product f , g 2 := M f gdμ and the corresponding norm · 2 . The Laplace-Beltrami operator Δ is self-adjoint, positive semi-definite in L 2 (M), and admits the Laplacian orthonormal eigensystem (λ n , φ n ) +∞ n=0 [73], Δφ n = λ n φ n , with λ 0 = 0 and λ n ≤ λ n+1 . Let H := L 2 (R) ∩ C 0 (R) be the filter space and ϕ : R → R be a positive, continuous, and square-integrable filter in H. Assuming that the filter admits the power-series representation ϕ(t) := +∞ n=0 α n t n , we define the Laplacian spectral operator ϕ(Δ) := +∞ n=0 α n Δ n i.e., the action of ϕ(Δ) on f is equal to the scalar product between the spectral kernel K ϕ (p, q) := +∞ n=0 ϕ(λ n )φ n (p)φ n (q) and the input signal. Signal analysis in the spectral frequency domain applies graph filters, by selectively modifying the Fourier coefficients [11,24,41]. Spectral signal processing [24,40,68] represents an input signal in terms of the eigenvectors of a graph operator (e.g., the graph Laplacian, a kernel matrix) to define the Fourier transform and the convolution with another signal [22,81,82]. Then, the FT decomposes signals into frequencies and maps signals localised in the time domain to signals spread out across the frequency domain, and vice versa. According to this simple approach, whose discrete counterpart is based mainly on numerical linear algebra, spectral graph processing has been applied successfully to characterise geometric and topological properties of graphs and to dimensionality reduction [6], through the projection of the input data on low-dimensional subspaces generated by a small set of Laplacian [12] or kernel [34] eigenfunctions. However, the FT of a signal on arbitrary data is not a continuous function as in the 1D case but a discrete infinite real sequence that converges to zero. Spectral signal processing on graphs [22,81,82] includes signal interpolation [37,56] and diffusion wavelets [12,40,50,80]. Further applications are graph embedding [11], the definition of uncertainty principles [68] and random walks [74], data representation [95] and classification [60]. The modelling and training of convolutional neural networks [75] and the design of fast localised convolutional filters on graphs [24] have been recently addressed in the spectral domain and generalised to 3D data through geometric deep learning [4,66]. Convolutional filters [11,24] extract local features that are shared across the input domain and reduce the number of free parameters for generic deep architectures. Localised filters or kernels [10,24] are applied to identify similar features of the input data, independently of their size. To reduce the computational complexity of the evaluation of these filters, previous work [40,88] has proposed specific filter representations in terms of the Chebyshev polynomial basis, which avoids the evaluation of the Laplacian spectrum. A generalisation of CNNs to arbitrary data (e.g., graphs) is not straightforward as translation operators are not trivially defined. Graph Neural Networks [13,16,24,93,94] apply either filters in the spectral domain or spatial-based graph convolution, which combines information over neighbouring. Spectral wavelets Subdivision wavelets on the sphere have been defined through lifting schemes [86] and have been extended to arbitrary surfaces through multi-resolution analysis [47]. Further generalisations of these approaches include the B-spline wavelets [7], the subdivision [59,92] and spherical Haar [48] wavelets. Applying the properties of the Laplacian spectrum as a generalisation of the Fourier analysis, the diffusion wavelets [19] define a sequence of nested subspaces for multi-resolution signal representation. The spectral wavelets [19,58,83] are achieved as localised and scaled versions of ϕ(tΔ) f (c.f., Eq. (2)), by selecting f := δ p and scaling the filtered eigenvalues by a factor t. According to the selection of the filter function, we define different classes of wavelets, such as the Haar and diffusion wavelets [19,71]. The Mexican hat wavelets [38] is defined as the negative first-order derivatives of the heat kernel with respect to time; in particular, these wavelets are localised in both space and frequencies. Generalising the multi-resolution analysis for 1D and 2D signals [55], wavelet analysis on arbitrary data (e.g., graphs, 3D domains) is constructed by defining polynomials of a differential operator (e.g., the Laplace-Beltrami operator, the heat diffusion operator) [18,57], or applying a filter function to the Laplacian spectrum [21], or encoding energy spectral densities in the design of tight frames or graphs [14], or approximating smooth functions with tight wavelet frames [28].

LT and Meshless Methods for PDEs
Meshless methods compute the solution that satisfies the input equation at a set of collocation points, which have a non-ordered spatial point distribution and no pre-defined connectivity or spatial relationships. The Kansa method [43,44] lacks a symmetric collocation matrix due to the boundary collocation of mixed boundary conditions. Combining the differential and boundary operators of the input PDE with their self-adjoint operators, the Hermite collocation method [29] introduces a symmetric collocation matrix. To improve the accuracy of the solution in regions close to the boundary of Ω, the modified Kansa method [31] applies an additional set of collocation points inside and outside Ω, and close to ∂Ω. In [65], a novel collocation method with Radial Basis Functions has been applied to the solution of inhomogeneous parabolic equations by rewriting the solution in terms of the exponential operator that is approximated through the Padè-Chebyshev approximation of the 1D Gaussian function. The resulting meshless solver uniformly converges to the ground-truth solution, as the degree of the rational polynomial increases, and is independent of the evaluation of the spectrum of the input operator.
A time-dependent PDE is reduced to an inhomogeneous modified Helmholtz equation through the LT or a time-stepping approach. More precisely, the LT (Lg)(s) = +∞ 0 g(t) exp(−st)dt removes the time variable and transforms the input PDE in a Helmholtz-type boundary value problem, which is split into the solution to a homogeneous equation and the computation of a particular solution. The particular solution is approximated as a linear combination of RBFs and the method of fundamental solutions is applied to solve the homogeneous PDE. The Stiehfest's method [87] inverts the LT and provides the solution to the input and time-dependent PDE. However, the numerical inversion of the LT is ill-posed, as small truncation errors can be increased by the inversion process. Indeed, time-stepping methods (e.g., θ -method, time stepping with local refinement in space and time [1])) avoid the numerical instability of the inverse LT. In particular, the method of lines [39] applies a spatial discretisation of the input domain Ω as a set of points, evaluates the unknown solution at these points, and transforms the input PDE into a system of ordinary differential equations. Specialisations of these techniques include meshless solvers for transient heat conduction analysis [33], for diffusion problems and image halftoning [96].
For the numerical solution of time-fractional diffusion equations, the Laplace-transformed boundary particle method [30] applies an LT-based technique to obtain the corresponding time-independent inhomogeneous equation, which is then solved through a meshless boundary particle method that does not require any inner nodes. For the solution of time fractional derivative equations, previous work typically combines the finite difference method for temporal discretization with other numerical methods for spatial discretization, such as the Fourier method [20], spectral method [53], finite element method [32,42,54,69], boundary element method [45], and radial basis function meshless collocation methods [9,23,36,49]. These methods can reduce the computational cost for large computational domain problems but encounter an exponential increase of computing costs with advancing time, as a matter of the underlying finite difference schemes for temporal discretization.
In [46], a novel space-time meshless method has been applied to the solution of the backward heat conduction problem through a numerical approximation based on the Trefftz basis functions of the heat equation. Instead of applying collocation methods based on unstruc- tured points, boundary points are collocated in the space-time coordinate system such that the initial and boundary conditions are treated as boundary conditions on the space-time domain boundary. Then, the backward heat conduction problem is converted into an inverse boundary value problem. In [90], the solution of a time-fractional diffusion equation is computed through a local meshless method based on the LTs. Instead of a large global collocation matrix, differentiation matrices are constructed by solving small systems over small local domains and combined with LT for a temporal variable.

Spectral LT
Filtering the Laplacian spectrum, we introduce the Laplacian spectral kernels and the corresponding LTs (Sect. 3.1). Then, we define the spectral LT (Sect. 3.2) and study the relations between the spectral LT and FT on arbitrary domains. Considering the eigensystem of the 1D second-order derivatives, the definition of the spectral LT boils down to the 1D case. Furthermore, the spectral LT satisfies the main properties of the 1D and 2D LTs, such as dilation, translation, scaling, derivation, and localisation. Given a function on a 3D domain or its LT transform F, we visualise the corresponding level sets γ α := {p : F(p) = α} and colour-map, which varies from blue (global minimum) to red (global maximum). Given a function on a graph or its LT transform, we visualise the level sets of its underlying approximation and colour map, which begins with red, passes through yellow, green, cyan, blue, and magenta, and returns to red.

LT of Laplacian Spectral Wavelets and Lernels
The LT of the spectral wavelet ψ t,p in Eq. (3) is In particular, the LT of where Δ † is the pseudo-inverse of the Laplace-Beltrami operator and ϕ is the LT of the input filter. Indeed, L f t solves the harmonic equation Noting , as s → 0 + (Fig. 1), the evaluation of its LT reduces to the solution of the regularised Laplace equation (Δ+sI)g = f . Similarly, the LT of the "modulated" diffusion wavelet H t := t K t is equal to LH t (p, q) = +∞ n=0 (s+λ n ) −2 φ n (p)φ n (q). In particular, we get the following limits: lim s→+∞ (LH t )(s, p, q) = 0 and lim s→+∞ (LH t )(s, p, q) = 0.

Spectral LT: Definition and Properties
According to previous work on spectral filtering (Sect. 2.2), we define the spectral operator where Φ ϕ is the spectral function induced by ϕ (Figs. 2, 3a). The continuity of the filter function allows us to evaluate (ϕ(λ n )) +∞ n=0 and the L 2 (R)-integrability of the filter ensures the well-posedness of the spectral operator. In fact, the spectral operator is linear and continuous, according to the upper bound The properties of the spectral operator depend only on the behaviour of the filter in the spectral domain and on the Laplacian spectrum; increasing or decreasing the filter decay to zero encodes global or local properties of the input data, respectively. Furthermore (Fig. 4), a generally low number of Laplacian eigenfunctions in the sum in Eq. (6) guarantees a high approximation accuracy of the spectral functions. Generalising the LT of the spectral wavelets in Eqs. (4), (5) and recalling the definitions of the spectral function (6), we define the spectral LT of Φ ϕ and ϕ(Δ) as (Figs. 2, 3b).

Examples of spectral LTs
• For the Heasvide function

Polynomial and Rational Approximations of the Spectral LT
Approximating the input filter ϕ with a polynomial, or a rational polynomial, ρ (r ) , we show that a low approximation error between ϕ and ρ (r ) corresponds to a low approximation error between the corresponding spectral LTs Φ ϕ , Φ ρ (r) (Sect. 4.1) and the corresponding operators ϕ(Δ) and ρ (r ) (Δ). In the first approach, we approximate the input filter with a polynomial and show that the spectral LT can be iteratively computed (Sect. 4.2). The second approach approximates the 1D LT ϕ with a rational polynomial and maps this approximation on the corresponding spectral LTs (Sect. 4.3). Both methods allow us to compute the spectral LT without evaluating the Laplacian spectrum, which is generally computationally expensive and numerically unstable in the case of large data.

Approximation of the Spectral LT
First of all, we recall that the 1D LT is continuous. To this end, we notice that Let us assume that ρ : R → R is a function that approximates the filter ϕ. Applying the upper bound (9) and the linearity of the LT, we get that i.e., the approximation of Φ ϕ reduces to the computation of a function ρ that well approximates the filter ϕ in the L 2 (R)-norm. Approximating ϕ with a polynomial or a rational polynomial ρ (r ) , the sequence ( ρ (r ) (Δ) f ) +∞ r =0 , ρ (r ) (Δ) f := +∞ n=0 ρ (r ) (λ n ) f , φ n 2 φ n , induced by the rational polynomial approximation ρ (r ) of ϕ, converges to ϕ(Δ) f ; in fact, In particular, the error between the input spectral LT ϕ(Δ) and its approximation ρ (r ) (Δ) is controlled by the approximation error ρ (r ) −ϕ 2 (Fig. 7). The representation of ρ must guarantee that the evaluation of Φ ρ is performed efficiently and is independent of the evaluation of the Laplacian spectrum, which is computationally infeasible for large data and numerically unstable. To this end, we propose a fast polynomial approximation of the input filters based on a discrete scalar product induced by Chebyshev nodes (Sect. 4.2). For the evaluation of ρ(Δ) f , we apply a rational approximation (Sect. 4.3), which is efficiently evaluated through recursive relations, in analogy to Chebyshev polynomials.

Polynomial Approximation of the Spectral LT
For the computation of the spectral LT, we introduce a polynomial approximation of the input filter with canonical or Chebyshev polynomials. This last choice defines a recursive and fast computation of a polynomial approximation of the input filter and the corresponding spectral LT (Fig. 2).  6 Color map and level sets of the spectral function Φ ϕ , its LT F Φ ϕ and inverse LT F −1 Φ ϕ Canonical polynomials for the computation of the spectral LT Recalling the power series representation ϕ(t) = +∞ n=0 α n t n , introduced in Sect. 3.1, and that the LT of the function ψ(t) := t n is ψ(s) = n!/s n+1 , we get that the corresponding LT is ϕ(t) = +∞ n=0 α n n! s n+1 (Fig. 7), i.e., Chebyshev polynomials for the computation of the spectral LT Since the approximation of the input filter in terms of the canonical basis is generally expensive as the polynomial degree increases, we represent the filter in terms of the Chebyshev polynomials and apply a recursive formulation of the computation. The Chebyshev polynomials of the first kind (T n ) +∞ n=0 are defined by the identity cos(nθ) = T n (cos θ), e.g., T 0 (t) = 1, T 1 (t) = t, T 2 (t) = 2t 2 − 1, . . . Equivalently, these polynomials are defined through the recursive relation The Chebyshev polynomials of the first kind are orthogonal with respect to the weighted L 2 (R) scalar product induced by (1 − t 2 ) −1/2 on the interval [−1, 1]. Indeed, a filter ϕ : [−1, 1] → R is represented in terms of the Chebyshev polynomials as ϕ(t) = +∞ n=0 a n T n (t), whose coefficients are a 0 = π −1 ϕ, 1 w , a n = (2π) −1 ϕ, T n w , According to the relations the coefficients are evaluated in linear time according to a 0 = 1 Once the representation ϕ(t) = +∞ n=0 a n T n (t) has been computed, we apply the relation ϕ(Δ) = +∞ n=0 a n T n (Δ), where T n (Δ) is evaluated as done for the canonical basis in Eq. (10).

Rational Approximation of the Spectral LT
For the approximation and evaluation of the rational polynomial ρ(·) that approximates ϕ, we consider (i) the canonical polynomial basis and the Chebyshev polynomial basis and (ii) the Chebyshev rational polynomials. Among these options, the recursive representation of the Chebyshev rational polynomials of the first kind allows us to reduce the computation cost for the evaluation of the approximation and resembles the polynomial case (Fig. 3).
Rational approximation with canonical polynomials Instead of approximating the input filter ϕ before applying the LT, we approximate the LT ϕ with a rational polynomial of degree where P n , Q m are polynomials of degree n and m respectively. Through the inverse 1D LT, we get the approximation of both the spectral LT and the rational approximation of the input filter. The Padè approximant is uniquely defined if the constant term at the denominator has been set equal to 1 (c.f., Eq. 12); otherwise, the approximant is defined up to a multiplication by a constant. For the computation of the coefficients of the polynomials P n , Q m , we can apply Wylm's epsilon algorithm or the extended Euclidean algorithm [8], and the sequence transformation [15]. Through the canonical basis and the representation (12), we evaluate the right-side and left-side terms of If Q has distinct roots (α k ) n k=0 , Q(s) := k n=1 (s − α k ), then Rational approximation with Chebyshev polynomials Instead of expanding the numerator and denominator in terms of the monomial basis, we use the Chebyshev polynomials T k , In the case of the Chebyshev polynomials (14), we evaluate the right-side and left-side terms of through the recursive relations (11).

Rational approximation with rational Chebyshev polynomials
The Chebyshev rational polynomials are defined as on the interval [0, +∞). These rational polynomials are orthogonal with respect to the weighted scalar product induced by w(s) := s −1/2 (1 + s) −1 , according to the relation For an arbitrary function ϕ ∈ L 2 (R), the orthogonality of the Chebyshev rational polynomials allows us to apply the relation ϕ(s) = +∞ n=0 F n R n (s), F n := ϕ, R n w . Through the recursive relation (16), the spectral LT of ϕ(Δ) f is efficiently evaluated as where g n−1 is the solution to (Δ + id)g n−1 = (Δ − id) f n−1 and g n−1 is uniquely defined, as the operator (Δ + id) is positive-definite. Computing the polynomial and the rational polynomial approximations, a generally low number of rational polynomials (Fig. 8) allows us to achieve a high approximation accuracy (e.g., 10 −15 , n ≤ 15) and high numerical accuracy and robustness (y-axis), measured as the conditioning number of the coefficient Fig. 11 Robustness of the spectral LT to increasing noise. The colour map and level sets of the input function Φ ϕ and of its LT ϕ(Δ) are computed on the noise surface but visualised on the smooth surface to improve the analysis of the results. To generate a noisy surface. we perturb the coordinates of the mesh vertices with a Gaussian noise, whose magnitude is set equal to the 10% of the bounding box of the input surface. The maximum variation between the smooth (a) and noisy spectral LTs (b-d) remains lower than 5% matrix. In all the paper examples, the spectral LT is computed with a rational polynomial of degree r := 7 (if not differently stated).

Discrete Spectral LT and Experimental Results
Firstly, we introduce the discretisation of the Laplace-Beltrami operator and spectral operators; then, we discuss the discrete LT and its main properties (Sect. 5.1). For our tests, we consider the solution of the diffusion equation with the filtered Laplace operator (Sect. 5.2) and the Laplace-Beltrami operator (Sect. 5.3) through the LT.

Spectral Laplacian Operator: Discretisation and Properties
Graph Laplacian For the discretisation of the LT, let us consider a domain M and a signal f : M → R represented as the vector f := ( f (p i )) n i=1 of f -values at the nodes of the input graph or the vertices of a triangle mesh M. For unweighted graphs, the domain of the continuous FT and the continuous convolution is taken to be the graph connectivity (i.e., the binary adjacency matrix) and for a weighted graph, we consider the continuous FT induced by the weighted Laplacian matrix. On meshes, we select the linear FEM Laplacian matrix [78] and its spectrum for the evaluation of the continuous FT. The graph Laplacian [17,89] is defined asL := D −1 L, where L := W − D, and the normalised graph Laplacian is L := D −1/2 LD −1/2 . Here, W is the weight matrix whose entry (i, j) is the weight associated with the corresponding edge, and D is the diagonal matrix whose entries are the sum of the rows of L. The Laplace-Beltrami operator on a triangle mesh with n vertices is discretised as the n × n matrix L := D −1 L, where the mass matrix D encodes the variation of the areas of the Voronoi-regions [26] or triangles [78], and the stiffness matrix L encodes the cotangent of the triangle angles [67]. In all these cases, the spectral decomposition is LX = DXΛ, X DX = I, where X := [x 1 , . . . , x n ] is the eigenvectors' matrix and Λ is the diagonal matrix of the eigenvalues (λ i ) n i=1 [35,52,84]. Discretization According to Eq. (8), the discretisation of the spectral LT is ϕ( The LT of ϕ is computed analytically and then approximated with a rational polynomial to evaluate ϕ( L) without computing the Laplacian spectrum, which is computationally infeasible for large data sets. For the discretisation and computation of the spectral LT, we replace the Laplace-Beltrami operator with the Laplacian matrix; in this way, the equations (13), (15), (17) reduce to sparse equations that are efficiently solved with iterative sparse solvers. The computational cost for the continuous FT and convolution with an optimal number k of Laplacian eigenpairs is O(kn log n) − O(kn 2 ), according to the sparsity of the Laplacian matrix. In Fig. 9, we report the computation time for the continuous FT of a signal or surface with a different sampling density and with a different optimal number of Laplacian eigenpairs (i.e., a different target approximation accuracy). The evaluation of the LT is generally robust to domain sampling and connectivity (Fig. 10), and geometric noise (Fig. 11). Selecting different 1D filters and the corresponding LTs, these functions are easily "extended" to 3D domains through the Laplacian spectrum (Fig. 2).
The computational cost for the approximation of the spectral LT (c.f., Eq. 8) depends on the sparsity degree of the Laplacian matrix [35] and takes from O(kn log n) to O(kn 2 ) time, where k is the number of selected eigenpairs. Choosing a polynomial or rational approximation of the input filter of degree (r , l), the evaluation of the spectral LT is reduced to r linear systems whose coefficient matrix is L. Through iterative solvers, the computational cost is O(r τ (n)), where τ (n) is the cost for the solution of a sparse linear system, which varies from O(n) to O(n 2 ), according to the sparsity of the coefficient matrix, and it is O(n log n) in the average case. Indeed, the polynomial and rational polynomial approximations have the same order of computational complexity. Fig. 13 Concerning the tests in Fig. 12 (1st, 2nd rows), approximation error ∞ (y-axis) for the solution to the heat equation on 5 samplings of the input domain (Fig. 12, first row) with a different separation distance (xaxis). As a solver, we apply the linear FEM [62] (blue line) and meshless [90] (green line) P.-C. approximation (best approximation accuracy), and the LT-based solver (red line). The separation distance of the point set P := {p i } n i=1 is defined as the radius S(P) := 1 2 min 1≤i≤ j≤n { p i − p j 2 } of the smallest ball anywhere without any point in its interior, but with at least two points of P on its boundary (Color figure online)
Applying the LT to both sides of the new diffusion equation in Eq. (19) and noting that v(p, 0) = 0, we get the homogeneous filtered Laplace equation where w is the LT of w. If u 0 does not belong to the kernel of ϕ(Δ), then we consider the non-homogeneous filtered Laplace equation Representing the solution u(·, t) = +∞ n=0 α n (t)φ n in terms of the Laplacian eigenfunctions and imposing that it satisfies Eq. (20), we get the relation +∞ n=0 α n (s) ϕ(λ n ) − s κ(p) φ n = u 0 (p) κ(p) . Rewriting the initial condition in spectral form as u 0 (p) = +∞ n=0 α n (0)φ n , the pre-

Spectral LT Solver and Heat Equation
In the following experimental tests, we show that the LT is particularly useful for the solution of the heat equation at small scales, where it is generally difficult to achieve a high approximation accuracy and a numerically stable computation, due to the discretisation of the temporal variable (e.g., θ -method [1,70]) and heuristic criteria for the selection of the number of Laplacian eigenpairs with respect to time and target accuracy. To this end, we consider the homogeneous 2D heat equation ∂ t u − Δu = 0, Ω := [0, 1] × [0, 1], with (i) initial condition δ p , at t := 0 (Figs., 12,13) or (ii) the inhomogeneous heat equation with analytic solutions (Fig. 14). The approximation error is evaluated in terms of the difference between the computed and ground-truth solutions. The solution is computed with a linear FEM solver, the θ -method (e.g., θ := 1/2), the truncates spectral approximation, and the LT-based solver. Comparing (Figs. 12, 13) the linear FEM [62] (blue line) and meshless [90] (green line) P.-C. approximation (best approximation accuracy), and the LT-based solver (red line) for the solution to the homogeneous heat equation on 5 samplings of the input domain with a different separation distance (x-axis), the LT-based solver provides the highest approximation accuracy of the ground-truth solution. Furthermore, the spectral truncated approximation of the solution to the heat equation (Fig. 12, 100 eigenfunctions) with initial condition u(0, ·) = δ p is generally affected by undulations of the solution as we move far from the source point, especially at small scales (i.e., t → 0 + ). The shape and distribution of the level sets confirm the high accuracy of the spectral LT-based solver at small and large scales. According to these tests, a smaller scale generally requires a higher number (i) of eigenfunctions with larger frequencies for the spectral approximation and (ii) of iterations for the θ -method, which introduce oscillations in the approximate solution. In contrast, the LT-based solver provides the highest approximation accuracy of the ground-truth solution. The LT-based solver generally provides an approximation accuracy higher than the truncated spectral approximation, and the θ -method (θ := 1/2, 1) at all scales (i.e., t := 0.1, 0.01, 0.001). At small scales, the FEM solution to the heat equation generally shows small undulations (e.g., due to the Gibbs phenomenon), which are removed by the LT-based solver. A local discontinuity of the   (Fig., 14), our approach achieves an error lower than 10 −4 . Comparing the ground-truth solution to the homogeneous heat equation with its P.-C. approximation (r := 7) induced by the Gaussian kernel, whose radius is selected according to [27] and [72], and the distribution of the approximation error (Fig. 15), the spectral solver achieves an error lower than 1%, which is localised in a neighbourhood of the boundary of Ω, as a matter of the partial support of the RBFs whose centres are located close to ∂Ω.
On a 3D domain represented as a surface (Figs. 16, 17), we compute the solution to the homogeneous heat equation with initial condition δ p , at t := 0, whose solution at small scales is numerically unstable due to the discontinuity of the initial condition. These instabilities result in the Gibbs phenomenon for the θ -method, FEM solvers, and truncated spectral approximation. In contrast, the spectral LT solver is oblivious to the Gibbs phenomenon at all scales. The spectral LT removes the temporal scale, which is the main source of numerical instabilities of previous work, as t → 0. Furthermore, we require only to discretise the LB operator on the input domain Ω, without extracting a tetrahedral mesh of Ω as required by FEM solvers. An analogous approach applies to the solution of the heat equation on a graph (Fig. 18), where it is not possible to apply FEM solvers; in this case, the Laplace-Beltrami operator is discretised by the graph Laplacian [63,66].

Conclusions and Future Work
This paper has addressed the definition and properties of the spectral LT on arbitrary data (e.g., curved volumes, surfaces), as an extension of its 1D and 2D counterparts. The definition and computation of the LT on arbitrary domains are challenging for several reasons. On the one hand, the solution to PDEs on curved volumes and surfaces is still an open research problem, due to the complexity of partitioning the input domain and constructing proper functional spaces. On the other hand, the definition of the LT on arbitrary domains generalises the properties of the LT from the 1D case to the nD case, also taking into account the geometry and topology of the input domain. The spectral LT satisfies standard properties of the 1D/2D transform, such as dilation, translation, scaling, derivation, localisation, and Parseval's equality, and is l enough to be applied to different representations of arbitrary domains, such as surfaces (Figs. 16,17) or graphs (Fig. 18). To this end, we have defined the spectral LT in terms of the spectrum of the Laplace-Beltrami operator, which is intrinsic t = 10 −3 t = 10 −2 t = 10 −1 Fig. 18 Multi-scale behaviour of the spectral LT of a delayed impulse ϕ(x) := δ(x − t) at different time shifts t of unit impulse and centred at a seed node of the input graph. Here, the 1D LT is ϕ(s) = exp(−st). See also Fig. 17. Here, the colour map begins with red, passes through yellow, green, cyan, blue, and magenta, and returns to red (Color figure online)

Fig. 19
Correspondences on non-isometric shapes of the TOSCA data set [5] evaluated with the functional map framework [61], which uses diffusion kernels computed with the spectral LT-based solver to the input data and signals and allows us to encode their geometric and topological properties at different levels of complexity in the spectral LT. Through polynomial and rational polynomial approximations, the spectral LT is efficiently computed for signals defined on (i) structured (e.g., regular/irregular) data, (ii) time-depending data, where the model is not constructed once from a set of "static" data but is updated and adapted to the input data, and (iii) sparse or dense data. In the experimental analysis, we have discussed the application of the generalisation of the LT to the solution of the heat diffusion equation on arbitrary domains (e.g., 3D domains) and graph signal processing. The LT-based solution to the heat equation is also useful for the computation of diffusion functions that are applied to the evaluation of correspondences [61] between arbitrary 3D domains discretised by a triangle mesh with different geometry and connectivity. For shape correspondence, we select a low number k of ground-truth landmarks and consider the corresponding k spectral functions or sk diffusion functions at s scales computed with the spectral LT (Fig. 19). The ground-truth landmarks initialise the functional map and are used to compute the functions that generate the subspace on which the shape descriptors will be projected to define the functional map. For the rigid and articulated shapes, the computed correspondences correctly map local and global features on the source and target shape. In future work, we plan to focus on the link between the Laplace and Fourier transforms and extend the definition of convolution operators on arbitrary domains.