Accurate and Fast Evaluation of Elementary Symmetric Functions Stef Graillat LIP6/PEQUAN - Université Pierre et Marie Curie (Paris 6) - CNRS Joint work with Hao Jiang and Roberto Barrio 21st IEEE International Symposium on Computer Arithmetic Austin, Texas, USA, April 7-10, 2013 S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 1 / 28
Motivations (1/2) Polynomials play a central role in computational and applied mathematics The determination of the zeros of polynomials is a classical problem of computational mathematics Inverse problem : given the zeros, determine the coefficients of the polynomial S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 2 / 28
Motivations (2/2) Characteristic polynomial of a n × n matrix A det( λ I − A ) = λ n + c 1 λ n − 1 +···+ c n − 1 λ + c n c 1 = trace( A ) c n = det( A ) Eigenvalues: ( λ i ) for i = 1,..., n n n � � c 1 = λ i c n = λ i i = 1 i = 1 → the c i are elementary symmetric functions of the λ i S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 3 / 28
Outline of the talk Motivations Classical Summation Algorithm Error-free transformations Compensated Summation Algorithm Conclusion and future work S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 4 / 28
Elementary Symmetric Functions (ESF) Definition 1 The k-th Elementary Symmetric Function (ESF) associated with a vector of n numbers X = ( x 1 ,..., x n ) is defined by � S ( n ) k ( X ) = x π 1 x π 2 ... x π k with 1 ≤ k ≤ n 1 ≤ π 1 < ... < π k ≤ n For k = 0, S ( n ) = 1 0 � n � The k -th function S ( n ) k ( X ) consists of summands k → straightforward computation is very expensive S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 5 / 28
Applications of computing ESF The ESFs appear when expanding a linear factorization of a polynomial n n n � � � c i x i = ( − 1) n − i S ( n ) n − i ( x 1 ,..., x n ) x i ( x − x i ) = i = 1 i = 0 i = 0 One can evaluate polynomial’s coefficients { c i } n i = 0 from its zeros { x i } n i = 1 , specially compute characteristic polynomials from eigenvalues Part of conditional maximum likelihood estimation (CMLE) of item parameters under the Rasch model in psychological measurement. Accurate evaluation allows much more items to be calibrated Thermodynamic properties of systems of fermions S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 6 / 28
Condition number of ESF Condition numbers measure the sensitivity of the solution of a problem to perturbation in the data Definition 2 (Condition number of the k -th ESF) � | S ( n ) k ( X +△ X ) − S ( n ) � k ( X ) | cond ( S ( n ) k ( X )) = lim ε → 0 sup : |△ X | < ε | X | ε | S ( n ) k ( X ) | A direct calculation yields kS ( n ) k ( | X | ) cond ( S ( n ) k ( X )) = | S ( n ) k ( X ) | n ( X )) = cond ( � n In particular, cond ( S ( n ) i = 1 x i ) = n and � n 1 ( X )) = cond ( � n i = 1 | x i | cond ( S ( n ) i = 1 x i ) = i = 1 x i | . | � n S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 7 / 28
Classic Summation Algorithm Algorithm 1 Input: X = ( x 1 ,..., x n ) and k Output: k-th ESF S ( n ) k ( X ) = S ( n ) k function S ( n ) k = SumESF ( X , k ) S ( i ) S ( i ) S (1) 0 = 1, 1 ≤ i ≤ n − 1; j = 0, j > i ; 1 = x 1 ; for i = 2 : n for j = max{1, i + k − n } : min{ i , k } S ( i ) j = S ( i − 1) + x i S ( i − 1) j − 1 ; j end end j ( x 1 ,..., x i ) = � S ( i ) j = S ( i ) 1 ≤ π 1 < ... < π j ≤ i x π 1 x π 2 ... x π j Substitution of j = 1 : i for j = max{1, i + k − n } : min{ i , k } makes it possible to compute all ESF simultaneously → Algorithm used in MATLAB poly function S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 8 / 28
Standard model of floating-point arithmetic Assume floating point arithmetic adhering IEEE 754 with rounding to nearest with rounding unit u (no underflow nor overflow) Let x , y ∈ F and ◦ ∈ { + , − , · ,/}. The result x ◦ y is not in general a floating-point number fl( x ◦ y ) = ( x ◦ y )(1 + δ ), | δ | ≤ u IEEE 754 standard (2008) Type Size Mantissa Exponent Unit rounding Interval u = 2 1 − 24 ≈ 1,92 × 10 − 7 ≈ 10 ± 38 binary32 32 bits 23+1 bits 8 bits u = 2 1 − 53 ≈ 2,22 × 10 − 16 ≈ 10 ± 308 binary64 64 bits 52+1 bits 11 bits We denote n u γ n : = 1 − n u S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 9 / 28
Rounding error analysis Theorem 1 (Rehman, Ipsen (2011)) If X = ( x 1 ,..., x n ) is a vector of floating-point numbers, the computed S ( n ) S ( n ) k-th elementary symmetric function � = � k ( X ) by Algorithm 1 in k floating-point arithmetic verifies � S ( n ) k − S ( n ) � � � ≤ 1 � � k k γ 2( n − 1) cond ( S ( n ) � k ), 2 ≤ k ≤ n − 1, S ( n ) k � n � S ( n ) 1 − S ( n ) � � i = 1 | x i | � � 1 � ≤ γ n − 1 cond ( S ( n ) � 1 ) = γ n − 1 | � n i = 1 x i | , k = 1, S ( n ) 1 � � S ( n ) n − S ( n ) � � ≤ 1 � � n n γ n − 1 cond ( S ( n ) n ) = γ n − 1 , k = n . � S ( n ) n S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 10 / 28
Getting more accuracy with compensated algorithms Error-free transformations are properties and algorithms to compute the generated elementary rounding errors, a , b entries ∈ F , a ◦ b = fl( a ◦ b ) + e , with e ∈ F Key tools for accurate computation fixed length expansions libraries: double-double (Briggs, Bailey, Hida, Li), quad-double (Bailey, Hida, Li) arbitrary length expansions libraries: Priest, Shewchuk compensated algorithms (Kahan, Priest, Ogita-Rump-Oishi, Graillat-Langlois-Louvet, etc.) S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 11 / 28
EFT for the summation x = fl( a ± b ) ⇒ a ± b = x + y with y ∈ F , Algorithms of Dekker (1971) and Knuth (1974) Algorithm 2 (EFT of the sum of 2 floating point numbers with | a | ≥ | b | ) function [ x , y ] = FastTwoSum ( a , b ) x = a ⊕ b y = ( a ⊖ x ) ⊕ b Algorithm 3 (EFT of the sum of 2 floating point numbers) function [ x , y ] = TwoSum ( a , b ) x = a ⊕ b z = x ⊖ a y = ( a ⊖ ( x ⊖ z )) ⊕ ( b ⊖ z ) S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 12 / 28
EFT for the product (1/3) x = fl( a · b ) a · b = x + y with y ∈ F , ⇒ Algorithm TwoProduct by Veltkamp and Dekker (1971) a = x + y and x and y non overlapping with | y | ≤ | x | . Algorithm 4 (Error-free split of a floating point number into two parts) function [ x , y ] = Split ( a ) factor = 2 s + 1 % u = 2 − p , s = ⌈ p /2 ⌉ c = factor ⊗ a x = c ⊖ ( c ⊖ a ) y = a ⊖ x S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 13 / 28
EFT for the product (2/3) Algorithm 5 (EFT of the product of 2 floating point numbers) function [ x , y ] = TwoProduct ( a , b ) x = a ⊗ b [ a 1 , a 2 ] = Split ( a ) [ b 1 , b 2 ] = Split ( b ) y = a 2 ⊗ b 2 ⊖ ((( x ⊖ a 1 ⊗ b 1 ) ⊖ a 2 ⊗ b 1 ) ⊖ a 1 ⊗ b 2 ) Theorem 2 Let a , b ∈ F and let x , y ∈ F such that [ x , y ] = TwoProduct ( a , b ) . Then, a · b = x + y , x = fl( a · b ), | y | ≤ u | x | , | y | ≤ u | a · b | , The algorithm TwoProduct requires 17 flops. S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 14 / 28
EFT for the product (3/3) Given a , b , c ∈ F , FMA ( a , b , c ) is the nearest floating point number a · b + c ∈ F Algorithm 6 (EFT of the product of 2 floating point numbers) function [ x , y ] = TwoProductFMA ( a , b ) x = a ⊗ b y = FMA ( a , b , − x ) The FMA is available for example on PowerPC, Itanium, Cell, Xeon Phi processors. S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 15 / 28
Compensated Summation Algorithm Algorithm 7 Input: X = ( x 1 ,..., x n ) and k ( n ) ( n ) Output: k -th ESF S k ( X ) = S k ( n ) function S k = CompSumESF ( X , k ) ( i ) S ( i ) S ( i ) S (1) � 0 = 1, 1 ≤ i ≤ n − 1; � j = 0, j > i ; � 1 = x 1 ; � ǫ S j = 0, ∀ i , j for i = 2 : n for j = max{1, i + k − n } : min{ i , k } [ p , β ( i ) S ( i − 1) % S ( i ) j = S ( i − 1) + x i S ( i − 1) j ] = TwoProd ( x i , � j − 1 ); j j − 1 S ( i ) j , σ ( i ) S ( i − 1) [ � j ] = TwoSum ( � , p ); j ( i ) ( i − 1) ( i − 1) ⊕ ( β ( i ) j ⊕ σ ( i ) � j = � j ) ⊕ x i ⊗ � ǫ S ǫ S ǫ S j j − 1 end end ( n ) ( n ) S ( n ) = � k ⊕ � S ǫ S k k S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 16 / 28
Error bound on the Compensated Summation Algorithm Theorem 3 For a vector of n floating-point numbers X = ( x 1 ,..., x n ) , the relative forward error bound in Algorithm satisfies ( n ) k − S ( n ) � � S � ≤ u + 1 � � k 2( n − 1) cond ( S ( n ) k γ 2 � k ( X )), S ( n ) k � S ( n ) 1 − S ( n ) � � � � 1 n − 1 cond ( S ( n ) � ≤ u + γ 2 � 1 ), S ( n ) 1 � � S ( n ) n − S ( n ) � � ≤ u + 1 � � n n γ n γ 2 n cond ( S ( n ) n ), � S ( n ) n with 2 ≤ k ≤ n − 1 , k = 1 , k = n, respectively. S. Graillat (Univ. Paris 6) Accurate and Fast Evaluation of ESF 17 / 28
Recommend
More recommend
