RotEqNet: Rotation-Equivariant Network for Fluid Systems with Symmetric High-Order Tensors

In the recent application of scientific modeling, machine learning models are largely applied to facilitate computational simulations of fluid systems. Rotation symmetry is a general property for most symmetric fluid systems. However, in general, current machine learning methods have no theoretical way to guarantee rotational symmetry. By observing an important property of contraction and rotation operation on high-order symmetric tensors, we prove that the rotation operation is preserved via tensor contraction. Based on this theoretical justification, in this paper, we introduce Rotation-Equivariant Network (RotEqNet) to guarantee the property of rotation-equivariance for high-order tensors in fluid systems. We implement RotEqNet and evaluate our claims through four case studies on various fluid systems. The property of error reduction and rotation-equivariance is verified in these case studies. Results from the comparative study show that our method outperforms conventional methods, which rely on data augmentation.


Introduction
With recent developments in data science and computational tools, machine learning algorithms have been increasingly applied in different engineering and science areas to model physical phenomena. The data from physical experiments and numerical simulations are a source of knowledge about the physical world, on which data-driven methods could be performed to extract new physical laws [1,2,3,4,5,6]. For example, in turbulence RANS modeling in fluid mechanics, traditional modeling methods have failed in many flow scenarios. A unified RANS model that can successfully describe complex flows, including boundary layer, a strong rotation, separation still does not exist according to the author's knowledge [7,8]. On the other hand, advanced measurement and direct numerical simulations provide plenty of data that could be utilized to establish and validate new models.
With the above argument, data-driven methods are particularly suitable for turbulence modeling and some other areas in physics and engineering. There have been many attempts to discover new turbulence models using machine learning methods. Milano and Koumoutsakos [9] reconstruct near-wall flow applying neural networks and compared their results with linear methods (POD). Zhang and Duraisamy [10] used Gaussian process regression combined with an artificial neural network to predict turbulent channel flow and bypass transition. Beck, Flad, and Munz [11] applied residual neural network for Large Eddy Simulation. Chen et. al. proposed an ODE network to generally learn differential equations [12].
The physical laws often appear in the form of tensorial equalities which inherently obey certain types of symmetry. For example, the constitution laws in fluid and solid mechanics should obey translation and rotation invariance [13]. The turbulence RANS model is local tensorial equality between mean velocity gradient and Reynolds stress. The turbulence RANS models should also be rotation invariant [14,15]. However, machine learning methods for RANS modeling do not automatically guarantee rotation invariance, if we use Cartesian components of tensors as input and output of training data. This problem has been addressed by [16,5]. In [16,17], Reynolds stress is expressed as a general expansion of nonlinear integrity basis multiplied by scalar functions of invariants of strain rate and rotation rate tensors. Machine learning is performed to find these scalar functions of tensor invariants of strain rate and rotation rate tensors. Mathematically this expansion comes from an application of the Caylay-Hamilton theory. The special case used in [16,17] is derived by S.B.Pope in [15].
Although such construction is general and possible for higher-order tensors and tensor tuples containing multiple tensors, the number of this basis and the derivation complexity will grow exponentially and become prohibitive for real applications [18,19].
Why would this problem of rotation-equivariance be hard to solve? At first glance, if a system has the property of rotation-equivariance, one has more information for this system. Therefore, this added property of rotationequivariance would lower the performance of a learner. More specifically, adding this new rule of rotation symmetry in a system will require the machine learning algorithm to extract more rules from existing data [20]. In this case, the property of rotation-equivariance could be considered as a continuous group action. There is limited research in the field of deep learning that considers the preservation of symmetries under continuous group actions for physical systems. To address our second point, continuous information is hard to be absorbed. If we consider a machine learning algorithm as an information compression process from input to output [21], a continuous transformation as rotation will be difficult for learning algorithms to absorb.
Given the universal approximation theorem by [22], it would seem that the application of neural networks, especially deep neural networks could solve any problem. As formulated by [23,24,25,26], advanced machine learning methods, especially deep neural networks [27], seem to provide a new opportunity for physical equations approximation. However, in this case of rotation symmetry, if we use a multiple layer perceptron M to learn the relation f , then most likely M does not preserve rotation-equivariance. Generally, the neural network function classes do not satisfy rotation equivariance.
There have been previous works considering group-equivariance with convolutional neural networks in image recognition. A general method has been proposed using group convolution [28,29,30]. Based on the idea of using convolution, several methods composed a steerable filter for rotation-equivariance in convolutional neural networks [31,32,33,34]. However, these works cannot be applied in physical systems as well. One of the most important reasons is that the rotation operation on the image is different from rotation operation on physical systems. Consider a rotation operation on a specific image. We are thinking of a transformation from polar coordinates centering at a certain point [35]. This kind of transformation is different from rotation operation on tensors. Additionally, these methods have a strong restriction that this model must be built on convolutional neural networks. Yet, considering physical systems, convolutional neural networks might not be the best choice since they are designed for image processing.
The problem of rotation-equivariance is also quite impossible to be simply solved by data augmentation and preprocessing. Mentioned by previous works [16], a typical solution is to apply the technique of data augmentation. However, the method of data augmentation fails to have a theoretical guarantee of obtaining the property of rotation-equivariance with finite sample set. Data augmentation method has a theoretical foundation that at infinite sample limit it will asymptotically reach rotation equivariance. However, such a dataset is not only difficult to obtain but also requires much higher computation power while training the model. In the case of using naive preprocessing methods, the problem is that there are limited theoretical tools to deal with high-order tensors, and only limited methods to use for low order tensors. It is hard to apply specific techniques, such as diagonalization, in the case of high-order tensors. Since naive data preprocessing methods are impossible to apply, a more complex method with a theoretical guarantee should be proposed in order to solve this problem.
In this paper, we establish Rotation-Equivariant Network (RotEqNet), a new data-driven framework, which guarantees rotation-equivariance at a theoretical level. Different from previous methods, we first find a method to preserve rotation operation via tensor contraction. In our proposed position standardization algorithm, it could properly link a high-order tensor to a low order tensor with the same rotation operation. By applying mathematical tools for low order matrices (diagonalization and QR factorization), a desired standard position could be derived by the rotation matrix from the previous step. Standard position algorithm is proven to be rotation-invariant in Theorem 3.1, i.e. two tensors differ by a rotation would have the same standard position. Therefore, the learning rules based on standard position are forming a quotient space of the original rules in random rotated plural position [31,36]. In this way, RotEqNet lowers the training difficulty of a randomly positioned dataset. Further, RotEqNet is also proven to be rotation-equivariant, as we have shown in Theorem 3.2. These advantages of RotEqNet would result in an observable error reduction compared to previously introduced data-driven methods. We applied RotE-qNet into four different case studies ranging from second-order, third-order, and fourth-order. These case studies are designed based on Newtonian fluids, Large-eddy simulations, and Electrostriction. Improved performances could be observed for using RotEqNet. The error is reduced for 99.6%, 15.62%, and 54.63% for second, third, forth-order case studies, respectively. Our contribution in this paper is three-fold: 1. We showed an important property of contraction operation on tensors. Contraction operation will preserve rotation operation on tensor with arbitrary order. This is stated in Lemma 2.3.

