A Fast Algorithm for the Moments of Bingham Distribution
AA Fast Algorithm for the Moments of Bingham Distribution
Yixiang Luo , Jie Xu and Pingwen Zhang Courant Institute of Mathematical Sciences, New York University, New York 10012, USA Department of Mathematics, Purdue University, West Lafayette 47907, USA LMAM & School of Mathematical Sciences, Peking University, Beijing 100871, ChinaEmail: [email protected], [email protected], [email protected]
December 6, 2016
Abstract
We propose a fast algorithm for evaluating the moments of Bingham distribution.The calculation is done by piecewise rational approximation, where interpolation andGaussian integrals are utilized. Numerical test shows that the algorithm reaches themaximal absolute error less than 5 × − remarkably faster than adaptive numericalquadrature. We apply the algorithm to a model for liquid crystals with the Binghamdistribution to examine the defect patterns of rod-like molecules confined in a sphere,and find a different pattern from the Landau-de Gennes theory. Keywords : Bingham distribution, directional data, piecewise rational approxima-tion, liquid crystals.
The Bingham distribution is an important antipodally symmetric distribution on the unitsphere S . Although introduced from a statistical perspective [4], it has found applicationsin liquid crystals [7, 6, 3, 8], palaeomagnetism [18, 12, 11], and various other fields involvingdata on the sphere [5, 1, 15, 17, 21].The density function of the Bingham distribution is given by f ( x | B ) = exp (cid:88) i,j =1 B ij x i x j (cid:44)(cid:90) S exp (cid:88) i,j =1 B ij x i x j d x , x ∈ S , (1)where B is a 3 × (cid:104) x n x n x n (cid:105) = (cid:90) S f ( x | B ) x n x n x n d x . (2)Denote Z n n n ( B ) = (cid:90) S x n x n x n exp (cid:88) i,j =1 B ij x i x j d x . (3)1 a r X i v : . [ m a t h . NA ] D ec hen the moments can be expressed as (cid:104) x n x n x n (cid:105) = Z n n n ( B ) /Z ( B ).Even when solving a single problem, the evaluation of moments (2) may need to be donerepeatedly. This is a typical case in the simulations of liquid crystals. In each iterationor time step, (2) is computed at each grid point. Generally speaking, the number of spacediscretization is O ( N ). If we calculate (2) by direct numerical quadrature, it costs O ( N )operations for every single calculations, leading to a total cost of O ( N ). On the otherhand, it should be noted that the density function (1) is determined only by B , not relevantto parameters (and domains, etc.) specified by the problem to be solved. Therefore, it isdesirable to have a fast algorithm for the evaluation of (2).The existing approximations of (2) are designed only for special cases, and are not accu-rate enough to meet the demand of simulations in many problems. Kent [10] proposed simpleexpansions for the zeroth and second moments. The relative error is about 0.1%. Kume andWood [14, 13] developed a method to compute the Z ( B ) by using saddle-point approx-imation. It is accurate for the final estimation result when applying this method in doingmaximum likelihood estimation, but not accurate enough for evaluating Z ( B ). Moreover,the approximation cannot be easily extended to general Z n n n ( B ). Wang et. al. [20] usedpiecewise linear interpolation to compute B from Z n n n /Z where n + n + n = 2. Thisapproach works well for B not far from zero matrix, but is inaccurate when it is not the case.We also mention that in [7] the fourth-order moments Z n n n /Z, ( n + n + n = 4), areapproximated by polynomials of the second-order moments Z n n n , ( n + n + n = 2), witha relative error of 5 × − . This approach is restricted to the cases where B is not involvedexplicitly.In this paper, we introduce a fast and accurate algorithm for evaluating Z n n n ( B ). Wedivide B into three cases and use different approximation method for each case. The maintechniques we utilize are interpolation and Gaussian integrals. We have implemented themethod for n + n + n ≤ BinghamMoments . It is freely availableonline [16], in which pre-calculations are done and saved as constants in the routine to raisethe real-time efficiency. The cost of evaluating Z n n n is reduced to O (1) compared with O ( N ) in numerical integration. Numerical experiments show that the absolute error is lessthan 5 × − in the routine, while 10 times faster than adaptive numerical quadrature withthe same accuracy. We apply the method to a liquid crystal model proposed in [3, 8]. Themodel substitutes the polynomial bulk energy in the widely-used Landau-de Gennes theorywith the entropy term expressed by the Bingham distribution. By this substitution the orderparameters are confined in the physical range, and it is shown in [8] that this model canbe derived from molecular theory. We examine the defect patterns for rod-like moleculesconfined in a sphere, and find a different structure from the Landau-de Gennes theory. Therest of paper is organized as follows. In Sec. 2, we present the approximation method. Thenumerical accuracy is examined in Sec. 3. An application to liquid crystals is given in Sec.4. Concluding remarks are stated in Sec. 5. We diagonalize B using an orthogonal matrix T with det T = 1, B = T diag( b , b , b ) T T . f ( x | B ) = exp (cid:32) (cid:88) i =1 b i ( T T x ) i (cid:33)(cid:44)(cid:90) S exp (cid:32) (cid:88) i =1 b i ( T T x ) i (cid:33) d x . (4)Thus, by the transformation x −→ T T x , Z n n n ( B ) = (cid:90) S x n x n x n exp (cid:32) (cid:88) i =1 b i ( T T x ) i (cid:33) d x = (cid:90) S ( T x ) n ( T x ) n ( T x ) n exp (cid:32) (cid:88) i =1 b i x i (cid:33) d x (5)becomes a linear combination of Z m m m (diag( b , b , b )). Furthermore, the distribution f ( x | diag( b , b , b )) is invariant under changes ( b , b , b ) → ( b + h, b + h, b + h ) for anyreal number h . Without loss of generality, we assume that b ≤ b ≤ b = 0. Denote Z n n n ( b , b ) = Z n n n (diag( b , b , Z n n n ( b , b ) is nonzeroonly if n i are even numbers. Then by x = 1 − x − x , we can express Z n n n ( b , b )linearly by Z nm ( b , b ). Hence it suffices to compute Z nm ( b , b ), denoted in abbreviate by Z nm ( b , b ).Choosing a parameter d >
0, we divide ( b , b ) ∈ ( −∞ , into three regions,( −∞ , − d ] , ( −∞ , − d ] × ( − d, ∪ ( − d, × ( −∞ , − d ] , ( − d, . and use different approximation method for each region. The following Gaussian integral isused in the approximation, (cid:90) R x n exp( − αx )d x = (cid:114) πα (2 n − α ) n , α > . (6) b , b ≤ − d We transform the integral domain into the unit circle, Z nm ( b , b ) =2 (cid:90) (cid:90) x + x < x n x m · exp (cid:0) b x + b x (cid:1) · (cid:112) − x − x d x d x , =2 (cid:88) j,k ≥ (cid:18) j + kj (cid:19) (2 j + 2 k − j + 2 k )!! (cid:90) (cid:90) x + x < x j + n x k + m exp (cid:0) b x + b x (cid:1) d x d x . (7)The series converges because b , b <
0. We truncate the series at j + k ≤ N . Moreover, if d is large, then x j + n x k + m increases with polynomial rate, while exp (cid:0) b x + b x (cid:1) decreaseswith exponential rate. Thus we expand the integral domain to R in the truncated series,which yields the following approximation formula,ˆ Z nm ( b , b ) =2 (cid:88) j + k ≤ N (cid:18) j + kj (cid:19) (2 j + 2 k − j + 2 k )!! (cid:90) (cid:90) R x j + n x k + m exp (cid:0) b x + b x (cid:1) d x d x . = (cid:88) j + k ≤ N (cid:18) j + kj (cid:19) (2 j + 2 k − j + 2 k )!! (cid:115) π b b (2 j + n − k + m − b ) j + n/ (2 b ) k + m/ . (8)3 .2 b > − d , b ≤ − d or b ≤ − d , b > − d We explain our approximation method by the case b ≤ − d , b > − d . Rewrite Z nm ( b , b ) as Z nm ( b , b ) = 4 (cid:90) − x n · exp (cid:0) b x (cid:1) · g m ( b , x ) d x , (9)where g m ( b , x ) = (cid:90) √ − x x m · exp (cid:0) b x (cid:1) · (cid:112) − x − x d x . (10)Denote a = 1 − x and r = x /a , then we have g m = (cid:90) √ a x m · exp (cid:0) b x (cid:1) · (cid:112) a − x d x = a m/ (cid:90) r m · exp (cid:0) b ar (cid:1) · √ − r d r = a m/ · √ π · Γ[( m + 1) / m + 2) / · F (cid:18) m + 12 ; m + 22 ; b a (cid:19) , where Γ( t ) = (cid:90) ∞ x t − exp( − x )d x is the gamma function, and F denotes the confluent hypergeometric function.Note that F ( m +12 ; m +22 ; b a ) is an entire function about a ∈ C . Therefore g m ( b , x )equals to its Taylor’s series at x = 0 for x ∈ ( − , g m ( b , x ) = (cid:88) j ≥ j )! (cid:32) ∂ j ∂x j g m ( b , (cid:33) x j . Similar to the case b , b ≤ − d , we truncate the series at j ≤ N . Again noticing b ≤ − d ,we expand the integral interval in (9) to R , leading to the approximation formulaˆ Z nm ( b , b ) =4 (cid:88) j ≤ N j )! (cid:32) ∂ j ∂x j g m ( b , (cid:33) (cid:90) R exp (cid:0) b x (cid:1) x j + n d x =4 (cid:88) j ≤ N j )! (cid:32) ∂ j ∂x j g m ( b , (cid:33) · (cid:114) π − b (2 j + n − − b ) j + n/ . (11)Next, we explain how to calculate the derivatives ∂ j g m ( b , /∂x j . Denote h ( a ) = a m/ , h ( a ) = F (cid:18) m + 12 ; m + 22 ; b a (cid:19) . Then we have ∂ j g∂a j = 12 √ π · Γ[( m + 1) / m + 2) / · j (cid:88) k =0 (cid:18) jk (cid:19) ∂ ka h · ∂ j − ka h , (12)4ith ∂ ka h = ( m/ m/ − k )! a m − k , k ≤ m , ∂ ka h = 0 , k > m , (13)and ∂ ka h = b k (cid:18) m + 12 (cid:19) ( k ) (cid:46) (cid:18) m + 22 (cid:19) ( k ) · F (cid:18) m + 12 + k ; m + 22 + k ; b a (cid:19) (14)where x (0) = 1 , x ( k ) = x ( x + 1)( x + 2) · · · ( x + k − ∂ a∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 = − , ∂ i a∂x i (cid:12)(cid:12)(cid:12)(cid:12) x =0 = 0 , i (cid:54) = 2 , and the chain rule, we arrive at ∂ j ∂x j g ( x m | b , x ) (cid:12)(cid:12)(cid:12)(cid:12) x =0 = ( − j · (2 j )! j ! · ∂ j g∂a j (cid:12)(cid:12)(cid:12)(cid:12) a =1 . (15)The derivatives ∂ j g m ( b , /∂x j are functions of b . In the routine BinghamMoments , weprecompute the values on grid points b = 0 . k , and compute the values between the gridpoints by linear interpolation. b > − d , b > − d In this bounded region of ( b , b ), we use interpolation for Z and Z mn /Z . We computethem and their derivatives about b , b , ∂Z ∂b = Z , ∂ ( Z nm /Z ) ∂b = Z n +2 ,m Z − Z nm Z Z , on the grid ( b , ) j = − j ∆ b, ≤ j ≤ − d/ ∆ b . These values are computed in advance andsaved as constants in the routine BinghamMoments . For Z nm not on the grid points, wecalculate with the interpolation described below. Suppose we already know f ( x i , y j ) , f x ( x i , y j ) , f y ( x i , y j ) , j = 1 , . To obtain the approxiamte value f ( x, y ) on ( x, y ) ∈ [ x , x ] × [ y , y ], we first calculate f ( x, y ) , f ( x, y ) , f ( x , y ) , f ( x , y )with third-order Hermite interpolation, f ( x, y ) = f ( x , y ) · (1 + 2 x − xx − x )( x − x x − x ) + f ( x , y ) · (1 + 2 x − xx − x )( x − x x − x ) + f x ( x , y ) · ( x − x )( x − x x − x ) + f x ( x , y ) · ( x − x )( x − x x − x ) . Next we calculate f y ( x, y ) , f y ( x, y ) , f x ( x , y ) , f x ( x , y )5ith linear interpolation, f y ( x, y ) = f y ( x , y ) x − xx − x + f y ( x , y ) x − x x − x . Then we can calculate f ( x, y ) with third order Hermite interpolation by f ( x, y ) = f ( x , y ) · (1 + 2 x − xx − x )( x − x x − x ) + f ( x , y ) · (1 + 2 x − xx − x )( x − x x − x ) + f x ( x , y ) · ( x − x )( x − x x − x ) + f x ( x , y ) · ( x − x )( x − x x − x ) , (16)or f ( x, y ) = f ( x, y ) · (1 + 2 y − yy − y )( y − y y − y ) + f ( x, y ) · (1 + 2 y − yy − y )( y − y y − y ) + f y ( x, y ) · ( y − y )( y − y y − y ) + f y ( x, y ) · ( y − y )( y − y y − y ) . (17)We compute f ( x, y ) as the average of (16) and (17). We have introduced four parameters in the above: the size d for dividing the domain, the orderof truncation N and N , and the grid size for the interpolation ∆ b . We choose parametersas d = 30, N = 5, N = 5, ∆ b = 0 .
025 for Z , and ∆ b = 0 . Z nm /Z in theroutine BinghamMoments , achieving maximal absolute error less than 5 × − for Z and (cid:104) x n x m (cid:105) , n + m ≤
4. We will verify this in Sec. 3.2. With these parameters, the memoryneeded for loading precomputed values (including ∂ j g m ( b , /∂x j in the case 2.2, and thevalues on the grid points in the case 2.3) is about 75MB, which is available for commoncomputers. We give an error estimate for the case 2.1 with some special functions. Denote F ( x ) = e − x (cid:90) x e t d t as the Dawson function, γ ( n, x ) = (cid:90) x t n − e − t d t as the lower incomplete gamma function, and α n ( z ) = E − n ( z ) = n ! z − n − e − z (cid:18) z + z
2! + · · · + z n n ! (cid:19) as the exponential integral function. 6 heorem 3.1. Let ˆ Z nm be defined in (8) and denote N = N . For b , b ≤ − d , it holds | Z nm − ˆ Z nm | ≤ π F ( √ d ) √ d − π N (cid:88) j =0 (2 j − j )!! · d − j − γ ( j + 1 , d )+2 π N +max( n,m ) (cid:88) j =0 (2 j − j )!! α j ( d ) . (18) Proof.
We can divide the error into two parts: e = Z nm ( b , b ) − (cid:90) (cid:90) B (0 , x n x m exp (cid:0) b x + b x (cid:1) N (cid:88) j =0 (2 j − j )!! · ( x + x ) j d x =2 (cid:90) (cid:90) B (0 , x n x m exp (cid:0) b x + b x (cid:1) (cid:88) j>N (2 j − j )!! ( x + x ) j d x , (19) e =2 (cid:90) (cid:90) R \ B (0 , x n x m exp (cid:0) b x + b x (cid:1) N (cid:88) j =0 (2 j − j )!! ( x + x ) j d x . (20)For e , we have e ≤ (cid:90) (cid:90) B (0 , exp (cid:0) − d ( x + x ) (cid:1) (cid:88) j>N (2 j − j )!! ( x + x ) j d x = 2 (cid:90) (cid:90) B (0 , exp (cid:0) − d ( x + x ) (cid:1) (cid:112) − x − x − (cid:88) j ≤ N (2 j − j )!! ( x + x ) j d x ≤ π (cid:90) e − dr r − r d r − π N (cid:88) j =0 (2 j − j )!! (cid:90) r j +1 e − dr d r = 4 π F ( √ d ) √ d − π N (cid:88) j =0 (2 j − j )!! d − j − γ ( j + 1 , d ) . (21)In the above, we use the polar coordinate transformation x = r cos θ, x = r sin θ . For e ,denote M = max { n, m } , then we have e ≤ π N + M (cid:88) j =0 (2 j − j )!! (cid:90) ∞ r j +1 e − dr d r = 2 π N + M (cid:88) j =0 (2 j − j )!! α n ( d ) . (22)Combining (21) and (22), we get (18).For our chosen parameters d = 30 and N = 5, the upper bound given by (18) is 6 . × − for n + m ≤
4. We also give the upper bound calculated from (18) for a few d and N inTable 1. The estimate (18) is also helpful to choosing parameters under different demand ofaccuracy, which will be shown in Table 3. 7
13 16 20 26 N . × − . × − . × − . × − Table 1: Absolute error bound by (18) under different values of d and N for n + m ≤ Z Z /Z Z /Z Maximal error 6 . × − . × − . × − Moment Z /Z Z /Z Z /Z Maximal error 4 . × − . × − . × − Table 2: Maximal absolute error for the 30 ,
000 pairs of ( b , b ). We compare the results calculated by our method and the results calculated by numericalintegration to testify the accuracy of our method numrically. The parameters in our methodare chosen as N = 5, N = 5 and d = 30. For numerical integration, we use adaptiveSimpson’s method to control the absolute error less than 10 − . We select randomly 10 , b i for each of the three cases: b , b ≤ − d , max( b , b ) > − d & min( b , b ) ≤ − d ,and b , b > − d , and calculate Z and the moments Z nm /Z where n + m = 2 ,
4. Table 2shows the maximal absolute errors of Z and Z nm /Z among the 30000 samples, which areunder the magnitude of 10 − . In particular, the errors of Z nm /Z are less than 5 × − . Wealso examine the distribution of the absolute errors of Z (Figure 1(a)), Z /Z (Figure 1(b))and Z /Z (Figure 1(c)) for the 10000 samples in each of three cases respectively, and findthat for most b i the absolute errors are less than 10 − . Moreover, the numerical test alsoshows our method is very fast. Calculating all these 30 ,
000 examples, the adaptive Simpson’smethod with the target accuracy 5 × − spend 3117 .
761 seconds while our method only0 .
193 seconds. Both routines are written in C and run in the same computer with a CPUclock speed 2 . d , N and N in Table 3 for different demandedaccuracy for Z nm /Z , which are also testified numerically with 30 ,
000 random samples. Bycomparing with the errors in Table 3 and Table 1, we find that the upper bound given by(18) are indicative for the choice of parameters.Demanded maximal absolute error 5 × − × − × − × − d
13 16 20 26 N N d , N and N under different demanded absoluteerror. 8 E E E E (a) absolute errors of Z E E E E E (b) absolute errors of Z /Z E E E E E (c) absolute errors of Z /Z Figure 1: Distribution of error. Blue bars: b , b ≤ − d ; Green bars: max( b , b ) > − d &min( b , b ) ≤ − d ; Yellow bars: b , b > − d . E i represents the interval [10 − i − , − i ) for i = 7 , , E = [0 , − ) and E = [10 − , ∞ ). In this section, we apply our algorithm to a Q -tensor model for rod-like liquid crystals.Compared with the original Landau-de Gennes Q -tensor theory, the model is able to constrainthe tensor within the physical range [3], and is closely connected to molecular theory [8]. Butthe Bingham distribution in the model brings difficulty in numerical simulations. We willexplain how our fast algorithm accelerates the computation.Suppose that the rod-like molecules are confined inside the unit sphere. Then the anchoringeffect on the spherical surface will induce defects for the alignment of the molecules. Weconsider the following simplified free energy, F = (cid:90) Ω d x d y d z (cid:104) ( B : ( Q + I − log Z ) − α | Q | + 12 α |∇ Q | (cid:105) + F p , (23)where the region Ω is chosen as the unit sphere, I is the identity matrix, and Q ij ( x ) = (cid:90) S ( x i x j − δ ij ) f ( x | B )d S is a symmetric traceless matrix describing the orientational distribution of rod-like moleculesat each spatial point, with f ( x | B ) and Z = Z ( B ) defined in (1) and (2). Here δ ij isthe Kronecker notation. The first two terms in the integral are the bulk energy describingthe nematic phase in equilibrium. This bulk energy is the only terms distinct from thephenomenological Landau-de Gennes theory, where the bulk energy is given as a polynomial a tr( Q ) − a tr( Q ) + a (tr( Q )) . The gradient term is the energy contribution of the spatial inhomogeneity. The boundarypenalty term F p = (cid:90) ∂ Ω d S [ Q xy − Q ( x −
13 )] +[ Q z − Q y ] +[ Q xy − Q ( y −
13 )] +[ Q z − Q x ]
9s added to enforce the value of Q on the sphere to be approximately Q = λ x − xy xzxy y − yzxz yz z − . In fact, if Q is given as above, then F p = 0. Our aim is to find local minimizers of the energyfunctional (23) that describe metastable states.Express B as B = T diag( b , b , T T , where T is orthogonal with det T = 1 and can beexpressed by Euler angles, T = cos α cos γ − cos β sin α sin γ cos γ sin α + cos α cos β sin γ sin β sin γ − cos β cos γ sin α − cos α sin γ cos α cos β cos γ − sin α sin γ cos γ sin β sin α sin β − cos α sin β cos β . In this case, Q = T diag( q , q , q ) T T , where the eigenvalues are given by q = Z ( b , b ) /Z ( b , b ), q = Z ( b , b ) /Z ( b , b ), and q = 1 − q − q .We use the spherical coordinates ( r, θ, φ ) to represent the position, i.e., x = r sin θ cos φ, y = r sin θ sin φ, z = r cos θ. (24)The integral becomes (cid:82) ( · )d x d y d z = (cid:82) ( · ) r sin θ d r d θ d φ , and the gradient term becomes |∇ Q | = | ∂ r Q | + 1 r | ∂ θ Q | + 1 r sin θ | ∂ φ Q | . (25)The free energy is discretized at N × N × N = 32 Gaussian quadrature nodes ( r j , θ k , φ l ) in[0 , × [0 , π ] × [0 , π ]. At each node ( b , b , α, β, γ ) jkl act as the basic variables, from which Q jkl is computed. The gradient term is computed using the spectral-collocation method.From the value of Q at the discretized nodes, a polynomial Q ( r, θ, φ ) = N − (cid:88) j =0 M − (cid:88) k =0 L − (cid:88) l =0 c jklQ r j θ k φ l is constructed through interpolation. The derivatives about ( r, θ, φ ), as well as the valueson the boundary, are then computed from the above polynomial. We refer to [19] where thedetails about the spectral-collocation method are illustrated. The free energy is minimizedusing the BFGS method (see, for instance, [2]). In the iteration we need to compute thederivatives of F about ( b i ) jkl , where fourth moments are involved. For instance, ∂∂b Q = T diag( ∂q ∂b , ∂q ∂b , ∂ ( − q − q ) ∂b ) T T , where ∂q ∂b = Z Z − Z Z . It is worth pointing out that at each point, the value of Q and Z are computed from B .Therefore, our algorithm is executed O ( N ) times in each BFGS iteration step, which greatlyaccelerates the simulation. Another thing is that the Bingham distribution remains the same10 a) Radial hedgehog (b) Ring disclination (c) Sphere ring band Figure 2: Three axisymmetric defect patterns, shown by the slice of x - x plane, where x is the axis of symmetry. White rods represent principal eigenvectors. The background colordescribes the biaxiality µ , with red indicates biaxial and blue indicates uniaxial. In all threecases α = 0 .
04, and α are chosen as: (a) α = 11; (b) α = 16; (c) α = 22.when we alter the parameters α , , the domain (from sphere to cylinder or ellipse, etc.), andadd some terms like in [8]. Thus our algorithm is suitable for all these cases.Before looking at the results, we first define the biaxiality. When Q (cid:54) = 0, we say Q isuniaxial if it has two identical eigenvalues, and is biaxial if it has distinct eigenvalues. Notethat tr Q = 0. The biaxiality is measured by µ = 1 − Q ) (tr Q ) . For uniaxial Q , we have µ = 0; for biaxial Q , we have 0 < µ ≤
1. We examine thedefect pattern under different α and α . At each point, the favored direction of the rod-likemolecules is the principal unit eigenvector n of Q . While Q is continuous in the unit sphere, n might be discontinuous at the points where Q = 0 or Q has two identical positive eigenvalues.Defect patterns are classified by the configuration of these points.We fix α = 0 .
04 and let α vary. Three defect patterns are observed and drawn inFigure 2: radial hedgehog (Figure 2(a)), when α = 11; ring disclination (Figure 2(b)), when α = 16; sphere ring band (Figure 2(c)), when α = 22. In the radial hedgehog pattern, Q is uniaxial everywhere with the principal eigenvector along the radial direction. The spherecenter, where Q = 0, is the only point defect. In the ring disclination pattern, the pointswhere Q has two identical positive eigenvalues form a circle in the x - y plane, round whichis a torus of biaxial region. In the sphere ring band pattern, the points where Q = 0 formtwo rings on the spherical surface. In the band between these two rings on the sphericalsurface, Q has two identical positive eigenvalues. A strong biaxial region is observed insidethe sphere near the band. The last pattern is not found in the Landau-de Gennes theory [9].We believe that this novel pattern come from the term B : ( Q + I/ − log Z , since it is theonly term different from the Landau-de Gennes theory. Hence, it is necessary for this modelto be further examined. We develop a fast and accurate algorithm to evaluate the moments of Bingham distribution.Numerical test shows that it is remarkbly faster than direct numerical quadrature, while11aintaining high accuracy. We apply the algorithm to the liquid crystal model that containsthe Bingham distribution, which is able to constrain the order parameters within the physicalrange. We examine the defect patterns of liquid cystals confined inside a sphere and find anovel pattern, suggesting that the model be examined thoroughly and compared with theLandau-de Gennes theory in future studies. Armed with our algorithm, these studies willbecome much less expensive computationally.
Acknowledgment
Pingwen Zhang is supported by National Natural Science Foundationsof China (Grant No. 11421101 and No. 11421110001).
References [1] V Alastru´e, P S´aez, MA Mart´ınez, and M Doblar´e. On the use of the Bingham statisti-cal distribution in microsphere-based constitutive models for arterial tissue.
MechanicsResearch Communications , 37(8):700–706, 2010.[2] Mordecai Avriel.
Nonlinear programming: analysis and methods . Courier Corporation,2003.[3] John M Ball and Apala Majumdar. Nematic liquid crystals: from Maier-Saupe to acontinuum theory.
Molecular crystals and liquid crystals , 525(1):1–11, 2010.[4] Christopher Bingham. An antipodally symmetric distribution on the sphere.
The Annalsof Statistics , pages 1201–1225, 1974.[5] Maxime Descoteaux, Rachid Deriche, Thomas R Kn¨osche, and Alfred Anwander. Deter-ministic and probabilistic tractography based on complex fibre orientation distributions.
Medical Imaging, IEEE Transactions on , 28(2):269–286, 2009.[6] J Feng, CV Chaubal, and LG Leal. Closure approximations for the Doi theory: Whichto use in simulating complex flows of liquid-crystalline polymers?
Journal of Rheology(1978-present) , 42(5):1095–1119, 1998.[7] Massimiliano Grosso, Pier L Maffettone, and Francois Dupret. A closure approxima-tion for nematic liquid crystals based on the canonical distribution subspace theory.
Rheologica acta , 39(3):301–310, 2000.[8] Jiequn Han, Yi Luo, Wei Wang, Pingwen Zhang, and Zhifei Zhang. From microscopictheory to macroscopic theory: a systematic study on modeling for liquid crystals.
Archivefor Rational Mechanics and Analysis , 215(3):741–809, 2014.[9] Yucheng Hu, Yang Qu, and Pingwen Zhang. On the disclination lines of nematic liquidcrystals.
Communications in Computational Physics , 19:354–379, 2016.[10] John T Kent. Asymptotic expansions for the Bingham distribution.
Applied statistics ,pages 139–144, 1987.[11] JT Kent, JC Briden, and KV Mardia. Linear and planar structure in ordered multi-variate data as applied to progressive demagnetization of palaeomagnetic remanence.
Geophysical Journal International , 75(3):593–621, 1983.1212] JL Kirschvink. The least-squares line and plane and the analysis of palaeomagnetic data.
Geophysical Journal International , 62(3):699–718, 1980.[13] Alfred Kume, SP Preston, and Andrew TA Wood. Saddlepoint approximations for thenormalizing constant of Fisher–Bingham distributions on products of spheres and Stiefelmanifolds.
Biometrika , 100(4):971–984, 2013.[14] Alfred Kume and Andrew TA Wood. Saddlepoint approximations for the Bingham andFisher–Bingham normalising constants.
Biometrika , 92(2):465–476, 2005.[15] Karsten Kunze and Helmut Schaeben. The bingham distribution of quaternions andits spherical radon transform in texture analysis.
Mathematical Geology , 36(8):917–943,2004.[16] Yixiang Luo and Jie Xu. https://github.com/yixiangLuo/Bingham-moment-function/,2016.[17] Thomas P Minka. Automatic choice of dimensionality for PCA. In
NIPS , volume 13,pages 598–604, 2000.[18] Tullis C Onstott. Application of the Bingham distribution function in paleomagneticstudies.
Journal of Geophysical Research: Solid Earth , 85(B3):1500–1510, 1980.[19] Jie Shen, Tao Tang, and Li-Lian Wang.
Spectral methods: algorithms, analysis andapplications , volume 41. Springer Science & Business Media, 2011.[20] Han Wang, Kun Li, and Pingwen Zhang. Crucial properties of the moment closure modelFENE-QE.
Journal of Non-Newtonian Fluid Mechanics , 150(2):80–92, 2008.[21] Feng Zhao, Jian Peng, and Jinbo Xu. Fragment-free approach to protein folding usingconditional neural fields.