Real Harmonic Analysis on the Special Orthogonal Group

This paper presents theoretical analysis and software implementation for real harmonics analysis on the special orthogonal group. Noncommutative harmonic analysis for complex-valued functions on the special orthogonal group has been studied extensively. However, it is customary to treat real harmonic analysis as a special case of complex harmonic analysis, and there have been limited results developed specifically for real-valued functions. Here, we develop a set of explicit formulas for real-valued irreducible unitary representations on the special orthogonal group, and provide several operational properties, such as derivatives, sampling, and Clebsch-Gordon coefficients. Furthermore, we implement both of complex and real harmonics analysis on the special orthogonal group into an open source software package that utilizes parallel processing through the OpenMP library. The efficacy of the presented results are illustrated by benchmark studies and an application to spherical shape matching.

representations weighted by Fourier parameters. In particular, harmonic analysis on the special orthogonal group, or the rotation group, has been utilized in quantum physics [21,2]. Recently, it also has been applied to various engineering problems in robotics, controls, and machine learning [4,14,5]. Computationally efficient numerical implementation and fast Fourier transform algorithms on the special orthogonal group are considered in [17,13].
Despite extensive prior works in broad areas of science and engineering, all of the aforementioned references deal with complex harmonic analysis, where matrix representations and Fourier-parameters are complex-valued. There have been limited efforts in noncommutative harmonic analysis for real-valued functions on the special orthogonal group. One of the earlier studies in real harmonic analysis include [18], where real-valued matrix representations for several atomic orbitals are evaluated up to a certain order. Later, in [11,12], recurrence relations to evaluate real matrix representations are presented in terms of nine elements of a rotation matrix. In [3], real-valued matrix representations are constructed by transforming complex-valued matrix representations, namely wigner D-matrices.
In this paper, we develop a new form of matrix representations for real harmonic analysis on the special orthogonal group. Compared with [11], these are formulated in terms of Euler angles so that a real fast Fourier transform on the special orthogonal group can be developed utilizing various discrete fast Fourier transform algorithms in the Euclidean space. The presented form is based on the wigner d-function of the second Euler angle, which can be computed by a recurrence relation. But, once the wigner d-function is evaluated, it is explicit with respect to the remaining two other Euler angles. Compared with [3], the presented real-valued matrix representation can be constructed without need for evaluating complex-valued matrix representation, or the wigner D-matrices.
Furthermore, we present several operational properties of real harmonic analysis. First, a sample theorem is presented such that Fourier transform of a bandlimited function with a bandwidth B can be exactly computed with (2B) 3 function evaluations. As such, there is no need to approximate an integration over the special orthogonal group with a quadrature rule, when computing Fourier transforms. This is utilized to construct a fast Fourier transform for real-valued functions on the special orthogonal group. Next, we develop Clebsch-Gordon coefficients for real harmonics so that a product of two real matrix representations is written as a linear combination of other real matrix representations. Finally, we show an explicit expression for the derivatives of real matrix representations to formulate representations for lie algebra.
Beyond these theoretical contributions, an open source software package has been developed for real harmonic analysis on the special orthogonal group. This library includes fast Fourier transform, inverse Fourier transform, and evaluation of Clebsch-Gordon coefficients and derivatives. While this paper focuses on real harmonic analysis on the special orthogonal group, this software library also provides complex harmonic analysis, including wigner D-functions, and spherical harmonics on the unit-sphere as well. There are developed in c++ while utilizing the OpenMP library [6] to support accelerated parallel computing for multithread processors. These are verified by various software unit-testing. A benchmark study and a particular application to spherical shape matching for Earth topological data are presented as well.
In short, the main contribution of this paper is theoretical and numerical analysis for the foundation of real harmonic analysis on the special orthogonal group. The presented form of real matrix representations and its derivatives, and real Clebsch-Gordon coefficients have been unprecedented. Also, the proposed software library can be utilized in various engineering application of noncommutative harmonic analysis on the special orthogonal group.

Complex Harmonic Analysis on SO(3)
The three-dimensional special orthogonal group is composed of 3 × 3 orthogonal matrices with determinant one, i.e., (1) Harmonic analysis for complex-valued functions on SO(3) has been studied extensively, for example in [4,21,15]. In this section, we summarize selected results that are used for the subsequent developments of real harmonics on SO(3).