2.
We propose a properly designed RotEqNet with a position standardization algorithm to guarantee the property of rotation-equivariant. We proved the property of rotation-invariant of position standardization algorithm in Theorem 3.1 and the property of rotation-equivariant of RotEqNet rigorously in Theorem 3.2.
3. We implement our proposed algorithm and the architecture of RotEqNet. We further conduct case studies to show its credibility in design and superiority compared to baseline methods.
To provide a general architecture of our paper, in Section 2 we introduce basic definitions of rotation for arbitrary order tensor (tuples) and related concepts. In Section 2.3 we formulate rotation invariance (equivariance) on supervised learning methods. The RotEqNet and main algorithm is presented in Section 3, and numerical results are shown in Section 4.

Tensor and its operations
In this section, we first introduce an abstract way of defining tensor. One reason for us to introduce the more abstract way to think about tensors is that it provides a convenient formalism for the operations we will do on the tonsorial data discussed in the previous section. The operations are 1. Linear transformation

Contraction
The formalism helps us to provethat these two operations commute which lays theoretical ground for the computation of a representative of rotationally-relatated tensors. We will call this representative standard position

Abstract definition of tensors
Following [37], fix a vector space V of dimension n over R. A tensor product V ⊗ V is a vector space with the property that R-bilinear maps V ×V → R are in natural one-to-one correspondence with R-linear maps V ⊗V → R.
The tensor product V ⊗ V can be constructed as the quotient vector space V × V/C, where C is generated by vectors of the following types where x and y are vectors in V and a is a scalar in R. This means any element in C can be written as a linear combination of vectors of the above form. C is not necessarily a vector space of finite dimension. But the quotient space V ⊗ V is. Let g : V × V → V ⊗ V be the natural projection map, then we use x ⊗ y to denote the image of (x, y) under g.
Let e 1 , · · · , e n be a basis of V, then e i ⊗ e j for i = 1, ..., n and j = 1, ..., n form a basis of V ⊗ V. This means any vector p ∈ V ⊗ V can be written as for some a i j ∈ R.
Here are some relations of tensors which come directly as a consequence of the relations generating C: a(e i ⊗ e j ) = ae i ⊗ e j = e i ⊗ ae j (2.3) (a i e i + a j e j ) ⊗ (a k e k ) = a i a k (e i ⊗ e j ) + a j a k (e j ⊗ e k ) (2.4) The representation of a tensor in V ⊗ V is similar to the representation of a linear map V → V, i.e. a matrix.
In fact, there is a natural way to think of a tensor as a linear map: For each element e i ⊗ e j in the basis of V ⊗ V, we can think of it as a linear map V → V by defining e i ⊗ e j (v) = e i < e j , v >, where <, > is the natural inner product on V. Extend the definition linearly to every element in V ⊗ V, we obtain a way to identify V ⊗ V as the space of linear map V → V. In fact, the tensor i, j a i, j e i ⊗ e j corresponds to the linear map represented by the matrix [a i j ].
We have defined the tensor product V ⊗ V over V. The definition/construction of order k tensor k V ⊗ · · · ⊗ V follows the same course. We will denote order k tensor by ⊗ k V.
The basis of ⊗ k V is given by e i1 ⊗ · · · ⊗ e ik , where i = 1, ..., n and j = 1, ..., k. With respect to this basis, any order k tensor can be written as i1,··· ,ik a i1,...,ik e i1 ⊗ · · · ⊗ e ik . Analogous to the order 2 case, we can think of an order k tensor as a k-dimensional matrix, the typical way a tensor in physical experiments are represented.
We will use T k to denote a tensor of order k, i.e. a vector in ⊗ k V. k is called the rank of the tensor.

