Numerical exploration of the range of shape functionals using neural networks

Technical Summary: Numerical Exploration of the Range of Shape Functionals Using Neural Networks

Problem Statement
The paper addresses the numerical exploration of Blaschke–Santaló diagrams, which are geometric representations of the range of vector-valued shape functionals. These diagrams characterize the possible inequalities relating prescribed shape functionals (e.g., volume, perimeter, eigenvalues) within a specific class of sets, typically convex bodies. While theoretical results exist for purely geometric functionals, obtaining explicit descriptions or qualitative properties for diagrams involving spectral or PDE-related quantities remains highly challenging. Previous numerical approaches, such as random sampling of convex polygons or methods relying on Voronoi tessellation, suffer from non-uniform sampling, failure to capture boundary extremals, or reliance on a priori knowledge of diagram convexity (e.g., vertical convexity) that is difficult to establish for general cases.

Methodology
The authors propose a novel numerical framework combining a specific neural network architecture for shape parametrization with an interacting particle system for uniform sampling.

1. Invertible Neural Network Parametrization:

   • Sublinear Architecture: The method parametrizes convex bodies in $\mathbb{R}^d$ as the image of the unit ball $B$ under a homeomorphism $\phi_\theta$. This map is constructed using a gauge function $p_\theta$, which is represented by a neural network.
   • Convexity Conservation: To ensure the resulting sets are inherently convex, the network architecture is designed to satisfy sublinearity constraints. Specifically, the gauge function is defined as $p_\theta(x) = \beta \|x\| \text{LSE}(W^T x / \|x\|)$, where LSE is the LogSumExp function. This smooth approximation of the maximum of linear forms ensures $p_\theta$ is continuous, sublinear, and $C^\infty$ away from the origin.
   • Universal Approximation: The authors rely on a "Universal Approximation Theorem" (proven in prior work [44]) stating that the set of convex bodies representable by this architecture is dense in the space of all convex bodies with respect to the Hausdorff distance.
   • Symmetry Encoding: The framework allows for the imposition of symmetries (e.g., reflectional symmetry) by averaging the gauge function over a finite group of isometries, ensuring the resulting shape is invariant under that group.

2. Uniform Sampling via Interacting Particles:

   • To generate a uniform distribution of points within the Blaschke–Santaló diagram, the authors minimize a Riesz energy functional. This is formulated as an optimization problem where $N$ shapes $\Omega_1, \dots, \Omega_N$ are optimized to maximize their pairwise distances in the functional space:
     $$\min_{\theta_1, \dots, \theta_N} \sum_{1 \le i \neq j \le N} |F(\Omega_{\theta_i}) - F(\Omega_{\theta_j})|^{-s}$$
   • This approach is motivated by the "Poppy-seed Bagel Theorem," which guarantees that minimizing such energy leads to a uniform distribution on the manifold as $N$ increases.
   • The optimization is performed using automatic differentiation in PyTorch (via l-BFGS), allowing gradients to flow through the entire computational graph, including the solution of PDEs and integral evaluations.

3. Computation of Shape Functionals:

   • Integral Quantities: Volume, perimeter, and moment of inertia are computed by pulling back integrals from the shape domain $\Omega_\theta$ to the reference unit ball $B$ using change-of-variable formulas involving the Jacobian of the neural map.
   • PDE-Related Quantities:
     • Torsional Rigidity: Solved using a Method of Fundamental Solutions (MFS) approach, where the solution is approximated by a linear combination of fundamental solutions. The boundary condition is fitted via least squares, a process fully differentiable in PyTorch.
     • Neumann Eigenvalues: Computed using a mesh-free RBF-Galerkin method. The weak formulation is pulled back to the reference domain, discretized using Radial Basis Functions, and solved as a generalized eigenvalue problem. To enable differentiation, the problem is transformed into a standard eigenvalue problem via Cholesky decomposition of the mass matrix.
   • Regularization: To prevent the optimization from collapsing into flat, degenerate shapes (which causes numerical instability), a regularizer based on the condition number of the Jacobian is added to the loss function.

Key Contributions

• Unconstrained Convex Optimization: The introduction of a smooth, invertible neural network parametrization that intrinsically preserves convexity without requiring explicit inequality constraints or penalization terms.
• Uniform Diagram Sampling: The application of interacting particle systems to minimize Riesz energy, providing a method to uniformly sample Blaschke–Santaló diagrams without relying on prior knowledge of the diagram's convexity properties.
• Differentiable PDE Solvers: The successful integration of automatic differentiation with mesh-free methods (MFS and RBF-Galerkin) to compute and differentiate spectral and torsional functionals with respect to shape parameters.
• Open-Source Implementation: The release of a comprehensive PyTorch package (ConvexShapeOpt) and notebooks to facilitate reproducibility.

Results
The method was tested on various diagrams in $\mathbb{R}^2$ and $\mathbb{R}^3$, demonstrating effectiveness across geometric and PDE-based functionals:

• Geometric Diagrams (Volume, Perimeter, Moment of Inertia): The method successfully recovered known theoretical bounds and extremal shapes (e.g., stadiums, triangles, rhombi) for both general and symmetric convex bodies in 2D and 3D.
• Torsional Rigidity Diagrams: The approach accurately mapped the relationship between volume, perimeter, and torsional rigidity, recovering known inequalities (e.g., Saint-Venant, Polya, Makai) and visualizing the corresponding optimal shapes.
• Spectral Diagrams (Volume, First Two Neumann Eigenvalues): The method explored the diagram for $\mu_1$ and $\mu_2$. While it recovered the general shape of the domain, the authors note limitations in precision for the upper boundary due to memory constraints associated with differentiating through the eigenvalue assembly.
• Willmore Energy Diagrams: The framework successfully handled higher-order boundary terms, visualizing the trade-offs between volume, perimeter, and Willmore energy for planar and 3D convex sets.
• Comparison: Comparisons with naive random sampling of convex polygons showed that the proposed method captures regions of the diagram that random sampling misses, particularly near the boundaries.

Significance and Claims
The paper claims to provide a "novel numerical framework" that overcomes the limitations of previous methods by offering a "fully unconstrained" approach to shape optimization. The significance lies in:

1. Flexibility: The ability to handle diverse functionals (geometric, spectral, PDE-based) within a single, unified differentiable framework.
2. Robustness: The method provides a "satisfying description" of diagrams in both 2D and 3D without requiring a priori assumptions about the diagram's geometry (such as vertical convexity).
3. Scalability: While computational costs increase in 3D due to quadrature requirements, the method remains robust and does not require additional implementation complexities compared to 2D.
4. Accessibility: By leveraging PyTorch's automatic differentiation, the framework lowers the barrier for exploring new functionals and is made available to the research community.

The authors maintain a modest tone regarding limitations, acknowledging that the method's memory efficiency is constrained by the size of the computational graph required for high-precision differentiation, particularly for spectral problems. They also note that while the gauge function parametrization is powerful, it makes certain geometric quantities (like minimal width) harder to compute than with support functions, though the two are theoretically interchangeable.