Euler Angles
Any R ∈ SO(3) can be parameterized by three angles α, γ ∈ [0, 2π) and where e 2 = (0, 1, 0), e 3 = (0, 0, 1) ∈ R 3 , and the hat map· : R 3 → so (3) is defined such that x × y =xy for any x, y ∈ R 3 . The Lie algebra is defined as These are referred to as 3-2-3 Euler angles. While this parameterization involves complicated combinations of trigonometric functions and inherent singularities in the resulting kinematics equation, there are several advantages in harmonic analysis, such as convenience in applying fast Fourier transform techniques.

Irreducible Unitary Representation: Wigner D-Matrix
As a transformation group, SO(3) acts on R 3 while preserving its inherent structures of R 3 . For example, Rx ∈ R 3 for any R ∈ SO(3) and x ∈ R 3 via the matrix multiplication. Let L 2 (R 3 ) be the set of complex-valued, square-integrable functions on R 3 . The left regular representation on SO(3), namely D(R) for R ∈ SO (3) is an operator D(R) : for any f ∈ L 2 (R 3 ) and x ∈ R 3 , i.e., it transforms a function such that it becomes equivalent to rotating the input argument by R T . It is straightforward to show that Also, it is a linear operator on L 2 (R 3 ). Consequently, by selecting a basis on L 2 (R 3 ), D(R) can be represented with a matrix, thereby resulting in matrix representations.
The matrix representation on SO(3) is often induced from that of SU(2) = {U ∈ C 2×2 | U * U = I 2×2 , det[U] = 1} utilizing one-to-two correspondence between SO(3) and SU (2). More specifically, consider representation on SU(2) as an operator on the analytic functions from C 2 to C. By selecting the set of homogeneous polynomials as the basis of the analytic functions, one can derive matrix representation of SU(2) as which is referred to as the wigner D-matrix that is common in quantum mechanics [21]. Here the range of the index l is {0, 1 2 , 1, 3 2 . . .}, and for each l, the indices m, n vary in {−l, −l + 1, . . . , l − 1, l}. In (4), the real-valued d l m,n is called wigner d-matrix. An explicit form of d l m,n is tabulated up to l = 5 in [21], and a recursive algorithm to evaluate it for arbitrary order is presented in [3]. The above expression for matrix representation SU(2) results in a matrix representation of SO(3) by taking the integer values of l only. Throughout this paper, D l (R) ∈ C (2l+1)×(2l+1) is considered as a square matrix where the row index and the column index are m and n, respectively, which vary from −l to l in ascending order.
From (3), it is trivial to show D l (I 3×3 ) = I 2l+1×2l+1 , or D l m,n (I 3×3 ) = δm,n. This representation is irreducible, i.e., cannot be block diagonalized consistently with a similarity transform, and it is unitary, i.e., where the last equality is from the homomorphism property. When α = γ = 0, these also imply While this paper follows the convention of 3-2-3 Euler angles, using other types of Euler angles in fact does not matter as it would yield an equivalent form of matrix representations that can be constructed by a similarity transform of (4).

(8)
Therefore, the Fourier parameters in (6) can be discovered by the inner product, Equations (9) and (6) are considered as the Fourier transform of f (R), and its inverse transform, respectively. A function f ∈ L(SO(3)) is band-limited with the band B, if its Fourier parameters vanish, i.e., F l m,n = 0 for any l ≥ B. The classical sampling theorem states that the Fourier transform of a band-limited function can be recovered from the sample values that are chosen at a uniform grid with a frequency that is at least twice of the band-limit. Using this, a fast Fourier transform technique is presented in [13].

Clebsch-Gordon coefficients
The product of two wigner D-matrices can be rewritten as a finite linear combination of other wigner D-matrices as follows.
Interestingly, the coupling coefficients are split into two parts as written above, and they are referred to as Clebsch-Gordon coefficients. There are several properties of Clebsch-Gordon coefficients listed in [21], including Using the last property, the double summation at (14) reduces to for l = max{|l 1 − l 2 |, |m 1 + m 2 |, |n 1 + n 2 |}. While (14) commonly appears in the literature, (19) is not available at least in the cited references. A computational scheme to evaluate Clebsch-Gordon coefficients is proposed in [20].
In [15], it is proposed to rearranged the coefficients into a matrix C l1,l2 ∈ R (2l1+1)(2l2+1)×(2l1+1)(2l2+1) according to the following ordering rules: which begins from 0 following the convention of the C programming language that is adopted for software implementation in this paper. Under this matrix formulation, (14) is rearranged into where ⊗ denotes the Kronecker product.