Rotation on tensors: a linear transformation
A linear transformation on higher-order tensor is a generalization of a linear transformation on the first-order tensor, i.e. a vector.
Let g : V → V be a linear transformation. Use the basis e 1 , · · · , e n of V, we can represent this expression with the equation Let M(g) denote the matrix representation of g with respect to the basis e 1 , · · · , e n . Then i.e. the transpose of the matrix [a i j ] The map g naturally induces a map ⊗ k g on ⊗ k V. On the basis element e i1 ⊗ · · · ⊗ e ik , the action of ⊗ k g is defined as e i1 ⊗ · · · ⊗ e ik → g(e i1 ) ⊗ · · · ⊗ g(e ik ) (2.7) For any tensor T ∈ ⊗ k V, we will use g(T ) to denote the extension of g on ⊗ k V There is a convenient way to represent a linear transformation of 2-tensor as matrix multiplication.
For a 2-tensor T = i, j b i j e i ⊗ e j , use M(T ) be the matrix whose (i, j) term is b i j .
Lemma 2.1. Rotation operation by matrix R on second-order tensor (matrix) is a change of basis operation.

8)
where × here means the usual matrix multiplication. Lemma 2.1 and remark 2.2 will be used in the proof of Theorem 3.1. As we have shown in this subsection, one could use a matrix form of rotation operation with certain rules of matrix multiplication to perform a rotation on the tensor. In the following proofs of this paper, we applied this idea to perform rotation operation on tensors via matrix multiplication.

Contraction on tensors: reduction of order
Let , be the standard inner product on V. Using this inner product, we can define the contraction of a tensor. It "merges" vectors on the specified axes using the inner product and reduces the rank of the tensor by 2.
Formally, let C(a, b) denote the contraction along axis-a and axis-b. Here, the axis means the ordinal of V in ⊗ k V.
For example, axis-1 refers to the first copy of V in ⊗ k V.
On the element ⊗ k j=1 v i j , C(a, b) acts on it by pairing v ia and v ib via the inner product , , i.e.
wherev means v is not present.
We can then define C(a, b) on ⊗ k V by extending linearly. When k = 2, contraction is nothing other than taking the trace of the corresponding matrix.
Lemma 2.3 shows an interesting connection between rotation operation and contraction. To understand this lemma, it represents that the contraction of a tensor is compatible with a linear transformation if this linear transformation is a rotation. This is an important lemma which is the foundation of the entire analysis in this paper. We would further utilize this lemma for extracting its rotation operation from higher (arbitrary) orders. We show the proof in 5.

Supervised learning setup
In our problem, given data set D = {X i ; y i } i=1,...,N . The data set contains N input-output pairs (X i ; y i ). The input here is a tensor tuple: N x is the length of X i . Normally, we only have one output.
Generally speaking, following the definition of [38,39], parametric supervised learning can be viewed as a type of a model composed from two parts. The first part is a predictor. Given parameter θ, we have: , whereŷ is the prediction output of learning model M, θ is the parameter of M. As stated, it predicts value based on input X i .

6
The second part is an optimizer, which updates the parameter θ based on a loss function. For a regression model, a typical loss function would be defined as: where · represents 2-norm.
We usually hopes to minimize this loss function. It is formulated by: where M is a learning model and Mˆθ is the optimal solution. Specifically, in this work, we applied Neural Networks [40] and Random Forests [41] in the case studies.

Obtaining rotation-equivariance properties in systems using supervised learning
Group equivariance is an important property for most physical systems. Typical examples of group equivariance could be rotation group equivariance, scaling group equivariance, and translation group equivariance.
Mathematically, group equivariance is a property of a mapping f : X → Y to commute from X to Y under rotation group actions. Specifically, let R ∈ S O(n) be a rotation action. f : As a special case of rotation-equivariant, a function f : Since supervised learning models could be considered as functions, name a machine learning model as M θ .
For a rotation operation R, we hope to obtain the property that: For analysis below in Sec. 3.3, we prove the rotation-equivariance property following the definition stated here in Equ. 2.18. In other words, if a system would satisfy the property in Equ. 2.18, then this system is rotationequivariant.

Modeling symmetric fluid systems via supervised learning
The machine learning approach to the fluid dynamics modeling involves training a supervised learning model M using X i as features and Y as label.
In our case, the underlying space S of the fluid dynamic system is complete with respect to rotation. This means for all rotation R : R n → R n , R(p) ∈ F for all p ∈ S . The objects we want to model via machine learning are rotation-equivariant tensorial fields on S .
Let X be a tensorial field on S , for any point p ∈ S , we use X(p) to denote the tensor at p (for example, pressure at a particular point in a fluid dynamics system). X is said to be rotation-equivariant if for all point p ∈ S and all rotation R X(R(p)) = R(X(p))  Suppose one has tensorial fields X 1 , · · · , X n , Y on S such that X i and Y are related by some unknown physical Supervised machine learning methods can be used here to learn a function M that approximates f such that M generalizes well on new data.
Suppose those tensorial fields are rotation-equivariant, then naturally the model f as well its proxy M

