Paper deep dive

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

Ruihan Liu, Yu Ji, Jianbo Yu, Shifu Yan, Qingchao Jiang

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 100

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%

Last extracted: 8/29/2026, 4:36:08 AM

Summary

This paper introduces a framework for E(3)-equivariant uncertainty quantification (UQ) in symmetric rank-2 tensor prediction. The method ensures strictly positive-definite (SPD) covariance matrices by mapping from the flat Lie algebra sym(6) to the SPD manifold via matrix exponentiation, preserving rotational symmetry. It utilizes irreducible representation decomposition (Sym^2(rho_c)) and a Log-Euclidean Equivariant Scoring Objective (LE-ESO) based on the Multivariate Laplace distribution for robust optimization. The approach is validated on ModelNet40 inertia tensors and Materials Project dielectric tensors, demonstrating competitive performance and physically consistent uncertainty estimates.

Entities (9)

Equivariant Covariance Tensors · method · 98%Materials Project · dataset · 96%ModelNet40 · dataset · 96%E(3)-equivariant neural networks · model · 95%LE-ESO · lossfunction · 94%SPD manifold · mathematicalconcept · 93%Symmetric rank-2 tensor · datastructure · 92%Kelvin-Mandel coordinates · representation · 90%e3nn · library · 88%

Relation Signals (6)

Equivariant Covariance Tensors → validateson → ModelNet40

confidence 96% · Validation on ModelNet40 inertia tensors and Materials Project dielectric tensors demonstrates that our method achieves competitive performance

Equivariant Covariance Tensors → validateson → Materials Project

confidence 96% · Validation on ModelNet40 inertia tensors and Materials Project dielectric tensors demonstrates that our method achieves competitive performance

Equivariant Covariance Tensors → uses → matrix exponentiation

confidence 95% · By mapping from the flat Lie algebra sym(6) to the curved SPD manifold via matrix exponentiation, we strictly ensure positive-definite covariances

Equivariant Covariance Tensors → uses → LE-ESO

confidence 94% · we formulate a Log-Euclidean Equivariant Scoring Objective (LE-ESO)—a robust surrogate loss based on the Multivariate Laplace distribution

Equivariant Covariance Tensors → decomposes → Sym^2(rho_c)

confidence 92% · Our approach decomposes the covariance into irreducible representations Sym^2(rho_c) ≅ 2×(l=0) ⊕ 2×(l=2) ⊕ 1×(l=4)

E(3)-equivariant neural networks → lacks → rigorous confidence measures

confidence 90% · While E(3)-equivariant neural networks excel at point estimates, they lack rigorous confidence measures.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

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 a $6\times6$ covariance 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 representations $\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 $\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.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

96,504 characters extracted from source content.

Expand or collapse full text