Relation to Spherical Harmonics
Consider the unit-sphere, with the associated Legendre polynomials P m l for l ∈ {0, 1, . . .} and −l ≤ m ≤ l. Let dx = 1 4π sin θdφdθ be the measure of S 2 , normalized such that S 2 dx = 1. We define an inner product on the square-integrable functions, namely L 2 (S 2 ) as Spherical harmonics satisfies the following orthogonality with respect to the above inner product, yielding Spherical harmonics are closely related to the wigner D-function. Using the homomorphism property of the group representation [4], we obtain Therefore, using (24), the wigner D-function can be rediscovered from the spherical harmonics as Let Y l be the (2l + 1) × 1 column vector whose elements are composed of Y l m for m ∈ {−l, . . . l} in ascending order. The above equation is rewritten in a matrix form as

Real Harmonic Analysis on SO(3)
The objective of this paper is to develop the counterparts of Section 2 for realvalued functions on SO (3), resulting in real harmonic analysis on SO(3). Instead of formulating as a specialized form of wigner D-matrices, real matrix representations are directly constructed in terms of Euler angles, and they are utilized for fast Fourier transform of real-valued functions on SO(3). Further, Clebsch-Gordon coefficients and derivatives are formulated as well.

Real Irreducible Unitary Representations
We follow the approaches presented in [3], where real harmonics on SO(3) is constructed from real spherical harmonics. Orthogonal basis for real-valued functions on S 2 , namely S l (x) ∈ R 2l+1 is constructed by the following transform, where x ∈ S 2 and the matrix T l ∈ C (2l+1)×(2l+1) is defined as or in a matrix form, It is straightforward to show that the columns of T l are mutually orthonormal, i.e., T l is unitary so that Let U l ∈ R (2l+1)×(2l+1) be the matrix for the l-th real harmonics on SO(3). Motivated by (27), it is defined as Substituting (28) and rearranging, This is a homomorphism, i.e., where the last equality is from the homomorphism property and U l (I 3×3 ) = I (2l+1)×(2l+1) . Also, similar with (25) and (8), While U l (R) can be evaluated by transforming D l (R) according to (33) as in [3], the procedure will involve unnecessary steps with complex variables. Here, we present an alternative, explicit formulation as follows.
Theorem 1 Real harmonics on SO(3) defined in (33) is equivalent to the following formulations.
satisfying Ψ l m,n = Ψ l m,−n . Or in a matrix formulation, From (29), the expression in the summation vanishes when |p| = |m| or |q| = |n|. Consequently, the above double summation reduces to the following cases, depending on whether any sub index is zero or not, This theorem states that real matrix representation on the special orthogonal group can be constructed directly in terms of Euler angles without need for evaluating complex, wigner D-matrices. The expression presented in (36) is composed of sine and cosine terms for multiples of α, γ, that are similar to those appear in real spherical harmonics, and Ψ terms defined by the wigner d-matrices. When written in a matrix form as (38), it is given by a product of three terms, where each term depends on one of Euler-angles. In contrast to the recursive formulation written in terms of elements of a rotation matrix [11], this structure is useful for a fast Fourier transform algorithm.
In particular, when l = 1, (38) results in In other words, the real matrix representation of the order l = 1 is similar to the rotation matrix itself.

Fourier Transform on SO(3)
From Peter-Weyl theorem, real harmonics yield a complete, orthogonal basis on the square-integrable, real-valued functions on SO(3). More explicitly, any realvalued f ∈ L 2 (SO (3)) is expanded as for real-valued Fourier parameters F l m,n ∈ R. From (35), the Fourier parameters are obtained by Next, we present a sampling theorem to evaluate (42) exactly with a finite number of samples.  w k U l m,n (R(α j1 , β k , γ j2 ))f(R(α j1 , β k , γ j2 )), where the weighting parameter w k ∈ R is defined as Proof Define a sampling distribution as s(R(α, β, γ)) = 2B−1 j1,k,j2=0 which is the linear combination of grid points weighted by the parameter w k . Since w k U l m,n (R(α j1 , β k , γ j2 )).
For the selected grid, it is straightforward to show 2B−1 j1=0 sin mα j1 = 0, 2B−1 j1=0 cos mα j1 = 2Bδ m,0 . Therefore, from (37), S l m,n = 4B 2 δ m,0 δ n,0 In [7], it is shown that the selected weight satisfies Next, define g(R) = f (R)s(R). From (46) and using the property of the delta function, it is straightforward to show that the Fourier parameters for g(R) is given by On the other hand, using (47), Now, we show that the above expression for g(R) yields the Fourier parameters that are identical to those of f (R). Since f (R) can be expanded as a linear combination of U l1 for 0 ≤ l 1 ≤ B − 1. The last term of the above equation is expanded by the product U l1 U l2 with 0 ≤ l 1 ≤ B − 1 and 2B ≤ l 2 . According to the Clebsch-Gordon theorem, U l1 U l2 is a linear combination of U l3 for |l 1 − l 2 | ≤ l 3 ≤ l 1 + l 2 . We have min |l 1 − l 2 | = 2B − 1 − B = B + 1. As such, f (R) and g(R) share the Fourier coefficients in the given band limit, i.e., F l m,n = G l m,n for l ∈ {0, . . . B − 1}. This shows (44).
Therefore, Fourier parameters of any band-limited function with the bandwidth B can be computed exactly with (2B) 3 samples evaluated at the given grid points. Utilizing this, we present a fast Fourier transform.