Rotation Equivariant Network
In this section, we would like to propose Rotation Equivariant Network (RotEqNet) to solve rotation problems for high order tensors in fluid systems. RotEqNet is based on the position standardization algorithm, as we would further discuss in Section 3.2. We first provide a general description of the whole architecture in 3.1.

Architecture
As shown in Figure1, RotEqNet generally goes through three important steps: position standardization, prediction of kernel predictor, and position resetting. To be specific, the position standardization is an algorithm to transfer incoming tensor to its standard position. In Figure1, the 'even order standardization' and 'odd order standardization' sections denote this algorithm in position standardization. Then, X s is considered as a standard position of input tensor X, and R is an extracted rotation operation to transfer between standard position and original position. The output of kernel predictor is only dealing with standard positions. This will result the output y s in its standard position as well. Finally, apply R −1 to output y s will be our final prediction. A general mathematical description of this process could be described as: How would this process help to solve rotation problems for high order tensors?
An important reason is related to a reduced function space for learning. When a learning model is only training with the standard position, it would no longer still have to deal with the entire group action causing a group-equivariant, but only need to focus on the pattern by the related physical equation. Name the rotation group as G, and consider a full function space C(X, Y). As mentioned in [31], instead of performing regression on C(X, Y), RotEqNet is essentially exploring a much smaller space C(X/G, Y/G). The reduction of input-output dimensionality makes the training easier. With the same number of samples, the pattern for learning requires a far smaller space to explore. The second reason is RotEqNet could provide a theoretical guarantee of the property of rotation-equivariant. Utilizing rotation symmetry as a strong prior for most physical systems, RotEqNet have a better generalized result learning from limited amount of data.
The following subsections will introduce position standardization algorithm in a complete manner with proof on its property of rotation-invariant in Theorem 3.1. Then, we will demonstrate the proof of showing RotEqNet is rotation-equivariant in Theorem 3.2.

Position standardization algorithm for High Order Tensors
Let D denote our data set. The first stage of RotEqNet is to find a good representative of all tensors that are related to each other by rotation. We will call this representative the sample in "standard position," and we will denote it by (X s , Y s ). We will use S to denote the position standardization algorithm and S (X, Y) = (X s , Y s ) to mean reducing (X, Y) to its standard position.
S has the following property that ∀(X, Y) ∈ D and all rotation operation R ∈ S O(n), This means, S produces exactly the same output no matter how (X, Y) is rotated, i.e. it is rotation-invariant.
Intuitively, for a tensor T , we are selecting a representative on the orbit , as the rotation invariant of a T [42]. In our algorithm, we initially perform a tensor contraction to higherorder tensors, reducing the dimension to obtain a lower order tensor. Then using diagonalization for even cases and QR factorization for odd cases, the algorithm could obtain a rotation operation acting on T . Finally, it could get a tensor in standard position by rotating T the original tensor with the inverse of the obtained rotation matrix.
This operation is compatible with the theoretical result shown in Lemma 2.3.

Tensor of even order
Given a symmetric tensor of even order T n ∈ ⊗ n V (n is even). Let C denote a sequence of contraction along the first two axes until we reach a second-order tensor. Applying C to T n we get: Then we find the orthonormal eigen-vectors of T 2 and use them to form the orthonormal matrix R that diagonalize T 2 Since R is an orthonormal matrix, we have We will call D the standard position of T 2 . We write R(T 2 ) = D to shorten the notation Since contraction and rotation are compatible by Theorem 2.3. We can apply R to T n before we apply contraction, and we will have For the even tensor T n , we define

Tensor of odd order
Why would it be different for even order and odd order? Since odd order cannot directly reduce its dimension to 2 by contraction. Due to the fact that each contraction will reduce the dimensionality by 2, the reduced dimension will also be an odd number, which cannot be 2. Involving in this problem, this would further be impossible to extract the rotation matrix, which is impossible to rotate the tensor into a standard position. The following described is the method that we use to solve the problem.
Given a symmetric tensor of an odd order tensor T n ∈ ⊗ n V (n is odd). Let C denote a sequence of contraction along the first two axes until we reach a third-order tensor. Applying C to T n we get: After we obtain T 3 , we could obtain three different order one tensors by contracting it. Name the contracted results, which are first-order tensors i.e. vectors, V 1 , V 2 , V 3 . We could get an order two tensor by concatenating them.
Then, we could perform the a similar process as before. We perform QR factorization to obtain rotation matrix R. For odd tensor, we define: The pseudocode of our proposed algorithm is documented in Algorithm 1. We evaluate our method of position standardization algorithm in Section 4.

