Equivariant Covariance Tensors: Guaranteed SPD Uncertainty for Tensor-Valued Geometric Learning

# Guaranteed SPD Uncertainty for Tensor-Valued Geometric Learning
Source: [https://arxiv.org/html/2608.24386](https://arxiv.org/html/2608.24386)
## Equivariant Covariance Tensors: Guaranteed SPD Uncertainty for Tensor\-Valued Geometric Learning

Ruihan LiuAffiliation:College of Intelligent Robotics and Advanced Manufacturing, Fudan University, Shanghai 200433, ChinaYu JiAffiliation:College of Intelligent Robotics and Advanced Manufacturing, Fudan University, Shanghai 200433, ChinaJianbo YuAffiliation:School of Microelectronics, Fudan University, Shanghai 200433, ChinaCorrespondence to:[jb\_yu@fudan\.edu\.cn](mailto:[email protected])Qingchao JiangAffiliation:School of Information Science and Engineering, East China University of Science and Technology, Shanghai 200237, China

###### Abstract

Tensor\-valued prediction is fundamental to geometric deep learning, yet uncertainty quantification \(UQ\) for such outputs remains an open challenge\. While E\(3\)\-equivariant neural networks excel at point estimates, they lack rigorous confidence measures\. We focus on symmetric rank\-2 tensor prediction, where the target has six Kelvin–Mandel coordinates and full uncertainty is represented by a6×66\\times 6covariance matrix\. We introduce a framework for E\(3\)\-equivariant UQ, modeling the full predictive distribution where both mean and covariance preserve rotational symmetry\. Our approach decomposes the covariance into irreducible representationsSym2​\(ρc\)≅2×\(l=0\)⊕2×\(l=2\)⊕1×\(l=4\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\cong 2\\times\(l=0\)\\oplus 2\\times\(l=2\)\\oplus 1\\times\(l=4\)\. By mapping from the flat Lie algebra𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)to the curved SPD manifold via matrix exponentiation, we strictly ensure positive\-definite covariances while maintaining exact equivariance\. Furthermore, we formulate a Log\-Euclidean Equivariant Scoring Objective \(LE\-ESO\)—a robust surrogate loss based on the Multivariate Laplace distribution—providing robustness to heavy\-tailed errors and stable optimization\. Validation on ModelNet40 inertia tensors and Materials Project dielectric tensors demonstrates that our method achieves competitive performance and provides physically consistent, symmetry\-preserving uncertainty estimates with useful risk and OOD sensitivity\.

###### Keywords:

Machine Learning, Equivariant Neural Networks, Uncertainty Quantification, Geometric Deep Learning, Tensor\-Valued Prediction

## 1Introduction

Tensor\-valued predictions are fundamental to scientific computing and geometric deep learning, with applications spanning material properties \(elasticity tensors, dielectric response\), biomedical imaging \(diffusion MRI\), and computational fluid dynamics\. While E\(3\)\-equivariant neural networks \(ENNs\) have achieved remarkable success in predicting tensorial properties like dielectric constants\([Batatia et al\., 2022](https://arxiv.org/html/2608.24386#bib.bib12);[Heilman et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib3)\), they are inherently deterministic—outputting single point estimates without any confidence measure\. This is a critical limitation: overconfident tensor predictions in scientific decisions can lead to costly experimental failures\.

In this paper, the primary output is a symmetric rank\-2 tensorC∈ℝsym3×3C\\in\\mathbb\{R\}^\{3\\times 3\}\_\{\\text\{sym\}\}, such as a dielectric tensor\. AlthoughCCis a3×33\\times 3matrix, symmetry leaves six independent degrees of freedom\. We represent it by the Kelvin\-Mandel vector𝐜∈ℝ6\\mathbf\{c\}\\in\\mathbb\{R\}^\{6\}, for which rotations act by an orthogonal six\-dimensional representationρc​\(R\)\\rho\_\{c\}\(R\)\. Therefore, uncertainty overCCis not a3×33\\times 3covariance, but a6×66\\times 6covariance over the Kelvin\-Mandel coordinates\.

Extending ENNs with uncertainty quantification poses a fundamental challenge: uncertainty estimates must themselves transform equivariantly\. Mathematically, the covarianceΣ∈ℝ6×6\\Sigma\\in\\mathbb\{R\}^\{6\\times 6\}must satisfy:

Σ⁡\(R⋅X\)=ρc​\(R\)​Σ​\(X\)​ρc​\(R\)⊤,∀R∈O⁡\(3\),\\Sigma\(R\\\!\\cdot\\\!X\)=\\rho\_\{c\}\(R\)\\Sigma\(X\)\\rho\_\{c\}\(R\)^\{\\top\},\\quad\\forall R\\in O\(3\),whereρc​\(R\)\\rho\_\{c\}\(R\)is the rotation matrix in Kelvin\-Mandel notation\. This creates a parameterization challenge: Cholesky\-type covariance heads guarantee SPD but are not equivariant in Kelvin–Mandel covariance coordinates, while direct equivariant regression preserves the transformation law but does not guarantee SPD\.

In this work, we introduce a principled framework for E\(3\)\-equivariant full\-covariance uncertainty quantification in symmetric tensor\-valued prediction\. Our contributions address this equivariant SPD parameterization challenge through: \(1\) an equivariant matrix\-exponential head parameterizing the SPD manifold via irreducible decompositionSym2​\(ρc\)≅2×\(ℓ=0\)⊕2×\(ℓ=2\)⊕1×\(ℓ=4\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\cong 2\\times\(\\ell=0\)\\oplus 2\\times\(\\ell=2\)\\oplus 1\\times\(\\ell=4\); \(2\) a stable joint optimization strategy via a Log\-Euclidean scoring objective that integrates uncertainty calibration with geometric feature extraction; and \(3\) rigorous validation on ModelNet40 inertia tensors with competitive performance on Materials Project dielectric prediction \(MAE 1\.55\)\.

## 2Related Work

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

Figure 1:E\(3\)\-equivariant tensor uncertainty framework\.\(1\) Feature Decomposition: covariance head models symmetry via irreducible representations \(ℓ=0,2,4\\ell=0,2,4\)\. \(2\) Manifold Mapping: network predicts unconstrainedA∈𝔰​𝔶​𝔪​\(6\)A\\in\\mathfrak\{sym\}\(6\)\. \(3\) Validity Constraint:Σ=exp⁡\(A\)\\Sigma=\\exp\(A\)projects to SPD manifold𝒫6\\mathcal\{P\}\_\{6\}\.#### Equivariant Tensor Prediction\.

E\(3\)\-equivariant neural networks \(ENNs\) achieve state\-of\-the\-art tensor prediction across materials science and 3D geometry\([Batatia et al\., 2022](https://arxiv.org/html/2608.24386#bib.bib12);[Heilman et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib3);[Hua et al\., 2026](https://arxiv.org/html/2608.24386#bib.bib29);[Pakornchote et al\., 2023](https://arxiv.org/html/2608.24386#bib.bib4);[Fung et al\., 2021](https://arxiv.org/html/2608.24386#bib.bib1);[Reiser et al\., 2022](https://arxiv.org/html/2608.24386#bib.bib5);[Du et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib21);[Equer et al\., 2023](https://arxiv.org/html/2608.24386#bib.bib19)\)\. However, all existing ENNs are inherently deterministic—they cannot quantify confidence in their predictions\. Our framework preserves equivariance guarantees while adding rigorous uncertainty quantification, addressing this critical gap\. While our approach utilizes the spherical harmonic basis provided by e3nn\([Geiger and Smidt, 2022](https://arxiv.org/html/2608.24386#bib.bib26)\), recent advances in representation theory, such as the High\-Rank Irreducible Cartesian Tensor \(ICT\) decomposition framework\([Shao et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib30)\), offer alternative analytical paths for constructing equivariant bases directly in Cartesian space\. Such methods could potentially simplify the implementation for higher\-rank tensor properties\.

#### SPD Constraints in Neural Networks\.

Ensuring symmetric positive\-definite \(SPD\) covariances is essential for valid uncertainty\. Prior work uses Cholesky decomposition, eigendecomposition, or Riemannian optimization\([Jekel et al\., 2022](https://arxiv.org/html/2608.24386#bib.bib16);[Pouliquen et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib7);[Zhao et al\., 2023](https://arxiv.org/html/2608.24386#bib.bib14)\)\. These methods guarantee SPD in a fixed coordinate parameterization, but they do not by themselves provide an equivariant mapX↦Σ⁡\(X\)X\\mapsto\\Sigma\(X\)under the covariance representationρc\\rho\_\{c\}\. Orthogonal conjugation preserves SPD; thus the difficulty is not a conflict between the SPD property and rotation itself\. Rather, the challenge is to design a neural parameterization that is simultaneously equivariant and constrained to the SPD cone\. Cholesky\-type heads enforce SPD but are not equivariant in Kelvin\-Mandel covariance coordinates, whereas direct equivariant regression can preserve the transformation law but does not guarantee SPD\. Our matrix\-exponential head resolves this parameterization problem by predicting an equivariant symmetric operatorA⁡\(X\)A\(X\)and settingΣ⁡\(X\)=exp⁡\(A⁡\(X\)\)\\Sigma\(X\)=\\exp\(A\(X\)\), which satisfiesexp⁡\(ρc​\(R\)​A​ρc​\(R\)⊤\)=ρc​\(R\)​exp⁡\(A\)​ρc​\(R\)⊤\\exp\(\\rho\_\{c\}\(R\)A\\rho\_\{c\}\(R\)^\{\\top\}\)=\\rho\_\{c\}\(R\)\\exp\(A\)\\rho\_\{c\}\(R\)^\{\\top\}and therefore is jointly SPD and equivariant in all rotated frames\.

#### Probabilistic Methods and Equivariant GPs\.

Probabilistic extensions for ENNs are severely underdeveloped\. While ensemble methods can provide heuristic uncertainty estimates\([Rudner et al\., 2022](https://arxiv.org/html/2608.24386#bib.bib15)\)and even exhibit emergent equivariance\([Gerken and Kessel, 2024](https://arxiv.org/html/2608.24386#bib.bib31)\), our framework explicitly learns the full 21\-parameter aleatoric covariance tensor, which is essential for modeling the inherent anisotropic noise in physical properties\. Bayesian neural networks face computational challenges at scale\([Rensmeyer et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib8);[Olivier et al\., 2021](https://arxiv.org/html/2608.24386#bib.bib6);[Doan et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib23);[Sheinkman and Wade, 2025](https://arxiv.org/html/2608.24386#bib.bib24)\)\. Existing equivariant Bayesian approaches focus on scalar or vector quantities\([Zhou et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib9)\), missing tensor\-valued uncertainty\. E\(3\)\-equivariant Gaussian processes provide theoretically sound uncertainties but are typically restricted to isotropic or diagonal covariance approximations\([Steinert et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib10);[Bevanda et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib13)\)—scaling them to the full 21\-parameter covariance structure of rank\-2 tensors remains computationally prohibitive\. Our neural network approach offers computational scalability while learning complete equivariant tensor correlations\. Recent work on uncertainty calibration and stable tensor operations provides theoretical grounding for our design\([Berman et al\., 2026](https://arxiv.org/html/2608.24386#bib.bib11);[Gruber and Buettner, 2022](https://arxiv.org/html/2608.24386#bib.bib17);[Fakour et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib18);[Newman et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib20)\)\.

## 3Methods

### 3\.1Problem Formulation

We address the fundamental challenge of predicting tensor\-valued quantities while quantifying predictive uncertainty\. We present the construction for symmetric rank\-2 tensorsC∈ℝsym3×3C\\in\\mathbb\{R\}^\{3\\times 3\}\_\{\\text\{sym\}\}, which already require a full6×66\\times 6covariance representation\. The SPD construction and scoring objective are representation\-agnostic once an equivariant symmetric operatorA⁡\(X\)A\(X\)is available\. However, the parameterization ofA⁡\(X\)A\(X\)is representation\-specific and must be constructed separately for each tensor order and symmetry group\. We predictCCfrom a 3D structureX=\{\(𝐱i,fi\)\}i=1NX=\\\{\(\\mathbf\{x\}\_\{i\},f\_\{i\}\)\\\}\_\{i=1\}^\{N\}, where𝐱i∈ℝ3\\mathbf\{x\}\_\{i\}\\in\\mathbb\{R\}^\{3\}denotes spatial coordinates andfif\_\{i\}represents associated features \(such as atomic species in materials\)\. This formulation encompasses important applications such as predicting dielectric tensors in crystal structures\. Although we demonstrate this framework on materials science, the approach is applicable to any 3D point cloud with rank\-2 tensorial attributes, ranging from biological molecules to geometric shapes\.

Rather than producing deterministic point estimates, we model the full predictive distribution

p⁡\(C∣X\)=Laplace​\(C∣μ⁡\(X\),Σ⁡\(X\)\),p\(C\\mid X\)=\\text\{Laplace\}\(C\\mid\\mu\(X\),\\Sigma\(X\)\),\(1\)whereμ⁡\(X\)∈ℝsym3×3\\mu\(X\)\\in\\mathbb\{R\}^\{3\\times 3\}\_\{\\text\{sym\}\}is the predicted mean tensor andΣ⁡\(X\)\\Sigma\(X\)captures the predictive uncertainty through a covariance structure\.

To maintain geometric consistency across all domains, both the mean predictionμ\\muand covarianceΣ\\Sigmamust satisfy equivariance constraints with respect to rotations and reflections\. Since we predict global tensor properties \(e\.g\., dielectric tensor of a unit cell\), the predictions are invariant to translations but equivariant to orthogonal transformations:

μ⁡\(R⋅X\)\\displaystyle\\mu\(R\\cdot X\)=R​μ​\(X\)​R⊤,\\displaystyle=R\\mu\(X\)R^\{\\top\},Σ⁡\(R⋅X\)\\displaystyle\\Sigma\(R\\cdot X\)=ρc​\(R\)​Σ​\(X\)​ρc​\(R\)⊤\.\\displaystyle=\\rho\_\{c\}\(R\)\\,\\Sigma\(X\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\.\(2\)Our construction ensures fullO⁡\(3\)O\(3\)equivariance\. For symmetric rank\-2 tensors, the transformation under reflection \(detR=−1\\det R=\-1\) is handled by the even\-parity representation of the Kelvin\-Mandel basisρc​\(R\)=R⊗sR\\rho\_\{c\}\(R\)=R\\otimes\_\{s\}R, where\(detR\)2=1\(\\det R\)^\{2\}=1guarantees consistent behavior for both chiral and achiral structures\. We numerically verifyO⁡\(3\)O\(3\)equivariance in Section[4\.5](https://arxiv.org/html/2608.24386#S4.SS5)\(Table[4](https://arxiv.org/html/2608.24386#S4.T4)\), with detailed analysis in Appendix[E\.3](https://arxiv.org/html/2608.24386#A5.SS3)\. Predictions are translation\-invariant since global tensor properties depend only on relative atomic positions\. This universal equivariance constraint forms the foundation of our domain\-agnostic framework\. We emphasize that while Eq\.[1](https://arxiv.org/html/2608.24386#S3.E1)presents the distribution in standard predictive form, the training optimizes the robustified LE\-ESO objective \(Section[3\.5](https://arxiv.org/html/2608.24386#S3.SS5.SSS0.Px2)\) which generalizes the Multivariate Laplace negative log\-likelihood with enhanced stability against outliers\.

### 3\.2Voigt Representation and Covariance Structure

To facilitate neural network implementation while preserving the tensor’s geometric structure, we employ the Kelvin\-Mandel notation to flatten the symmetric tensorCCinto a66\-dimensional vector\. Unlike standard Voigt representation, Kelvin\-Mandel notation maintains the isometry property between tensor and vector spaces:

𝐜KM=\[C11,C22,C33,2​C23,2​C13,2​C12\]⊤,\\mathbf\{c\}\_\{\\text\{KM\}\}=\[C\_\{11\},C\_\{22\},C\_\{33\},\\sqrt\{2\}C\_\{23\},\\sqrt\{2\}C\_\{13\},\\sqrt\{2\}C\_\{12\}\]^\{\\top\},\(3\)
which preserves the Frobenius norm:‖𝐜KM‖2=‖C‖F\\\|\\mathbf\{c\}\_\{\\text\{KM\}\}\\\|\_\{2\}=\\\|C\\\|\_\{F\}\. Under rotationRR,𝐜KM′=ρc​\(R\)​𝐜KM\\mathbf\{c\}\_\{\\text\{KM\}\}^\{\\prime\}=\\rho\_\{c\}\(R\)\\mathbf\{c\}\_\{\\text\{KM\}\}whereρc​\(R\)\\rho\_\{c\}\(R\)is an orthogonal6×66\\times 6matrix\. The covariance transforms asΣ′=ρc​\(R\)​Σ​ρc​\(R\)⊤\\Sigma^\{\\prime\}=\\rho\_\{c\}\(R\)\\Sigma\\rho\_\{c\}\(R\)^\{\\top\}, maintaining coordinate invariance of the physical uncertainty\.

### 3\.3Irreducible Representation Decomposition

The mathematical structure of equivariant uncertainty quantification becomes clear through representation theory\. In group theory, irreducible representations are fundamental building blocks that cannot be further decomposed into smaller invariant subspaces\. The 6D representationρc\\rho\_\{c\}for symmetric3×33\\times 3tensors decomposes into irreducible representations ofSO⁡\(3\)\\mathrm\{SO\}\(3\)as

ρc≅l=0⊕l=2,\\rho\_\{c\}\\cong l=0\\oplus l=2,\(4\)corresponding respectively to the isotropic \(trace\) and deviatoric \(traceless\) components of the tensor\. Here≅\\congdenotes an isomorphism of representations, not equality of matrices: after a fixed change of basis, the six Kelvin\-Mandel coordinates split into a one\-dimensional isotropic trace component and a five\-dimensional traceless deviatoric component\.

SinceΣ\\Sigmatransforms asρc⊗ρc\\rho\_\{c\}\\otimes\\rho\_\{c\}, its representation decomposes as:

ρc⊗ρc\\displaystyle\\rho\_\{c\}\\otimes\\rho\_\{c\}=\(l=0⊕l=2\)⊗\(l=0⊕l=2\)\\displaystyle=\(l=0\\oplus l=2\)\\otimes\(l=0\\oplus l=2\)=\(l=0\)⊕2​\(l=2\)⊕\(l=4\)⊕\(l=1,3\)antisym\.\\displaystyle=\(l=0\)\\oplus 2\(l=2\)\\oplus\(l=4\)\\oplus\(l=1,3\)\_\{\\text\{antisym\}\}\.\(5\)Theℓ=1\\ell=1andℓ=3\\ell=3components lie in the antisymmetric part ofρc⊗ρc\\rho\_\{c\}\\otimes\\rho\_\{c\}, corresponding to operators that change sign under exchange of the two covariance indices\. Since covariance matrices satisfyΣ=Σ⊤\\Sigma=\\Sigma^\{\\top\}, only the symmetric squareSym2​\(ρc\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)remains, yielding:

Sym2​\(ρc\)≅2×\(l=0\)⊕2×\(l=2\)⊕1×\(l=4\),\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\cong 2\\times\(l=0\)\\oplus 2\\times\(l=2\)\\oplus 1\\times\(l=4\),\(6\)This provides 21 independent degrees of freedom for symmetric6×66\\times 6SPD covariance matrices\.

### 3\.4Equivariant Neural Architecture

Figure[1](https://arxiv.org/html/2608.24386#S2.F1)summarizes the data flow\. A shared E\(3\)\-equivariant encoder first maps the input structureXXto latent irreducible features\. The mean head selects theℓ=0⊕ℓ=2\\ell=0\\oplus\\ell=2components and outputs the Kelvin\-Mandel mean vector𝝁KM​\(X\)∈ℝ6\\boldsymbol\{\\mu\}\_\{\\text\{KM\}\}\(X\)\\in\\mathbb\{R\}^\{6\}, which is mapped back to a symmetric3×33\\times 3tensor\. The covariance head outputs coefficients in2×\(ℓ=0\)⊕2×\(ℓ=2\)⊕1×\(ℓ=4\)2\\times\(\\ell=0\)\\oplus 2\\times\(\\ell=2\)\\oplus 1\\times\(\\ell=4\), which are assembled into an equivariant symmetric operatorA⁡\(X\)∈𝔰​𝔶​𝔪​\(6\)A\(X\)\\in\\mathfrak\{sym\}\(6\)\. Finally, the predictive covariance isΣ⁡\(X\)=exp⁡\(A⁡\(X\)\)\\Sigma\(X\)=\\exp\(A\(X\)\), and the loss is computed from the Kelvin\-Mandel residualΔ​𝐜=𝐜KM−𝝁KM​\(X\)\\Delta\\mathbf\{c\}=\\mathbf\{c\}\_\{\\text\{KM\}\}\-\\boldsymbol\{\\mu\}\_\{\\text\{KM\}\}\(X\)\.

The covariance headfΣ:X↦A⁡\(X\)f\_\{\\Sigma\}:X\\mapsto A\(X\)must satisfy the transformation property

A⁡\(R⋅X\)=ρc​\(R\)​A​\(X\)​ρc​\(R\)⊤\.A\(R\\\!\\cdot\\\!X\)=\\rho\_\{c\}\(R\)\\,A\(X\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\.\(7\)
Following the irreducible representation decomposition in Eq\.[6](https://arxiv.org/html/2608.24386#S3.E6), we constructA⁡\(X\)A\(X\)through structured Clebsch\-Gordan combinations\. Letϕ\(L\)​\(X\)∈ℝFL×\(2​L\+1\)\\phi^\{\(L\)\}\(X\)\\in\\mathbb\{R\}^\{F\_\{L\}\\times\(2L\+1\)\}denote spherical tensor features of orderLLproduced by the equivariant backbone\. We index the irreps appearing inSym2​\(ρc\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)by

ℐ=\{\(0,1\),\(0,2\),\(2,1\),\(2,2\),\(4,1\)\},\\mathcal\{I\}=\\\{\(0,1\),\\,\(0,2\),\\,\(2,1\),\\,\(2,2\),\\,\(4,1\)\\\},where\(L,r\)\(L,r\)refers to therr\-th copy of theLL\-irrep, and write

A⁡\(X\)=∑\(L,r\)∈ℐ∑m=−LLaL,m\(r\)​\(X\)​BL,m\(r\),A\(X\)=\\sum\_\{\(L,r\)\\in\\mathcal\{I\}\}\\sum\_\{m=\-L\}^\{L\}a^\{\(r\)\}\_\{L,m\}\(X\)\\,B^\{\(r\)\}\_\{L,m\},\(8\)with the equivariant coefficients

aL,m\(r\)​\(X\)=∑f=1FLwL,r,f​ϕf,m\(L\)​\(X\)\.a^\{\(r\)\}\_\{L,m\}\(X\)=\\sum\_\{f=1\}^\{F\_\{L\}\}w\_\{L,r,f\}\\,\\phi^\{\(L\)\}\_\{f,m\}\(X\)\.HereBL,m\(r\)∈ℝsym6×6B^\{\(r\)\}\_\{L,m\}\\in\\mathbb\{R\}^\{6\\times 6\}\_\{\\text\{sym\}\}are*fixed, input\-independent*basis matrices for therr\-th copy of theLL\-irrep inSym2​\(ρc\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\), computed once from the Clebsch\-Gordan decomposition\. All input dependence resides in the equivariant coefficientsaL,m\(r\)​\(X\)a^\{\(r\)\}\_\{L,m\}\(X\)\. Under rotation, the coefficient vector\(aL,m\(r\)\)m=−LL\(a^\{\(r\)\}\_\{L,m\}\)\_\{m=\-L\}^\{L\}and the basis\(BL,m\(r\)\)m=−LL\(B^\{\(r\)\}\_\{L,m\}\)\_\{m=\-L\}^\{L\}transform by the same Wigner\-D\(L\)D^\{\(L\)\}representation, so their contraction yieldsA⁡\(R⋅X\)=ρc​\(R\)​A​\(X\)​ρc​\(R\)⊤A\(R\\cdot X\)=\\rho\_\{c\}\(R\)\\,A\(X\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\. This explicit tensor basis construction guarantees that any choice of weightswL,r,fw\_\{L,r,f\}preserves the equivariance property, providing hard\-constrained geometric consistency rather than soft\-regularized approximation\. We note that while we rely on spherical tensor products, the orthogonal ICT decomposition matrices\([Shao et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib30)\)provide an equivalent and highly efficient basis for higher\-order Cartesian tensors, which may offer computational advantages for future extensions to rank\-4 tensors like elasticity\.

Our architecture implements an E\(3\)\-equivariant neural network using thee3nnlibrary, following the established paradigm of equivariant message passing\. The implementation proceeds through two key stages\. First, an equivariant message passing backbone processes the atomic structure to produce latent features transforming under mixed irreps up toℓm​a​x=4\\ell\_\{max\}=4\. Second, an equivariant linear layer—implemented via a fourth\-order Cartesian tensor with symmetryi​j​k​l=j​i​k​l=i​j​l​k=k​l​i​jijkl=jikl=ijlk=klij—assembles these latent features into the symmetric block structure ofAA, ensuring consistency with theSym2​\(ρc\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)representation while automatically filtering to thel=0,2,4l=0,2,4components required for rank\-4 covariance output\. The covariance head uses a residual connectionA=Abase⋅I\+Δ​AA=A\_\{\\text\{base\}\}\\cdot I\+\\Delta A, whereAbaseA\_\{\\text\{base\}\}provides an isotropic baseline andΔ​A\\Delta Acaptures the anisotropic uncertainty structure learned from the data\.

We employ an end\-to\-end joint optimization strategy, where both the mean and covariance heads are trained simultaneously\. The inherent stability of our matrix\-exponential mapping and the LE\-ESO loss eliminates the need for gradient detachment, allowing the backbone to learn geometric features that are mutually informative for both point prediction and uncertainty quantification\. This design guarantees that the raw network outputAAnaturally possesses the correct equivariance properties, setting the stage for positive\-definite covariance construction\.

### 3\.5Positive\-Definite Covariance Construction and Training Objective

![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/riemannian_diagram_icml_4k.png)Figure 2:Geometric interpretation of equivariant covariance construction\.The network operates in the Lie algebra𝔤≅ℝsym6×6\\mathfrak\{g\}\\cong\\mathbb\{R\}\_\{\\mathrm\{sym\}\}^\{6\\times 6\}\(left\), producingA⁡\(X\)A\(X\)\. The exponential mapΣ=exp⁡\(A\)\\Sigma=\\exp\(A\)\(center\) projects onto the SPD manifold𝒫6\\mathcal\{P\}\_\{6\}\(right\), guaranteeing valid uncertainties while preservingE⁡\(3\)E\(3\)\-equivariance\.#### From Curved Manifold to Flat Tangent Space\.

To strictly enforce the SPD constraint while maintaining equivariance, we leverage the Log\-Euclidean framework\([Arsigny et al\., 2006](https://arxiv.org/html/2608.24386#bib.bib2)\)and parameterizeΣ⁡\(X\)\\Sigma\(X\)via the matrix exponential mapping:

Σ⁡\(X\)=exp⁡\(A⁡\(X\)\),\\Sigma\(X\)=\\exp\(A\(X\)\),\(9\)We thus optimize within the tangent space𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)—a flat Euclidean vector space at the identity—and the matrix exponential liftsA⁡\(X\)A\(X\)to the SPD manifold𝒫6\\mathcal\{P\}\_\{6\}while preservingE⁡\(3\)E\(3\)\-equivariance, sinceρc​\(R\)\\rho\_\{c\}\(R\)acts by orthogonal conjugation \(Figure[2](https://arxiv.org/html/2608.24386#S3.F2)\)\.

#### Log\-Euclidean Equivariant Scoring Objective \(LE\-ESO\)\.

Given the Kelvin\-Mandel residualΔ​𝐜=𝐜KM−𝝁KM​\(X\)\\Delta\\mathbf\{c\}=\\mathbf\{c\}\_\{\\text\{KM\}\}\-\\boldsymbol\{\\mu\}\_\{\\text\{KM\}\}\(X\), we define the Log\-Euclidean Mahalanobis distance

DM​\(A,Δ​𝐜\)=Δ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜\.D\_\{M\}\(A,\\Delta\\mathbf\{c\}\)=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\}\.\(10\)To control the influence of extreme outliers, we use a robustified distanceD~M\\tilde\{D\}\_\{M\}with transition thresholdτ\\tau:

D~M=\{DM,DM<τ,τ\+log⁡\(1\+DM−τ\),DM≥τ\.\\tilde\{D\}\_\{M\}=\\begin\{cases\}D\_\{M\},&D\_\{M\}<\\tau,\\\\ \\tau\+\\log\(1\+D\_\{M\}\-\\tau\),&D\_\{M\}\\geq\\tau\.\\end\{cases\}\(11\)The final LE\-ESO objective combines uncertainty volume regularization with the \(robustified\) data\-fit term:

ℒLE\-ESO=α​Tr⁡\(A\)\+D~M,\\mathcal\{L\}\_\{\\text\{LE\-ESO\}\}=\\alpha\\operatorname\{Tr\}\(A\)\+\\tilde\{D\}\_\{M\},\(12\)whereα\>0\\alpha\>0controls the trade\-off between uncertainty volume and data fit\.

Whenα=1\\alpha=1andD~M=DM\\tilde\{D\}\_\{M\}=D\_\{M\}, this is the multivariate Laplace negative log\-likelihood in Log\-Euclidean form, a strictly proper scoring rule on𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)\([Gneiting and Raftery, 2007](https://arxiv.org/html/2608.24386#bib.bib32)\)\. With log\-tail robustification orα≠1\\alpha\\neq 1, the objective is arobust surrogate scoring objective, trading strict propriety for outlier stability\. We setα=1\\alpha=1in our primary experiments; Appendix[D\.4](https://arxiv.org/html/2608.24386#A4.SS4)reports a validation sweep overα∈\{0\.03,0\.1,0\.3,1\.0\}\\alpha\\in\\\{0\.03,0\.1,0\.3,1\.0\\\}showing stable training, withα=1\\alpha=1also giving the lowest validation MAE\. Detailed derivation and gradient analysis are in Appendix[A\.4](https://arxiv.org/html/2608.24386#A1.SS4)\.

#### Training Stability via Joint Optimization\.

In our experiments, the Lie algebra parameterization combined with eigenvalue clamping and the log\-tail robustification in Eq\.[11](https://arxiv.org/html/2608.24386#S3.E11)keeps gradients and the matrix exponential numerically well\-behaved during training, enabling stable end\-to\-end joint optimization of both the mean and covariance heads without gradient detachment, so the backbone can receive informative gradients from the uncertainty quantification objective\. This joint training paradigm allows the learned geometric features to be optimized for both prediction accuracy and uncertainty calibration, without sacrificing numerical stability or equivariance guarantees\.

## 4Experiments

![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/uncertainty_3d_separated.png)Figure 3:3D uncertainty visualization on ModelNet40\.\(Top\) Point clouds with uncertainty ellipsoids \(blue\) and principal axes: predicted \(red\) vs ground truth \(green\)\. \(Bottom\) Anisotropic covariance matrices\. Samples show varying errors \(0\.293–0\.472\) and Mahalanobis distances \(0\.37–0\.50\), demonstrating geometric\-adaptive uncertainty that rotates with the object shape\.We evaluate the proposed framework in two main settings: \(i\) controlled geometric validation on ModelNet40 inertia tensors\([Wu et al\., 2015](https://arxiv.org/html/2608.24386#bib.bib27)\), and \(ii\) real\-data dielectric tensor prediction on the Materials Project\([Barroso\-Luque et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib22);[Jain et al\., 2013](https://arxiv.org/html/2608.24386#bib.bib25)\)\. These main experiments are complemented by additional studies in Appendix[D](https://arxiv.org/html/2608.24386#A4), including ModelNet40 shape\-covariance prediction \(Appendix[D\.1](https://arxiv.org/html/2608.24386#A4.SS1)\), rank\-4 elasticity tensor prediction \(Appendix[D\.2](https://arxiv.org/html/2608.24386#A4.SS2)\), runtime profiling \(Appendix[D\.3](https://arxiv.org/html/2608.24386#A4.SS3)\), and sensitivity to the LE\-ESO weightα\\alpha\(Appendix[D\.4](https://arxiv.org/html/2608.24386#A4.SS4)\)\. Dataset statistics are detailed in Appendix[C\.1](https://arxiv.org/html/2608.24386#A3.SS1)\.

We compare our equivariant full\-covariance model against two primary baselines: \(i\) a deterministic model trained with MSE, and \(ii\) a diagonal UQ model assuming independent components\. For ablation analysis, we evaluate non\-equivariant and non\-SPD variants \(Section[4\.5](https://arxiv.org/html/2608.24386#S4.SS5)\)\. All estimators utilize an E\(3\)\-equivariant backbone withℓm​a​x=4\\ell\_\{max\}=4to support rank\-4 covariance output; detailed hyperparameters and training protocols are provided in Appendix[B](https://arxiv.org/html/2608.24386#A2)\.

Our experiments verify three central hypotheses: \(1\) geometric validation—exact equivariance preservation on complex 3D shapes; \(2\) accuracy—maintained point\-prediction performance while modeling full covariance structure; and \(3\) calibration—covariance matrices that reflect true error distributions while respecting underlying geometry\.

### 4\.1Controlled Geometric Validation: Inertia Tensor Prediction

The ModelNet40 experiments are intended as controlled geometric validation rather than as replacements for closed\-form tensor estimators\. For inertia tensors, analytic formulas exist; our goal is to isolate whether the learned covariance remains equivariant, SPD, and geometrically meaningful under controlled point\-cloud perturbations\. The inertia tensorℐ\\mathcal\{I\}transforms equivariantly under rotation \(ℐ′=R​ℐ​R⊤\\mathcal\{I\}^\{\\prime\}=R\\mathcal\{I\}R^\{\\top\}\), making it a clean target for verifying geometric consistency\. We use the official split of 12,311 CAD models and introduce aleatoric uncertainty via Gaussian jitter applied to point positions \(detailed preprocessing in Appendix[C\.1](https://arxiv.org/html/2608.24386#A3.SS1)\)\. For this validation task, the mean and covariance heads are trained jointly without gradient detachment, as the synthetic noise is well\-behaved\. We provide a second controlled rank\-2 tensor validation on ModelNet40 shape covariance in Appendix[D\.1](https://arxiv.org/html/2608.24386#A4.SS1), showing that the construction is not tied to the inertia target\.

We quantify equivariance error using the relative Frobenius normEequiv=‖Σ⁡\(R⋅X\)−ρc​\(R\)​Σ​\(X\)​ρc​\(R\)⊤‖F/‖Σ⁡\(X\)‖FE\_\{\\text\{equiv\}\}=\\\|\\Sigma\(R\\\!\\cdot\\\!X\)\-\\rho\_\{c\}\(R\)\\Sigma\(X\)\\rho\_\{c\}\(R\)^\{\\top\}\\\|\_\{F\}/\\\|\\Sigma\(X\)\\\|\_\{F\}\. Prediction accuracy and uncertainty scores are reported in Table[1](https://arxiv.org/html/2608.24386#S4.T1); equivariance and SPD\-validity results are analyzed separately in Table[4](https://arxiv.org/html/2608.24386#S4.T4), where the equivariance errors remain at the level of10−710^\{\-7\}, confirming near\-machine\-precision symmetry preservation\. Our full\-covariance model reduces MAE by 15% relative to the diagonal baseline while maintaining perfect SPD properties \(\>\>99\.9% validity, median condition number 6\.8\)\. Detailed SPD analysis is provided in Appendix[E\.5](https://arxiv.org/html/2608.24386#A5.SS5)\.

Table 1:ModelNet40 inertia tensor prediction\. Our equivariant full\-covariance framework achieves strong performance with exact geometric consistency and physical constraints\.Figure[3](https://arxiv.org/html/2608.24386#S4.F3)provides visual validation that our uncertainty estimates are physically meaningful: uncertainty ellipsoids align with principal shape axes \(demonstrating E\(3\)\-equivariance\), expand in regions with sparse point density \(capturing sampling ambiguity\), and preserve tensorial correlations across components\.

### 4\.2Application: Dielectric Tensor Prediction

![Refer to caption](https://arxiv.org/html/2608.24386v1/paper_figure_accuracy.png)Figure 4:Performance on Materials Project\.\(a\) Parity plot of diagonal components \(R2=0\.659R^\{2\}=0\.659\)\. \(b\) Multivariate Laplace reliability diagram \(MACE=0\.0489=0\.0489\)\.Table 2:Performance on Materials Project dielectric dataset\. GoeCTP\([Hua et al\., 2026](https://arxiv.org/html/2608.24386#bib.bib29)\)employs a scalable equivariant architecture\. DTNet\([Mao et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib28)\)uses universal potential embeddings\. MACE\-Ens\. is a 5\-model deep ensemble\. Ours \(Calib\.\) applies temperature scaling \(T≈0\.05T\\approx 0\.05\)\. ES generalizes CRPS to multivariate settings\.Figure 5:Uncertainty diagnostics and utility analysis\.\(a\) Sharpness distribution of 95% confidence volumes, demonstrating the model’s ability to assign heteroscedastic uncertainties\. \(b\) Risk\-coverage analysis comparing ranking by directional uncertainty \(λmax​\(Σ\)\\lambda\_\{\\max\}\(\\Sigma\)\) and total uncertainty \(Tr⁡\(Σ\)\\mathrm\{Tr\}\(\\Sigma\)\)\. Directional uncertainty provides a modest but consistent advantage in identifying high\-error samples\.#### Dataset and Preprocessing\.

We utilize the Materials Project dielectric tensor dataset\([Barroso\-Luque et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib22);[Jain et al\., 2013](https://arxiv.org/html/2608.24386#bib.bib25)\)with static dielectric tensors computed via DFPT\. To ensure data quality and consistency, we apply systematic filtering criteria \(detailed in Appendix[C\.1](https://arxiv.org/html/2608.24386#A3.SS1)\), including structure size constraints \(3≤atoms≤303\\leq\\text\{atoms\}\\leq 30\), SPD positive\-definiteness verification, and value range constraints\. We apply Matrix Log\-Normalization to reduce the dynamic range of dielectric tensors before converting them to Kelvin–Mandel coordinates, ensuring that the six independent components are scaled consistently across diverse crystal structures\. We process atomic structures into graphs with5\.05\.0Å cutoff\. Architecture and training details are provided in Appendix[B](https://arxiv.org/html/2608.24386#A2)\.

### 4\.3Prediction Accuracy

Our goal is not state\-of\-the\-art point prediction alone, but competitive accuracy with full\-covariance, symmetry\-preserving uncertainty\. Table[2](https://arxiv.org/html/2608.24386#S4.T2)shows our method achieves competitive MAE among UQ models \(1\.55, vs\. 1\.96 for the MACE deep ensemble and 2\.25 for diagonal UQ\), and remains close to deterministic point predictors such as DTNet\([Mao et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib28)\)\(1\.91\) and GoeCTP\([Hua et al\., 2026](https://arxiv.org/html/2608.24386#bib.bib29)\)\(1\.41\), which do not provide calibrated equivariant covariance estimates\. The primary contribution is the principled, backbone\-agnostic UQ mechanism, not a new point estimator\. Modeling the full covariance manifold does not hinder mean estimation; the gains are most visible on off\-diagonal components, which capture anisotropic directional dependencies\.

Crucially, the parity plot in Figure[4](https://arxiv.org/html/2608.24386#S4.F4)a demonstrates that the highR2R^\{2\}in Kelvin\-Mandel log\-space suggests that the E\(3\)\-equivariant backbone effectively captures the underlying physics of dielectric properties while the UQ branch provides necessary aleatoric regularization\. The predictive covariance is SPD by construction due to the matrix exponential\. Separately, we also check whether the predicted mean dielectric tensors satisfy the expected positive\-definiteness of the physical tensor; all test\-set mean predictions satisfy this constraint in our run\.

#### Extensibility to Higher\-Order Tensors\.

To test extensibility beyond rank\-2, Appendix[D\.2](https://arxiv.org/html/2608.24386#A4.SS2)reports a real\-data rank\-4 elasticity experiment, where the mean target is a rank\-4 elasticity tensor with 21 independent components under standard minor/major symmetries\. This is supporting evidence rather than a comprehensive rank\-4 benchmark; the model achieves competitive MAE, improves empirical coverage from∼\\sim35% for a naive UQ baseline to∼\\sim52%, and preserves numerical equivariance and SPD validity\.

### 4\.4Uncertainty Calibration

A primary contribution of our work is the calibration of tensor\-valued uncertainty\. Our model is trained using the Multivariate Laplace negative log\-likelihood via the LE\-ESO objective \(Eq\.[12](https://arxiv.org/html/2608.24386#S3.E12)\), which naturally accounts for the heavy\-tailed error distributions common in materials data\. Figure[4](https://arxiv.org/html/2608.24386#S4.F4)b shows the reliability diagram evaluated against this same Multivariate Laplace distribution, confirming that our training objective and calibration assessment are properly matched\.

The model achieves a MACE of 0\.0489, indicating good agreement between predicted and empirical confidence levels under the multivariate Laplace evaluation protocol\. Compared with the Gaussian\-NLL objective in Table[3](https://arxiv.org/html/2608.24386#S4.T3), LE\-ESO gives lower calibration error and better accuracy on this dataset\. While the curve remains slightly below the diagonal in the high\-confidence regime, this indicates that the model isconservative\(under\-confident\), which is preferable for high\-stakes materials screening as it avoids over\-optimistic predictions\.

#### Utility for Risk\-Informed Decision Making\.

To evaluate the utility of our equivariant covariance tensors for risk\-informed decision making, we perform a risk\-coverage analysis comparing two ranking metrics: the total uncertainty \(Trace​\(Σ\)\\text\{Trace\}\(\\Sigma\)\) and the directional uncertainty \(λmax​\(Σ\)\\lambda\_\{\\max\}\(\\Sigma\)\)\. As illustrated in Figure[5](https://arxiv.org/html/2608.24386#S4.F5)b, the maximum eigenvalueλmax\\lambda\_\{\\max\}—representing the variance along the most uncertain principal axis—serves as a more targeted indicator of directional prediction risk, capturing failure modes that are partially obscured by scalar total uncertainty\.

At 90% coverage, ranking byλmax​\(Σ\)\\lambda\_\{\\max\}\(\\Sigma\)improves retained\-set MAE by 3\.1% relative to the full test set, compared with a 0\.8% improvement when ranking byTrace​\(Σ\)\\text\{Trace\}\(\\Sigma\)\. When compared directly against Trace\-based ranking under the same retained\-set protocol, the advantage ofλmax\\lambda\_\{\\max\}is smaller but consistent, and persists at lower coverage levels where Trace\-based ranking can fall below the full\-dataset baseline\. Appendix[D\.5](https://arxiv.org/html/2608.24386#A4.SS5)provides the full retained\-set comparison, including the diagonal\-UQ baseline\. These results indicate that for anisotropic physical properties like dielectric tensors, capturing the directional components of uncertainty is informative for identifying potential failure modes beyond what isotropic or diagonal approximations expose\.

#### Ablation of Training Objectives\.

This ablation isolates the effect of the scoring objective while keeping the equivariant backbone and matrix\-exponential covariance head fixed\. We compare LE\-ESO against the standard Gaussian NLL and Multivariate Energy Score \(ES\) training\. As summarized in Table[3](https://arxiv.org/html/2608.24386#S4.T3), the Laplace\-based LE\-ESO achieves the lowest MAE \(1\.55\) and calibration error \(MACE 0\.049, ES 0\.66\)\. The performance gap relative to Gaussian NLL \(MAE 1\.78\) is consistent with the heavy\-tailed nature of materials residuals, where the linear Mahalanobis penalty of the Laplace formulation is less sensitive to extreme deviations than the quadratic Gaussian penalty\. Energy Score training is competitive on calibration but exhibits higher MAE \(1\.64\)\. Overall, the results suggest that the Laplace\-based LE\-ESO objective provides a better accuracy\-calibration trade\-off than Gaussian NLL or Energy Score training in this dataset\.

Table 3:Ablation of training objectives on Materials Project\.All variants use the same E\(3\)\-equivariant backbone and matrix\-exponential head\. LE\-ESO \(Ours\) demonstrates superior robustness to heavy\-tailed noise compared to Gaussian NLL\.Table 4:Equivariance and SPD\-validity analysis\. Baseline B′isolates the Cholesky covariance failure mode \(equivariantμ\\mu, non\-equivariantΣ\\Sigma\)\. Only the matrix\-exponential head achieves both covariance equivariance and strict SPD validity\.
#### Chemical Out\-of\-Distribution Analysis\.

Beyond internal calibration, a key utility of symmetry\-preserving UQ is identifying Out\-of\-Distribution \(OOD\) samples during materials screening\. We perform aChemical Substitution Analysisby replacing common atoms in the test set with unseen elements \(e\.g\., Actinides U, Pu; Rare Earths Gd, Sm\)\. As shown in Figure[6](https://arxiv.org/html/2608.24386#S4.F6), the predicted directional uncertaintyλmax\\lambda\_\{\\max\}rises monotonically from 2\.71 to 5\.23 \(\+93\.2%\) as the substitution ratio reaches 100%, suggesting that the covariance head captures patterns correlated with chemical distribution shift; we do not claim to disentangle aleatoric and epistemic uncertainty in this analysis\.

![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/ood_chemical_test.png)Figure 6:Chemical OOD sensitivity analysis\.The predicted uncertainty increases as the crystal lattice is populated with unseen chemical species, suggesting that the model provides a useful distribution\-shift risk signal\.

### 4\.5Equivariance and Validity Verification

Table[4](https://arxiv.org/html/2608.24386#S4.T4)isolates equivariance and SPD validity across four baselines: A \(non\-equivariant GNN \+ Cholesky\), B \(E3NN \+ coordinate\-wise heads\), B′\(our equivariant mean head \+ Cholesky covariance\), and C \(direct equivariant regression ofA⁡\(X\)A\(X\)without the matrix exponential\)\. Implementation details are in Appendix[C\.2](https://arxiv.org/html/2608.24386#A3.SS2)\.

Baseline B′reaches near\-machine\-precisionℰμ≈1\.4×10−6\\mathcal\{E\}\_\{\\mu\}\\\!\\approx\\\!1\.4\\times 10^\{\-6\}butℰΣ≈0\.43\\mathcal\{E\}\_\{\\Sigma\}\\\!\\approx\\\!0\.43, confirming that the failure is specific to Cholesky: its lower\-triangular structure is not preserved under orthogonal conjugation byρc​\(R\)\\rho\_\{c\}\(R\)\. Only the matrix\-exponential head achievesO⁡\(10−7\)O\(10^\{\-7\}\)equivariance error with strict SPD validity\.

#### Computational Overhead\.

Appendix[D\.3](https://arxiv.org/html/2608.24386#A4.SS3)reports runtime profiling\. The full\-covariance model is more expensive than the deterministic baseline because of the covariance branch and its backpropagation, but it provides the full anisotropic covariance in a single forward pass without ensembling\.

## 5Discussion

Our primary validated setting is E\(3\)\-equivariant UQ for symmetric rank\-2 tensors\. The additional shape\-covariance experiment in Appendix[D\.1](https://arxiv.org/html/2608.24386#A4.SS1)shows that the construction is not specific to the inertia target, while the rank\-4 elasticity experiment in Appendix[D\.2](https://arxiv.org/html/2608.24386#A4.SS2)provides supporting evidence for higher\-order tensor targets\. However, the higher\-order empirical study is not yet exhaustive, and broader validation across additional equivariant backbones, tensor orders, and symmetry groups remains future work\. The construction separates a group\-agnostic SPD/UQ core—namely, the matrix\-exponential mapping and the Log\-Euclidean scoring objective—from a representation\-specific equivariant parameterization of the symmetric operatorA⁡\(X\)A\(X\)\. Extending the framework to other groups therefore depends on the availability of suitable representation\-theoretic bases and implementation tools\.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China \(Grant No\. 62573132\) and the Fundamental Research Funds for the Central Universities \(Grant No\. 2025SMECP012\)\.

## Impact Statement

This paper introduces a framework for equivariant uncertainty quantification in tensor\-valued geometric learning\. Its primary potential impact lies in scientific machine learning, particularly in materials discovery tasks where tensor\-valued properties such as dielectric or elastic responses are important\. By providing symmetry\-preserving uncertainty estimates, the method may support more reliable and risk\-aware computational screening, helping researchers prioritize candidates for further validation\. At the same time, predictions from such models should not be treated as substitutes for experimental or domain\-expert verification\. Responsible use of the framework requires human\-in\-the\-loop assessment, careful calibration checks, and validation on the target scientific domain\.

## References

- Arsignyet al\.\(2006\)V\. Arsigny, P\. Fillard, X\. Pennec, and N\. AyacheLog\-euclidean metrics for fast and simple calculus on diffusion tensors\.Magnetic Resonance in Medicine56\(2\),pp\. 411–421\.Cited by:[§3\.5](https://arxiv.org/html/2608.24386#S3.SS5.SSS0.Px1.p1.1)\.
- Barroso\-Luqueet al\.\(2024\)L\. Barroso\-Luque, M\. Shuaibi, X\. Fu, B\. M\. Wood, M\. Dzamba, M\. Gao, A\. Rizvi, C\. L\. Zitnick, and Z\. W\. UlissiOpen materials 2024 \(omat24\) inorganic materials dataset and models\.arXiv preprint arXiv:2410\.12771\.Cited by:[§C\.1](https://arxiv.org/html/2608.24386#A3.SS1.SSS0.Px2.p1.1),[§4\.2](https://arxiv.org/html/2608.24386#S4.SS2.SSS0.Px1.p1.1),[§4](https://arxiv.org/html/2608.24386#S4.p1.1)\.
- Batatiaet al\.\(2022\)I\. Batatia, D\. P\. Kovacs, G\. Simm, C\. Ortner, and G\. CsányiMACE: higher order equivariant message passing neural networks for fast and accurate force fields\.Advances in neural information processing systems35,pp\. 11423–11436\.Cited by:[§1](https://arxiv.org/html/2608.24386#S1.p1.1),[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Bermanet al\.\(2026\)E\. Berman, J\. Ginesin, M\. Pacini, and R\. WaltersOn Uncertainty Calibration for Equivariant Functions\.Transactions on Machine Learning Research\.External Links:[Link](https://mlanthology.org/tmlr/2026/berman2026tmlr-uncertainty/)Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Bevandaet al\.\(2025\)P\. Bevanda, M\. Beier, A\. Capone, S\. G\. Sosnowski, S\. Hirche, and A\. LedererKoopman\-equivariant gaussian processes\.InProceedings of The 28th International Conference on Artificial Intelligence and Statistics,Y\. Li, S\. Mandt, S\. Agrawal, and E\. Khan \(Eds\.\),Proceedings of Machine Learning Research, Vol\.258,pp\. 3151–3159\.External Links:[Link](https://proceedings.mlr.press/v258/bevanda25a.html)Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Doanet al\.\(2025\)B\. G\. Doan, A\. Shamsi, X\. Guo, A\. Mohammadi, H\. Alinejad\-Rokny, D\. Sejdinovic, D\. Teney, D\. C\. Ranasinghe, and E\. AbbasnejadBayesian low\-rank learning \(Bella\): a practical approach to Bayesian neural networks\.Proceedings of the AAAI Conference on Artificial Intelligence39\(15\),pp\. 16298–16307\.External Links:[Document](https://dx.doi.org/10.1609/aaai.v39i15.33790),[Link](https://ojs.aaai.org/index.php/AAAI/article/view/33790)Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Duet al\.\(2024\)H\. Du, J\. Wang, J\. Hui, L\. Zhang, and H\. WangDenseGNN: universal and scalable deeper graph neural networks for high\-performance property prediction in crystals and molecules\.npj Computational Materials10\(1\),pp\. 292\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Equeret al\.\(2023\)L\. Equer, T\. K\. Rusch, and S\. MishraMulti\-scale message passing neural pde solvers\.arXiv preprint arXiv:2302\.03580\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Fakouret al\.\(2024\)F\. Fakour, A\. Mosleh, and R\. RamezaniA structured review of literature on uncertainty in machine learning & deep learning\.arXiv preprint arXiv:2406\.00332\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Funget al\.\(2021\)V\. Fung, J\. Zhang, E\. Juarez, and B\. G\. SumpterBenchmarking graph neural networks for materials chemistry\.npj Computational Materials7\(1\),pp\. 84\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Geiger and Smidt \(2022\)M\. Geiger and T\. SmidtE3nn: euclidean neural networks\.arXiv preprint arXiv:2207\.09453\.Cited by:[§B\.2](https://arxiv.org/html/2608.24386#A2.SS2.p1.1),[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Gerken and Kessel \(2024\)J\. E\. Gerken and P\. KesselEmergent equivariance in deep ensembles\.arXiv preprint arXiv:2403\.03103\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Gneiting and Raftery \(2007\)T\. Gneiting and A\. E\. RafteryStrictly proper scoring rules, prediction, and estimation\.Journal of the American Statistical Association102\(477\),pp\. 359–378\.Cited by:[§3\.5](https://arxiv.org/html/2608.24386#S3.SS5.SSS0.Px2.p2.1)\.
- Gruber and Buettner \(2022\)S\. Gruber and F\. BuettnerBetter uncertainty calibration via proper scores for classification and beyond\.Advances in Neural Information Processing Systems35,pp\. 8618–8632\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Heilmanet al\.\(2024\)A\. Heilman, C\. Schlesinger, and Q\. YanEquivariant graph neural networks for prediction of tensor material properties of crystals\.arXiv preprint arXiv:2406\.03563\.Cited by:[§1](https://arxiv.org/html/2608.24386#S1.p1.1),[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1),[Table 2](https://arxiv.org/html/2608.24386#S4.T2.5.2.1.1)\.
- Huaet al\.\(2026\)H\. Hua, J\. Yang, W\. Lin, and P\. ZhouRevisiting the Canonicalization for Fast and Accurate Crystal Tensor Property Prediction\.Proceedings of the AAAI Conference on Artificial Intelligence40\(1\),pp\. 417–425\.External Links:[Document](https://dx.doi.org/10.1609/aaai.v40i1.37004),[Link](https://ojs.aaai.org/index.php/AAAI/article/view/37004)Cited by:[§E\.6](https://arxiv.org/html/2608.24386#A5.SS6.SSS0.Px1.p1.1),[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1),[§4\.3](https://arxiv.org/html/2608.24386#S4.SS3.p1.1),[Table 2](https://arxiv.org/html/2608.24386#S4.T2),[Table 2](https://arxiv.org/html/2608.24386#S4.T2.4),[Table 2](https://arxiv.org/html/2608.24386#S4.T2.5.4.1.1)\.
- Jainet al\.\(2013\)A\. Jain, S\. P\. Ong, G\. Hautier, W\. Chen, W\. D\. Richards, S\. Dacek, S\. Cholia, D\. Gunter, D\. Skinner, G\. Ceder,et al\.Commentary: the materials project: a materials genome approach to accelerating materials innovation\.APL materials1\(1\)\.Cited by:[§C\.1](https://arxiv.org/html/2608.24386#A3.SS1.SSS0.Px2.p1.1),[§D\.2](https://arxiv.org/html/2608.24386#A4.SS2.p1.1),[§4\.2](https://arxiv.org/html/2608.24386#S4.SS2.SSS0.Px1.p1.1),[§4](https://arxiv.org/html/2608.24386#S4.p1.1)\.
- Jekelet al\.\(2022\)C\. F\. Jekel, K\. E\. Swartz, D\. A\. White, D\. A\. Tortorelli, and S\. E\. WattsNeural network layers for prediction of positive definite elastic stiffness tensors\.arXiv preprint arXiv:2203\.13938\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px2.p1.1)\.
- Kuleshovet al\.\(2018\)V\. Kuleshov, N\. Fenner, and S\. ErmonAccurate uncertainties for deep learning using calibrated regression\.InInternational Conference on Machine Learning,pp\. 2796–2804\.Cited by:[§B\.3](https://arxiv.org/html/2608.24386#A2.SS3.SSS0.Px6.p1.1)\.
- Maoet al\.\(2024\)Z\. Mao, W\. Li, and J\. TanDielectric tensor prediction for inorganic materials using latent information from preferred potential\.npj Computational Materials10\(1\),pp\. 265\.Cited by:[§4\.3](https://arxiv.org/html/2608.24386#S4.SS3.p1.1),[Table 2](https://arxiv.org/html/2608.24386#S4.T2),[Table 2](https://arxiv.org/html/2608.24386#S4.T2.4),[Table 2](https://arxiv.org/html/2608.24386#S4.T2.5.3.1.1)\.
- Newmanet al\.\(2024\)E\. Newman, L\. Horesh, H\. Avron, and M\. E\. KilmerStable tensor neural networks for efficient deep learning\.Frontiers in Big Data7,pp\. 1363978\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Olivieret al\.\(2021\)A\. Olivier, M\. D\. Shields, and L\. Graham\-BradyBayesian neural networks for uncertainty quantification in data\-driven materials modeling\.Computer methods in applied mechanics and engineering386,pp\. 114079\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Pakornchoteet al\.\(2023\)T\. Pakornchote, A\. Ektarawong, and T\. ChotibutStrainTensorNet: predicting crystal structure elastic properties using se\(3\)\-equivariant graph neural networks\.Physical Review Research5\(4\),pp\. 043198\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Pouliquenet al\.\(2025\)C\. Pouliquen, M\. Massias, and T\. VayerSchur’s positive\-definite network: deep learning in the spd cone with structure\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 71401–71416\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px2.p1.1)\.
- Reiseret al\.\(2022\)P\. Reiser, M\. Neubert, A\. Eberhard, L\. Torresi, C\. Zhou, C\. Shao, H\. Metni, C\. van Hoesel, H\. Schopmans, T\. Sommer,et al\.Graph neural networks for materials science and chemistry\.Communications Materials3\(1\),pp\. 93\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1)\.
- Rensmeyeret al\.\(2024\)T\. Rensmeyer, B\. Craig, D\. Kramer, and O\. NiggemannHigh accuracy uncertainty\-aware interatomic force modeling with equivariant bayesian neural networks\.Digital Discovery3\(11\),pp\. 2356–2366\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Rudneret al\.\(2022\)T\. G\. Rudner, Z\. Chen, Y\. W\. Teh, and Y\. GalTractable function\-space variational inference in bayesian neural networks\.Advances in Neural Information Processing Systems35,pp\. 22686–22698\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Shaoet al\.\(2025\)S\. Shao, Y\. Li, Z\. Lin, and Q\. CuiHigh\-rank irreducible cartesian tensor decomposition and bases of equivariant spaces\.Journal of Machine Learning Research26\(175\),pp\. 1–53\.External Links:[Link](http://jmlr.org/papers/v26/25-0134.html)Cited by:[§E\.6](https://arxiv.org/html/2608.24386#A5.SS6.p1.1),[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px1.p1.1),[§3\.4](https://arxiv.org/html/2608.24386#S3.SS4.p3.4)\.
- Sheinkman and Wade \(2025\)A\. Sheinkman and S\. WadeThe architecture and evaluation of bayesian neural networks\.arXiv e\-prints,pp\. arXiv–2503\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Steinertet al\.\(2025\)T\. Steinert, D\. Ginsbourger, A\. Lykke\-Møller, O\. Christiansen, and H\. MossIntegration\-free kernels for equivariant gaussian process modelling\.InForty\-second International Conference on Machine Learning,Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.
- Wuet al\.\(2015\)Z\. Wu, S\. Song, A\. Khosla, F\. Yu, L\. Zhang, X\. Tang, and J\. Xiao3d shapenets: a deep representation for volumetric shapes\.InProceedings of the IEEE conference on computer vision and pattern recognition,pp\. 1912–1920\.Cited by:[§D\.1](https://arxiv.org/html/2608.24386#A4.SS1.p1.1),[§4](https://arxiv.org/html/2608.24386#S4.p1.1)\.
- Zhaoet al\.\(2023\)W\. Zhao, F\. Lopez, J\. M\. Riestenberg, M\. Strube, D\. Taha, and S\. TrettelModeling graphs beyond hyperbolic: graph neural networks in symmetric positive definite matrices\.InJoint European Conference on Machine Learning and Knowledge Discovery in Databases,pp\. 122–139\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px2.p1.1)\.
- Zhouet al\.\(2024\)X\. Zhou, Z\. Liu, and H\. XiaoBi\-eqno: generalized approximate bayesian inference with an equivariant neural operator framework\.arXiv preprint arXiv:2410\.16420\.Cited by:[§2](https://arxiv.org/html/2608.24386#S2.SS0.SSS0.Px3.p1.1)\.

## Appendix ATheoretical Proofs and Derivations

In this section, we provide formal statements and proofs establishing the mathematical validity of our equivariant uncertainty formulation\. We establish the theoretical foundations for \(i\) the representation\-theoretic decomposition of covariance tensors, \(ii\) the matrix exponential construction ensuring both positive\-definiteness and equivariance, and \(iii\) the numerical stability of our loss formulation\.

### A\.1Rotation Matrices in Kelvin\-Mandel Space

For a rotation matrixR∈S​O​\(3\)R\\in SO\(3\), the corresponding6×66\\times 6transformation matrixρc​\(R\)\\rho\_\{c\}\(R\)in Kelvin\-Mandel space can be derived from the Kronecker product structure\. The vectorization operationvec​\(C\)\\text\{vec\}\(C\)maps a symmetric tensorCCto a 9\-dimensional vector, and under rotation:

vec​\(C′\)=\(R⊗R\)​vec​\(C\),\\text\{vec\}\(C^\{\\prime\}\)=\(R\\otimes R\)\\text\{vec\}\(C\),\(13\)where⊗\\otimesdenotes the Kronecker product\.

The matrixρc​\(R\)\\rho\_\{c\}\(R\)is obtained by projectingR⊗RR\\otimes Ronto the 6\-dimensional symmetric subspace and applying the Kelvin\-Mandel scaling matrix𝐏\\mathbf\{P\}:

ρc​\(R\)=𝐏⋅𝒮⋅\(R⊗R\)⋅𝒮T⋅𝐏−1,\\rho\_\{c\}\(R\)=\\mathbf\{P\}\\cdot\\mathcal\{S\}\\cdot\(R\\otimes R\)\\cdot\\mathcal\{S\}^\{T\}\\cdot\\mathbf\{P\}^\{\-1\},\(14\)where𝒮\\mathcal\{S\}is the6×96\\times 9selection matrix\.

For practical computation,ρc​\(R\)\\rho\_\{c\}\(R\)has the explicit block structure:

ρc​\(R\)=\[R112R122R1322​R12​R132​R11​R132​R11​R12R212R222R2322​R22​R232​R21​R232​R21​R22R312R322R3322​R32​R332​R31​R332​R31​R322​R21​R312​R22​R322​R23​R33R22​R33\+R23​R32R21​R33\+R23​R31R21​R32\+R22​R312​R11​R312​R12​R322​R13​R33R12​R33\+R13​R32R11​R33\+R13​R31R11​R32\+R12​R312​R11​R212​R12​R222​R13​R23R12​R23\+R13​R22R11​R23\+R13​R21R11​R22\+R12​R21\]\.\\footnotesize\\rho\_\{c\}\(R\)=\\begin\{bmatrix\}R\_\{11\}^\{2\}&R\_\{12\}^\{2\}&R\_\{13\}^\{2\}&\\sqrt\{2\}R\_\{12\}R\_\{13\}&\\sqrt\{2\}R\_\{11\}R\_\{13\}&\\sqrt\{2\}R\_\{11\}R\_\{12\}\\\\ R\_\{21\}^\{2\}&R\_\{22\}^\{2\}&R\_\{23\}^\{2\}&\\sqrt\{2\}R\_\{22\}R\_\{23\}&\\sqrt\{2\}R\_\{21\}R\_\{23\}&\\sqrt\{2\}R\_\{21\}R\_\{22\}\\\\ R\_\{31\}^\{2\}&R\_\{32\}^\{2\}&R\_\{33\}^\{2\}&\\sqrt\{2\}R\_\{32\}R\_\{33\}&\\sqrt\{2\}R\_\{31\}R\_\{33\}&\\sqrt\{2\}R\_\{31\}R\_\{32\}\\\\ \\sqrt\{2\}R\_\{21\}R\_\{31\}&\\sqrt\{2\}R\_\{22\}R\_\{32\}&\\sqrt\{2\}R\_\{23\}R\_\{33\}&R\_\{22\}R\_\{33\}\+R\_\{23\}R\_\{32\}&R\_\{21\}R\_\{33\}\+R\_\{23\}R\_\{31\}&R\_\{21\}R\_\{32\}\+R\_\{22\}R\_\{31\}\\\\ \\sqrt\{2\}R\_\{11\}R\_\{31\}&\\sqrt\{2\}R\_\{12\}R\_\{32\}&\\sqrt\{2\}R\_\{13\}R\_\{33\}&R\_\{12\}R\_\{33\}\+R\_\{13\}R\_\{32\}&R\_\{11\}R\_\{33\}\+R\_\{13\}R\_\{31\}&R\_\{11\}R\_\{32\}\+R\_\{12\}R\_\{31\}\\\\ \\sqrt\{2\}R\_\{11\}R\_\{21\}&\\sqrt\{2\}R\_\{12\}R\_\{22\}&\\sqrt\{2\}R\_\{13\}R\_\{23\}&R\_\{12\}R\_\{23\}\+R\_\{13\}R\_\{22\}&R\_\{11\}R\_\{23\}\+R\_\{13\}R\_\{21\}&R\_\{11\}R\_\{22\}\+R\_\{12\}R\_\{21\}\\end\{bmatrix\}\.\(15\)This explicit form ensures thatρc​\(R\)\\rho\_\{c\}\(R\)maintains orthogonality in Kelvin\-Mandel space:ρc​\(R\)T​ρc​\(R\)=I6\\rho\_\{c\}\(R\)^\{T\}\\rho\_\{c\}\(R\)=I\_\{6\}\.

#### Note on Voigt vs\. Kelvin\-Mandel Notation\.

While standard Voigt notation mapsCi​jC\_\{ij\}to\[C11,C22,C33,C23,C13,C12\]T\[C\_\{11\},C\_\{22\},C\_\{33\},C\_\{23\},C\_\{13\},C\_\{12\}\]^\{T\}, it does not preserve the Frobenius norm\. Kelvin\-Mandel notation applies2\\sqrt\{2\}scaling to shear components, ensuring‖𝐜KM‖2=‖C‖F\\\|\\mathbf\{c\}\_\{\\text\{KM\}\}\\\|\_\{2\}=\\\|C\\\|\_\{F\}\. This isometric property is crucial for maintaining geometric consistency in uncertainty quantification\.

### A\.2Irreducible Representation Decomposition

###### Proposition A\.1\(Irreducible decomposition of the covariance representation\)\.

Letρc\\rho\_\{c\}denote the 6\-dimensional real representation ofSO⁡\(3\)\\mathrm\{SO\}\(3\)corresponding to symmetric rank\-2 tensors, i\.e\.

ρc≅l=0⊕l=2\.\\rho\_\{c\}\\cong l=0\\oplus l=2\.Then the symmetric tensor product representation ofρc\\rho\_\{c\}decomposes as

Sym2​\(ρc\)≅2×\(l=0\)⊕2×\(l=2\)⊕1×\(l=4\),\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\;\\cong\\;2\\times\(l=0\)\\;\\oplus\\;2\\times\(l=2\)\\;\\oplus\\;1\\times\(l=4\),which possesses2121independent degrees of freedom—equal to that of a symmetric6×66\\times 6covariance matrix\.

###### Proof\.

Sinceρc≅l=0⊕l=2\\rho\_\{c\}\\cong l=0\\oplus l=2, the symmetric square decomposes as

Sym2​\(ρc\)≅Sym2​\(l=0\)⊕\(l=0⊗l=2\)⊕Sym2​\(l=2\)\.\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\cong\\mathrm\{Sym\}^\{2\}\(l=0\)\\oplus\(l=0\\otimes l=2\)\\oplus\\mathrm\{Sym\}^\{2\}\(l=2\)\.We haveSym2​\(l=0\)=l=0\\mathrm\{Sym\}^\{2\}\(l=0\)=l=0,l=0⊗l=2=l=2l=0\\otimes l=2=l=2, andSym2​\(l=2\)=l=0⊕l=2⊕l=4\\mathrm\{Sym\}^\{2\}\(l=2\)=l=0\\oplus l=2\\oplus l=4\. Combining these yields

Sym2​\(ρc\)≅2×\(l=0\)⊕2×\(l=2\)⊕1×\(l=4\),\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)\\cong 2\\times\(l=0\)\\oplus 2\\times\(l=2\)\\oplus 1\\times\(l=4\),which has dimension2⋅1\+2⋅5\+1⋅9=212\\cdot 1\+2\\cdot 5\+1\\cdot 9=21\. ∎

### A\.3Matrix Exponential Properties

###### Proposition A\.2\(Positive\-definiteness and equivariance of the exponential map\)\.

LetA⁡\(X\)∈ℝsym6×6A\(X\)\\in\\mathbb\{R\}^\{6\\times 6\}\_\{\\mathrm\{sym\}\}satisfy the equivariance condition

A⁡\(R⋅X\)=ρc​\(R\)​A​\(X\)​ρc​\(R\)⊤∀R∈O⁡\(3\)\.A\(R\\\!\\cdot\\\!X\)=\\rho\_\{c\}\(R\)\\,A\(X\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\\quad\\forall R\\in O\(3\)\.Then the matrix exponential

Σ⁡\(X\)=exp⁡\(A⁡\(X\)\)\\Sigma\(X\)=\\exp\(A\(X\)\)is \(i\) symmetric positive\-definite for allXX, and \(ii\) equivariant under the same group action:

Σ⁡\(R⋅X\)=ρc​\(R\)​Σ​\(X\)​ρc​\(R\)⊤\.\\Sigma\(R\\\!\\cdot\\\!X\)=\\rho\_\{c\}\(R\)\\,\\Sigma\(X\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\.

###### Proof\.

For any real symmetricAA, there exists an orthogonalQQand real diagonalΛ\\Lambdasuch thatA=Q​Λ​Q⊤A=Q\\Lambda Q^\{\\top\}\. Then

exp⁡\(A\)=Q​exp⁡\(Λ\)​Q⊤,\\exp\(A\)=Q\\,\\exp\(\\Lambda\)\\,Q^\{\\top\},whereexp⁡\(Λ\)\\exp\(\\Lambda\)has strictly positive diagonal entriesexp⁡\(λi\)\>0\\exp\(\\lambda\_\{i\}\)\>0\. Thusexp⁡\(A\)\\exp\(A\)is symmetric positive\-definite\.

For equivariance, note thatρc​\(R\)\\rho\_\{c\}\(R\)is orthogonal\. The matrix exponential satisfiesexp⁡\(S​A​S−1\)=S​exp⁡\(A\)​S−1\\exp\(SAS^\{\-1\}\)=S\\exp\(A\)S^\{\-1\}for any invertibleSS\. TakingS=ρc​\(R\)S=\\rho\_\{c\}\(R\)gives

exp⁡\(ρc​\(R\)​A​ρc​\(R\)⊤\)=ρc​\(R\)​exp⁡\(A\)​ρc​\(R\)⊤,\\exp\(\\rho\_\{c\}\(R\)A\\rho\_\{c\}\(R\)^\{\\top\}\)=\\rho\_\{c\}\(R\)\\exp\(A\)\\rho\_\{c\}\(R\)^\{\\top\},which proves equivariance\. ∎

###### Proposition A\.3\(Equivariance of spectral functions\)\.

Letf:ℝ→ℝf:\\mathbb\{R\}\\to\\mathbb\{R\}be a scalar function\. For a symmetric matrixAAwith eigenvalue decompositionA=Q​Λ​Q⊤A=Q\\Lambda Q^\{\\top\}, define the spectral functionF⁡\(A\)=Q​diag​\(f⁡\(λ1\),…,f⁡\(λn\)\)​Q⊤F\(A\)=Q\\,\\text\{diag\}\(f\(\\lambda\_\{1\}\),\\dots,f\(\\lambda\_\{n\}\)\)\\,Q^\{\\top\}\. Sinceρc​\(R\)\\rho\_\{c\}\(R\)is orthogonal for allR∈O⁡\(3\)R\\in O\(3\), it follows that

F⁡\(ρc​\(R\)​A​ρc​\(R\)⊤\)=ρc​\(R\)​F​\(A\)​ρc​\(R\)⊤\.F\(\\rho\_\{c\}\(R\)\\,A\\,\\rho\_\{c\}\(R\)^\{\\top\}\)=\\rho\_\{c\}\(R\)\\,F\(A\)\\,\\rho\_\{c\}\(R\)^\{\\top\}\.Thus, eigenvalue clamping and anisotropic jitter \(as defined in Eq\.[18](https://arxiv.org/html/2608.24386#A2.E18)\) preserveO⁡\(3\)O\(3\)equivariance\.

###### Proof\.

For any orthogonal matrixUU, the spectral function commutes with orthogonal similarity transformations:F⁡\(U​A​U⊤\)=U​F​\(A\)​U⊤F\(UAU^\{\\top\}\)=UF\(A\)U^\{\\top\}\. This follows from the fact thatU​A​U⊤=Q​Λ​Q⊤UAU^\{\\top\}=Q\\Lambda Q^\{\\top\}whereQ=U​Q0Q=UQ\_\{0\}for the original eigenvectorsQ0Q\_\{0\}ofAA\. Applying the definition ofFF:

F⁡\(U​A​U⊤\)=\(U​Q0\)​diag​\(f⁡\(λi\)\)​\(U​Q0\)⊤=U⁡\(Q0​diag​\(f⁡\(λi\)\)​Q0⊤\)​U⊤=U​F​\(A\)​U⊤\.F\(UAU^\{\\top\}\)=\(UQ\_\{0\}\)\\,\\text\{diag\}\(f\(\\lambda\_\{i\}\)\)\\,\(UQ\_\{0\}\)^\{\\top\}=U\(Q\_\{0\}\\,\\text\{diag\}\(f\(\\lambda\_\{i\}\)\)\\,Q\_\{0\}^\{\\top\}\)U^\{\\top\}=UF\(A\)U^\{\\top\}\.TakingU=ρc​\(R\)U=\\rho\_\{c\}\(R\)completes the proof\. ∎

### A\.4Numerically Stable Loss Function

We present two formulations: the standard Gaussian NLL \(for comparison\) and the Multivariate Laplace NLL used in our implementation\.

###### Proposition A\.4\(Gaussian NLL in Log\-Euclidean form\)\.

LetAAbe a symmetric matrix,Σ=exp⁡\(A\)\\Sigma=\\exp\(A\), andΔ​𝐜=𝐜true−μ\\Delta\\mathbf\{c\}=\\mathbf\{c\}\_\{\\text\{true\}\}\-\\mu\. The standard Gaussian negative log\-likelihood is

ℒGauss=12​log​detΣ\+12​Δ​𝐜⊤​Σ−1​Δ​𝐜\.\\mathcal\{L\}\_\{\\text\{Gauss\}\}=\\frac\{1\}\{2\}\\log\\det\\Sigma\+\\frac\{1\}\{2\}\\Delta\\mathbf\{c\}^\{\\top\}\\Sigma^\{\-1\}\\Delta\\mathbf\{c\}\.Then the following loss is algebraically equivalent and numerically stable:

ℒGauss=12​Tr⁡\(A\)\+12​Δ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜\.\\boxed\{\\mathcal\{L\}\_\{\\text\{Gauss\}\}=\\frac\{1\}\{2\}\\operatorname\{Tr\}\(A\)\+\\frac\{1\}\{2\}\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\.\}

###### Proof\.

UsingΣ=exp⁡\(A\)\\Sigma=\\exp\(A\)and the identitydet\(exp⁡\(A\)\)=exp⁡\(Tr⁡\(A\)\)\\det\(\\exp\(A\)\)=\\exp\(\\operatorname\{Tr\}\(A\)\), we obtainlogdetΣ=Tr\(A\)\\log\\det\\Sigma=\\operatorname\{Tr\}\(A\)\. Since\(exp⁡\(A\)\)−1=exp⁡\(−A\)\(\\exp\(A\)\)^\{\-1\}=\\exp\(\-A\), the Mahalanobis term becomesΔ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\. ∎

###### Proposition A\.5\(Multivariate Laplace NLL in Log\-Euclidean form\)\.

LetAAbe a symmetric matrix,Σ=exp⁡\(A\)\\Sigma=\\exp\(A\), andΔ​𝐜=𝐜true−μ\\Delta\\mathbf\{c\}=\\mathbf\{c\}\_\{\\text\{true\}\}\-\\mu\. The Multivariate Laplace negative log\-likelihood \(with unit scale\) is

ℒLaplace=logdetΣ\+Δ​𝐜⊤​Σ−1​Δ​𝐜\.\\mathcal\{L\}\_\{\\text\{Laplace\}\}=\\log\\det\\Sigma\+\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\Sigma^\{\-1\}\\Delta\\mathbf\{c\}\}\.The numerically stable form in the Lie algebra𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)is:

ℒLaplace=Tr⁡\(A\)\+DM,\\boxed\{\\mathcal\{L\}\_\{\\text\{Laplace\}\}=\\operatorname\{Tr\}\(A\)\+D\_\{M\},\}whereDM=Δ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜D\_\{M\}=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\}is the Mahalanobis distance in the Log\-Euclidean metric\.

###### Proof\.

The log\-determinant term follows identically:logdetΣ=Tr\(A\)\\log\\det\\Sigma=\\operatorname\{Tr\}\(A\)\. For the Mahalanobis distance term, note thatΣ−1=exp⁡\(−A\)\\Sigma^\{\-1\}=\\exp\(\-A\), so:

DM=Δ​𝐜⊤​Σ−1​Δ​𝐜=Δ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜\.D\_\{M\}=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\Sigma^\{\-1\}\\Delta\\mathbf\{c\}\}=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\}\.The key difference from the Gaussian case is thesquare root: the Laplace NLL is linear inDMD\_\{M\}rather than quadratic inDM2D\_\{M\}^\{2\}\. This provides robustness to outliers, as large residuals contribute linearly rather than quadratically to the loss\. ∎

#### Laplacian\-Huber Compound Robust Loss\.

To handle extreme outliers in materials data, we introduce aLaplacian\-Huber schemethat combines two complementary robustness mechanisms\. For the Mahalanobis distanceDM=Δ​𝐜⊤​Σ−1​Δ​𝐜D\_\{M\}=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\Sigma^\{\-1\}\\Delta\\mathbf\{c\}\}, our robust loss is:

D~M=\{DM,DM<τ\(linear region\)τ\+log⁡\(1\+DM−τ\),DM≥τ\(log\-tail region\)\\tilde\{D\}\_\{M\}=\\begin\{cases\}D\_\{M\},&D\_\{M\}<\\tau\\quad\\text\{\(linear region\)\}\\\\ \\tau\+\\log\(1\+D\_\{M\}\-\\tau\),&D\_\{M\}\\geq\\tau\\quad\\text\{\(log\-tail region\)\}\\end\{cases\}\(16\)This design has a clear statistical interpretation: \(1\) thelinear region\(DM<τD\_\{M\}<\\tau\) preserves the core Laplace distribution assumption, providing natural robustness through a linear rather than quadratic penalty on residuals; \(2\) thelog\-tail region\(DM≥τD\_\{M\}\\geq\\tau\) smoothly compresses the contribution of extreme residuals so that the residual penalty grows logarithmically rather than linearly, mitigating large updates from rare outliers in practice\. We setτ=5\.0\\tau=5\.0based on validation analysis—approximately 99% of well\-predicted samples haveDM<5D\_\{M\}<5, while extreme outliers beyond this threshold are smoothly bounded without affecting the majority of the data distribution\.

#### Gradient Behavior of the Laplace NLL\.

Using eigenvalue decompositionA=Q​Λ​Q⊤A=Q\\Lambda Q^\{\\top\}, define the whitened residual𝐳=Q⊤​Δ​𝐜\\mathbf\{z\}=Q^\{\\top\}\\Delta\\mathbf\{c\}\. The Mahalanobis distance becomes:

DM=∑i=16zi2​exp⁡\(−λi\)\.D\_\{M\}=\\sqrt\{\\sum\_\{i=1\}^\{6\}z\_\{i\}^\{2\}\\exp\(\-\\lambda\_\{i\}\)\}\.The gradient with respect to eigenvaluesΛ\\Lambdain the linear region \(DM<τD\_\{M\}<\\tau\) is:

∂ℒLaplace∂λk=1−zk2​exp⁡\(−λk\)2​DM,\\frac\{\\partial\\mathcal\{L\}\_\{\\text\{Laplace\}\}\}\{\\partial\\lambda\_\{k\}\}=1\-\\frac\{z\_\{k\}^\{2\}\\exp\(\-\\lambda\_\{k\}\)\}\{2D\_\{M\}\},\(17\)Compared to the Gaussian gradient12​\(1−zk2​exp⁡\(−λk\)\)\\frac\{1\}\{2\}\(1\-z\_\{k\}^\{2\}\\exp\(\-\\lambda\_\{k\}\)\), the Laplace gradient carries an additional1/\(2​DM\)1/\(2D\_\{M\}\)factor that reduces the growth rate of the data\-fit term, but it does not by itself yield a uniform bound: when a single directionkkdominatesDMD\_\{M\}\(soDM≈\|zk\|exp\(−λk/2\)D\_\{M\}\\approx\|z\_\{k\}\|\\exp\(\-\\lambda\_\{k\}/2\)\), the contribution behaves like\|zk\|exp\(−λk/2\)/2\|z\_\{k\}\|\\exp\(\-\\lambda\_\{k\}/2\)/2and is unbounded asλk→−∞\\lambda\_\{k\}\\to\-\\infty\. In the logarithmic tail region \(DM≥τD\_\{M\}\\geq\\tau\), the gradient becomes:

∂D~M∂λk=−zk2​exp⁡\(−λk\)2​\(1\+DM−τ\)​DM,\\frac\{\\partial\\tilde\{D\}\_\{M\}\}\{\\partial\\lambda\_\{k\}\}=\-\\frac\{z\_\{k\}^\{2\}\\exp\(\-\\lambda\_\{k\}\)\}\{2\(1\+D\_\{M\}\-\\tau\)D\_\{M\}\},which approaches a finite constant \(→−1/2\\to\-1/2in a single dominant direction\) rather than diverging linearly withzk2​exp⁡\(−λk\)z\_\{k\}^\{2\}\\exp\(\-\\lambda\_\{k\}\)as in the unmodified Laplace case\. We therefore do not claim that the Lie algebra parameterization alone guarantees bounded gradients near the SPD boundary; rather, in our implementation, gradient and matrix\-exponential overflow are controlled jointly by \(i\) eigenvalue clampingλk∈\[λmin,λmax\]\\lambda\_\{k\}\\in\[\\lambda\_\{\\min\},\\lambda\_\{\\max\}\]before the exponential map, and \(ii\) the log\-tail robustification in Eq\.[16](https://arxiv.org/html/2608.24386#A1.E16)\.

#### Comparison: Gaussian vs\. Laplace NLL\.

The key distinctions are summarized below:

Table 5:Comparison of Gaussian and Laplace negative log\-likelihood formulations\.

## Appendix BModel Architecture and Implementation Details

### B\.1Network Architecture and Data Flow

Our architecture implements an E\(3\)\-equivariant neural network using standard message passing layers\. The input atomic numbers are first projected into 119\-dimensional Magpie feature embeddings, which then pass throughL=2L=2interaction layers with a hidden dimension of 64\. The backbone employs SiLU activations for scalar features and gated\-tanh for higher\-order tensors to preserve equivariance throughout computation\. To support rank\-4 covariance output, we set the maximum rotation orderℓm​a​x=4\\ell\_\{max\}=4, enabling the fullSym2​\(ρc\)\\mathrm\{Sym\}^\{2\}\(\\rho\_\{c\}\)representation required by our theoretical decomposition\.

The network branches into two distinct heads that operate in parallel\. The mean head predicts Voigt components throughℓ=0⊕ℓ=2\\ell=0\\oplus\\ell=2irreducible representations, directly outputting the tensor mean prediction\. The covariance head outputs the symmetric tensor basis defined as2×\(ℓ=0\)⊕2×\(ℓ=2\)⊕1×\(ℓ=4\)2\\times\(\\ell=0\)\\oplus 2\\times\(\\ell=2\)\\oplus 1\\times\(\\ell=4\), which is then linearly projected to the Lie algebra elementA⁡\(X\)∈ℝs​y​m6×6A\(X\)\\in\\mathbb\{R\}^\{6\\times 6\}\_\{sym\}\. This dual\-head architecture ensures that both mean and uncertainty predictions respect the underlying geometric symmetries\.

#### Joint Training Stability\.

The Lie algebra parametrization combined with the LE\-ESO loss provides inherent numerical stability that, in our experiments, enables end\-to\-end joint optimization of the mean and covariance heads without gradient detachment\. To validate this robustness, we conducted an ablation study comparing joint training with gradient\-detached training \(where UQ gradients are blocked from flowing back to the backbone\)\. Both approaches achieved comparable performance \(MAE difference<0\.02<0\.02\), with joint training showing marginally better uncertainty calibration\. This empirical finding is consistent with the geometric advantage of the Log\-Euclidean framework: by operating in the flat tangent space𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)rather than on the curved SPD manifold directly, gradients remain well\-conditioned even when the covariance head receives informative error signals from the scoring objective\. We did not observe variance collapse or shortcut learning in this setting\.

### B\.2Implementation Details of the Equivariant Covariance Head

To strictly enforce the symmetry properties of the covariance tensor, we employ theCartesianTensorformalism from thee3nnlibrary\([Geiger and Smidt, 2022](https://arxiv.org/html/2608.24386#bib.bib26)\)\. The covariance of a symmetric rank\-2 tensor is mathematically a rank\-4 tensor𝒞i​j​k​l\\mathcal\{C\}\_\{ijkl\}with specific permutation symmetries\. First, the covariance exhibits symmetry of the first tensor argument such that𝒞i​j​k​l=𝒞j​i​k​l\\mathcal\{C\}\_\{ijkl\}=\\mathcal\{C\}\_\{jikl\}\. Second, it maintains symmetry of the second tensor argument with𝒞i​j​k​l=𝒞i​j​l​k\\mathcal\{C\}\_\{ijkl\}=\\mathcal\{C\}\_\{ijlk\}\. Third, the covariance itself is symmetric, satisfying𝒞i​j​k​l=𝒞k​l​i​j\\mathcal\{C\}\_\{ijkl\}=\\mathcal\{C\}\_\{klij\}\.

In our implementation, we define the output space using the formula"ijkl=jikl=ijlk=klij", which restricts the learnable basis to the subspace ofℝ3×3×3×3\\mathbb\{R\}^\{3\\times 3\\times 3\\times 3\}satisfying these symmetries\. Thee3nnlibrary automatically computes the change\-of\-basis matrix from the irreducible representations \(irreps\) ofS​O​\(3\)SO\(3\)to this symmetric Cartesian basis\. The projection to the6×66\\times 6Kelvin\-Mandel matrixA⁡\(X\)A\(X\)proceeds in two systematic steps\.

First, in the irreps to Cartesian mapping, the features are mapped to the rank\-4 Cartesian tensor𝒞i​j​k​l\\mathcal\{C\}\_\{ijkl\}using the precomputed equivariant basis:

𝒞i​j​k​l=∑L,mwL,m​Yi​j​k​lL,\\mathcal\{C\}\_\{ijkl\}=\\sum\_\{L,m\}w\_\{L,m\}Y^\{L\}\_\{ijkl\},whereYi​j​k​lLY^\{L\}\_\{ijkl\}are the Clebsch\-Gordan coefficients projecting the spherical harmonics onto the Cartesian tensor components\. Second, in the Cartesian to Kelvin\-Mandel transformation, the3×3×3×33\\times 3\\times 3\\times 3tensor is flattened into a6×66\\times 6matrixAK​MA\_\{KM\}using the Kelvin\-Mandel isometry\. This mapping preserves the Frobenius norm \(i\.e\.,‖𝒞‖F=‖AK​M‖F\\\|\\mathcal\{C\}\\\|\_\{F\}=\\\|A\_\{KM\}\\\|\_\{F\}\) by scaling the off\-diagonal shear components by2\\sqrt\{2\}\. For indices mappingi​j→αij\\to\\alphaandk​l→βkl\\to\\beta\(whereα,β∈\{6\}\\alpha,\\beta\\in\\\{1\.\.\\\!6\\\}\), the entryAα​βA\_\{\\alpha\\beta\}is given by:

Aα​β=ηα​ηβ​𝒞i​j​k​l,A\_\{\\alpha\\beta\}=\\eta\_\{\\alpha\}\\eta\_\{\\beta\}\\mathcal\{C\}\_\{ijkl\},whereη=\{1,1,1,2,2,2\}\\eta=\\\{1,1,1,\\sqrt\{2\},\\sqrt\{2\},\\sqrt\{2\}\\\}corresponds to the indices\{x​x,y​y,z​z,y​z,x​z,x​y\}\\\{xx,yy,zz,yz,xz,xy\\\}\. This construction guarantees that the predicted matrixA⁡\(X\)A\(X\)strictly lies in the symmetric subspace𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)and transforms exactly according toρc⊗ρc\\rho\_\{c\}\\otimes\\rho\_\{c\}\.

### B\.3Training Protocol and Stability Measures

All models were optimized using AdamW with hyperparametersβ1=0\.9,β2=0\.999\\beta\_\{1\}=0\.9,\\beta\_\{2\}=0\.999and a weight decay of10−410^\{\-4\}\. We employed a OneCycleLR scheduler with a peak learning rate of10−310^\{\-3\}, warming up for 20% of the total 50 epochs before gradually decaying\. To ensure training stability during the critical early phases of covariance learning, we implemented several key strategies\.

We introduced acritical loss annealing strategywhere the auxiliary MSE warmup weightλMSE\\lambda\_\{\\mathrm\{MSE\}\}gradually decays from0\.90\.9to full LE\-ESO optimization by epoch 5 \(note thatλMSE\\lambda\_\{\\mathrm\{MSE\}\}is distinct from the LE\-ESO weightα\\alphain Eq\.[12](https://arxiv.org/html/2608.24386#S3.E12)\)\. This gradual transition is essential for preventing early training instability \- it allows the network to first learn reasonable mean predictions before tackling the more complex uncertainty quantification task\. Furthermore, we applied eigenvalue clamping within\[λmin,λmax\]=\[−4,3\]\[\\lambda\_\{\\min\},\\lambda\_\{\\max\}\]=\[\-4,3\]to prevent numerical overflow in the matrix exponential computation\. This constraint ensures that the resulting covariance eigenvalues remain in\[e−4,e3\]≈\[0\.018,20\.1\]\[e^\{\-4\},e^\{3\}\]\\approx\[0\.018,20\.1\], preventing both variance collapse and explosion during early training\.

#### Anisotropic Jitter as Numerical Gradient Stabilizer\.

A subtle numerical issue arises in automatic differentiation of eigenvalue decompositions: when the covariance matrix has degenerate eigenvalues \(λi=λj\\lambda\_\{i\}=\\lambda\_\{j\}\), the Jacobian contains singular terms\(λi−λj\)−1\(\\lambda\_\{i\}\-\\lambda\_\{j\}\)^\{\-1\}\. This is particularly problematic for high\-symmetry crystals \(e\.g\., cubic systems\) where physical symmetry can cause eigenvalue degeneracy\. To ensure differentiability, we introduce anumerical gradient stabilizer—a tiny anisotropic perturbation applied only during the spectral decomposition step:

λ~i=λi\+ϵ⋅i,i=0,…,5,\\tilde\{\\lambda\}\_\{i\}=\\lambda\_\{i\}\+\\epsilon\\cdot i,\\quad i=0,\\ldots,5,\(18\)withϵ≈10−6\\epsilon\\approx 10^\{\-6\}for float64 precision\. This design is crucial: anisotropicshift \(ϵ⋅I\\epsilon\\cdot I\) would preserve degeneracy and fail to resolve the singularity, whereas the anisotropic pattern guaranteesλi−λj≠0\\lambda\_\{i\}\-\\lambda\_\{j\}\\neq 0for alli≠ji\\neq j\. Importantly, this jitter isnotan architectural choice—it is a numerical safeguard with magnitudeO⁡\(10−6\)O\(10^\{\-6\}\)that is negligible compared to typical eigenvalue scales \(∼1\\sim 1\)\. Empirically, our equivariance verification \(Table[4](https://arxiv.org/html/2608.24386#S4.T4)\) shows errors on the order of10−710^\{\-7\}, confirming that this minimal perturbation does not compromise the geometric fidelity of the learned representations\. The jitter operates entirely within the numerical solver and is invisible to the upstream equivariant architecture\.

#### Discussion on Jitter and Calibration Impact\.

The introduced jitter \(ϵ≈10−6\\epsilon\\approx 10^\{\-6\}\) is several orders of magnitude smaller than the predicted eigenvalues \(∼1\.0\\sim 1\.0\)\. We observe that this perturbation is essential for maintaining stable gradients during joint training but has a negligible impact on both calibration \(MACE change<10−4<10^\{\-4\}\) and equivariance \(errors remain at the level of10−710^\{\-7\}as reported in Table[4](https://arxiv.org/html/2608.24386#S4.T4)\)\.

#### Hyperparameter Details for Numerical Stability\.

We provide the specific hyperparameter values used in our implementation\. For eigenvalue clamping, we use\[λmin,λmax\]=\[−4,3\]\[\\lambda\_\{\\min\},\\lambda\_\{\\max\}\]=\[\-4,3\], which constrains the covariance eigenvalues to\[e−4,e3\]≈\[0\.018,20\.1\]\[e^\{\-4\},e^\{3\}\]\\approx\[0\.018,20\.1\]\. This range was chosen to prevent both variance collapse \(eigenvalues≪1\\ll 1\) and explosion \(eigenvalues≫1\\gg 1\) during early training\. For the Huber robustification, the threshold is set toτ=5\.0\\tau=5\.0, meaning that Mahalanobis distances above 5\.0 transition to logarithmic scaling\. This threshold was selected based on validation set analysis to be significantly above typical well\-predicted samples \(DM≈2D\_\{M\}\\approx 2–33\) while effectively capping the influence of extreme outliers\.

#### Log\-Euclidean Framework and Information Geometry\.

The space of SPD matrices𝒫6\\mathcal\{P\}\_\{6\}is not a vector space but a Riemannian manifold with non\-Euclidean geometry\. Direct optimization on this manifold introduces path\-dependent gradients and numerical instabilities near the boundary\. By working in the tangent space𝔰​𝔶​𝔪​\(6\)\\mathfrak\{sym\}\(6\)—the Lie algebra of symmetric matrices—we obtain a flat Euclidean vector space where standard optimization is geometrically well\-defined\. The matrix exponential serves as the Riemannian exponential map, lifting points from the tangent space to the curved manifold while preserving the geometric structure\.

#### Geometric Interpretation ofα\\alpha\.

The parameterα\\alphacontrols the tightness of the equivariant confidence hull, a process analogous to entropy regularization in information\-theoretic learning\. The termlogdetΣ\\log\\det\\Sigmarepresents the infinitesimal volume element of the uncertainty manifold in the Riemannian geometry of SPD matrices\. By adjustingα\\alpha, we effectively control the trade\-off between: \(i\)Information\-theoretic volume:α​log​detΣ=α​Tr​\(A\)\\alpha\\log\\det\\Sigma=\\alpha\\operatorname\{Tr\}\(A\)penalizes excessive uncertainty spread; and \(ii\)Geometric fit:DMD\_\{M\}measures the normalized prediction error in the metric induced byΣ\\Sigma\. Theoretically, for the standard Multivariate Laplace distribution,α=1\\alpha=1\(no1/21/2coefficient as in the Gaussian case\)\. However, we treatα\\alphaas a tunable hyperparameter to balance model confidence with coverage: largerα\\alphaencourages tighter confidence regions \(lower uncertainty volume\), while smallerα\\alphaallows more conservative uncertainty estimates\. This flexibility is valuable for materials science applications where the true noise level may vary across different datasets and measurement modalities\.

#### Temperature Scaling for Calibration\.

To ensure the predicted covariance tensorsΣ\\Sigmareflect the empirical error distribution, we apply post\-hoc temperature scaling\([Kuleshov et al\., 2018](https://arxiv.org/html/2608.24386#bib.bib33)\)\. The optimal temperatureT≈0\.05T\\approx 0\.05was determined on the validation set via a robust median\-matching strategy\. The small value ofTTreflects the heavy\-tailed nature of the initial residuals, requiring the model to significantly contract its uncertainty hulls after training with the robustified LE\-ESO\. This adjustment yields a calibrated covarianceΣ′=T⋅Σ\\Sigma^\{\\prime\}=T\\cdot\\Sigma, which is equivalent to an additive shiftA′=A\+ln⁡\(T\)​IA^\{\\prime\}=A\+\\ln\(T\)Iin the Lie algebra\. This scaling effectively aligns the predictive distribution with the requirements of scoring rules \(e\.g\., Energy Score\) without affecting the mean prediction or the exact E\(3\)\-equivariance\.

#### Spectral Bounding for Manifold Consistency\.

To maintain numerical consistency with the Riemannian structure of𝒫6\\mathcal\{P\}\_\{6\}, we constrain the Lie algebra eigenvalues to a bounded interval before computingexp⁡\(Λ\)\\exp\(\\Lambda\)\. This spectral bounding ensures that the resulting covariance eigenvalues remain in a geometrically valid range, preventing both variance collapse \(near\-zero eigenvalues\) and explosion \(excessively large eigenvalues\) during early training when predictions may be far from the data manifold\. The bounds\[λmin,λmax\]\[\\lambda\_\{\\min\},\\lambda\_\{\\max\}\]are chosen to map to a physically meaningful covariance spectrum\[eλmin,eλmax\]\[e^\{\\lambda\_\{\\min\}\},e^\{\\lambda\_\{\\max\}\}\]under the matrix exponential\.

#### Laplacian\-Huber Compound Robust Loss\.

To handle extreme outliers in material property data, we introduce aLaplacian\-Huber schemewith two regimes: when the Mahalanobis distanceDMD\_\{M\}is below thresholdτ=5\.0\\tau=5\.0, we apply alinear penalty\(the Laplace distribution core\); whenDMD\_\{M\}exceedsτ\\tau, we switch to alogarithmic penalty\(τ\+log⁡\(1\+DM−τ\)\\tau\+\\log\(1\+D\_\{M\}\-\\tau\)\) that compresses the residual contribution for very largeDMD\_\{M\}\. This design preserves the statistical interpretation of the Multivariate Laplace distribution for normal samples while limiting the influence of extreme outliers, so that the residual term grows logarithmically instead of linearly in the tail\.

#### Gradient Analysis: Practical Control on the Lie Algebra\.

We do not claim that the Lie algebra parameterization by itself yields a uniform gradient bound near the SPD\-cone boundary; rather, in practice, gradient and matrix\-exponential overflow are controlled jointly by eigenvalue clamping and the log\-tail robustification\. Let𝐳=Q⊤​\(𝐜true−μ\)\\mathbf\{z\}=Q^\{\\top\}\(\\mathbf\{c\}\_\{\\text\{true\}\}\-\\mu\)be the rotated residuals in the eigenbasis\. The loss gradient with respect to eigenvaluesΛ=diag​\(λ1,…,λ6\)\\Lambda=\\text\{diag\}\(\\lambda\_\{1\},\\ldots,\\lambda\_\{6\}\)is:

∂ℒLE\-ESO∂Λ=α​I−∂D~M∂Λ\.\\frac\{\\partial\\mathcal\{L\}\_\{\\text\{LE\-ESO\}\}\}\{\\partial\\Lambda\}=\\alpha I\-\\frac\{\\partial\\tilde\{D\}\_\{M\}\}\{\\partial\\Lambda\}\.\(19\)As derived in Proposition[A\.5](https://arxiv.org/html/2608.24386#A1.Thmtheorem5), the Laplace residual term carries an extra1/\(2​DM\)1/\(2D\_\{M\}\)factor relative to the Gaussian case, which reduces the growth rate of∂DM/∂λk\\partial D\_\{M\}/\\partial\\lambda\_\{k\}but does not by itself produce a uniform bound: when a single direction dominatesDMD\_\{M\}, the contribution to∂DM/∂λk\\partial D\_\{M\}/\\partial\\lambda\_\{k\}still grows as\|zk\|exp\(−λk/2\)/2\|z\_\{k\}\|\\exp\(\-\\lambda\_\{k\}/2\)/2asλk→−∞\\lambda\_\{k\}\\to\-\\infty\. We therefore enforce the constraintλk∈\[λmin,λmax\]\\lambda\_\{k\}\\in\[\\lambda\_\{\\min\},\\lambda\_\{\\max\}\]before the matrix exponential, which preventsexp⁡\(−λk\)\\exp\(\-\\lambda\_\{k\}\)from diverging and keeps∂DM/∂λk\\partial D\_\{M\}/\\partial\\lambda\_\{k\}finite over the optimization trajectory\. The log\-tail region \(DM≥τD\_\{M\}\\geq\\tau\) further compresses the residual gradient:∂D~M/∂λk\\partial\\tilde\{D\}\_\{M\}/\\partial\\lambda\_\{k\}approaches a finite constant \(→−1/2\\to\-1/2in a dominant direction\) instead of growing withzk2​exp⁡\(−λk\)z\_\{k\}^\{2\}\\exp\(\-\\lambda\_\{k\}\), mitigating the influence of rare extreme outliers\. Empirically, this combination keeps training stable enough to support end\-to\-end joint optimization without gradient detachment\. The loss maintainsO⁡\(3\)O\(3\)\-invariance since both the trace and matrix exponential preserve equivariance under orthogonal transformations\.

The numerically stable loss in Eq\.[12](https://arxiv.org/html/2608.24386#S3.E12)follows directly from algebraic identities proven in Proposition[A\.5](https://arxiv.org/html/2608.24386#A1.Thmtheorem5)\. Importantly, both the trace termTr⁡\(A\)\\operatorname\{Tr\}\(A\)and the Mahalanobis distanceDM=Δ​𝐜⊤​exp⁡\(−A\)​Δ​𝐜D\_\{M\}=\\sqrt\{\\Delta\\mathbf\{c\}^\{\\top\}\\exp\(\-A\)\\Delta\\mathbf\{c\}\}are invariant under any orthogonal transformationρc​\(R\)\\rho\_\{c\}\(R\):

Tr⁡\(ρc​\(R\)​A​ρc​\(R\)⊤\)=Tr⁡\(A\),DM​\(ρc​\(R\)​Δ​𝐜,ρc​\(R\)​A​ρc​\(R\)⊤\)=DM​\(Δ​𝐜,A\)\.\\operatorname\{Tr\}\(\\rho\_\{c\}\(R\)A\\rho\_\{c\}\(R\)^\{\\top\}\)=\\operatorname\{Tr\}\(A\),\\qquad D\_\{M\}\(\\rho\_\{c\}\(R\)\\Delta\\mathbf\{c\},\\rho\_\{c\}\(R\)A\\rho\_\{c\}\(R\)^\{\\top\}\)=D\_\{M\}\(\\Delta\\mathbf\{c\},A\)\.Consequently, the loss function provides an exact symmetry\-preserving training objective, in contrast to approximate equivariant regularizations or data augmentation\-based approaches\.

Training was conducted on a single NVIDIA RTX 4060 Ti with batch sizes of 32 for ModelNet40 and 16 for the Materials Project, requiring approximately 10 hours for complete convergence\. Additional implementation details include processing atomic structures into graphs with a5\.05\.0Å cutoff distance and using Magpie feature embeddings of dimension 119 for atomic number representations\.

## Appendix CExperimental Setup and Analysis

### C\.1Dataset Configuration and Preprocessing

We evaluate our framework on two distinct datasets that provide complementary validation of our equivariant uncertainty quantification approach\. ModelNet40 serves for geometric validation with physically defined tensor properties, while the Materials Project provides a real\-world materials science application with experimentally relevant predictions\.

#### ModelNet40\.

The dataset comprises 12,311 CAD models across 40 categories\. We adhere to the official split, utilizing 9,843 models for training and 2,468 for testing\. To simulate measurement uncertainty and validate our probabilistic framework, we sampleN=2048N=2048points uniformly from mesh surfaces and apply Gaussian jitter withσn​o​i​s​e=0\.01\\sigma\_\{noise\}=0\.01\. This noise injection creates the aleatoric uncertainty necessary for testing our framework’s ability to capture geometric ambiguity arising from point cloud sampling\.

#### Materials Project Dielectric Dataset\.

We source precomputed dielectric tensor predictions from the Materials Project database\([Barroso\-Luque et al\., 2024](https://arxiv.org/html/2608.24386#bib.bib22);[Jain et al\., 2013](https://arxiv.org/html/2608.24386#bib.bib25)\)\. To ensure data quality and consistency, we apply systematic filtering criteria: \(1\)structure size—we exclude crystals with fewer than 3 atoms or more than 30 atoms to balance computational efficiency and representation learning; \(2\)positive\-definiteness—we verify that all dielectric tensors have eigenvalues strictly greater than10−410^\{\-4\}, excluding numerically singular matrices; \(3\)value range—we remove samples with dielectric constants outside\[−10,50\]\[\-10,50\]or with diagonal entries below 1\.0\. After filtering, the dataset comprises 5,002 crystalline structures, partitioned into 4,236 for training, 485 for validation, and 281 for testing\. We apply Matrix Log\-Normalization with parametersμl​o​g=1\.24\\mu\_\{log\}=1\.24andσl​o​g=0\.86\\sigma\_\{log\}=0\.86to handle the wide dynamic range while preserving the relationships between different crystal structures\.

### C\.2Equivariance Ablation Study Design

To systematically isolate the contributions of equivariance and SPD constraints, we designed four baseline variants that progressively incorporate different architectural components\. Baseline A employs a standard non\-equivariant GNN with a Cholesky covariance head to test the necessity of equivariant message passing\. Baseline B upgrades the backbone to an equivariant neural network \(ENN\), but keeps coordinate\-wise scalar MLP outputs for both the mean and Cholesky covariance heads, evaluating whether equivariant features alone are sufficient to guarantee equivariant tensor outputs\. Baseline B′further replaces the scalar mean head with the same equivariant mean construction used in our model while retaining the Cholesky covariance head\. This baseline isolates the covariance\-parameterization failure mode: the mean branch is equivariant by construction, whereas the Cholesky covariance is SPD but not equivariant under the Kelvin–Mandel covariance representation\. Baseline C uses an ENN backbone with direct equivariant regression of the symmetric operatorA⁡\(X\)A\(X\)but omits the matrix exponential, testing the importance of the SPD projection\. Finally, our full method combines the ENN backbone with the matrix\-exponential covariance head to simultaneously guarantee covariance equivariance and SPD validity\.

## Appendix DAdditional Experimental Results

This appendix collects supplementary experiments that complement the two main experiments in the body of the paper\. They are intended as supporting evidence for the scope and robustness of the proposed equivariant SPD/UQ construction rather than as a comprehensive benchmarking study\.

### D\.1ModelNet40 Shape\-Covariance Validation

The shape\-covariance experiment is a controlled geometric validation benchmark on the ModelNet40 dataset\([Wu et al\., 2015](https://arxiv.org/html/2608.24386#bib.bib27)\)\. Like inertia tensor prediction, the target admits a closed\-form estimator from the point cloud\. We therefore do not present this task as a real\-world setting where neural prediction is necessary\. Instead, it tests whether the proposed equivariant SPD/UQ construction remains valid on a second symmetric rank\-2 tensor target beyond inertia, demonstrating that the framework is not specific to the inertia formulation\.

Table 6:ModelNet40 shape\-covariance validation\. The goal is controlled geometric validation rather than replacing the closed\-form estimator\. “Mean Tensor PSD Rate” refers to the fraction of predicted mean shape\-covariance tensors that are positive semi\-definite, distinct from the SPD validity of the predictive covarianceΣ⁡\(X\)\\Sigma\(X\)\.The point\-prediction MAE/RMSE of the UQ model is comparable to the deterministic baseline, while the predictive covarianceΣ⁡\(X\)\\Sigma\(X\)is exactly E\(3\)\-equivariant and SPD by construction\. As in the inertia setting, we observe near\-machine\-precision equivariance \(errors on the order of10−710^\{\-7\}\) and strict SPD validity for the predictive covariance\.

### D\.2Rank\-4 Elasticity Tensor Prediction

We evaluate the framework on a real\-data elasticity tensor prediction task from the Materials Project\([Jain et al\., 2013](https://arxiv.org/html/2608.24386#bib.bib25)\)\. Unlike the rank\-2 dielectric setting, the mean target here is directly a rank\-4 elasticity tensor\. Under the standard minor and major symmetries, the elasticity tensor has 21 independent components\. This experiment is intended as supporting evidence that the proposed structured equivariant SPD/UQ construction can be extended beyond the six\-dimensional symmetric rank\-2 setting; it is not intended as a comprehensive study of all higher\-order tensor parameterizations\.

On this benchmark, the model achieves a test MAE of approximately 5\.0 GPa, which is essentially on par with a deterministic baseline and noticeably better than a naive UQ baseline\. The structured UQ model also improves uncertainty quality over the naive baseline: empirical coverage rises from approximately35%35\\%to52%52\\%, and the uncertainty–error correlation rises from approximately−0\.15\-0\.15to0\.310\.31\. At the same time, both numerical equivariance/prediction consistency and predictive covariance SPD validity remain at100%100\\%, matching the structural behavior observed in the rank\-2 dielectric setting\. We emphasize that this single experiment is supporting evidence that the proposed structured equivariant SPD/UQ construction remains feasible on a higher\-order tensor target, rather than an exhaustive higher\-order benchmark\.

### D\.3Computational Overhead

We profile per\-batch wall\-clock time on the Materials Project dielectric task using a single NVIDIA RTX 4060 Ti, with batch size 16 and identical input pipelines\.

Table 7:Runtime profiling on Materials Project \(RTX 4060 Ti, batch size 16\)\. The full\-covariance model is more expensive primarily because of the covariance branch and its backpropagation, rather than the matrix exponential alone\.The diagonal\-UQ overhead is negligible \(1\.5%\), confirming that anisotropy modeling, not uncertainty quantification per se, dominates cost\. For inference, a single forward pass yields the full anisotropic covariance, in contrast to ensemble methods that requireNNforward passes\.

### D\.4Sensitivity to the LE\-ESO Weight

The weightα\\alphacontrols the trade\-off between the log\-volume termTr\(A\)=logdetΣ\\mathrm\{Tr\}\(A\)=\\log\\det\\Sigmaand the geometric data\-fit term in LE\-ESO\. We useα=1\\alpha=1in the main experiments because it corresponds to the canonical coefficient in the multivariate Laplace objective motivating LE\-ESO\.

To evaluate sensitivity, we run a short validation sweep overα∈\{0\.03,0\.10,0\.30,1\.00\}\\alpha\\in\\\{0\.03,0\.10,0\.30,1\.00\\\}on the Materials Project dielectric task\. Table[8](https://arxiv.org/html/2608.24386#A4.T8)reports the best validation MAE in log\-Kelvin–Mandel space\. Across the tested values, the validation MAE remains in a moderate range \(0\.3520\.352–0\.4570\.457\), indicating that performance is not tied to a narrow value ofα\\alpha\. The canonical choiceα=1\\alpha=1also gives the lowest validation MAE in this sweep\.

We also report the best validation LE\-ESO value for completeness\. Importantly, this value is evaluated using the sameα\\alphaas the corresponding training run, and therefore should be interpreted as the optimized objective for that setting rather than as a fixed cross\-α\\alphanegative log\-likelihood\. Since changingα\\alphachanges the scoring objective itself, these LE\-ESO values are not directly comparable as absolute NLL values across differentα\\alpha\.

Table 8:Sensitivity to the LE\-ESO weightα\\alphaon the Materials Project dielectric task\. MAE is measured in log\-Kelvin–Mandel space and is comparable across rows\. The LE\-ESO value is evaluated with the sameα\\alphaused for training, so it reflects the optimized objective for each setting rather than a fixed cross\-α\\alphaNLL\.
### D\.5Additional Risk\-Coverage Analysis

To complement the risk\-coverage discussion in the main paper, we report the full retained\-set comparison betweenλmax​\(Σ\)\\lambda\_\{\\max\}\(\\Sigma\)ranking,Trace⁡\(Σ\)\\operatorname\{Trace\}\(\\Sigma\)ranking, and a diagonal\-UQ baseline that ignores off\-diagonal correlations\. At 90% coverage, ranking byλmax\\lambda\_\{\\max\}improves retained\-set MAE by 3\.1% relative to the full test set\. The improvement ofλmax\\lambda\_\{\\max\}overTrace\\operatorname\{Trace\}under the same retained\-set protocol is approximately 1\.5%—smaller than the headline 3\.1% number but consistent across coverage levels\. At 80% coverage,λmax\\lambda\_\{\\max\}continues to retain a positive improvement, whileTrace\\operatorname\{Trace\}\-based ranking can fall slightly below the full\-dataset baseline\. The diagonal\-UQ baseline, which lacks off\-diagonal covariance information, ranks the test set less informatively than eitherλmax\\lambda\_\{\\max\}orTrace\\operatorname\{Trace\}from the full\-covariance model\. These results support the interpretation that directional uncertainty captures failure modes that scalar total uncertainty partially obscures, while clarifying that the practical advantage overTrace\\operatorname\{Trace\}is moderate rather than dramatic\.

## Appendix EAdditional Results and Analysis

### E\.1Training Dynamics and Loss Analysis

Figure[7](https://arxiv.org/html/2608.24386#A5.F7)illustrates the complete training dynamics of our equivariant uncertainty framework\. To ensure a stable optimization landscape, we employ a two\-stage curriculum: the model is initially warmed up with a combined MSE\-LE\-ESO objective for 5 epochs to establish a reliable mean prediction baseline before transitioning to heavy\-tailed LE\-ESO optimization\.

Panel \(a\) reveals that the loss stabilizes rapidly upon transition, with no numerical spikes despite the non\-linear nature of the matrix exponential map\. Panel \(b\) demonstrates that the addition of the uncertainty branch does not compromise the underlying point\-prediction accuracy; instead, the MAE for both diagonal \(εi​i\\varepsilon\_\{ii\}\) and off\-diagonal \(εi​j\\varepsilon\_\{ij\}\) components plateaus at a state\-of\-the\-art level, benefiting from the robust regularization provided by the UQ branch\.

Most importantly, panel \(c\) highlights the sophisticated trade\-off mechanism inherent in our loss formulation\. As the validation epoch progresses, the network balances the data fit term \(Mahalanobis distance\) against the uncertainty regularization term \(logdetΣ\\log\\det\\Sigma\)\. The joint optimization allows both branches to benefit from shared geometric representations, preventing “shortcut learning” where the model might collapse its uncertainty to minimize the scoring rule\. The eventual convergence of the Mahalanobis distance toward a steady value confirms that the model has effectively learned to characterize the aleatoric noise in the dielectric property space\.

Figure 7:Training dynamics\.Two\-stage optimization: warmup \(5 epochs\) then LE\-ESO\. \(a\) LE\-ESO convergence\. \(b\) MAE stability for diagonal/off\-diagonal components\. \(c\) Balance between data fit \(𝔼⁡\[DM\]\\mathbb\{E\}\[D\_\{M\}\]\) and regularization \(logdetΣ\\log\\det\\Sigma\)\.
### E\.2Empirical Verification of Theoretical Guarantees

To validate the theoretical guarantees established in Appendix[A](https://arxiv.org/html/2608.24386#A1), we performed rigorous numerical checks throughout training that confirm both the mathematical correctness and practical stability of our implementation\.

The numerical stability of our approach stems from the eigenvalue decompositionA=Q​Λ​Q⊤A=Q\\Lambda Q^\{\\top\}used to compute the loss without explicitly formingΣ=exp⁡\(A\)\\Sigma=\\exp\(A\)\. SinceAAis symmetric, all eigenvaluesλi\\lambda\_\{i\}are real andexp⁡\(λi\)\\exp\(\\lambda\_\{i\}\)remains positive\. In practice, eigenvalue clamping keeps these exponentials bounded, preventing numerical overflow and improving gradient conditioning\. Together with the log\-tail robustification, this enables joint end\-to\-end training without explicit regularization on the covariance spectrum, addressing a critical limitation of direct covariance optimization approaches\.

For equivariance verification, we continuously monitored the relative Frobenius\-norm difference between rotated predictions and transformed predictions:

Eequiv=‖Σ⁡\(R⋅X\)−ρc​\(R\)​Σ​\(X\)​ρc​\(R\)⊤‖F‖Σ⁡\(X\)‖F\.E\_\{\\text\{equiv\}\}=\\frac\{\\\|\\Sigma\(R\\\!\\cdot\\\!X\)\-\\rho\_\{c\}\(R\)\\Sigma\(X\)\\rho\_\{c\}\(R\)^\{\\top\}\\\|\_\{F\}\}\{\\\|\\Sigma\(X\)\\\|\_\{F\}\}\.Across all random rotations tested during training, this error consistently remained at the level of10−710^\{\-7\}, confirming that our implementation achieves near\-machine\-precision equivariance rather than approximate symmetry preservation\.

For SPD validation, we monitored the spectrum of predicted covariance matrices throughout training\. The minimum eigenvalue ofΣ⁡\(X\)\\Sigma\(X\)across all batches remained strictly positive \(\>10−5\>10^\{\-5\}\), with no numerical violations of the SPD constraint observed \(see Figure[8](https://arxiv.org/html/2608.24386#A5.F8)a\)\.

Beyond the geometric validation on ModelNet40, we further analyzed the conditioning of the predicted covariances for the dielectric tensor task \(Materials Project\)\. As shown in Figure[8](https://arxiv.org/html/2608.24386#A5.F8)b, the distribution of condition numbersκ⁡\(Σ\)\\kappa\(\\Sigma\)remains numerically well\-conditioned for the final model, with a mean of 3\.80 and a maximum of 16\.4\.

This result is particularly significant because, unlike the synthetic jitter in ModelNet40, the uncertainty in dielectric tensors arises from complex physical and DFT approximation errors\. The low condition numbers indicate that our matrix exponential mapping naturally induces numerically stable, non\-degenerate uncertainty estimates without requiring auxiliary regularization terms \(e\.g\., hinge loss penalties on eigenvalues\)\. This confirms that the optimization landscape remains well\-behaved even for high\-dimensional material representations\.

### E\.3Reflection Symmetry and Chirality Handling

Our framework explicitly accounts for improper rotations \(reflections\) by ensuring that the representationρc\\rho\_\{c\}correctly tracks the parity of the tensorial outputs\. For the symmetric rank\-2 tensors considered here—such as dielectric or inertia tensors—the physical quantities areeven tensorsunder parity, meaning they are invariant to inversion\. The Kelvin\-Mandel representationρc​\(R\)\\rho\_\{c\}\(R\)used throughout this paper is defined by projectingR⊗RR\\otimes Ronto the symmetric subspace, as given in Eq\.[14](https://arxiv.org/html/2608.24386#A1.E14)of Appendix[A\.1](https://arxiv.org/html/2608.24386#A1.SS1); we use that construction directly here rather than introducing a separate definition\. Since this transformation is built from two factors ofRR, the determinant contribution\(detR\)2=1\(\\det R\)^\{2\}=1ensures that the framework handles chiral structures and their mirror images with consistent physical semantics\. We numerically verified fullO⁡\(3\)O\(3\)equivariance by testing improper rotations, achieving errors at the level of10−710^\{\-7\}consistent with theS​O​\(3\)SO\(3\)results reported in Table[4](https://arxiv.org/html/2608.24386#S4.T4)\. Consequently, our uncertainty quantification remains valid regardless of the handedness of the coordinate system, a critical requirement for modeling both chiral and achiral materials\.

### E\.4Spectral Analysis and Sharpness Distribution

To further investigate the UQ quality, we provide detailed spectral analysis in Figure[8](https://arxiv.org/html/2608.24386#A5.F8)\. The eigenvalue distribution \(Panel a\) confirms that all predicted covariances maintain strict positive\-definiteness with a minimum eigenvalueλmin≈0\.449\\lambda\_\{\\min\}\\approx 0\.449, safely avoiding variance collapse\. The condition number distribution in Figure[8](https://arxiv.org/html/2608.24386#A5.F8)b shows that the predicted covariance matrices remain numerically well\-conditioned, consistent with the verification in Appendix E\.2\. Complementary to the risk\-coverage analysis in Section[4\.4](https://arxiv.org/html/2608.24386#S4.SS4), the sharpness distribution in Figure[5](https://arxiv.org/html/2608.24386#S4.F5)a reveals that the model effectively differentiates between “simple” and “complex” atomic environments by assigning confidence volumes spanning several orders of magnitude\.

![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/paper_figure_spd.png)\(a\)Spectrum Validity
![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/analysis_condition_number.png)\(b\)Conditioning

Figure 8:Numerical stability analysis\.\(a\) Positive eigenvalues ensure SPD validity\. \(b\) Moderate condition numbers indicate numerical stability\.
### E\.5ModelNet40 SPD Analysis

The 3D uncertainty visualization in Figure[3](https://arxiv.org/html/2608.24386#S4.F3)demonstrates that our framework produces physically meaningful uncertainty estimates where uncertainty ellipsoids align with principal shape axes \(demonstrating E\(3\)\-equivariance\), expand in regions with sparse point density \(capturing sampling ambiguity\), and preserve tensorial correlations across components\.

We verify physical consistency through systematic validation of SPD properties \(Figure[9](https://arxiv.org/html/2608.24386#A5.F9)\)\. Our predictions maintain strict SPD requirements \(\>\>99\.9% validity\) with well\-conditioned covariance structures \(median condition number 6\.8\), contrasting sharply with unconstrained baselines that frequently violate physical constraints\.

![Refer to caption](https://arxiv.org/html/2608.24386v1/figures/spd_analysis_test.png)Figure 9:SPD validity on inertia\-tensor task\.\(a\) Minimum eigenvalue distribution \(all positive\)\. \(b\)log10\\log\_\{10\}condition numbers \(well\-conditioned\)\. \(c\) Total uncertaintyTr⁡\(Σ\)\\mathrm\{Tr\}\(\\Sigma\)\.
### E\.6Limitations and Future Work

The SPD construction and scoring objective are representation\-agnostic once an equivariant symmetric operatorA⁡\(X\)A\(X\)is available\. However, extending the full parameterization to higher\-order tensor predictions requires group\- and representation\-specific basis construction\. Our main implementation and empirical validation focus on symmetric rank\-2 tensors, with the rank\-4 elasticity experiment in Appendix[D\.2](https://arxiv.org/html/2608.24386#A4.SS2)serving as preliminary supporting evidence\. Extending to fourth\-order tensors \(e\.g\., elasticity tensors\) and beyond introduces two key computational challenges\. First, the tensor basis construction scales asO⁡\(dℓ\)O\(d^\{\\ell\}\)whereℓ\\ellis the tensor rank, making the basis enumeration for rank\-4 and higher tensors substantially more expensive\. Second, the covariance matrix dimension grows combinatorially—for a rank\-kksymmetric tensor in 3D, the Kelvin\-Mandel representation has dimension\(k\+1\)​\(k\+2\)/2\(k\+1\)\(k\+2\)/2, leading to covariance matrices of sizeO⁡\(k4\)O\(k^\{4\}\)\. This scaling necessitates careful memory management and may require approximations such as low\-rank covariance factorization or hierarchical uncertainty modeling\. Future work should explore more efficient equivariant basis constructions for higher\-order tensors\. In particular, integrating path\-matrix based ICT decompositions\([Shao et al\., 2025](https://arxiv.org/html/2608.24386#bib.bib30)\)could significantly reduce the overhead of basis enumeration for rank\-4 and higher tensors, enabling the extension of our uncertainty framework to complex properties like the full elasticity tensor\.

#### Modularity and Backbone Extensibility\.

A key strength of our framework is its modular design: the matrix\-exponential UQ head is completely backbone\-agnostic and can be integrated with any E\(3\)\-equivariant architecture\. While this study utilizes a standard message\-passing backbone to validate the UQ mechanism, future work will explore pairing our UQ head with higher\-accuracy architectures such as GoeCTP\([Hua et al\., 2026](https://arxiv.org/html/2608.24386#bib.bib29)\)to combine state\-of\-the\-art point prediction with calibrated, symmetry\-preserving uncertainty estimates\. This plug\-and\-play capability allows practitioners to add rigorous uncertainty quantification to existing equivariant models without architectural reengineering\.