Substituting (36) into (44), the Fourier transform represented by (44) can be executed in the following sequence. For each
f (R(α j1 , β k , γ j2 )) cos(mα j1 + nγ j2 ), which can be computed by a real-valued fast Fourier transform algorithm developed in R, such as [19]. The Fourier parameters defined in (36) are evaluated by

Derivatives of Representation
Similar with (10), the derivatives of U l (R) at R = I 3×3 results in the real matrix representation of so(3). Here we use the same notation as the complex valued case as for η ∈ R 3 . From the linearity, u l (η) for any η can be evaluated by the results of the following theorem.

Clebsch-Gordon Coefficients
Here we find the Clebsch-Gordon coefficients for real harmonics. The objective is to write a product of two real harminics as a linear combination of other harmonics. Repeatedly using the property of Kronecker product, namely (AC) ⊗ (BD) = (A ⊗ B)(C ⊗ D) for arbitrary compatible matrices A, B, C, D, (33) results in

Substituting (22), this is rearranged into
where the Clebsch-Gordon matrix for the real harmonics, namely c l1,l2 ∈ C (2l1+1)(2l2+1)×(2l1+1)(2l2+1) is defined as Interestingly, while the Clebsch-Gordon coefficients C l1,l2 for the complex harmonics are real-valued, those for real harmonics can be complex-valued as the matrix T l is composed of real or imaginary elements. In the element-wise form, The following theorem presents two properties of the real Clebsch-Gordon coefficients, and provides an alternative, simpler method to evaluate (54).  Table 1.
(58) where As stated above, the Clebsch-Gordon coefficients for real harmonics is complex in general. However, by investigating Table 1, one can show that the product c l,m l1,m1,l2,m2 c l,n l1,n1,l2,n2 in (58) is real always. This is expected as U l m,n (R) is realvalued always.

Software Implementation
The presented real harmonic analysis on SO(3) has been implemented by C++, and the resulting software package has been disclosed as an Open Source Library at https://github.com/fdcl-gwu/FFTSO3. This provides the following functionality and features:

Validation
The software package is validated by software unit-testing implemented via Google Test [8]. There are several tests to verify that numerical results are consistent with theoretical analysis. Here we show the results of a particular unit-testing designed to ensure that the composition of the inverse transform and the forward transform yields the identity map on the space of Fourier parameters [13].
Specifically, we define a function f (R) : SO(3) → R as the inverse Fourier transform presented in (41). We assume that the function is band-limited so that new set of Fourier parameters, namely G l m,n . Theoretically, G l m,n should be identical to the original Fourier parameters F l m,n of f (R), as the composition of the inverse transform and the forward transform is the identity map on the space of Fourier parameters. Numerical errors of the presented software library are measured by the discrepancy between F l m,n and G l m,n calculated as the sum of the matrix norm, Table 2 summarizes the above error for varying band-limits. The error increases as the band-limit, since the error measure is defined as the sum over the index l. However, all of the errors are under an acceptable level relative to the machine precision of double-type variables that the presented software uses.  Figure 1 illustrates the computation time with respect to the band-with for several number of threads. The computation time increases as B is increased, but it reduces as more threads are available. More specifically, the next subfigure shows the speed-up factor, the ratio of the computation time with a single thread to that of multiple threads. Ideally, the speed-up factor is identical to the number of threads, as indicated by the dotted line. However, due to the overhead in distributing the tasks in multiple threads and collecting results, as well as delay in accessing a shared memory, the speed-up factor is often lower than the number of threads in practice. It is illustrated that the speed-up factor increases as B is increased, i.e., parallel processing is more beneficial if the bandwidth is greater. Figure 2 is for the inverse transform. Compared with the forward transform, only a fraction of time is required. As such, the speed-up factor is relatively low even for higher bandwidth. At least, they are greater than one, indicating that the computation time is reduced by parallel processing.