Theoretical Analysis of Rotation-equivariant property
In our analysis, we aim to show the rotation-equivariant property of RotEqNet. As an important first step, we need to analyze the property of rotation invariant (standard position) derived by the position standardization algoriTheorem We hope to show Equ. 3.2 is true. Once Equ. 3.2 is proved true, RotEqNet will automatically be rotation-equivariant based on its architecture.
To sketch an outline about theorems below, Lemma 2.3 would serve as an important fact for preserving rotation information after contraction. Our algorithm has been analyzed by Theorem 3.1. T n = T 3: if T is in even then 4: # For even cases 5: while order(T n ) ë 2 do 6: T n*2 = contraction(T n ô T n*2 , 1, 2) # For odd cases 14: while order(T n ) ë 3 do 15: T n*2 = contraction(T n ô T n*2 , 1, 2, 3)

Theoretical Analysis of Rotation-equivariant property
In our analysis, we aim to show the rotation-equivariant property of RotEqNet. As an important first step, we need to analyze the property of rotation invariant (standard position) derived by the position standardization algorithm. We hope to show Equ. 21 is true. Once Equ. 21 is proved true, RotEqNet will automatically be rotation-equivariant based on its architecture.
To sketch an outline about theorems below, Thm. 1 would serve as an important fact for preserving rotation information after contraction. Our algorithm has been analyzed for even and odd order tensors respectively by Thm.
??, and Thm. ??. It is quite di erent between even and odd order tensors in our algorithm.

Main theorems and proofs
We provide a proof in Appendix D. We call S(T ) the standard position of S(T ) Using Thm. 2 we can automatically have the result on tensor with arbitrary order that the standard position, derived by position standardization algorithm, is rotation invariant.

Theorem 3. RotEqNet is rotation-equivariant.
Proof. Name RotEqNet as M R , kernel classifier as M k . Consider a input pair (X, y) À D. Suppose the result of standardize position algorithm S would have S(X) = P 1 (X), where P 1 denote a rotation operation.
First, consider M R (X) the process of RotEqNet is defined as: Proof. Name RotEqNet as M R , kernel classifier as M k . Consider a input pair (X, y) ∈ D. Suppose the result of standardize position algorithm S would have S (X) = P 1 (X), where P 1 denote a rotation operation.
First, consider M R (X) the process of RotEqNet is defined as: Consider another rotation operation P 2 in the matrix form acting on X, using Theorem. 3.1 we know that: where I is an identity matrix.
Then, consider M R (P 2 (X)) the process of RotEqNet is defined as: We know that S (P 2 (X)) = S (X) from Equ. 3.16. Therefore, we have M R (S (X)) = M R (S (P 2 (X))). Substitute M R (X) back into previous equation, To simplify, for rotation operation P 2 on input X, we have This is showing that M R is rotation-equivariant by definition. Therefore, RotEqNet is rotation-equivariant.
it has shown that Algorithm 1 is able to preserve rotation information for low dimension, and further extract using diagonalization for matrices. This part is a theoretical guarantee of our position standardization algoriTheorem

Case studies
In this section, a series of cases are provided to show the performance of RotEqNet. In the following subsections, cases are included from second-order, third-order, and fourth-order. We specifically investigate second-order cases with detailed studies on linear, and nonlinear test equations since, in current applications, second-order physical systems are widely used. Generally, we reported two properties of RotEqNet in every case study. The first property is a loss reduction property. We apply RotEqNet in each test physical equation and compared it to the baseline models (Neural Networks and Random Forests). Another one is the rotation-invariant property. We examine this property by letting RotEqNet and baseline models to perform prediction on rotations of randomly selected data. We report detailed information for these case studies in every subsection below. The interpretation of experimental results is also included in each subsection. Newtonian fluid is a type of fluid such that its viscous stress changes based on its flow. In this experiment, we aim to use simulation data to demonstrate this rule of Newtonian fluid. This would serve as a case study with second-order linear equations. Let σ ∈ R 3×3 be stress tensor, p ∈ R pressure and S ∈ R 3×3 strain rate. The rule of Newtonian fluid is an second-order physical equation which satisfies the following condition [43]: Another definition of the equation for Newtonian fluid would use the velocity vector field, defined as ∇v. ∇v could be expressed as a 3 × 3 matrix. Using this definition, the equation of Newtonian fluid could also be written as: This could also be considered as the definition of strain rate Based on this definition, we could observe that S = ∇v + ∇v T , and it is symmetric since S = S T . Since S is symmetric and I is an identity matrix, σ is also symmetric. Therefore, defining an arbitrary rotation matrix R, this system is equipped with the property of rotation-equivariant that R(σ) = R(−pI + µS ).
To quantify the stress for Newtonian fluid simulation, it would be useful to be able to predict the Newtonian fluid stress, given the simulation of pressure and velocity vector field. Based on this scenario, in this subsection, we provide a case study for the machine learning model on inputting the shear of Newtonian fluid and prediction of the stress.
To form a proper dataset D with N elements for a machine learning model for Newtonian fluid, the input x is set up to be a vector where x ∈ R 10 . Specifically, x is composed by p and flattened S in Eqn.4.1. The output y ∈ R 9 is a vector which is the flattened result of matrix σ derived by p and S . The dataset D would set up in the description above. To compare the difference of our method to the baseline method, we trained two models with the same hyper-parameter using different amounts of training data, ranging from 10, 000, 20, 000, ..., 100, 000. 85% of generated data is used for training and 15% of data is used for testing. A rotation set with 10,000 random rotation matrices is also generated for evaluating the property of . The machine learning model we apply here is neural networks and random forests because of the ability of these two models to approximate arbitrary functions. For Neural Networks, in our implementation, the logistic activation function is used as an activation function for every neuron. The number of neurons for two layers is 512 and 4, respectively. Adam optimizer [44] is applied as the algorithm for optimization, and the learning is set up to be 1 × 10 −3 . We also set the batch size to be 64. For random forests, 100 estimators are set up with mean squared error as the criterion. The maximum depth of random forests is 3 to lower the chance of overfitting. We used Sklearn for implementation [45]. The computation environment of this experiment is CPU.

