用于谱优化的神经网络

# Neural networks for spectral optimization
Source: [https://arxiv.org/html/2609.36047](https://arxiv.org/html/2609.36047)
###### Abstract\.

Given a functional dependent on the spectrum of a differential operator, we address the problem of finding a domain which optimizes this functional\. PDE solvers might be used to tackle this optimization\. It is however computationally expensive\. We propose two neural network models which learn the spectrum directly from the geometry of the domain and can be used to optimize the domain from one or more eigenvalues\. We investigate two representations\. The first encodes the domain through Fourier coefficients and a light MLP, which is efficient on star\-shaped geometries, achieving a precision of 0\.2%\. Through a rescaling of the coefficients the designed models satisfy the scaling law of the eigenvalues\. Additionally, averaging the outputs of the trained surrogates over rotations and reflections induces invariance for these transformations\. The second is a model that takes the landscape function, the indicator function and the gradient of the landscape function\. A Gram\-Schmidt process produces orthogonal eigenfunctions as output of the model along with the associated eigenvalues\. The landscape model reaches 1% mean relative error on the first ten eigenvalues, compared with 4% for an FNO model\. Replacing the landscape by an SDF worsened both prediction and optimization errors\. The trained model also generalizes from synthetic shapes to domains given as classical image dataset\. The resulting surrogates of both approaches recover classical spectral optima such as the disk for the first eigenvalue or the conjectured minima of higher eigenvalues\. This confirms that our models produce accurate differentiable estimates of eigenvalues, which can be used in shape optimization problems involving spectral quantities\.

## 1\.Introduction

Spectral geometry, which explores the interplay between the geometric properties of a domain and the spectrum of the differential operators defined on it, is a fundamental research field in mathematics\. Its questions, one of the most famous being "Why are all drums round?", find applications in numerous domains ranging from physics and biology to computer science\. A key problem in this field is shape optimization, which corresponds to the search for geometric shapes that optimize one or more eigenvalues of the Laplace operator\.

Understanding the spectrum of the Laplacian leverages many applications in geometry processing \(denoising, shape matching, shape retrieval\)\. Computing this spectrum with finite elements is accurate, but each new shape requires a new computationally expensive solve\. Neural surrogates aim to replace that by a single and fast forward pass\. However their accuracy depends on how the domain is represented\.

One of the major state\-of\-the\-art representation is theSigned Distance Function,S​D​FΩSDF\_\{\\Omega\}, withΩ\\Omegathe domain\. Briefly, the Signed Distance Function is a functionS​D​FΩ​\(𝐱\)SDF\_\{\\Omega\}\(\\mathbf\{x\}\)that gives the shortest distance from a point𝐱\\mathbf\{x\}to boundary∂Ω\\partial\\Omegaof the shape, with the sign indicating whether the point is inside or outside\. It is used in recent works in geometry processing and computer vision\[[6](https://arxiv.org/html/2609.36047#bib.bib10)\]\. A major issue of this representation is its topological sensitivity\. To illustrate this, let us consider a domain with an internal hole of radiusϵ\\epsilon\. Asϵ→0\\epsilon\\to 0\(capacity vanishing perturbation of radiusϵ\\epsilon\), the Dirichlet spectrum converges continuously to that of the domain without hole\. However theS​D​FSDFwill keep its zero value near the hole\.

Current state\-of\-the\-art approaches for predicting spectral properties also often rely onFourier Neural Operators\(FNO\)\[[11](https://arxiv.org/html/2609.36047#bib.bib4)\]orDeepONets\[[12](https://arxiv.org/html/2609.36047#bib.bib8)\]\.

DeepONets: DeepONets are based on the universal approximation theorem for operators\. The architecture consists of two sub\-networks: a Branch net, which encodes the input function \(e\.g\.,the indicator function of the domainΩ\\Omega\), and a trunk net, which encodes the coordinatesxxwhere the output function is evaluated\. The operator is represented as a dot product between the outputs of these two networks\.

Fourier Neural Operators \(FNO\)on the other hand perform spectral convolutions in the frequency domain\. The kernel integral operator is computed by applying a Fast Fourier Transform \(FFT\) to the input, multiplying by a learned weight matrix in the latent frequency space, and then applying an Inverse FFT\. Since this approach relies on global frequency transformations, FNOs often struggle to capture sharp, localized geometric features or complex boundary conditions \(Dirichlet/Neumann\)\. FNOs often suffer from artifacts near the edges and struggle to accurately resolve high\-frequency geometric features that are crucial for the upper part of the Laplacian spectrum\.

Finally, current approaches useVariational Physics\-Informed Neural Networks \(V\-PINNs\)for Eigenproblems\. While traditional PINNs solve for the eigenfunctionu⁡\(x\)u\(x\)given a shape\[[7](https://arxiv.org/html/2609.36047#bib.bib6),[13](https://arxiv.org/html/2609.36047#bib.bib5)\], recent Variational PINNs\[[8](https://arxiv.org/html/2609.36047#bib.bib7)\]use the Rayleigh Quotient:

\(1\.1\)ℒR​a​y​l​e​i​g​h=∫Ω\|∇u\|2​𝑑x∫Ω\|u\|2​𝑑x\\mathcal\{L\}\_\{Rayleigh\}=\\frac\{\\int\_\{\\Omega\}\|\\nabla u\|^\{2\}dx\}\{\\int\_\{\\Omega\}\|u\|^\{2\}dx\}as a loss function\. By minimizing this functional, the model can learn the eigenvaluesλk\\lambda\_\{k\}as the stationary values of the optimization process\. However, V\-PINNs are limited to finding the first few eigenvalues and require a new optimization for every new shape\.

### 1\.1\.Our contributions

In this paper, we propose differentiable surrogates for the first Dirichlet eigenvalues and eigenfunctions, with the explicit goal of shape optimization\. We compare two representations\. The first uses Fourier coefficients to represent the domain and has a light MLP structure that exploits translation and rotation invariance of specific geometries, its purpose being to quickly generate prototypes of optimal shapes for general functions\. The second approach is a topology\-aware encoder\-decoder that takes the landscape functionwΩw\_\{\\Omega\}\[[4](https://arxiv.org/html/2609.36047#bib.bib2)\], the indicator𝟏Ω\\mathbf\{1\}\_\{\\Omega\}and∇wΩ\\nabla w\_\{\\Omega\}, and predicts\(λk,ϕk\)\(\\lambda\_\{k\},\\phi\_\{k\}\)with a differentiable Gram\-Schmidt step\. We show that the landscape representation outperforms an SDF encoding and other state\-of\-the\-art models, recovers classical optimization outcomes such as the disk forλ1\\lambda\_\{1\}, and applies to domains given as images\.

## 2\.Background and related work

On every bounded open setΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}one considers the Laplace operator with Dirichlet boundary conditions\. The spectrum of the operators consists on a sequence of eigenvalues

0<λ1​\(Ω\)≤λ2​\(Ω\)≤⋯→\+∞,0<\\lambda\_\{1\}\(\\Omega\)\\leq\\lambda\_\{2\}\(\\Omega\)\\leq\\cdots\\to\+\\infty,multiplicity being counted\. These eigenvalues are the critical points of the Rayleigh quotient, namely

λk​\(Ω\)=minE∈𝒮k⁡maxu∈E∖\{0\}​∫Ω\|∇u\|2​𝑑x∫Ω\|u\|2​𝑑x,\\lambda\_\{k\}\(\\Omega\)=\\min\_\{E\\in\\mathcal\{S\}\_\{k\}\}\\max\_\{u\\in E\\setminus\\\{0\\\}\}\\frac\{\\int\_\{\\Omega\}\|\\nabla u\|^\{2\}dx\}\{\\int\_\{\\Omega\}\|u\|^\{2\}dx\},where𝒮k\\mathcal\{S\}\_\{k\}denotes the family of subspaces of dimensionkkin the Sobolev spaceH01​\(Ω\)H^\{1\}\_\{0\}\(\\Omega\)\. If the boundary ofΩ\\Omegais smooth, there exists a functionuk∈H01​\(Ω\)∖\{0\}u\_\{k\}\\in H^\{1\}\_\{0\}\(\\Omega\)\\setminus\\\{0\\\}such that

−Δ​uk=λk​\(Ω\)​uk​in​Ω,uk=0​on​∂Ω\.\-\\Delta u\_\{k\}=\\lambda\_\{k\}\(\\Omega\)u\_\{k\}\\mbox\{ in \}\\Omega,\\;u\_\{k\}=0\\mbox\{ on \}\\partial\\Omega\.If the boundary ofΩ\\Omegais nonsmooth the equation above is satisfied in a weak, variational, sense\.

In spectral geometry, mathematicians are interested in shape optimization problems which formally may be written as

minΩ∈𝒰a​d⁡F⁡\(λ1​\(Ω\),…,λk​\(Ω\)\),\\min\_\{\\Omega\\in\\mathcal\{U\}\_\{ad\}\}F\(\\lambda\_\{1\}\(\\Omega\),\\dots,\\lambda\_\{k\}\(\\Omega\)\),where𝒰a​d\\mathcal\{U\}\_\{ad\}is the family of admissible open sets andF:ℝk→ℝF:\\mathbb\{R\}^\{k\}\\to\\mathbb\{R\}\. A typical class𝒰a​d\\mathcal\{U\}\_\{ad\}of interest consists of all open subsets ofℝd\\mathbb\{R\}^\{d\}of prescribed measure\. The simplest such problem, forF=λ1​\(Ω\)F=\\lambda\_\{1\}\(\\Omega\), has a long history starting with the conjecture by lord Rayleigh from 1877 claiming that among all membranes with the same area the circular membrane \(disk\) has the lowest fundamental frequency \(this encodes the circular shape of all drums\)\. The conjecture was proved only in 1923 by Faber and Krahn, but opened the way to a series of deep conjectures relating geometry to the spectrum\. Nowadays, these problems are studied by the spectral geometry community\. While most of the questions are of pure mathematical nature, the theory of shape optimization initially developed by the engineering community provides effective tools to support or disprove these conjectures\. However, running a shape optimization algorithm requires fine knowledges from PDEs, calculus of variations, optimization and scientific computing, making such procedures inaccessible to the most of mathematicians working in spectral geometry\.

Several mathematical methods are available to approximate the eigenvalues of a given setΩ\\Omega: finite elements, finite differences, fundamental solutions, etc\. and several shape optimization algorithms can be used \(shape derivative, topological derivative, relaxation, level sets, etc\)\. For each specificFF, one has to readapt the algorithms taking into account its specific features\. The purpose of this paper is to propose a unified, simple method, based on neural networks which contains the full sequence \(computations of eigenvalues \- shape derivatives \- shape optimization\) and applies to a generalFF\.

A key question with a deep impact on the robustness of the computation, but also having its own mathematical interest, is the following: does a small geometric perturbation ofΩ\\Omegamay have a strong impact on the variation of the spectrum? What is the right tool to capture the fine properties of the geometric dependence of the spectrum? Learning the spectrum from geometry and making this procedure stable, relies precisely on this analysis\.

A first intuition is to take the characteristic function as an input for learning the spectrum\. This is mathematically inconsistent with the constraint in terms of volume, since characteristic functions which differ up to a set of Lebesgue measure00have the sameL1L^\{1\}\-norm\. While theL1L^\{1\}norms is stable for small variation of the characteristic function, the spectrum is very sensitive\. A second possible input is the signed distance function to the boundary\. This choice is again problematic, since a point in the middle of a disc would dramatically change the distance function, without any influence on the spectrum \(points have00capacity in dimension22\)\.

The choice we make is to use as main input the torsion \(or landscape\) function\. For everyΩ\\Omegathe landscape functionwΩw\_\{\\Omega\}solves \(in a weak sense if the boundary is not smooth,wΩ∈H01​\(Ω\)w\_\{\\Omega\}\\in H^\{1\}\_\{0\}\(\\Omega\)\)

−Δ​wΩ=1​in​Ω,wΩ=0​on​∂Ω\.\-\\Delta w\_\{\\Omega\}=1\\mbox\{ in \}\\Omega,w\_\{\\Omega\}=0\\mbox\{ on \}\\partial\\Omega\.It has been mathematically proved that the landscape function controls the spectrum, precisely ifΩ1⊂Ω2\\Omega\_\{1\}\\subset\\Omega\_\{2\}then

\|λk​\(Ω1\)−λk​\(Ω2\)\|≤C​∫wΩ2−wΩ1​𝑑x,\|\\lambda\_\{k\}\(\\Omega\_\{1\}\)\-\\lambda\_\{k\}\(\\Omega\_\{2\}\)\|\\leq C\\int w\_\{\\Omega\_\{2\}\}\-w\_\{\\Omega\_\{1\}\}dx,with a constantCCwhich is controlled\.

Equipped with this representation, several recent works in operator learning relate to our method\.\[[15](https://arxiv.org/html/2609.36047#bib.bib14)\]exploited the Rayleigh quotient inside a neural architecture to compute eigenvalues, but, like the V\-PINN approach discussed in the introduction, their network is re\-optimized for every new shape\.\[[10](https://arxiv.org/html/2609.36047#bib.bib12)\]proposed a hybrid finite\-element/neural framework \(FEENet\) that still requires a mesh of the domain, whereas the torsion\-based encoding only needs an indicator function\.\[[3](https://arxiv.org/html/2609.36047#bib.bib11)\]and\[[16](https://arxiv.org/html/2609.36047#bib.bib13)\]address related PDE\-on\-complex\-geometry and geometry\-processing problems, but do not target the Dirichlet spectrum or shape optimization\. On the shape\-optimization problem,\[[2](https://arxiv.org/html/2609.36047#bib.bib16)\]used a variational neural network to perform volume\-preserving optimization of the Dirichlet energy, which is methodologically close to our differentiable pipeline but targets a different \(non\-spectral\) functional\. Finally, the SDF representation discussed above is based on the representation presentated by\[[14](https://arxiv.org/html/2609.36047#bib.bib15)\]for shape reconstruction\. Overall, unlike these approaches, our two method predicts the firstK=10K=10eigenvalues in a single forward pass, for arbitrary domains given as a raster image, and is differentiable\.

## 3\.Fourier coefficients based model

### 3\.1\.Approach

When performing numerical shape optimization in the plane, one of the most natural class of admissible shapes is the class of bounded, star shaped, open sets\. These shapes can be described in polar coordinates by

Ω=\{\(ϕ,ρ\)∈\[0,2​π\)×ℝ\+\|ρ<ρΩ​\(ϕ\)\},\\Omega=\\\{\(\\phi,\\rho\)\\in\[0,2\\pi\)\\times\\mathbb\{R\}^\{\+\}~\|~\\rho<\\rho\_\{\\Omega\}\(\\phi\)\\\},where the functionρΩ\\rho\_\{\\Omega\}is a generic positive, continuous,2​π2\\pi\-periodic function onℝ\\mathbb\{R\}\. As such, it can be expressed as a Fourier series,

ρΩ​\(ϕ\)=a0\+∑i≥1\(ai​cos⁡\(i​ϕ\)\+bi​sin⁡\(i​ϕ\)\),\\rho\_\{\\Omega\}\(\\phi\)=a\_\{0\}\+\\sum\_\{i\\geq 1\}\(a\_\{i\}\\cos\(i\\phi\)\+b\_\{i\}\\sin\(i\\phi\)\),the setΩ\\Omegabeing uniquely defined by the sequences of coefficients\(ai\)i≥0\(a\_\{i\}\)\_\{i\\geq 0\}and\(bi\)i≥1\(b\_\{i\}\)\_\{i\\geq 1\}\. It is therefore natural to focus on sets defined by a fixed finite number of coefficients2​N\+12N\+1\. In this study we work in the class of shapes

𝒰15=\{Ω⊂ℝ2\|star shaped,ρΩ\(ϕ\)=a0\+∑i=115\(aicos\(iϕ\)\+bisin\(iϕ\)\)\}\.\\mathcal\{U\}\_\{15\}=\\Big\\\{\\Omega\\subset\\mathbb\{R\}^\{2\}~\\Big\|~\\text\{ star shaped\},\\ \\rho\_\{\\Omega\}\(\\phi\)=a\_\{0\}\+\\sum\_\{i=1\}^\{15\}\(a\_\{i\}\\cos\(i\\phi\)\+b\_\{i\}\\sin\(i\\phi\)\)\\Big\\\}\.A strong quality of this class is that even though it is of finite dimension each shape has an exact analytic formulation and is independent of any discretization used to numerically compute the eigenvalues\. Additionally, computing the area of a shapeΩ\\Omegain the class𝒰15\\mathcal\{U\}\_\{15\}is explicit, indeed

\(3\.1\)\|Ω\|=π​a02\+π2​∑i=115\(ai2\+bi2\)\.\|\\Omega\|=\\pi a\_\{0\}^\{2\}\+\\frac\{\\pi\}\{2\}\\sum\_\{i=1\}^\{15\}\(a\_\{i\}^\{2\}\+b\_\{i\}^\{2\}\)\.Yet, there is no trivial method to determine that a given truncated Fourier series stays positive on the whole interval\[0,2​π\)\[0,2\\pi\)\. To create the database needed to train our models, we add a decay assumption of the coefficients\(ai\)0≤i≤15\(a\_\{i\}\)\_\{0\\leq i\\leq 15\}and\(bi\)1≤i≤15\(b\_\{i\}\)\_\{1\\leq i\\leq 15\}to ensure the positivity of the border function, namely,

\(3\.2\)∀1≤i≤15,\(ai2\+bi2\)12≤0\.75​a0i2\.\\forall 1\\leq i\\leq 15,\\qquad\(a\_\{i\}^\{2\}\+b\_\{i\}^\{2\}\)^\{\\frac\{1\}\{2\}\}\\leq 0\.75\\,\\frac\{a\_\{0\}\}\{i^\{2\}\}\.
We generate a database of2000020000entries \(see Figure[1](https://arxiv.org/html/2609.36047#S3.F1)\) satisfying this decay assumption and computing the first1010Dirichlet eigenvalues of each entry using P2 elements and a finite element PDE solver written in FreeFEM\[[5](https://arxiv.org/html/2609.36047#bib.bib9)\]\.

Figure 1\.Sample from the Fourier database\.Our aim is to build a surrogate to the PDE solver that computes the list of eigenvaluesL⁡\(Ω\)=\{λ1​\(Ω\),…,λ10​\(Ω\)\}L\(\\Omega\)=\\\{\\lambda\_\{1\}\(\\Omega\),\\dots,\\lambda\_\{10\}\(\\Omega\)\\\}given the parameters\(ai\)0≤i≤15\(a\_\{i\}\)\_\{0\\leq i\\leq 15\}and\(bi\)1≤i≤15\(b\_\{i\}\)\_\{1\\leq i\\leq 15\}that we couple with the symbolic differentiation library SymPy to perform the optimization of a generic functionalF⁡\(λ1,…,λ10\)F\(\\lambda\_\{1\},\.\.\.,\\lambda\_\{10\}\)\. The surrogate needs not only to be accurate but also to be sufficiently stable to perform optimization from its partial derivatives\. Additionally, it needs to satisfy some qualitative properties of the eigenvalues: their scaling law and their invariance by rotations and reflections\. Note that we do not try to impose invariance under translations of the shape, indeed for a given shapeΩ\\Omegain𝒰15\\mathcal\{U\}\_\{15\}the translated shapeΩ′=Ω\+x\\Omega^\{\\prime\}=\\Omega\+xis in general not in𝒰15\\mathcal\{U\}\_\{15\}\.

Concerning the scaling law, for eachkkand eachΩ⊂ℝ2\\Omega\\subset\\mathbb\{R\}^\{2\}thekk\-th eigenvalue satisfies for anyt\>0t\>0the relationλk​\(t​Ω\)=t−2​λk​\(Ω\)\\lambda\_\{k\}\(t\\Omega\)=t^\{\-2\}\\lambda\_\{k\}\(\\Omega\)which implies in particular that the normalized eigenvaluejk​\(Ω\)=\|Ω\|​λk​\(Ω\)j\_\{k\}\(\\Omega\)=\|\\Omega\|\\lambda\_\{k\}\(\\Omega\)is invariant for dilations ofΩ\\Omega\. Noting that the area of a given shapeΩ\\Omegacan be explicitly computed from its coefficients by \([3\.1](https://arxiv.org/html/2609.36047#S3.E1)\), we decide to train models that do not learn directly the listL⁡\(Ω\)L\(\\Omega\)but rather the associated list of normalized eigenvaluesJ⁡\(Ω\)J\(\\Omega\)\. Because of this normalization, we also notice that it is not needed to keep track of the first coefficient to computeJJas it is uniquely determined by the reduced set of coefficientsq=\(a1/a0,…,a15/a0,b1/a0,…,b15/a0\)q=\(a\_\{1\}/a\_\{0\},\.\.\.,a\_\{15\}/a\_\{0\},b\_\{1\}/a\_\{0\},\.\.\.,b\_\{15\}/a\_\{0\}\)\. We then design the surrogate to computeJJgiven the data of a vectorq∈ℝ30q\\in\\mathbb\{R\}^\{30\}, and it satisfies the scaling law directly from its structure\.

To increase the stability of the derivatives of the surrogate we aggregate the results given by three independently trained MLPs\. Each MLP has the same structure :30→192→192→192→1030\\to 192\\to 192\\to 192\\to 10\. Both the width and the activation function have been tested for performance and SiLU has been retained \(see Section[3\.2](https://arxiv.org/html/2609.36047#S3.SS2)\)\. The MLP is trained to compute for a givenqqthe quantitylog⁡J\\log J\(thelog\\logis applied term by term on the listJJ\) and we denote bylog⁡J​\(q\)\(m\)\\log J\(q\)^\{\(m\)\}the output of a network\(m\)\(m\)\.

To account for the invariance by reflection and rotations we average the output of the network\(m\)\(m\)over1212fixed rotations and their associated reflections for any given shapeΩ\\Omega\. Doing so, we compute the function

log⁡J​\(q\)\(m\)¯=1\#​𝒯​∑τ∈𝒯log⁡J​\(τ⁡\(q\)\)\(m\),\\overline\{\\log J\(q\)^\{\(m\)\}\}=\\frac\{1\}\{\\\#\\mathcal\{T\}\}\\sum\_\{\\tau\\in\\mathcal\{T\}\}\\log J\(\\tau\(q\)\)^\{\(m\)\},
where𝒯\\mathcal\{T\}denotes the set of transformations corresponding to the reflections and rotations and\#​𝒯\\\#\\mathcal\{T\}its cardinal\. This averaged function,q↦log⁡J​\(q\)\(m\)¯q\\mapsto\\overline\{\\log J\(q\)^\{\(m\)\}\}, is by design invariant by reflections of the shape and by the twelve rotations\. This induces the needed reflection invariance and a close to rotational invariance to the final surrogate without the necessity of learning it in the weights of the networks\.

Finally, we aggregate the three networks computing the following average

J~​\(q\)=13​∑\(m\)exp⁡\(log⁡J​\(q\)\(m\)¯\)\.\\tilde\{J\}\(q\)=\\frac\{1\}\{3\}\\sum\_\{\(m\)\}\\exp\(\\overline\{\\log J\(q\)^\{\(m\)\}\}\)\.

### 3\.2\.Experiments

To test our surrogate, we use two different types of validation\. We first validate the model against test data unused during the training process for its accuracy\. We then validate the model through optimization experiments, we check that it produces the expected sets and values compared to the shape optimization literature\. As comparison we use the values obtained in\[[1](https://arxiv.org/html/2609.36047#bib.bib1)\]in a similar class of sets but with a more traditional shape gradient based optimization process initialized with a genetic algorithm\. For reference we give in Table[1](https://arxiv.org/html/2609.36047#S3.T1)the optimal values they obtained for\{λ5,…,λ10\}\\\{\\lambda\_\{5\},\.\.\.,\\lambda\_\{10\}\\\}\.

Table 1\.Optimal values for\{λ5,…,λ10\}\\\{\\lambda\_\{5\},\\dots,\\lambda\_\{10\}\\\}computed by Antunes and Freitasλ5=78\.20\\lambda\_\{5\}=78\.20λ6=88\.52\\lambda\_\{6\}=88\.52λ7=106\.14\\lambda\_\{7\}=106\.14λ8=118\.92\\lambda\_\{8\}=118\.92λ9=132\.68\\lambda\_\{9\}=132\.68λ10=142\.72\\lambda\_\{10\}=142\.72We first perform tests to determine the best structure for the MLP\. We test different activation functions and widths for the hidden layers\. We test for the activation function in \{ReLU, SiLU, Tanh\} and the width in\{90,192,300\}\\\{90,192,300\\\}for networks computing the list of the first1010rescaled Dirichlet eigenvalues\. For each couple \(activation, width\), we independently train three MLP networks of structure :30→witdh→witdh→witdh→1030\\to\\text\{witdh\}\\to\\text\{witdh\}\\to\\text\{witdh\}\\to 10for12001200epoch and the best epoch is kept\. We measure the total time needed to train the three networks and the mean relative error of the aggregated model \(see Table[2](https://arxiv.org/html/2609.36047#S3.T2)\)\. We observe that for each width, the most accurate model is the aggregate of networks activated by SiLU and the least accurate is the Tanh aggregate\. We also notice that although the models of width300300are more accurate than those of width192192the training time is much greater for a very marginal gain\. Our choice is then to use networks of width192192and hidden layers30→192→192→192→1030\\to 192\\to 192\\to 192\\to 10for further tests and compare the results obtained with ReLU and SiLU as activation functions\.

Table 2\.Mean relative error \(%\) by eigenvalue and training time \(s\) of the surrogate depending on the width of the hidden layers of the networks and their activation function\.Now performing tests to verify that it is possible to use our models to perform optimization we choose to compare their performance while optimizingλ5\\lambda\_\{5\}throughλ10\\lambda\_\{10\}among sets of fixed area\. During the optimization process, we follow the partial derivatives of the surrogate obtained through an automatic differentiation of the composing networks and we couple this gradient descent with a projection of the new shape on𝒰15\\mathcal\{U\}\_\{15\}\. We initialize the optimization on ten different shapes: the disk, the four best shapes of the training data for the given eigenvalue and five randomly chosen shapes among the rest of training shapes\. We then keep the best of the resulting optimized shapes\. We use the same starting shapes for the two models activated by ReLU and SiLU and compare their outputs by computing the eigenvalues of the two resulting optimized shapes with the FreeFEM solver \(Figure[2](https://arxiv.org/html/2609.36047#S3.F2)\)\. We observe that the results are similar for both activation function and that there is no clear better function\. Both models give shapes that are close to those of Antunes and Freitas in aspect and in eigenvalue and would be good candidates to initialize an optimization process involving a traditional PDE solver\.

λ5=79\.72\\lambda\_\{5\}=79\.72

λ6=89\.21\\lambda\_\{6\}=89\.21

λ7=109\.17\\lambda\_\{7\}=109\.17

λ8=120\.83\\lambda\_\{8\}=120\.83

λ9=137\.13\\lambda\_\{9\}=137\.13

λ10=144\.31\\lambda\_\{10\}=144\.31

λ5=79\.96\\lambda\_\{5\}=79\.96

λ6=89\.24\\lambda\_\{6\}=89\.24

λ7=108\.89\\lambda\_\{7\}=108\.89

λ8=120\.77\\lambda\_\{8\}=120\.77

λ9=138\.69\\lambda\_\{9\}=138\.69

λ10=144\.08\\lambda\_\{10\}=144\.08

Figure 2\.Optimized shapes for\{λ5,…,λ10\}\\\{\\lambda\_\{5\},\\dots,\\lambda\_\{10\}\\\}with activation function ReLU \(top\) and SiLU \(bottom\) and associated eigenvalues computed with the PDE solver\.We therefore retain the function SiLU as it gives a slightly better accuracy and since it is smooth, so is the resulting model\. Having a smooth model would allow for more involved optimization processes involving for instance second order derivatives\. With this retained model coupled with the library SymPy we are able to perform optimization of anyC1C^\{1\}function of the ten first eigenvalues among sets in𝒰15\\mathcal\{U\}\_\{15\}of area11\. We give in Figure[3](https://arxiv.org/html/2609.36047#S3.F3)some examples, note all the proposed functionals indeed admit optimizers as we work in the class of simply connected shapes\.

min⁡λ1\+λ2\\min\\lambda\_\{1\}\+\\lambda\_\{2\}

min∑i=110λi\\min\\sum\_\{i=1\}^\{10\}\\lambda\_\{i\}

max⁡λ6/λ1\\max\\lambda\_\{6\}/\\lambda\_\{1\}

min⁡λ2\+cos⁡\(λ5−λ1\)\\min\\lambda\_\{2\}\+\\cos\(\\lambda\_\{5\}\-\\lambda\_\{1\}\)

Figure 3\.Optimal shapes for different functionsFF\.

## 4\.Topological aware model with the landscape function

### 4\.1\.Approach

In this section, we consider as input representation the landscape functionwΩw\_\{\\Omega\}, introduced in Section[2](https://arxiv.org/html/2609.36047#S2), together with the indicator function

𝟏Ω​\(x\)=1​if​x∈Ω∪∂Ω,𝟏Ω​\(x\)=0​otherwise\.\\mathbf\{1\}\_\{\\Omega\}\(x\)=1\\text\{ if \}x\\in\\Omega\\cup\\partial\\Omega,\\qquad\\mathbf\{1\}\_\{\\Omega\}\(x\)=0\\text\{ otherwise\.\}In contrast to the SDF, the variations of the landscape functionwΩw\_\{\\Omega\}are consistent with variations of the Dirichlet spectrum, both are robust to capacity\-vanishing perturbations, as discussed in Section[2](https://arxiv.org/html/2609.36047#S2)\.

Since𝟏Ω\\mathbf\{1\}\_\{\\Omega\}is not differentiable, we replace it by a smoothedWΩ=σ⁡\(k​wΩ/wmax\)W\_\{\\Omega\}=\\sigma\(k\\,w\_\{\\Omega\}/w\_\{\\max\}\), whereσ\\sigmais the sigmoid function\. The network input is a 3\-channel image on a64×6464\\times 64grid\. The first channel is the smoothed indicatorWΩW\_\{\\Omega\}, the second one is the normalized landscapewΩ/wmax∈\[0,1\]w\_\{\\Omega\}/w\_\{\\max\}\\in\[0,1\], and the third one is the squared gradient\|∇wΩ/wmax\|2\|\\nabla w\_\{\\Omega\}/w\_\{\\max\}\|^\{2\}\.

Contrary to the Fourier coefficients of Section[3](https://arxiv.org/html/2609.36047#S3), our grid image has no rotation equivariance/invariance: the discretization and the resampling to64×6464\\times 64all yields non rotation invariance\. We thus rely on data augmentation during training, we simply apply rotations byπ2\\frac\{\\pi\}\{2\}together with a horizontal flip\.

We use a U\-Net encoder and decoder with skip connections, groupNorm and GELU activations\. The decoder outputsK=10K=10candidate functions on the grid\. The output is orthonormalized thanks to a differentiable Gram\-Schmidt process, this leads to eigenfunctions approximationsu^1,…,u^K\\hat\{u\}\_\{1\},\\dots,\\hat\{u\}\_\{K\}\. The ratiosλ^k/λ^1\\hat\{\\lambda\}\_\{k\}/\\hat\{\\lambda\}\_\{1\},k=2,…,Kk=2,\\dots,K, are computed with the Rayleigh quotient of eachu^k\\hat\{u\}\_\{k\}\. A scalar head applied to the pooled encoder bottleneck predictslog⁡\(λ^1⋅wmax\)\\log\(\\hat\{\\lambda\}\_\{1\}\\cdot w\_\{\\max\}\)\. We compare different architectures with the same representation and eigenvalues in the following subsection\.

The loss during training is the sum of three terms\. More precisely, aL1L^\{1\}loss on the ratiosλ^k/λ^1\\hat\{\\lambda\}\_\{k\}/\\hat\{\\lambda\}\_\{1\}, anL1L^\{1\}loss onlog⁡\(λ^1⋅wmax\)\\log\(\\hat\{\\lambda\}\_\{1\}\\cdot w\_\{\\max\}\), a Rayleigh\-quotient term\[[8](https://arxiv.org/html/2609.36047#bib.bib7)\]onu^k\\hat\{u\}\_\{k\}directly\.

Concerning the dataset, we generateN=100000N=100000domains including near\-circular, star\-shaped, non\-convex, multiply\-connected, several parts, polygons and equilateral triangles\. For some of the domains, we computeΦ=h−r2/4\\Phi=h\-r^\{2\}/4, wherehhis harmonic \(Δ​h=0\\Delta h=0\) andΔ⁡\(r2/4\)=1\\Delta\(r^\{2\}/4\)=1, this giveswΩ=Φ\|Ωw\_\{\\Omega\}=\\Phi\|\_\{\\Omega\}analytically, with no expensive solving\. For polygons,wΩw\_\{\\Omega\}is obtained by finite differences on a subset of aℤ2\\mathbb\{Z\}^\{2\}grid\. Note that this does not lead to large errors since the following subsection shows that polygons are the shapes associated to the smallest relative errors\. Across the dataset,λ1\\lambda\_\{1\}ranges over\[1\.88,751\.89\]\[1\.88,751\.89\]\. Furthermore, the input domains are approximated by a discretisation which induces some errors on the prediction of the spectrum\. On a grid, a finite difference used to compute the eigenvalues only needs the discrete boundary of the domain whereas a FEM method constructs a mesh that approximate more precisely the domain boundary\. We then trained our model on FreeFEM computed eigenvalues, leading to a mean relative error of1\.15%1\.15\\%\(median0\.53%0\.53\\%\) over allKKeigenvalues\. Figure[4](https://arxiv.org/html/2609.36047#S4.F4)illustrates some representative examples of our dataset\.

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_1.png)

holed

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_2.png)

several parts

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_3.png)

nonconvex

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_4.png)

polygon

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_5.png)

simple

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_6.png)

star

![Refer to caption](https://arxiv.org/html/2609.36047v1/torFunctdataset_7.png)

triangle

Figure 4\.One representative domain for each different shape\.
### 4\.2\.Experiments and ablation study

We evaluate our model and method after training\. Firstly, Table[3](https://arxiv.org/html/2609.36047#S4.T3)gives the mean relative error for each eigenvalue during the test of 15000 shapes not seen during training\.

Table 3\.Mean relative test error for each eigenvalue after validation\.The Table[4](https://arxiv.org/html/2609.36047#S4.T4)gives a comparison per shape type\.

Table 4\.Mean relative test error for each shape after validation\.Models comparisonWe compare our model, FNO and DeepONet after training on the same dataset ofN=100000N=100000domains\. We compare these architectures under two input representations, first the landscape\-based, and second an indicator function described previously\. All models are trained for300300epochs, and evaluated by the mean relative error over\{λ1,…,λ10\}\\\{\\lambda\_\{1\},\\dots,\\lambda\_\{10\}\\\}on a test set\.

Table 5\.Mean relative error of\{λ1,…,λ10\}\\\{\\lambda\_\{1\},\\dots,\\lambda\_\{10\}\\\}, by architecture and input representation\.Replacing the landscape\-based representation by the indicator alone increases the error by a factor of66to1010\. That shows that thanks to the landscape function, the model can recover the domain shape\. Furthermore, for FNO and DeepONet, the indicator function based models show a train/validation gap happening after the first2020epochs making the indicator function impossible to generalize from our dataset\.

If we replace the landscape function by the SDF function and keeping the same architecture and addinglog⁡\(SDFmax\)\\log\(\\text\{SDF\}\_\{\\max\}\)results in a test relative error of33\.06%33\.06\\%\. It is even worse for DeepONet model and FNO, with respectively45%45\\%and54\.68%54\.68\\%

Generalization to existing dataset \(MPEG\-7\)To show the generalization ability of our model, we evaluate our final trained model with14001400binary silhouettes from the MPEG\-7 CE\-Shape\-1 dataset\[[9](https://arxiv.org/html/2609.36047#bib.bib3)\]\. This dataset includes7070classes from animals and tools\. Ground truth is computed by finite differences on a200×200200\\times 200grid\. The global mean error to4\.83%4\.83\\%, please refer to Table[6](https://arxiv.org/html/2609.36047#S4.T6)for per eigenvalue error\.

Table 6\.Relative error per eigenvalue on MPEG\-7, mean relative error6\.88%6\.88\\%, median relative error4\.83%4\.83\\%\.Our model performs well both on convex and non\-convex domains\.

Shape optimization case studiesIn a similar fashion as the previous section, we experiment shape optimization by gradient descent through the surrogate and check the result against two classical spectral optimization theorems \(Faber\-Krahn and Hong\-Krahn\-Szego\) and one example from the previous section\.

Forλ1\\lambda\_\{1\}, we start from an ellipse of aspect ratioa/b=6a/b=6, the gradient descent converges toa/b=1\.02a/b=1\.02, withλ^1\\hat\{\\lambda\}\_\{1\}decreasing from104\.96104\.96to36\.8636\.86, only1\.43%1\.43\\%above the exact valuej012​π/A=36\.34j\_\{01\}^\{2\}\\pi/A=36\.34\. The relative error of the resulting landscape function has also a low value of2\.35%2\.35\\%compared to the theoretical landscape\. In a second experiment we start from a disk of fixed radius with two circular holes of radiusrr\.λ^1\\hat\{\\lambda\}\_\{1\}decreases monotonically asr→0r\\to 0in every one of2020testedrr\. \.

Forλ2\\lambda\_\{2\}, it is also well known that the minimizer over domains of constant area is the union of two equal disks\. We start from two eccentric ellipses with size ratioR1/R2=1\.35R\_\{1\}/R\_\{2\}=1\.35, aspect ratiosa1/b1=1\.82a\_\{1\}/b\_\{1\}=1\.82anda2/b2=0\.55a\_\{2\}/b\_\{2\}=0\.55and optimize their three shape parameters \(relative size, both eccentricities\) to minimizeλ^2\\hat\{\\lambda\}\_\{2\}\. The gradient descent converges toR1/R2=0\.99R\_\{1\}/R\_\{2\}=0\.99,a1/b1=1\.08a\_\{1\}/b\_\{1\}=1\.08,a2/b2=1\.05a\_\{2\}/b\_\{2\}=1\.05, whileλ^2\\hat\{\\lambda\}\_\{2\}decreases from111\.52111\.52to74\.3874\.38, only2\.35%2\.35\\%above the exact theoretical value2​π​j012/A=72\.672\\pi j\_\{01\}^\{2\}/A=72\.67, please refer to Figure[5](https://arxiv.org/html/2609.36047#S4.F5)\. Figure[5](https://arxiv.org/html/2609.36047#S4.F5)also shows the minimizer found forλ6\\lambda\_\{6\}presented by Antunes and Freitas\. We start from a near circle and it converges to a rounded triangle\. The resulting predicted eigenvalue has a2\.52\.5relative error compared to the eigenvalue computed with FreeFem\.

![Refer to caption](https://arxiv.org/html/2609.36047v1/two_ellipses_initial.png)

![Refer to caption](https://arxiv.org/html/2609.36047v1/two_ellipses_final.png)

![Refer to caption](https://arxiv.org/html/2609.36047v1/tor_lambda6_initial.png)

![Refer to caption](https://arxiv.org/html/2609.36047v1/tor_lambda6_final.png)

Figure 5\.On the left, we show the starting domain’s boundary, on the right, the boundary of the domain after the gradient descent\. On top,λ^2\\hat\{\\lambda\}\_\{2\}converges to the exact value for two equal disks \(2\.35%2\.35\\%of relative error\)\. After 100 iterationsa1/b1=1\.08a\_\{1\}/b\_\{1\}=1\.08,a2/b2=1\.05a\_\{2\}/b\_\{2\}=1\.05,R1/R2=0\.99R\_\{1\}/R\_\{2\}=0\.99\. On bottom,λ^6\\hat\{\\lambda\}\_\{6\}does not hallucinate with respect to the eigenvalue computed with FreeFem \(2\.5%2\.5\\%of relative error\)\.

## 5\.Conclusion and perspectives

In this study we proposed two very different approaches to tackle the problem of spectral optimization\. They differ in the methods used as well as in their scopes\.

The first approach, the Fourier based model is purposely build to answer the needs of mathematicians looking to support or disprove a given conjecture\. As such, the aim here is to obtain a very light surrogate to the PDE solver that may be used alongside the PDE solver and not just in its stead\. The design we propose is then constrained to use a shape formalism that is also efficient for traditional numerical shape optimization and this limits the possibility for the model to be generalized to larger classes of shapes\. The natural extension of this approach would be to learn other quantities such as the torsional rigidity, the Neumann eigenvalues or geometric quantities\.

The second approach aims to answer the question on a more fundamental level up to loosing some precision\. It uses the landscape function to represent the shapes which is the classical tool used in the analytical aspect of the field, for instance through the theory ofγ\\gamma\-convergence\. It showcases that this function is indeed the right notion of indicator function for sets when dealing with quantities linked to the Dirichlet laplacian\. It would be of interest to see if, using the same procedure, it is also possible to learn the spectrum of capacitary measures which are the natural relaxation of sets in this context and can also be represented with a notion of torsion function\.

### AI use statement

In this work, we used generative AI tools to expand our first codes as we explored new ideas\. We have not used generative AI tools for any other purposes such as drafting the present paper or proving mathematical claims\. We have reviewed all AI assisted work\. All AI produced code has been checked by at least two authors and most of it has been manually rewritten\. We take responsibility for the final content of this work, including text, claims or artifacts produced with the aid of generative AI\.

### Reproducibility statement

For reproducibility purposes, all random parameters and shapes have been generated with a fixed seed and every experiment can be reproduced\.

## References

- \[1\]P\. R\. S\. Antunes and P\. Freitas\(2012\)Numerical optimization of low eigenvalues of the Dirichlet and Neumann laplacians\.J\. Optim\. Theory Appl\.154\(1\),pp\. 235–257\.External Links:[Document](https://dx.doi.org/10.1007/s10957-011-9983-3)Cited by:[§3\.2](https://arxiv.org/html/2609.36047#S3.SS2.p1.1)\.
- \[2\]A\. Bélières, E\. Franck, V\. Michel\-Dansac, Y\. Privat,et al\.\(2025\)Volume\-preserving geometric shape optimization of the dirichlet energy using variational neural networks\.Neural Networks184,pp\. 106957\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.
- \[3\]F\. S\. Costabal, S\. Pezzuto, and P\. Perdikaris\(2024\)Δ\\Delta\-PINNs: physics\-informed neural networks on complex geometries\.Engineering Applications of Artificial Intelligence127,pp\. 107324\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.
- \[4\]G\. David, A\. Gloria, and S\. Mayboroda\(2025\)The landscape function onℝd\\mathbb\{R\}^\{d\}\.Annales de la Faculté des sciences de Toulouse : MathématiquesSer\. 6, 34\(2\),pp\. 315–337\(en\)\.External Links:[Document](https://dx.doi.org/10.5802/afst.1814),[Link](https://afst.centre-mersenne.org/articles/10.5802/afst.1814/)Cited by:[§1\.1](https://arxiv.org/html/2609.36047#S1.SS1.p1.1)\.
- \[5\]F\. Hecht\(2012\)New development in freefem\+\+\.Journal of numerical mathematics20\(3\-4\),pp\. 1–14\.Cited by:[§3\.1](https://arxiv.org/html/2609.36047#S3.SS1.p2.1)\.
- \[6\]K\. Ichimaru, D\. Thomas, T\. Iwaguchi, and H\. Kawasaki\(2025\)Neural sdf for shadow\-aware unsupervised structured light\.In2025 IEEE/CVF Winter Conference on Applications of Computer Vision \(WACV\),pp\. 287–296\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p3.1)\.
- \[7\]X\. Jin, S\. Cai, H\. Li, and G\. E\. Karniadakis\(2021\)NSFnets \(navier\-stokes flow nets\): physics\-informed neural networks for the incompressible navier\-stokes equations\.Journal of Computational Physics426,pp\. 109951\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p7.1)\.
- \[8\]E\. Kharazmi, Z\. Zhang, and G\. E\. Karniadakis\(2019\)Variational physics\-informed neural networks for solving partial differential equations\.arXiv preprint arXiv:1912\.00873\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p7.1),[§4\.1](https://arxiv.org/html/2609.36047#S4.SS1.p5.1)\.
- \[9\]L\.J\. Latecki, R\. Lakamper, and T\. Eckhardt\(2000\)Shape descriptors for non\-rigid shapes with a single closed contour\.InProceedings IEEE Conference on Computer Vision and Pattern Recognition\. CVPR 2000 \(Cat\. No\.PR00662\),Vol\.1,pp\. 424–429 vol\.1\.External Links:[Document](https://dx.doi.org/10.1109/CVPR.2000.855850)Cited by:[§4\.2](https://arxiv.org/html/2609.36047#S4.SS2.p6.1)\.
- \[10\]S\. Li and H\. Salahshoor\(2026\)Finite element eigenfunction network \(feenet\): a hybrid framework for solving pdes on complex geometries\.arXiv preprint arXiv:2602\.00870\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.
- \[11\]Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar\(2020\)Fourier neural operator for parametric partial differential equations\.arXiv preprint arXiv:2010\.08895\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p4.1)\.
- \[12\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis\(2021\)Learning nonlinear operators via deeponet based on the universal approximation theorem of operators\.Nature machine intelligence3\(3\),pp\. 218–229\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p4.1)\.
- \[13\]L\. Lu, X\. Meng, Z\. Mao, and G\. E\. Karniadakis\(2021\)DeepXDE: a deep learning library for solving differential equations\.SIAM review63\(1\),pp\. 208–228\.Cited by:[§1](https://arxiv.org/html/2609.36047#S1.p7.1)\.
- \[14\]J\. J\. Park, P\. Florence, J\. Straub, R\. Newcombe, and S\. Lovegrove\(2019\)Deepsdf: learning continuous signed distance functions for shape representation\.InProceedings of the IEEE/CVF conference on computer vision and pattern recognition,pp\. 165–174\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.
- \[15\]C\. Rowan, J\. Evans, K\. Maute, and A\. Doostan\(2025\)Solving engineering eigenvalue problems with neural networks using the rayleigh quotient\.arXiv preprint arXiv:2506\.04375\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.
- \[16\]R\. Williamson and N\. J\. Mitra\(2025\)Neural geometry processing via spherical neural surfaces\.InComputer Graphics Forum,Vol\.44,pp\. e70021\.Cited by:[§2](https://arxiv.org/html/2609.36047#S2.p7.1)\.