Spherical Shape Matching
We apply the above software package for spherical shape correlation and matching. Consider an object embedded in R 3 . Assuming that it has zero genus, i.e., no hole, it's shape can be completely described by deforming a sphere. Let f (x) : S 2 → R be the deformed radius along the radial direction specified by the unit-vector x ∈ S 2 . Suppose the object is rotated by a rotation described by R ∈ SO(3), and let g(x) : S 2 → R describe the rotated object. The spherical shape matching problem considered in this paper is to find the rotation matrix R ∈ SO(3) for given shape functions f (x) and g(x). This can be formulated as the following optimization problem on SO(3). For a rotation matrix R, let its cost be the discrepancy between g(x) and f (x) rotated by R, i.e., where we have used the fact f (R T x) = f (x) . The true rotation matrix R minimizes the cost function J (R), or equivalently, it maximizes the measure of correlation defined by C(R) = g(x), f (R T x) L(S 2 ) , i.e., R = arg max {C(R)}.
However, evaluating the above inner product requires an integration over S 2 , and it may be cumbersome to repeat it at every iteration of numerical optimization. Instead, the cost function can be evaluated without any integration, by utilizing the fundamental property of representation given by (3). Let the Fourier transform of f (x) and g(x) with a bandwidth B be The orthogonality of real spherical harmonics implies S l1 (x), (S l2 (x)) T L(S 2 ) = 1 4π δ l1,l2 I (2l1+1)×(2l1+1) . Therefore, the correlation function reduces to which is an algebraic equation of Fourier parameters that does not require any integration. Furthermore, its gradient can be evaluated by (49)-(51) as follow. For η = (η 1 , η 2 , η 3 ) ∈ R 3 , a directional derivatives of the correlation function is given by d dǫ ǫ=0 C(R exp(ǫη))) = ∇C(R) · η, where the gradient is considered as a 3 × 1 vector, i.e., ∇C(R) ∈ R 3 , whose i-th element is given by for i ∈ {1, 2, 3}. In the above expression, we have used the homomorphism property U l (R exp(ǫη)) = U l (R)U l (exp(ǫη)). Using these, one can apply any gradient-based numerical optimization algorithm, such as one summarized in Table 3. Table 3 Spherical shape matching algorithm 1: procedure Spherical shape matching

Earth Elevation Map
Here we consider a particular example with a worldwide elevation model, namely ETOPO5 [1]. The topological elevation model is interpolated such that for a given set of latitude and longitude the corresponding elevation is calculated. This function yields f (x), which is transformed by spherical harmonics with the bandwidth of B = 129 to obtain F l . We choose a particular rotation matrix with 3-2-3 Euler angle (α, β, γ) = ( π 6 , π 3 , π 4 ), and perform the Fourier transform of f (R T x) to obtain G l . See Figure 3 for the original elevation map, and the rotated one.
The optimization algorithm summarized in Table 3 is implemented with the initial guess (α, β, γ) = (0.3, 0.3, 0.3), a tolerance ǫ = 10 −6 and step-size δ = 5 × 10 −3 . The numerical optimization is terminated after 223 iterations, with the maximum absolute error of 7.98 × 10 −5 in terms of Euler angles. The evolution of the correlation function and its gradient magnitude over iterations, and the change of the rotation matrix in terms of Euler angles are illustrated in Figure 4.  We have presented real harmonic analysis on the special orthogonal group, including various operational properties such as fast Fourier transform, Clebsch-Gordon coefficients, and derivatives. There are further implemented into an open source software package supporting parallel processing for accelerated computing.
In particular, the presented form of real irreducible, unitary representations can be evaluated without constructing their counter parts in complex harmonic analysis, namely wigner D-matrices, and the given formulation in terms of Euler angles are suitable for fast Fourier transform. Future works include extension beyond the special orthogonal group, such as the special Euclidean group for various engineering applications regarding the coupled translational and rotational dynamics of a rigid body. Also, the intermediate term Ψ l l,m (β) in (37) may be directly evaluated without formulating real-valued wigner d-functions.