Results
There are two properties to evaluate, including error reduction and rotation-equivariant of RotEqNet. The effect of error reduction is evaluated for the first. A kernel predictor is trained by standard positions derived from training data. Then, the prediction algorithm is applied to both training and testing set to obtain the training and testing performances. The validation error E is defined as the Mean Squared Loss using the formulation that: In Eqn  Further, Fig. 4(b) shows the error reduction effect of RotEqNet using Random Forest as a kernel predictor.
Similarly, as shown in Fig. 4(b) with blue curves, it represents the performance of the baseline method (Random Forests). The second experiment group is RotEqNet with Random Forests as the kernel predictor. As shown in Fig. 4(b) in red curves, triangle curve represents training error and the circle curve represents testing error. We could observe a huge error reduction, 99.56% in training and 99.72% in testing, for RotEqNet compared to the error of using only the Random Forest predictor. For the last experiment group, as performances marked as black curves in the figure, trains the kernel predictor with standard position only. As stated before, this could also serve as a reason for the error reduction effect for RotEqNet on random forests.  second-order equations, the application of RotEqNet would result in better-generalized results in learning.
Another important property to evaluate for RotEqNet is rotation-equivariant. The experiment is designed on the definition of rotation-equivariant mentioned in Eqn. 2.18. First, we pick a randomly generated data (X 0 , y 0 ).
Then we apply the rotation set {R i } 10000 i=1 with 10,000 random rotation matrices to (X 0 , y 0 ). To fully investigate the property of rotation-equivariant, we apply an error evaluation method E D here to evaluate the error compared to real data, which is defined as: This error evaluation method (E D ) focuses more on the model's error on real data for all rotations. As shown in Tab. 2, for both baseline methods, using neural networks and random forests, there are large errors for E M and E D .
The baseline methods have no theoretical guarantee that it has the property of rotation-equivariant. However, there is an error reduction for both machine learning models when applying with RotEqNet's architecture. Especially for RotEqNet with Neural Networks as kernel predictor, we could observe that E D = 0.0013 with 99.8% of error reduction. This could guarantee the rotation-equivariant property of RotEqNet.

Problem statement
In this case, we consider the subgrid model of large eddy simulation (LES) of turbulent flow by Kosovic [46].
In this case study, as formulated previously in [47,48], we hope to obtain a learned model by simulation data from LES. This would serve as a case study with second-order non-linear equations. LES is defined as: Here τ i j is subgrid stress, which is a symmetric traceless 2nd order tensor. S i j and Ω i j are symmetric and antisymmetric parts of velocity gradient tensor G i j , where Tr(G) = 0. Further, C s , ∆, C 1 , C 2 are all constants. The configuration of constants above are reported in the next subsection.
In order to qualify the subgrid stress for LES, this study aims to predict the subgrid stress, given the simulation of velocity gradient tensor. This case study for the machine learning model on inputting the velocity gradient tensor.

Data generation and model description
Based on Eqn. 4.5, we first generate random data to obtain a simulated velocity gradient tensor G i j . The generation of random numbers follows a normal distribution from range (0, 1), and G i j is obtained from a random matrix G raw by subtracting 1 3 Tr(G raw ) from diagonal position. This would keep Tr(G) = 0. S i j and Ω i j could be obtained from G i j by getting its symmetric and anti-symmetric parts.  To form a proper dataset D with N elements for a machine learning model for Newtonian fluid, the input x is set up to be a vector where x ∈ R 9 . Specifically, x is composed by flattened G i j . The output y ∈ R 9 is a vector, which is the flattened result of matrix τ derived by G and other constants. To compare the difference of our method to the baseline method, we trained two models with the same hyper-parameter using different amounts of training data, ranging from 10, 000, 20, 000, ..., 100, 000. 85% of generated data is used for training, and 15% of data is used for testing. A rotation set with 10,000 random rotation matrices is also generated for evaluating the property of rotation-equivariant, denoted by {R i } 10000 i=1 . The model setup is the same compared to Sec. 4.1.2.