Equivariant Covariance Tensors: Guaranteed SPD Uncertainty for Tensor-Valued Geometric Learning Ruihan Liu Affiliation: College of Intelligent Robotics and Advanced Manufacturing, Fudan University, Shanghai 200433, China Yu Ji Affiliation: College of Intelligent Robotics and Advanced Manufacturing, Fudan University, Shanghai 200433, China Jianbo Yu Affiliation: School of Microelectronics, Fudan University, Shanghai 200433, China Correspondence to: jb_yu@fudan.edu.cn Shifu Yan Affiliation: ByteDance, Beijing, China Qingchao Jiang Affiliation: 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 a 6×66× 6 covariance 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 representations Sym2(ρc)≅2×(l=0)⊕2×(l=2)⊕1×(l=4)Sym^2( _c) 2×(l=0) 2×(l=2) 1×(l=4). By mapping from the flat Lie algebra (6) 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 1 Introduction 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; Heilman et al., 2024), 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 tensor C∈ℝsym3×3C ^3× 3_sym, such as a dielectric tensor. Although C is a 3×33× 3 matrix, symmetry leaves six independent degrees of freedom. We represent it by the Kelvin-Mandel vector ∈ℝ6c ^6, for which rotations act by an orthogonal six-dimensional representation ρc(R) _c(R). Therefore, uncertainty over C is not a 3×33× 3 covariance, but a 6×66× 6 covariance 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 ^6× 6 must satisfy: Σ⁡(R⋅X)=ρc(R)Σ(X)ρc(R)⊤,∀R∈O⁡(3), (R\!·\!X)= _c(R) (X) _c(R) , ∀ R∈ O(3), where ρc(R) _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 decomposition Sym2(ρc)≅2×(ℓ=0)⊕2×(ℓ=2)⊕1×(ℓ=4)Sym^2( _c) 2×( =0) 2×( =2) 1×( =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). 2 Related Work Figure 1: E(3)-equivariant tensor uncertainty framework. (1) Feature Decomposition: covariance head models symmetry via irreducible representations (ℓ=0,2,4 =0,2,4). (2) Manifold Mapping: network predicts unconstrained A∈(6)A∈ sym(6). (3) Validity Constraint: Σ=exp⁡(A) = (A) projects to SPD manifold 6P_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; Heilman et al., 2024; Hua et al., 2026; Pakornchote et al., 2023; Fung et al., 2021; Reiser et al., 2022; Du et al., 2024; Equer et al., 2023). 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 e3n (Geiger and Smidt, 2022), recent advances in representation theory, such as the High-Rank Irreducible Cartesian Tensor (ICT) decomposition framework (Shao et al., 2025), 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; Pouliquen et al., 2025; Zhao et al., 2023). These methods guarantee SPD in a fixed coordinate parameterization, but they do not by themselves provide an equivariant map X↦Σ⁡(X)X (X) under the covariance representation ρc _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 operator A⁡(X)A(X) and setting Σ⁡(X)=exp⁡(A⁡(X)) (X)= (A(X)), which satisfies exp⁡(ρc(R)Aρc(R)⊤)=ρc(R)exp⁡(A)ρc(R)⊤ ( _c(R)A _c(R) )= _c(R) (A) _c(R) 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) and even exhibit emergent equivariance (Gerken and Kessel, 2024), 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; Olivier et al., 2021; Doan et al., 2025; Sheinkman and Wade, 2025). Existing equivariant Bayesian approaches focus on scalar or vector quantities (Zhou et al., 2024), 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; Bevanda et al., 2025)—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; Gruber and Buettner, 2022; Fakour et al., 2024; Newman et al., 2024). 3 Methods 3.1 Problem Formulation We address the fundamental challenge of predicting tensor-valued quantities while quantifying predictive uncertainty. We present the construction for symmetric rank-2 tensors C∈ℝsym3×3C ^3× 3_sym, which already require a full 6×66× 6 covariance representation. The SPD construction and scoring objective are representation-agnostic once an equivariant symmetric operator A⁡(X)A(X) is available. However, the parameterization of A⁡(X)A(X) is representation-specific and must be constructed separately for each tensor order and symmetry group. We predict C from a 3D structure X=(i,fi)i=1NX=\(x_i,f_i)\_i=1^N, where i∈ℝ3x_i ^3 denotes spatial coordinates and fif_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 X)=Laplace(C μ(X), (X)), (1) where μ⁡(X)∈ℝsym3×3μ(X) ^3× 3_sym is the predicted mean tensor and Σ⁡(X) (X) captures the predictive uncertainty through a covariance structure. To maintain geometric consistency across all domains, both the mean prediction μ and covariance Σ must 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) μ(R· X) =Rμ(X)R⊤, =Rμ(X)R , Σ⁡(R⋅X) (R· X) =ρc(R)Σ(X)ρc(R)⊤. = _c(R)\, (X)\, _c(R) . (2) Our construction ensures full O⁡(3)O(3) equivariance. For symmetric rank-2 tensors, the transformation under reflection (detR=−1 R=-1) is handled by the even-parity representation of the Kelvin-Mandel basis ρc(R)=R⊗sR _c(R)=R _sR, where (detR)2=1( R)^2=1 guarantees consistent behavior for both chiral and achiral structures. We numerically verify O⁡(3)O(3) equivariance in Section 4.5 (Table 4), with detailed analysis in Appendix E.3. 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 presents the distribution in standard predictive form, the training optimizes the robustified LE-ESO objective (Section 3.5) which generalizes the Multivariate Laplace negative log-likelihood with enhanced stability against outliers. 3.2 Voigt 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 tensor C into a 66-dimensional vector. Unlike standard Voigt representation, Kelvin-Mandel notation maintains the isometry property between tensor and vector spaces: KM=[C11,C22,C33,2C23,2C13,2C12]⊤,c_KM=[C_11,C_22,C_33, 2C_23, 2C_13, 2C_12] , (3) which preserves the Frobenius norm: ‖KM‖2=‖C‖F\|c_KM\|_2=\|C\|_F. Under rotation R, KM′=ρc(R)KMc_KM = _c(R)c_KM where ρc(R) _c(R) is an orthogonal 6×66× 6 matrix. The covariance transforms as Σ′=ρc(R)Σρc(R)⊤ = _c(R) _c(R) , maintaining coordinate invariance of the physical uncertainty. 3.3 Irreducible 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 _c for symmetric 3×33× 3 tensors decomposes into irreducible representations of SO⁡(3)SO(3) as ρc≅l=0⊕l=2, _c l=0 l=2, (4) corresponding respectively to the isotropic (trace) and deviatoric (traceless) components of the tensor. Here ≅ denotes 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 Σ transforms as ρc⊗ρc _c _c, its representation decomposes as: ρc⊗ρc _c _c =(l=0⊕l=2)⊗(l=0⊕l=2) =(l=0 l=2) (l=0 l=2) =(l=0)⊕2(l=2)⊕(l=4)⊕(l=1,3)antisym. =(l=0) 2(l=2) (l=4) (l=1,3)_antisym. (5) The ℓ=1 =1 and ℓ=3 =3 components lie in the antisymmetric part of ρc⊗ρc _c _c, corresponding to operators that change sign under exchange of the two covariance indices. Since covariance matrices satisfy Σ=Σ⊤ = , only the symmetric square Sym2(ρc)Sym^2( _c) remains, yielding: Sym2(ρc)≅2×(l=0)⊕2×(l=2)⊕1×(l=4),Sym^2( _c) 2×(l=0) 2×(l=2) 1×(l=4), (6) This provides 21 independent degrees of freedom for symmetric 6×66× 6 SPD covariance matrices. 3.4 Equivariant Neural Architecture Figure 1 summarizes the data flow. A shared E(3)-equivariant encoder first maps the input structure X to latent irreducible features. The mean head selects the ℓ=0⊕ℓ=2 =0 =2 components and outputs the Kelvin-Mandel mean vector KM(X)∈ℝ6 μ_KM(X) ^6, which is mapped back to a symmetric 3×33× 3 tensor. The covariance head outputs coefficients in 2×(ℓ=0)⊕2×(ℓ=2)⊕1×(ℓ=4)2×( =0) 2×( =2) 1×( =4), which are assembled into an equivariant symmetric operator A⁡(X)∈(6)A(X)∈ sym(6). Finally, the predictive covariance is Σ⁡(X)=exp⁡(A⁡(X)) (X)= (A(X)), and the loss is computed from the Kelvin-Mandel residual Δ=KM−KM(X) =c_KM- μ_KM(X). The covariance head fΣ:X↦A⁡(X)f_ :X A(X) must satisfy the transformation property A⁡(R⋅X)=ρc(R)A(X)ρc(R)⊤.A(R\!·\!X)= _c(R)\,A(X)\, _c(R) . (7) Following the irreducible representation decomposition in Eq. 6, we construct A⁡(X)A(X) through structured Clebsch-Gordan combinations. Let ϕ(L)(X)∈ℝFL×(2L+1)φ^(L)(X) ^F_L×(2L+1) denote spherical tensor features of order L produced by the equivariant backbone. We index the irreps appearing in Sym2(ρc)Sym^2( _c) by ℐ=(0,1),(0,2),(2,1),(2,2),(4,1),I=\(0,1),\,(0,2),\,(2,1),\,(2,2),\,(4,1)\, where (L,r)(L,r) refers to the r-th copy of the L-irrep, and write A⁡(X)=∑(L,r)∈ℐ∑m=−LLaL,m(r)(X)BL,m(r),A(X)= _(L,r) _m=-L^La^(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)= _f=1^F_Lw_L,r,f\,φ^(L)_f,m(X). Here BL,m(r)∈ℝsym6×6B^(r)_L,m ^6× 6_sym are fixed, input-independent basis matrices for the r-th copy of the L-irrep in Sym2(ρc)Sym^2( _c), computed once from the Clebsch-Gordan decomposition. All input dependence resides in the equivariant coefficients aL,m(r)(X)a^(r)_L,m(X). Under rotation, the coefficient vector (aL,m(r))m=−L(a^(r)_L,m)_m=-L^L and the basis (BL,m(r))m=−L(B^(r)_L,m)_m=-L^L transform by the same Wigner-D(L)D^(L) representation, so their contraction yields A⁡(R⋅X)=ρc(R)A(X)ρc(R)⊤A(R· X)= _c(R)\,A(X)\, _c(R) . This explicit tensor basis construction guarantees that any choice of weights wL,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) 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 the e3n library, 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 ℓmax=4 _max=4. Second, an equivariant linear layer—implemented via a fourth-order Cartesian tensor with symmetry ijkl=jikl=ijlk=klijijkl=jikl=ijlk=klij—assembles these latent features into the symmetric block structure of A, ensuring consistency with the Sym2(ρc)Sym^2( _c) representation while automatically filtering to the l=0,2,4l=0,2,4 components required for rank-4 covariance output. The covariance head uses a residual connection A=Abase⋅I+ΔA=A_base· I+ A, where AbaseA_base provides an isotropic baseline and ΔA A captures 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 output A naturally possesses the correct equivariance properties, setting the stage for positive-definite covariance construction. 3.5 Positive-Definite Covariance Construction and Training Objective Figure 2: Geometric interpretation of equivariant covariance construction. The network operates in the Lie algebra ≅ℝsym6×6 g _sym^6× 6 (left), producing A⁡(X)A(X). The exponential map Σ=exp⁡(A) = (A) (center) projects onto the SPD manifold 6P_6 (right), guaranteeing valid uncertainties while preserving E⁡(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) and parameterize Σ⁡(X) (X) via the matrix exponential mapping: Σ⁡(X)=exp⁡(A⁡(X)), (X)= (A(X)), (9) We thus optimize within the tangent space (6) sym(6)—a flat Euclidean vector space at the identity—and the matrix exponential lifts A⁡(X)A(X) to the SPD manifold 6P_6 while preserving E⁡(3)E(3)-equivariance, since ρc(R) _c(R) acts by orthogonal conjugation (Figure 2). Log-Euclidean Equivariant Scoring Objective (LE-ESO). Given the Kelvin-Mandel residual Δ=KM−KM(X) =c_KM- μ_KM(X), we define the Log-Euclidean Mahalanobis distance DM(A,Δ)=Δ⊤exp⁡(−A)Δ.D_M(A, )= (-A) . (10) To control the influence of extreme outliers, we use a robustified distance D~M D_M with transition threshold τ: D~M=DM,DM<τ,τ+log⁡(1+DM−τ),DM≥τ. D_M= casesD_M,&D_M<τ,\\ τ+ (1+D_M-τ),&D_M≥τ. cases (11) The final LE-ESO objective combines uncertainty volume regularization with the (robustified) data-fit term: ℒLE-ESO=αTr⁡(A)+D~M,L_LE-ESO= (A)+ D_M, (12) where α>0α>0 controls the trade-off between uncertainty volume and data fit. When α=1α=1 and D~M=DM D_M=D_M, this is the multivariate Laplace negative log-likelihood in Log-Euclidean form, a strictly proper scoring rule on (6) sym(6) (Gneiting and Raftery, 2007). With log-tail robustification or α≠1α≠ 1, the objective is a robust surrogate scoring objective, trading strict propriety for outlier stability. We set α=1α=1 in our primary experiments; Appendix D.4 reports a validation sweep over α∈0.03,0.1,0.3,1.0α∈\0.03,0.1,0.3,1.0\ showing stable training, with α=1α=1 also giving the lowest validation MAE. Detailed derivation and gradient analysis are in Appendix A.4. Training Stability via Joint Optimization. In our experiments, the Lie algebra parameterization combined with eigenvalue clamping and the log-tail robustification in Eq. 11 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. 4 Experiments 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), and (i) real-data dielectric tensor prediction on the Materials Project (Barroso-Luque et al., 2024; Jain et al., 2013). These main experiments are complemented by additional studies in Appendix D, including ModelNet40 shape-covariance prediction (Appendix D.1), rank-4 elasticity tensor prediction (Appendix D.2), runtime profiling (Appendix D.3), and sensitivity to the LE-ESO weight α (Appendix D.4). Dataset statistics are detailed in Appendix C.1. We compare our equivariant full-covariance model against two primary baselines: (i) a deterministic model trained with MSE, and (i) a diagonal UQ model assuming independent components. For ablation analysis, we evaluate non-equivariant and non-SPD variants (Section 4.5). All estimators utilize an E(3)-equivariant backbone with ℓmax=4 _max=4 to support rank-4 covariance output; detailed hyperparameters and training protocols are provided in Appendix B. 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.1 Controlled 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 ℐI transforms equivariantly under rotation (ℐ′=RℐR⊤I =RIR ), 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). 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, showing that the construction is not tied to the inertia target. We quantify equivariance error using the relative Frobenius norm Eequiv=‖Σ⁡(R⋅X)−ρc(R)Σ(X)ρc(R)⊤‖F/‖Σ⁡(X)‖FE_equiv=\| (R\!·\!X)- _c(R) (X) _c(R) \|_F/\| (X)\|_F. Prediction accuracy and uncertainty scores are reported in Table 1; equivariance and SPD-validity results are analyzed separately in Table 4, where the equivariance errors remain at the level of 10−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. Table 1: ModelNet40 inertia tensor prediction. Our equivariant full-covariance framework achieves strong performance with exact geometric consistency and physical constraints. Method MAE RMSE LE-ESO SPD Ours (Full Cov) 0.0780.078 0.1280.128 10.6610.66 ✓ Diagonal Cov 0.092 0.156 12.84 ✓ Deterministic (MSE) 0.083 0.135 N/A ✗ Figure 3 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.2 Application: Dielectric Tensor Prediction 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) employs a scalable equivariant architecture. DTNet (Mao et al., 2024) uses universal potential embeddings. MACE-Ens. is a 5-model deep ensemble. Ours (Calib.) applies temperature scaling (T≈0.05T≈ 0.05). ES generalizes CRPS to multivariate settings. Method Type MAE (↓ ) LE-ESO (↓ ) ES (↓ ) SPD? SEConv (Heilman et al., 2024) Point 4.702 N/A N/A N/A DTNet (Mao et al., 2024) Point 1.91 N/A N/A N/A GoeCTP (Hua et al., 2026) Point 1.41 N/A N/A N/A Deterministic (MACE) Point 2.10 N/A N/A N/A MACE Ensemble (N=5N=5) UQ 1.96 -2.55 0.78 ≈ Diagonal UQ UQ 2.25 -1.85 0.87 ✓ Ours (Full Cov.) UQ 1.55 -2.43 0.87 ✓ Ours (Calibrated) UQ 1.55 -2.61 0.66 ✓ 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(Σ) _ ( )) and total uncertainty (Tr⁡(Σ)Tr( )). 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; Jain et al., 2013) with static dielectric tensors computed via DFPT. To ensure data quality and consistency, we apply systematic filtering criteria (detailed in Appendix C.1), including structure size constraints (3≤atoms≤303 ≤ 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 with 5.05.0Å cutoff. Architecture and training details are provided in Appendix B. 4.3 Prediction Accuracy Our goal is not state-of-the-art point prediction alone, but competitive accuracy with full-covariance, symmetry-preserving uncertainty. Table 2 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) (1.91) and GoeCTP (Hua et al., 2026) (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 4a demonstrates that the high R2R^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 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 ∼ 35% for a naive UQ baseline to ∼ 52%, and preserves numerical equivariance and SPD validity. 4.4 Uncertainty 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), which naturally accounts for the heavy-tailed error distributions common in materials data. Figure 4b 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, 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 is conservative (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(Σ)Trace( )) and the directional uncertainty (λmax(Σ) _ ( )). As illustrated in Figure 5b, the maximum eigenvalue λ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(Σ) _ ( ) improves retained-set MAE by 3.1% relative to the full test set, compared with a 0.8% improvement when ranking by Trace(Σ)Trace( ). When compared directly against Trace-based ranking under the same retained-set protocol, the advantage of λ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 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, 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. Training Objective MAE (↓ ) MACE (↓ ) ES (↓ ) Gaussian NLL 1.78 0.092 1.05 Energy Score (ES) 1.64 0.054 0.81 LE-ESO (Ours) 1.55 0.049 0.66 Table 4: Equivariance and SPD-validity analysis. Baseline B′ isolates the Cholesky covariance failure mode (equivariant μ, non-equivariant Σ ). Only the matrix-exponential head achieves both covariance equivariance and strict SPD validity. Method ℰμE_μ ℰΣE_ Mean equiv.? Cov. equiv.? SPD? Baseline A (non-equivariant GNN + Cholesky) 1.36×1001.36× 10^0 1.07×1001.07× 10^0 ✗ ✗ ✓ Baseline B (E3N + coordinate-wise heads) 1.28×1001.28× 10^0 1.03×1001.03× 10^0 ✗ ✗ ✓ Baseline B′ (equivariant mean + Cholesky cov.) 1.43×10−61.43× 10^-6 4.30×10−14.30× 10^-1 ✓ ✗ ✓ Baseline C (E3N + direct equivariant cov.) 7.59×10−77.59× 10^-7 4.76×10−74.76× 10^-7 ✓ ✓ ✗ Ours (E3N + matrix exp.) 2.39×−2.39× 10^-7 2.75×−2.75× 10^-7 ✓ ✓ ✓ 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 a Chemical Substitution Analysis by replacing common atoms in the test set with unseen elements (e.g., Actinides U, Pu; Rare Earths Gd, Sm). As shown in Figure 6, the predicted directional uncertainty λ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. 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.5 Equivariance and Validity Verification Table 4 isolates equivariance and SPD validity across four baselines: A (non-equivariant GNN + Cholesky), B (E3N + coordinate-wise heads), B′ (our equivariant mean head + Cholesky covariance), and C (direct equivariant regression of A⁡(X)A(X) without the matrix exponential). Implementation details are in Appendix C.2. Baseline B′ reaches near-machine-precision ℰμ≈1.4×10−6E_μ\!≈\!1.4× 10^-6 but ℰΣ≈0.43E_ \!≈\!0.43, confirming that the failure is specific to Cholesky: its lower-triangular structure is not preserved under orthogonal conjugation by ρc(R) _c(R). Only the matrix-exponential head achieves O⁡(10−7)O(10^-7) equivariance error with strict SPD validity. Computational Overhead. Appendix D.3 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. 5 Discussion Our primary validated setting is E(3)-equivariant UQ for symmetric rank-2 tensors. The additional shape-covariance experiment in Appendix D.1 shows that the construction is not specific to the inertia target, while the rank-4 elasticity experiment in Appendix D.2 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 operator A⁡(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 Arsigny et al. (2006) V. Arsigny, P. Fillard, X. Pennec, and N. Ayache Log-euclidean metrics for fast and simple calculus on diffusion tensors. Magnetic Resonance in Medicine 56 (2), p. 411–421. Cited by: §3.5. Barroso-Luque et al. (2024) L. Barroso-Luque, M. Shuaibi, X. Fu, B. M. Wood, M. Dzamba, M. Gao, A. Rizvi, C. L. Zitnick, and Z. W. Ulissi Open materials 2024 (omat24) inorganic materials dataset and models. arXiv preprint arXiv:2410.12771. Cited by: §C.1, §4.2, §4. Batatia et al. (2022) I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, and G. Csányi MACE: higher order equivariant message passing neural networks for fast and accurate force fields. Advances in neural information processing systems 35, p. 11423–11436. Cited by: §1, §2. Berman et al. (2026) E. Berman, J. Ginesin, M. Pacini, and R. Walters On Uncertainty Calibration for Equivariant Functions. Transactions on Machine Learning Research. External Links: Link Cited by: §2. Bevanda et al. (2025) P. Bevanda, M. Beier, A. Capone, S. G. Sosnowski, S. Hirche, and A. Lederer Koopman-equivariant gaussian processes. In Proceedings 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, p. 3151–3159. External Links: Link Cited by: §2. Doan et al. (2025) B. G. Doan, A. Shamsi, X. Guo, A. Mohammadi, H. Alinejad-Rokny, D. Sejdinovic, D. Teney, D. C. Ranasinghe, and E. Abbasnejad Bayesian low-rank learning (Bella): a practical approach to Bayesian neural networks. Proceedings of the AAAI Conference on Artificial Intelligence 39 (15), p. 16298–16307. External Links: Document, Link Cited by: §2. Du et al. (2024) H. Du, J. Wang, J. Hui, L. Zhang, and H. Wang DenseGNN: universal and scalable deeper graph neural networks for high-performance property prediction in crystals and molecules. npj Computational Materials 10 (1), p. 292. Cited by: §2. Equer et al. (2023) L. Equer, T. K. Rusch, and S. Mishra Multi-scale message passing neural pde solvers. arXiv preprint arXiv:2302.03580. Cited by: §2. Fakour et al. (2024) F. Fakour, A. Mosleh, and R. Ramezani A structured review of literature on uncertainty in machine learning & deep learning. arXiv preprint arXiv:2406.00332. Cited by: §2. Fung et al. (2021) V. Fung, J. Zhang, E. Juarez, and B. G. Sumpter Benchmarking graph neural networks for materials chemistry. npj Computational Materials 7 (1), p. 84. Cited by: §2. Geiger and Smidt (2022) M. Geiger and T. Smidt E3n: euclidean neural networks. arXiv preprint arXiv:2207.09453. Cited by: §B.2, §2. Gerken and Kessel (2024) J. E. Gerken and P. Kessel Emergent equivariance in deep ensembles. arXiv preprint arXiv:2403.03103. Cited by: §2. Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 102 (477), p. 359–378. Cited by: §3.5. Gruber and Buettner (2022) S. Gruber and F. Buettner Better uncertainty calibration via proper scores for classification and beyond. Advances in Neural Information Processing Systems 35, p. 8618–8632. Cited by: §2. Heilman et al. (2024) A. Heilman, C. Schlesinger, and Q. Yan Equivariant graph neural networks for prediction of tensor material properties of crystals. arXiv preprint arXiv:2406.03563. Cited by: §1, §2, Table 2. Hua et al. (2026) H. Hua, J. Yang, W. Lin, and P. Zhou Revisiting the Canonicalization for Fast and Accurate Crystal Tensor Property Prediction. Proceedings of the AAAI Conference on Artificial Intelligence 40 (1), p. 417–425. External Links: Document, Link Cited by: §E.6, §2, §4.3, Table 2, Table 2, Table 2. Jain et 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 materials 1 (1). Cited by: §C.1, §D.2, §4.2, §4. Jekel et al. (2022) C. F. Jekel, K. E. Swartz, D. A. White, D. A. Tortorelli, and S. E. Watts Neural network layers for prediction of positive definite elastic stiffness tensors. arXiv preprint arXiv:2203.13938. Cited by: §2. Kuleshov et al. (2018) V. Kuleshov, N. Fenner, and S. Ermon Accurate uncertainties for deep learning using calibrated regression. In International Conference on Machine Learning, p. 2796–2804. Cited by: §B.3. Mao et al. (2024) Z. Mao, W. Li, and J. Tan Dielectric tensor prediction for inorganic materials using latent information from preferred potential. npj Computational Materials 10 (1), p. 265. Cited by: §4.3, Table 2, Table 2, Table 2. Newman et al. (2024) E. Newman, L. Horesh, H. Avron, and M. E. Kilmer Stable tensor neural networks for efficient deep learning. Frontiers in Big Data 7, p. 1363978. Cited by: §2. Olivier et al. (2021) A. Olivier, M. D. Shields, and L. Graham-Brady Bayesian neural networks for uncertainty quantification in data-driven materials modeling. Computer methods in applied mechanics and engineering 386, p. 114079. Cited by: §2. Pakornchote et al. (2023) T. Pakornchote, A. Ektarawong, and T. Chotibut StrainTensorNet: predicting crystal structure elastic properties using se(3)-equivariant graph neural networks. Physical Review Research 5 (4), p. 043198. Cited by: §2. Pouliquen et al. (2025) C. Pouliquen, M. Massias, and T. Vayer Schur’s positive-definite network: deep learning in the spd cone with structure. In International Conference on Learning Representations, Vol. 2025, p. 71401–71416. Cited by: §2. Reiser et 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 Materials 3 (1), p. 93. Cited by: §2. Rensmeyer et al. (2024) T. Rensmeyer, B. Craig, D. Kramer, and O. Niggemann High accuracy uncertainty-aware interatomic force modeling with equivariant bayesian neural networks. Digital Discovery 3 (11), p. 2356–2366. Cited by: §2. Rudner et al. (2022) T. G. Rudner, Z. Chen, Y. W. Teh, and Y. Gal Tractable function-space variational inference in bayesian neural networks. Advances in Neural Information Processing Systems 35, p. 22686–22698. Cited by: §2. Shao et al. (2025) S. Shao, Y. Li, Z. Lin, and Q. Cui High-rank irreducible cartesian tensor decomposition and bases of equivariant spaces. Journal of Machine Learning Research 26 (175), p. 1–53. External Links: Link Cited by: §E.6, §2, §3.4. Sheinkman and Wade (2025) A. Sheinkman and S. Wade The architecture and evaluation of bayesian neural networks. arXiv e-prints, p. arXiv–2503. Cited by: §2. Steinert et al. (2025) T. Steinert, D. Ginsbourger, A. Lykke-Møller, O. Christiansen, and H. Moss Integration-free kernels for equivariant gaussian process modelling. In Forty-second International Conference on Machine Learning, Cited by: §2. Wu et al. (2015) Z. Wu, S. Song, A. Khosla, F. Yu, L. Zhang, X. Tang, and J. Xiao 3d shapenets: a deep representation for volumetric shapes. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 1912–1920. Cited by: §D.1, §4. Zhao et al. (2023) W. Zhao, F. Lopez, J. M. Riestenberg, M. Strube, D. Taha, and S. Trettel Modeling graphs beyond hyperbolic: graph neural networks in symmetric positive definite matrices. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, p. 122–139. Cited by: §2. Zhou et al. (2024) X. Zhou, Z. Liu, and H. Xiao Bi-eqno: generalized approximate bayesian inference with an equivariant neural operator framework. arXiv preprint arXiv:2410.16420. Cited by: §2. Appendix A Theoretical 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, (i) the matrix exponential construction ensuring both positive-definiteness and equivariance, and (i) the numerical stability of our loss formulation. A.1 Rotation Matrices in Kelvin-Mandel Space For a rotation matrix R∈SO(3)R∈ SO(3), the corresponding 6×66× 6 transformation matrix ρc(R) _c(R) in Kelvin-Mandel space can be derived from the Kronecker product structure. The vectorization operation vec(C)vec(C) maps a symmetric tensor C to a 9-dimensional vector, and under rotation: vec(C′)=(R⊗R)vec(C),vec(C )=(R R)vec(C), (13) where ⊗ denotes the Kronecker product. The matrix ρc(R) _c(R) is obtained by projecting R⊗R R onto the 6-dimensional symmetric subspace and applying the Kelvin-Mandel scaling matrix P: ρc(R)=⋅(R⊗R)⋅T⋅−1, _c(R)=P·S·(R R)·S^T·P^-1, (14) where S is the 6×96× 9 selection matrix. For practical computation, ρc(R) _c(R) has the explicit block structure: ρc(R)=[R112R122R1322R12R132R11R132R11R12R212R222R2322R22R232R21R232R21R22R312R322R3322R32R332R31R332R31R322R21R312R22R322R23R33R22R33+R23R32R21R33+R23R31R21R32+R22R312R11R312R12R322R13R33R12R33+R13R32R11R33+R13R31R11R32+R12R312R11R212R12R222R13R23R12R23+R13R22R11R23+R13R21R11R22+R12R21]. _c(R)= bmatrixR_11^2&R_12^2&R_13^2& 2R_12R_13& 2R_11R_13& 2R_11R_12\\ R_21^2&R_22^2&R_23^2& 2R_22R_23& 2R_21R_23& 2R_21R_22\\ R_31^2&R_32^2&R_33^2& 2R_32R_33& 2R_31R_33& 2R_31R_32\\ 2R_21R_31& 2R_22R_32& 2R_23R_33&R_22R_33+R_23R_32&R_21R_33+R_23R_31&R_21R_32+R_22R_31\\ 2R_11R_31& 2R_12R_32& 2R_13R_33&R_12R_33+R_13R_32&R_11R_33+R_13R_31&R_11R_32+R_12R_31\\ 2R_11R_21& 2R_12R_22& 2R_13R_23&R_12R_23+R_13R_22&R_11R_23+R_13R_21&R_11R_22+R_12R_21 bmatrix. (15) This explicit form ensures that ρc(R) _c(R) maintains orthogonality in Kelvin-Mandel space: ρc(R)Tρc(R)=I6 _c(R)^T _c(R)=I_6. Note on Voigt vs. Kelvin-Mandel Notation. While standard Voigt notation maps CijC_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 applies 2 2 scaling to shear components, ensuring ‖KM‖2=‖C‖F\|c_KM\|_2=\|C\|_F. This isometric property is crucial for maintaining geometric consistency in uncertainty quantification. A.2 Irreducible Representation Decomposition Proposition A.1 (Irreducible decomposition of the covariance representation). Let ρc _c denote the 6-dimensional real representation of SO⁡(3)SO(3) corresponding to symmetric rank-2 tensors, i.e. ρc≅l=0⊕l=2. _c l=0 l=2. Then the symmetric tensor product representation of ρc _c decomposes as Sym2(ρc)≅ 2×(l=0)⊕ 2×(l=2)⊕ 1×(l=4),Sym^2( _c)\; \;2×(l=0)\; \;2×(l=2)\; \;1×(l=4), which possesses 2121 independent degrees of freedom—equal to that of a symmetric 6×66× 6 covariance matrix. Proof. Since ρc≅l=0⊕l=2 _c l=0 l=2, the symmetric square decomposes as Sym2(ρc)≅Sym2(l=0)⊕(l=0⊗l=2)⊕Sym2(l=2).Sym^2( _c) ^2(l=0) (l=0 l=2) ^2(l=2). We have Sym2(l=0)=l=0Sym^2(l=0)=l=0, l=0⊗l=2=l=2l=0 l=2=l=2, and Sym2(l=2)=l=0⊕l=2⊕l=4Sym^2(l=2)=l=0 l=2 l=4. Combining these yields Sym2(ρc)≅2×(l=0)⊕2×(l=2)⊕1×(l=4),Sym^2( _c) 2×(l=0) 2×(l=2) 1×(l=4), which has dimension 2⋅1+2⋅5+1⋅9=212· 1+2· 5+1· 9=21. ∎ A.3 Matrix Exponential Properties Proposition A.2 (Positive-definiteness and equivariance of the exponential map). Let A⁡(X)∈ℝsym6×6A(X) ^6× 6_sym satisfy the equivariance condition A⁡(R⋅X)=ρc(R)A(X)ρc(R)⊤∀R∈O⁡(3).A(R\!·\!X)= _c(R)\,A(X)\, _c(R) ∀ R∈ O(3). Then the matrix exponential Σ⁡(X)=exp⁡(A⁡(X)) (X)= (A(X)) is (i) symmetric positive-definite for all X, and (i) equivariant under the same group action: Σ⁡(R⋅X)=ρc(R)Σ(X)ρc(R)⊤. (R\!·\!X)= _c(R)\, (X)\, _c(R) . Proof. For any real symmetric A, there exists an orthogonal Q and real diagonal Λ such that A=QΛQ⊤A=Q Q . Then exp⁡(A)=Qexp⁡(Λ)Q⊤, (A)=Q\, ( )\,Q , where exp⁡(Λ) ( ) has strictly positive diagonal entries exp⁡(λi)>0 ( _i)>0. Thus exp⁡(A) (A) is symmetric positive-definite. For equivariance, note that ρc(R) _c(R) is orthogonal. The matrix exponential satisfies exp⁡(SAS−1)=Sexp⁡(A)S−1 (SAS^-1)=S (A)S^-1 for any invertible S. Taking S=ρc(R)S= _c(R) gives exp⁡(ρc(R)Aρc(R)⊤)=ρc(R)exp⁡(A)ρc(R)⊤, ( _c(R)A _c(R) )= _c(R) (A) _c(R) , which proves equivariance. ∎ Proposition A.3 (Equivariance of spectral functions). Let f:ℝ→ℝf:R be a scalar function. For a symmetric matrix A with eigenvalue decomposition A=QΛQ⊤A=Q Q , define the spectral function F⁡(A)=Qdiag(f⁡(λ1),…,f⁡(λn))Q⊤F(A)=Q\,diag(f( _1),…,f( _n))\,Q . Since ρc(R) _c(R) is orthogonal for all R∈O⁡(3)R∈ O(3), it follows that F⁡(ρc(R)Aρc(R)⊤)=ρc(R)F(A)ρc(R)⊤.F( _c(R)\,A\, _c(R) )= _c(R)\,F(A)\, _c(R) . Thus, eigenvalue clamping and anisotropic jitter (as defined in Eq. 18) preserve O⁡(3)O(3) equivariance. Proof. For any orthogonal matrix U, the spectral function commutes with orthogonal similarity transformations: F⁡(UAU⊤)=UF(A)U⊤F(UAU )=UF(A)U . This follows from the fact that UAU⊤=QΛQ⊤UAU =Q Q where Q=UQ0Q=UQ_0 for the original eigenvectors Q0Q_0 of A. Applying the definition of F: F⁡(UAU⊤)=(UQ0)diag(f⁡(λi))(UQ0)⊤=U⁡(Q0diag(f⁡(λi))Q0⊤)U⊤=UF(A)U⊤.F(UAU )=(UQ_0)\,diag(f( _i))\,(UQ_0) =U(Q_0\,diag(f( _i))\,Q_0 )U =UF(A)U . Taking U=ρc(R)U= _c(R) completes the proof. ∎ A.4 Numerically 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). Let A be a symmetric matrix, Σ=exp⁡(A) = (A), and Δ=true−μ =c_true-μ. The standard Gaussian negative log-likelihood is ℒGauss=12logdetΣ+12Δ⊤Σ−1Δ.L_Gauss= 12 + 12 ^-1 . Then the following loss is algebraically equivalent and numerically stable: ℒGauss=12Tr⁡(A)+12Δ⊤exp⁡(−A)Δ. L_Gauss= 12Tr(A)+ 12 (-A) . Proof. Using Σ=exp⁡(A) = (A) and the identity det(exp⁡(A))=exp⁡(Tr⁡(A)) ( (A))= (Tr(A)), we obtain logdetΣ=Tr(A) =Tr(A). Since (exp⁡(A))−1=exp⁡(−A)( (A))^-1= (-A), the Mahalanobis term becomes Δ⊤exp⁡(−A)Δ (-A) . ∎ Proposition A.5 (Multivariate Laplace NLL in Log-Euclidean form). Let A be a symmetric matrix, Σ=exp⁡(A) = (A), and Δ=true−μ =c_true-μ. The Multivariate Laplace negative log-likelihood (with unit scale) is ℒLaplace=logdetΣ+Δ⊤Σ−1Δ.L_Laplace= + ^-1 . The numerically stable form in the Lie algebra (6) sym(6) is: ℒLaplace=Tr⁡(A)+DM, L_Laplace=Tr(A)+D_M, where DM=Δ⊤exp⁡(−A)ΔD_M= (-A) is the Mahalanobis distance in the Log-Euclidean metric. Proof. The log-determinant term follows identically: logdetΣ=Tr(A) =Tr(A). For the Mahalanobis distance term, note that Σ−1=exp⁡(−A) ^-1= (-A), so: DM=Δ⊤Σ−1Δ=Δ⊤exp⁡(−A)Δ.D_M= ^-1 = (-A) . The key difference from the Gaussian case is the square root: the Laplace NLL is linear in DMD_M rather than quadratic in DM2D_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 a Laplacian-Huber scheme that combines two complementary robustness mechanisms. For the Mahalanobis distance DM=Δ⊤Σ−1ΔD_M= ^-1 , our robust loss is: D~M=DM,DM<τ(linear region)τ+log⁡(1+DM−τ),DM≥τ(log-tail region) D_M= casesD_M,&D_M<τ (linear region)\\ τ+ (1+D_M-τ),&D_M≥τ (log-tail region) cases (16) This design has a clear statistical interpretation: (1) the linear region (DM<τD_M<τ) preserves the core Laplace distribution assumption, providing natural robustness through a linear rather than quadratic penalty on residuals; (2) the log-tail region (DM≥τD_M≥τ) 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τ=5.0 based on validation analysis—approximately 99% of well-predicted samples have DM<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 decomposition A=QΛQ⊤A=Q Q , define the whitened residual =Q⊤Δz=Q . The Mahalanobis distance becomes: DM=∑i=16zi2exp⁡(−λi).D_M= _i=1^6z_i^2 (- _i). The gradient with respect to eigenvalues Λ in the linear region (DM<τD_M<τ) is: ∂ℒLaplace∂λk=1−zk2exp⁡(−λk)2DM, _Laplace∂ _k=1- z_k^2 (- _k)2D_M, (17) Compared to the Gaussian gradient 12(1−zk2exp⁡(−λk)) 12(1-z_k^2 (- _k)), the Laplace gradient carries an additional 1/(2DM)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 direction k dominates DMD_M (so DM≈|zk|exp(−λk/2)D_M≈|z_k| (- _k/2)), the contribution behaves like |zk|exp(−λk/2)/2|z_k| (- _k/2)/2 and is unbounded as λk→−∞ _k→-∞. In the logarithmic tail region (DM≥τD_M≥τ), the gradient becomes: ∂D~M∂λk=−zk2exp⁡(−λk)2(1+DM−τ)DM, ∂ D_M∂ _k=- z_k^2 (- _k)2(1+D_M-τ)D_M, which approaches a finite constant (→−1/2→-1/2 in a single dominant direction) rather than diverging linearly with zk2exp⁡(−λk)z_k^2 (- _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] _k∈[ _ , _ ] before the exponential map, and (i) the log-tail robustification in Eq. 16. Comparison: Gaussian vs. Laplace NLL. The key distinctions are summarized below: Property Gaussian NLL Laplace NLL (LE-ESO) Mahalanobis term 12DM2 12D_M^2 DMD_M Penalty shape Quadratic Linear Log-det coefficient 12 12 11 (tunable as α) Outlier sensitivity High (quadratic growth) Low (linear growth) Gradient for large DMD_M ∝DM D_M ∝1 1 Statistical assumption Light-tailed errors Heavy-tailed errors Table 5: Comparison of Gaussian and Laplace negative log-likelihood formulations. Appendix B Model Architecture and Implementation Details B.1 Network 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 through L=2L=2 interaction 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 ℓmax=4 _max=4, enabling the full Sym2(ρc)Sym^2( _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 =0 =2 irreducible representations, directly outputting the tensor mean prediction. The covariance head outputs the symmetric tensor basis defined as 2×(ℓ=0)⊕2×(ℓ=2)⊕1×(ℓ=4)2×( =0) 2×( =2) 1×( =4), which is then linearly projected to the Lie algebra element A⁡(X)∈ℝsym6×6A(X) ^6× 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) 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.2 Implementation Details of the Equivariant Covariance Head To strictly enforce the symmetry properties of the covariance tensor, we employ the CartesianTensor formalism from the e3n library (Geiger and Smidt, 2022). The covariance of a symmetric rank-2 tensor is mathematically a rank-4 tensor ijklC_ijkl with specific permutation symmetries. First, the covariance exhibits symmetry of the first tensor argument such that ijkl=jiklC_ijkl=C_jikl. Second, it maintains symmetry of the second tensor argument with ijkl=ijlkC_ijkl=C_ijlk. Third, the covariance itself is symmetric, satisfying ijkl=klijC_ijkl=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×3R^3× 3× 3× 3 satisfying these symmetries. The e3n library automatically computes the change-of-basis matrix from the irreducible representations (irreps) of SO(3)SO(3) to this symmetric Cartesian basis. The projection to the 6×66× 6 Kelvin-Mandel matrix A⁡(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 ijklC_ijkl using the precomputed equivariant basis: ijkl=∑L,mwL,mYijklL,C_ijkl= _L,mw_L,mY^L_ijkl, where YijklLY^L_ijkl are the Clebsch-Gordan coefficients projecting the spherical harmonics onto the Cartesian tensor components. Second, in the Cartesian to Kelvin-Mandel transformation, the 3×3×3×33× 3× 3× 3 tensor is flattened into a 6×66× 6 matrix AKMA_KM using the Kelvin-Mandel isometry. This mapping preserves the Frobenius norm (i.e., ‖F=‖AKM‖F\|C\|_F=\|A_KM\|_F) by scaling the off-diagonal shear components by 2 2. For indices mapping ij→αij→α and kl→βkl→β (where α,β∈6α,β∈\1..\!6\), the entry AαβA_αβ is given by: Aαβ=ηαηβijkl,A_αβ= _α _βC_ijkl, where η=1,1,1,2,2,2η=\1,1,1, 2, 2, 2\ corresponds to the indices xx,yy,zz,yz,xz,xy\x,y,z,yz,xz,xy\. This construction guarantees that the predicted matrix A⁡(X)A(X) strictly lies in the symmetric subspace (6) sym(6) and transforms exactly according to ρc⊗ρc _c _c. B.3 Training Protocol and Stability Measures All models were optimized using AdamW with hyperparameters β1=0.9,β2=0.999 _1=0.9, _2=0.999 and a weight decay of 10−410^-4. We employed a OneCycleLR scheduler with a peak learning rate of 10−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 a critical loss annealing strategy where the auxiliary MSE warmup weight λMSE _MSE gradually decays from 0.90.9 to full LE-ESO optimization by epoch 5 (note that λMSE _MSE is distinct from the LE-ESO weight α in Eq. 12). 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][ _ , _ ]=[-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]≈[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 _i= _j), the Jacobian contains singular terms (λi−λj)−1( _i- _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 a numerical gradient stabilizer—a tiny anisotropic perturbation applied only during the spectral decomposition step: λ~i=λi+ϵ⋅i,i=0,…,5, λ_i= _i+ε· i, i=0,…,5, (18) with ϵ≈10−6ε≈ 10^-6 for float64 precision. This design is crucial: an isotropic shift (ϵ⋅Iε· I) would preserve degeneracy and fail to resolve the singularity, whereas the anisotropic pattern guarantees λi−λj≠0 _i- _j≠ 0 for all i≠ji≠ j. Importantly, this jitter is not an architectural choice—it is a numerical safeguard with magnitude O⁡(10−6)O(10^-6) that is negligible compared to typical eigenvalue scales (∼1 1). Empirically, our equivariance verification (Table 4) shows errors on the order of 10−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ε≈ 10^-6) is several orders of magnitude smaller than the predicted eigenvalues (∼1.0 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 of 10−710^-7 as reported in Table 4). Hyperparameter Details for Numerical Stability. We provide the specific hyperparameter values used in our implementation. For eigenvalue clamping, we use [λmin,λmax]=[−4,3][ _ , _ ]=[-4,3], which constrains the covariance eigenvalues to [e−4,e3]≈[0.018,20.1][e^-4,e^3]≈[0.018,20.1]. This range was chosen to prevent both variance collapse (eigenvalues ≪1 1) and explosion (eigenvalues ≫1 1) during early training. For the Huber robustification, the threshold is set to τ=5.0τ=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≈ 2–33) while effectively capping the influence of extreme outliers. Log-Euclidean Framework and Information Geometry. The space of SPD matrices 6P_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) 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 α. The parameter α controls the tightness of the equivariant confidence hull, a process analogous to entropy regularization in information-theoretic learning. The term logdetΣ represents the infinitesimal volume element of the uncertainty manifold in the Riemannian geometry of SPD matrices. By adjusting α, we effectively control the trade-off between: (i) Information-theoretic volume: αlogdetΣ=αTr(A)α = (A) penalizes excessive uncertainty spread; and (i) Geometric fit: DMD_M measures the normalized prediction error in the metric induced by Σ . Theoretically, for the standard Multivariate Laplace distribution, α=1α=1 (no 1/21/2 coefficient as in the Gaussian case). However, we treat α as a tunable hyperparameter to balance model confidence with coverage: larger α encourages tighter confidence regions (lower uncertainty volume), while smaller α allows 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 Σ reflect the empirical error distribution, we apply post-hoc temperature scaling (Kuleshov et al., 2018). The optimal temperature T≈0.05T≈ 0.05 was determined on the validation set via a robust median-matching strategy. The small value of T reflects 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⋅Σ =T· , which is equivalent to an additive shift A′=A+ln⁡(T)IA =A+ (T)I in 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 6P_6, we constrain the Lie algebra eigenvalues to a bounded interval before computing exp⁡(Λ) ( ). 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][ _ , _ ] are chosen to map to a physically meaningful covariance spectrum [eλmin,eλmax][e _ ,e _ ] under the matrix exponential. Laplacian-Huber Compound Robust Loss. To handle extreme outliers in material property data, we introduce a Laplacian-Huber scheme with two regimes: when the Mahalanobis distance DMD_M is below threshold τ=5.0τ=5.0, we apply a linear penalty (the Laplace distribution core); when DMD_M exceeds τ, we switch to a logarithmic penalty (τ+log⁡(1+DM−τ)τ+ (1+D_M-τ)) that compresses the residual contribution for very large DMD_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−μ)z=Q (c_true-μ) be the rotated residuals in the eigenbasis. The loss gradient with respect to eigenvalues Λ=diag(λ1,…,λ6) =diag( _1,…, _6) is: ∂ℒLE-ESO∂Λ=αI−∂D~M∂Λ. _LE-ESO∂ =α I- ∂ D_M∂ . (19) As derived in Proposition A.5, the Laplace residual term carries an extra 1/(2DM)1/(2D_M) factor relative to the Gaussian case, which reduces the growth rate of ∂DM/∂λk∂ D_M/∂ _k but does not by itself produce a uniform bound: when a single direction dominates DMD_M, the contribution to ∂DM/∂λk∂ D_M/∂ _k still grows as |zk|exp(−λk/2)/2|z_k| (- _k/2)/2 as λk→−∞ _k→-∞. We therefore enforce the constraint λk∈[λmin,λmax] _k∈[ _ , _ ] before the matrix exponential, which prevents exp⁡(−λk) (- _k) from diverging and keeps ∂DM/∂λk∂ D_M/∂ _k finite over the optimization trajectory. The log-tail region (DM≥τD_M≥τ) further compresses the residual gradient: ∂D~M/∂λk∂ D_M/∂ _k approaches a finite constant (→−1/2→-1/2 in a dominant direction) instead of growing with zk2exp⁡(−λk)z_k^2 (- _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 maintains O⁡(3)O(3)-invariance since both the trace and matrix exponential preserve equivariance under orthogonal transformations. The numerically stable loss in Eq. 12 follows directly from algebraic identities proven in Proposition A.5. Importantly, both the trace term Tr⁡(A)Tr(A) and the Mahalanobis distance DM=Δ⊤exp⁡(−A)ΔD_M= (-A) are invariant under any orthogonal transformation ρc(R) _c(R): Tr⁡(ρc(R)Aρc(R)⊤)=Tr⁡(A),DM(ρc(R)Δ,ρc(R)Aρc(R)⊤)=DM(Δ,A).Tr( _c(R)A _c(R) )=Tr(A), D_M( _c(R) , _c(R)A _c(R) )=D_M( ,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 a 5.05.0Å cutoff distance and using Magpie feature embeddings of dimension 119 for atomic number representations. Appendix C Experimental Setup and Analysis C.1 Dataset 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 sample N=2048N=2048 points uniformly from mesh surfaces and apply Gaussian jitter with σnoise=0.01 _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; Jain et al., 2013). 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 than 10−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 μlog=1.24 _log=1.24 and σlog=0.86 _log=0.86 to handle the wide dynamic range while preserving the relationships between different crystal structures. C.2 Equivariance 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 operator A⁡(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 D Additional 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.1 ModelNet40 Shape-Covariance Validation The shape-covariance experiment is a controlled geometric validation benchmark on the ModelNet40 dataset (Wu et al., 2015). 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) (X). Method MAE ↓ RMSE ↓ Mean Tensor PSD Rate Deterministic (MSE) 0.0333 0.0616 99.2% Full UQ (Ours) 0.0346 0.0633 99.5% The point-prediction MAE/RMSE of the UQ model is comparable to the deterministic baseline, while the predictive covariance Σ⁡(X) (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 of 10−710^-7) and strict SPD validity for the predictive covariance. D.2 Rank-4 Elasticity Tensor Prediction We evaluate the framework on a real-data elasticity tensor prediction task from the Materials Project (Jain et al., 2013). 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 approximately 35%35\% to 52%52\%, and the uncertainty–error correlation rises from approximately −0.15-0.15 to 0.310.31. At the same time, both numerical equivariance/prediction consistency and predictive covariance SPD validity remain at 100%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.3 Computational 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. Model Time / Batch Relative Cost Deterministic 130 ms 1.00× Diagonal UQ 132 ms 1.015× Full Covariance UQ 570 ms 4.4× 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 require N forward passes. D.4 Sensitivity to the LE-ESO Weight The weight α controls the trade-off between the log-volume term Tr(A)=logdetΣTr(A)= and the geometric data-fit term in LE-ESO. We use α=1α=1 in 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α∈\0.03,0.10,0.30,1.00\ on the Materials Project dielectric task. Table 8 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 α. The canonical choice α=1α=1 also 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 α as the corresponding training run, and therefore should be interpreted as the optimized objective for that setting rather than as a fixed cross-α negative log-likelihood. Since changing α changes the scoring objective itself, these LE-ESO values are not directly comparable as absolute NLL values across different α. Table 8: Sensitivity to the LE-ESO weight α on 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 α used for training, so it reflects the optimized objective for each setting rather than a fixed cross-α NLL. α Best val MAE ↓ Best val LE-ESO ↓ 0.03 0.3634 0.7911 0.10 0.4176 0.7165 0.30 0.4566 0.3469 1.00 0.3519 -2.9393 D.5 Additional Risk-Coverage Analysis To complement the risk-coverage discussion in the main paper, we report the full retained-set comparison between λmax(Σ) _ ( ) ranking, Trace⁡(Σ)Trace( ) ranking, and a diagonal-UQ baseline that ignores off-diagonal correlations. At 90% coverage, ranking by λmax _ improves retained-set MAE by 3.1% relative to the full test set. The improvement of λmax _ over TraceTrace 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 _ continues to retain a positive improvement, while TraceTrace-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 _ or TraceTrace 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 over TraceTrace is moderate rather than dramatic. Appendix E Additional Results and Analysis E.1 Training Dynamics and Loss Analysis Figure 7 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 (εii _i) and off-diagonal (εij _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Σ ). 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]E[D_M]) and regularization (logdetΣ ). E.2 Empirical Verification of Theoretical Guarantees To validate the theoretical guarantees established in Appendix A, 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 decomposition A=QΛQ⊤A=Q Q used to compute the loss without explicitly forming Σ=exp⁡(A) = (A). Since A is symmetric, all eigenvalues λi _i are real and exp⁡(λi) ( _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_equiv= \| (R\!·\!X)- _c(R) (X) _c(R) \|_F\| (X)\|_F. Across all random rotations tested during training, this error consistently remained at the level of 10−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) (X) across all batches remained strictly positive (>10−5>10^-5), with no numerical violations of the SPD constraint observed (see Figure 8a). 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 8b, the distribution of condition numbers κ⁡(Σ)κ( ) 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.3 Reflection Symmetry and Chirality Handling Our framework explicitly accounts for improper rotations (reflections) by ensuring that the representation ρc _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 are even tensors under parity, meaning they are invariant to inversion. The Kelvin-Mandel representation ρc(R) _c(R) used throughout this paper is defined by projecting R⊗R R onto the symmetric subspace, as given in Eq. 14 of Appendix A.1; we use that construction directly here rather than introducing a separate definition. Since this transformation is built from two factors of R, the determinant contribution (detR)2=1( R)^2=1 ensures that the framework handles chiral structures and their mirror images with consistent physical semantics. We numerically verified full O⁡(3)O(3) equivariance by testing improper rotations, achieving errors at the level of 10−710^-7 consistent with the SO(3)SO(3) results reported in Table 4. 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.4 Spectral Analysis and Sharpness Distribution To further investigate the UQ quality, we provide detailed spectral analysis in Figure 8. The eigenvalue distribution (Panel a) confirms that all predicted covariances maintain strict positive-definiteness with a minimum eigenvalue λmin≈0.449 _ ≈ 0.449, safely avoiding variance collapse. The condition number distribution in Figure 8b 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, the sharpness distribution in Figure 5a reveals that the model effectively differentiates between “simple” and “complex” atomic environments by assigning confidence volumes spanning several orders of magnitude. (a) Spectrum Validity (b) Conditioning Figure 8: Numerical stability analysis. (a) Positive eigenvalues ensure SPD validity. (b) Moderate condition numbers indicate numerical stability. E.5 ModelNet40 SPD Analysis The 3D uncertainty visualization in Figure 3 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). 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. Figure 9: SPD validity on inertia-tensor task. (a) Minimum eigenvalue distribution (all positive). (b) log10 _10 condition numbers (well-conditioned). (c) Total uncertainty Tr⁡(Σ)Tr( ). E.6 Limitations and Future Work The SPD construction and scoring objective are representation-agnostic once an equivariant symmetric operator A⁡(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 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 as O⁡(dℓ)O(d ) where ℓ is 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-k symmetric tensor in 3D, the Kelvin-Mandel representation has dimension (k+1)(k+2)/2(k+1)(k+2)/2, leading to covariance matrices of size O⁡(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) 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) 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.

Instant research discovery