Results
The effect of error reduction is evaluated for the first. The validation error E is defined as the Mean Squared Loss using the formulation in Eqn. 4.3. This evaluation E represents the expected error of model M with dataset D.    difficulty, and therefore RotEqNet could obtain better performance. Further, Fig. 5(b) shows the error reduction effect of RotEqNet using Random Forest as a kernel predictor. The general performance of using Random Forests as a kernel predictor is relatively worse compared to using Neural Networks as a kernel predictor. In Fig. 5 We name this testing equation as 'testing Newtonian fluid equation' for simplicity. The testing equation revised from Newtonian fluid equation can be described as: where σ ∈ R 3×3×3 is testing stress, p ∈ R is testing pressure, and v ∈ R 3×3×3 is testing velocity field. I ∈ R 3×3×3 is the identity third-order tensor.
Based on this testing equation, we could observe that (∇v + ∇v T ) T = ∇v + ∇v T . Since ∇v + ∇v T is symmetric, and I is symmetric, we have σ is also symmetric. Therefore, defining an arbitrary rotation matrix R, this system is equipped with the property of rotation-equivariant that R(σ) = R(−pI + µ(∇v + ∇v T )).
In order to qualify stress for testing the Newtonian fluid equation, this study aims to predict the stress, given the simulation of pressure and velocity gradient tensor. This case study for the machine learning model on inputting the pressure and velocity gradient tensor.

Data generation and model description
Based on Eqn. 4.6, we first generate random data to obtain ∇v and p. The generation of random numbers in ∇v follows a normal distribution from range (0, 1). σ could be obtained using the Eqn.4.6, derived from generated ∇v and p. Denote the dataset as , To form a proper dataset D with N elements for a machine learning model for Newtonian fluid, the input x is set up to be a vector where x ∈ R 28 . Specifically, x is composed by p and flattened (∇v + ∇v T ) in Eqn.4.6. The output y ∈ R 27 is a vector which is the flattened result of matrix σ. The dataset D would set up in the description above. To compare the difference of our method to the baseline method, we trained two models with the same hyper-parameter using different amounts of training data, ranging from 10, 000, 20, 000, ..., 100, 000. 85% of generated data is used for training and 15% of data is used for testing. A rotation set with 10,000 random rotation matrices is also generated for evaluating the property of rotation-equivariant, denoted by {R i } 10000 i=1 . The model setup is the same as Sec. 4.1.2.   therefore RotEqNet could obtain better performance. Moreover, RotEqNet has good generalization capability without overfitting. As shown in the blue curves of Fig. 6, if we apply raw data in learning directly on baseline models, the testing error is much higher compared to the training error. In this case study, introducing the archi- Further, Fig. 6(b) shows the error reduction effect of RotEqNet using Random Forest as a kernel predictor.

Results
The general performance of using Random Forests as a kernel predictor is relatively worse compared to using Neural Networks as a kernel predictor. In Fig. 5(b), blue curves represent the performance of training with raw data by Random Forests (baseline method); red curves represent the performance of RotEqNet; black curves represent the performance of Random Forest trained with standard positions. For the first point, we could observe an error reduction for 0.90% in training and 6.84% in testing for RotEqNet with Random Forests. As another point, RotEqNet is also obtaining a better-learned model for the model's capability in generalization. The testing error of the baseline method is observably higher than testing error, while the training and testing performance of RotEqNet is approximately the same. As suggested in Fig. 6(a), in second-order nonlinear cases, RotEqNet could reach a generalized learning result with remarkably lower error compared to baseline methods.
To evaluate the rotation-equivariant property of RotEqNet for this third-order case, we designed an experimental process that is close to the one stated in Sec. 4   This case study focuses on a linear relationship of fourth-order tensor. Nye [49] has introduced a fourthorder tensor in modeling elastic compliances and stiffnesses, which has been investigated using machine learning methods [50,51]. Generally, in the study of the properties of a crystalline and anisotropic elastic medium, a fourthorder tensor coefficient will typically be applied to model the relationship between two symmetric second-order tensors [52]. In this case, we study Electrostriction, a property causing all electrical non-conductors to change their shape under the application of an electric field. The relationship is described as: Here T i j ∈ R 3×3 is a symmetric traceless second-order strain tensor. S kl ∈ R 3×3 , S kl = S k S l where V k and V l are first-order electric polarization density. Note here V kl is symmetric. V i jkl ∈ R 3×3×3×3 is the electrostriction coefficient.
Based on the formulation above, this system is symmetric. Since S i j is symmetric, T T i j = (V i jkl S i j ) T = S kl V kli j = T i j . This could guarantee that T i j is also symmetric. Due to the face that the system is symmetric, applying a random rotation matrix R, R(T ) = R(VS ).
In order to qualify strain for study on Electrostriction, we aim to predict the strain, given the simulation of electrostriction coefficient and electric polarization density.

Data generation and model description
Based on Eqn. 4.7, we first generate random data to obtain simulated electrostriction coefficient tensor V i jkl and electric polarization density tensor S i j . The generation of random numbers follows a normal distribution. rotation-equivariant, denoted by {R i } 10000 i=1 . The model setup is the same compared to Sec. 4.1.2. We use Numpy to generate this simulated dataset by generating a random symmetric fourth-order tensor V, and second-order tensor S . T is computed from generated V and S by Eqn. 4.7.

Results
The effect of error reduction is evaluated for the first. The validation error E is defined as the Mean Squared Loss using the formulation in Eqn. 4.3. This evaluation E represents the expected error of model M with dataset D. Fig. 7 shows the performance of Neural Networks and Random Forests as kernel predictor separately. It is observable that in high-order cases, Neural Networks have huge superiority to Random Forests. Details will be demonstrated in the following paragraphs.
We are starting with Neural Networks, Fig. 7 To interpret the experimental results, firstly, RotEqNet could reach a better learning performance compared to simply applying Neural Networks (baseline method). A dataset with only standard positions has lower training difficulty compared to random positions. This claim is supported by black curves in Fig. 7(a)   Further, Fig. 7(b) shows the error reduction effect of RotEqNet using Random Forest as a kernel predictor.
At first glance, we could find that the curves for Random Forests are quite messy without certain patterns like Fig. 7(a). The general performance of using Random Forests as a kernel predictor is worse in both aspects of performance and generalization. In Tab. 7, we could observe a training error reduction for 0.58% and testing error reduction of 2.96%. Even if we could still see the general error of RotEqNet seems to be slightly lower than the baseline method. This result is not comparable to the error reduction performance with setting Neural Networks as a kernel predictor. As another point, selecting Random Forests as a kernel predictor fails to extract learning rules with the standard position. As we could observe the black curves in Fig. 7 Forest.
The large error reduction observed in case studies raised new opportunities in solving the problem of the physical system with rotation symmetry. Most physical systems have the property of rotation symmetry, and currently, there exist few works that could provide a theoretical guarantee to this property for machine learning methods. A key point in this problem is to design a properly defined algorithm to obtain rotation invariant for high-order tensors. This paper has shown RotEqNet with theoretical and experimental results aiming to solve the problem of rotation symmetry.
We first define a standard position as rotation invariant, which is compatible for high-order tensors. It allows us to extract the rotation invariant of high-order tensors using a contraction, diagonalization, and QR factorization.
The theoretical guarantee is shown in Thm. 3  Rotation-equivariant, this would cause RotEqNet could process this system with any rotation.
The two reasons above are the main reasons that are causing the error reduction for RotEqNet. There is also one point to mention is the selection of kernel predictor. The model selection of kernel predictor will affect the learning results significantly since the kernel predictor is essential in learning the physical system without the effect of rotation symmetry. Neural Networks is the best model in the design of the data-driven method for physical systems because of its flexibility to approximate arbitrary functions. We only reported the performance of Neural Networks and Random Forests as previous work by Ling [16]. As described in Sec. 4.4.3, the performance of Random Forests is limited compared to Neural Networks. Also, as a general trend in previous experiments, Neural Networks are usually reaching better performance compared to Random Forests. As a quick conclusion, we believe the application of Neural Networks as a kernel predictor has a series of advantages than other machine learning models.
We wish to further discuss about another error evaluation method of rotation-equivariance property that we do not mention in case studies. Consider a type of error evaluation, evaluating rotation error of model itself, the error 24 E M is defined as: The evaluation of this error is actually trivial since we have already rigorously provided a proof in Theorem 3.2 showing the rotation-equivariance property of RotEqNet. We applied this evaluation in first two case studies, and the estimated error is around 10 −16 for all these cases.
For future work, there are three directions to this paper: a better definition of standard position, application to other groups, and generalization to non-symmetric systems. For the first direction, for the current definition, the rotation invariant of odd-order tensors is not reaching equivalent performance as even-order tensors. It would be a good work for revising the definition of standard position for odd tensors. Second, besides rotation symmetries, there are also physical systems with other group-equivariant properties such as scaling and transaction. This work could provide a method in solving problems with other groups, but the detailed design of an algorithm should differ from case to case. Third, current work could only deal with the symmetric system. However, for a general case, if the system is not symmetric, there are certain methods to use RotEqNet in a symmetric system for solving a non-symmetric system. A good trick to consider, for example, is to deal with PP T , where P is a matrix. This is a great intuition to extend our current work into non-symmetric physical systems.  Proof. Since both C(a, b) and g are linear, we may assume that T is of the form v i1 ⊗ · · · ⊗ v in .

Proof of Theorem 1
Proof. Since the position standardization algorithm defines standard position differently for even and odd orders.
We show our proof on even and odd cases separately.
Suppose T has even order.
Let C be the sequence of contraction along the first two axes such that C(T ) = T 2 , where T 2 is a second-order tensor as described in the algorithm.
Given arbitrary even high order tensor T , we could perform contraction to a second order tensor T 2 via first two indices: For T 2 , using Lemma 2.1, there exists a rotation R such that: where R(T 2 s ) = RT 2 s R t . T 2 is diagonalizable because it is symmetric. Since R is represented by a orthonormal matrix, therefore R t = R −1 .
Based on Lemma 2.3, we know rotation commutes with contraction. Therefore, based on S the standard position S (T ) is defined as Consider a rotation operation P in its matrix form. When we act P on T we obtain a new tensor P(T ). For this new tensor, applying contraction we could have: To simplify, for a rotation operation P acting on an even high order tensor T ,

S (P(T )) = S (T ). (A-21)
This satisfy the definition of rotation invariant. Therefore, for even cases, the standard position S (T ) is a rotation invariant.
Suppose T has odd order.
Let C be the sequence of contraction along the first two axes such that C(T ) = T 3 , where T 3 is a third-order tensor as described in the algorithm. Based on S , we have In this case, Consider any rotation operation P acting on T n . We have, P(T 3 ) = P(C(T n )) (A-26) Using QR-factorization, The standard position of P(T n ) will be defined as: S (P(T n )) = R −1 2 (P(T n )) (A-28)