Evaluation of the matrix exponential function using finite elements in time

Ev aluation of the matrix exp onen tial function using finite elemen ts in time D H Gebremedhin, C A W eatherford, X Zhang, A Wynn I I I, G T anak a Department of Physics, Florida A & M University , T a lla hassee, FL 32307, USA E-mail: charl es.weather ford@famu.edu Abstract. The ev aluation of a matrix exponential function is a clas sic problem of computational linear algebr a. Many different metho ds have been e mployed for its nu merical ev aluation [Moler C and v an Loan C 1978 SIA M R eview 20 4], none of which pro duce a definitive algorithm whic h is broadly applicable a nd sufficiently acc urate, as well a s b eing reasonably fast. Herein, we e mploy a metho d which ev aulates a matrix exp onential as the solution to a first-order initia l v alue problem in a fictitious time v ariable. The new as pe ct of the pres ent implemen tatio n o f this metho d is to use finite elements in the fictitious time v ar iable. [W eatherford C A, Red E, and Wynn A 2002 Journal of Mole cular St ructur e 592 47 ] Then using an expansion in a prop er ly chosen time bas is, we a re a ble to make accurate c alculations of the exp onential of any given matrix as the s olution to a set of simultaneous equations . In t ro duction The ev aluation of the expo nen tial of a square matrix e A is a classic problem of computational linear alg ebra.[1] A large num b er of metho ds ha v e b een prop osed and used fo r its ev aluation. None of these metho ds hav e pro duced a generally applicable metho d whic h is suffi cien tly accurate as w ell as b eing reasonably fast. Thu s this problem migh t b e cons idered uns olv ed. It is clearly an imp o rtan t and p erv asiv e problem whic h arises in a wide v ariety of con texts.[2,3] A metho d whic h is capable of pro ducing a solution to esse n tially arbitrary precision w ould th us b e of great imp ortance. The presen t work uses a v ariant of a metho d described in Ref. 1, namely the in tro duction of an artificial time-parameter whic h pro duces an initial-v alue problem. Instead of calling a ‘canned’ s olv er, the pre sen t w o rk uses a metho d in tro duced by sev eral of the presen t authors [4 ] to solv e quantum mec hanical initial-v alue problems. In this metho d, a finite elemen t tec hnique is used to propagate from an initial condition at t = 0, whic h is the unit matrix, to the desired result at t = 1. The time axis is broken up into an arbitrar y n umber o f time elemen ts and t he solution is propagated from elemen t to elemen t, using a sp ecial basis in time introduced here for the first time. Evaluation of the matrix exp onential function using finite elements in time 2 The next section presen ts the analysis of the problem and describ es the solution algorithm. Then the algorithm is applied to ev aluate the exp onen tial of a num b er of test ma t r ices. Finally , the conclusions are presen ted. Analysis and Solution Algorit hm The problem at hand is the ev aluation of the exp onen tia l of a square (g enerally complex) matrix e A . The presen t metho d introduces an a r tificial time parameter so as to transform the ev aluation into the solution of an initial-v alue problem. F or a g iv en n × n square matrix A , consider the follo wing parametrized function definition: Ψ ( t ) ≡ e A t , (1) where Ψ is a lso an n × n square matrix and the desired solution is Ψ (1) = e A whic h ev olve s from the initial v alue giv en b y Ψ (0) = 1 ( 1 is the dia gonal unit matrix). This is a solution of the following linear ordinary differential equation, written for eac h matrix elemen t, ˙ Ψ ij ( t ) = n X k =1 A ik Ψ k j ( t ) , (2) where the ov er-dot stands f or the time-deriv ativ e. The time a xis extend s o v er the interv al [0 , 1]. No w break the time axis up into elemen ts t ha t extend b et w een no des t i and t i+1 , and define a lo cal time τ that spans [ − 1 , 1]. The lo cal time transformation is defined b y the relation, τ = q t − p, (3) where, q = 2 / ( t i+1 − t i ) and p = ( t i+1 + t i ) / ( t i+1 − t i ). Thus , for an arbitr a ry time elemen t e , Eq. (2) can b e written in terms of lo cal time τ as q ˙ Ψ ( e ) ij ( τ ) = n X k =1 A ( e ) ik Ψ ( e ) k j ( τ ) . (4) A t this p o in t , w e will use the follo wing ansatz for Ψ ( e ) to enforce con tin uit y b etw een t w o consecutiv e finite elemen ts Ψ ( e ) ij ( τ ) = f ( e ) ij ( τ ) + Ψ ( e − 1) ij (+1) , f ( e ) ij ( − 1) = 0 (5) and expand f ( e ) ij ( τ ) as f ( e ) ij ( τ ) = m − 1 X µ =0 B ij ( e ) µ s µ ( τ ) (6) in a basis we define b y s µ ( τ ) = Z τ − 1 T µ ( τ ) d τ (7) Evaluation of the matrix exp onential function using finite elements in time 3 where T µ ( τ ) are Cheb yshev Polynomials of the first kind.[5] Note t ha t these basis functions enforce the initial conditio n on the f ’s given in Eq. (5) since s µ ( − 1) = 0. The result f o r the decomp osition of f in m basis functions is Ψ ( e ) ij ( τ ) = m − 1 X µ =0 B ij ( e ) µ s µ ( τ ) + Ψ ( e − 1) ij (+1) (8) ˙ Ψ ( e ) ij ( τ ) = m − 1 X µ =0 B ij ( e ) µ T µ ( τ ) . (9) No w, insert (8) and (9) into Eq. (4), and m ultiply from the left b y w ( τ ) s µ ′ ( τ ) and in tegr a te f r o m − 1 to + 1 (note t hat w ( τ ) = (1 − τ 2 ) − 1 / 2 is the weigh ting function for Cheb yshev polynomials). Rearranging terms w e g et, q X µ [ Z 1 − 1 s µ ′ ( τ ) ω ( τ ) T µ ( τ ) d τ ]B ij(e) µ = X k µ A (e) ik [ Z 1 − 1 s µ ′ ( τ ) ω ( τ )s µ ( τ ) d τ ]B kj(e) µ + X k A ( e ) ik [ Z 1 − 1 s µ ′ ( τ ) ω ( τ ) T 0 ( τ ) d τ ]Ψ (e − 1) kj (+1) (10) where, T 0 ( τ ) = 1 . D efining, the inte grals in the ab ov e equation as C µ ′ µ ≡ Z 1 − 1 s µ ′ ( τ ) ω ( τ ) T µ ( τ ) d τ (11) D µ ′ µ ≡ Z 1 − 1 s µ ′ ( τ ) ω ( τ ) s µ ( τ ) d τ (12) g µ ′ ≡ Z 1 − 1 s µ ′ ( τ ) ω ( τ ) T 0 ( τ ) d τ (13) . and substituting Eqs. (11 - 13) into Eq. (10) giv es, q X µ C µ ′ µ B ij ( e ) µ = X k µ A ( e ) ik D µ ′ µ B k j ( e ) µ + g µ ′ X k A ( e ) ik Ψ ( e − 1) k j (+1) (14) or, r earr anging X µk ( q C µ ′ µ δ ik − A ( e ) ik D µ ′ µ ) B k j ( e ) µ = g µ ′ X k A ( e ) ik Ψ ( e − 1) k j (+1) (15) where δ ik is t he usual Kronec ker delta function. Then rewrite Eq. (15) as X µk Ω ( e ) ( µ ′ i )( µk ) B k j ( e ) µ = Γ ij ( e,e − 1) µ ′ (16) where Ω ( e ) ( µ ′ i )( µk ) ≡ ( q C µ ′ µ δ ik − A ( e ) ik D µ ′ µ ) (17) Γ ij ( e,e − 1) µ ′ ≡ g µ ′ X k A ( e ) ik Φ ( e − 1) k j (+1) . (18) Evaluation of the matrix exp onential function using finite elements in time 4 Equation (16) is a set of sim ultaneous equations of size ( n × m ), whic h can b e written in mat r ix for m as, Ω ( e ) B j ( e ) = Γ j ( e,e − 1) j = 1 , 2 , ..., n. (19) Here, Ω ( e ) is a (complex) matrix and for each j , Γ j ( e , e − 1 ) and B j ( e ) are v ectors. Eq . (1 9 ) applies for eac h time elemen t e . The solution is propagated from elemen t to elemen t, from t = 0 to t = 1. The ab ov e equation can be solv ed n umerically in man y w a ys, but we ha v e c hosen the metho d of LU decomp osition.[6] The prese n t metho d is ideally suited to high-p erfor ma nce computers w here the solv er of c hoice would pro bably b e iterative . In the presen t case, we apply LU decomp osition to Ω ( e ) and back substituting a ll of the Γ j ( e,e − 1) ’s, w e will hav e all the elemen ts for the matr ix (whic h can also b e view ed as three dimensional) B j ( e ) . This LU decomp osition only needs to b e done o nce since Ω ( e ) is indep enden t of time. Th us, the propaga tion just inv olv es a matrix v ector m ultiply . Then, w e emplo y Eq. (8 ) to solv e for Ψ ( e ) ( τ = 1) fo r the elemen t e, whic h, in turn, will b e used as Ψ ( e + 1 ) ( τ = − 1) for the ne xt elemen t e + 1. Starting off with a unit matrix for Ψ ( 1 ) ( t = 0), w e contin ue this pro cess till we calculate Ψ ( t = 1 ) at the last no de, whic h is the exp onen tial of the giv en matrix A . Results The calculations presen ted b elo w w ere done on a Macin tosh In tel laptop using G n u C++ whic h has mac hine accuracy limit of 2 . 22045 × 10 − 16 . As an illustration, let’s b or r ow a ’pathological’ matrix from [1], whic h w e hav e mo dified slightly to make it e v en w or se. Consider a matrix M1 give n b y , M1 = " − 73 36 − 96 47 # = " 1 3 2 4 # " − 1 0 0 − 25 # " 1 3 2 4 # − 1 . (20) The exp o nen t of M1 can b e easily calculated as, e M1 = " 1 3 2 4 # " e − 1 0 0 e − 25 # " − 2 3 / 2 1 − 1 / 2 # = " − 2 e − 1 + 3 e − 25 (3 / 2)( e − 1 − e − 25 ) − 4 e − 1 + 4 e − 25 3 e − 1 − 2 e − 25 # . The ab ov e matrix, exact to 16 decimal places, is giv en by e M1 ∼ = " − 0 . 73575888 2 3012208 0 . 5518191 6173633 1 6 − 1 . 47151776 4 6302175 1 . 1036383 2348655 1 1 # . (21) The r esult of our pro g ram is displa y ed b elo w and we run it b y using just 8 time steps and 8 basis f unctions. The result is accurate to 13 decimal places already . e M1 ∼ = " − 0 . 73575888 2 3012(181 ) 0 . 5 5181916 17363(35 8) − 1 . 47151776 4 6302(120 ) 1 . 1 0363832 34865(59 2) # . (22) Evaluation of the matrix exp onential function using finite elements in time 5 As an example of a non - diag onalizable ma t r ix, consider the fo llo wing matrix M2 , with complex eigenv a lues M2 = " 0 − 1 1 0 # (23) It can b e show n that, e M2 = " cos (1) − sin (1) sin (1) cos (1) # (24) W e are able to a c hiev e 14 decimal digit accuracy with 8 time steps and 8 ba sis functions. e M2 ∼ = " 0 . 5403023 0586814 − 0 . 841470984 80790 0 . 8414709 8480790 0 . 5403023 0586814 # (25) T able 1 sho ws the minim um n umber of basis functions, for a giv en n umber of t ime steps, whic h we re required to achiev e a precision of ± 1 × 10 − 14 on matrices whose exp onen t ial is know n exactly . The matrices c hosen are: the simplest po ssible matrix - a 2 × 2 real unit matrix, for whic h the r esult is the constant e on the diagonals, and the matrices M1 and M2 . Let’s c heck our prog r am on matr ices, which w e pic k ed randomly and for which we had no apriori kno wledge as to the result of their exp onentiation. W e fixed 8 time steps and/or 8 basis functions, and v aried the other corresp onding parameter fro m 5 to 4 0 and chec k ed how the results of the program v a ried in accuracy . F or the sak e of sa ving space, w e only displa yed the result of the last elemen t–the other elemen ts of the matrix exp o nen tia l behav ed similarly . The matrices c hosen are a 5 × 5 real matrix M3 , M3 =          − 0 . 1 − 0 . 2 − 0 . 3 − 0 . 4 − 0 . 5 − 0 . 6 − 0 . 7 − 0 . 8 − 0 . 9 − 1 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 1 2 3 4 0          (26) and a 3 × 3 complex matrix M4 M4 =     1 + i 1 − i i 1 2 i 0 1 + 2 i − 1 + i − 1 − i     . (27) F rom T able 2, one can see tha t the n umbers up to 12 decimal digits ha v e satura ted after 5 time steps a nd/ o r basis functions. Similarly , T able 3 sho ws 13 digits of accuracy as w e switc h the t w o parameters from 5 to 40, except for the cas e o f 5 basis func tions, whic h only sho ws 8 accurate significan t digits. This show s that f o r complex matrices, there is inheren tly more w ork for the progr a m to handle b ecause of the imaginary part of the mat rix elemen ts and there is apparently more sensitivit y to the num b er of basis functions used t ha n to the num b er of time steps. Evaluation of the matrix exp onential function using finite elements in time 6 T able 1. Minim um n um be r of basis functions and time s teps required for a precis ion of ± 1 × 10 − 14 for a 2 × 2 unit matrix, M1 and M 2 . 2 × 2 unit matrix M1 M2 Time steps Basis functions Time steps Basis functions Time steps Basis functions 1 11 05 - 1 11 2 9 08 7 2 9 4 8 16 6 4 8 8 7 50 5 8 7 16 6 256 4 15 6 58 5 - - 40 5 T able 2. Results of matrix e M 3 55 for typical runs of 8 time steps and 8 basis functions. Result Time s teps Basis functions 3.2103 0930 5973 118 5 8 3.2103 0930 5973 288 40 8 3.2103 0931 5373 377 8 5 3.2103 0930 5973 281 8 40 T able 3. Results of matrix e M 4 33 for typical runs of 8 time steps and 8 basis functions. Result Time steps Basis functions -0.511 9771 222980 63 - i 0 .0897 7 2811 3135 12 5 8 -0.511 9771 222980 81 - i 0.08 9772 8 1131 35 26 40 8 -0.511 9771 2 12646 60 - i 0 .0897 72810979965 8 5 -0.511 9771 222980 82 - i 0.08 9772 8 1131 35 26 8 40 Conclusion W e ha v e presen ted a robust, easily used, a nd accurate algorithm for the ev aluation of the exp onen tial of a mat r ix. W e did this by introducing an artificial time parameter and ev aluating the matrix exp o nential as the solutio n of an initial-v alue problem in this artificial time. W e solv ed the initial- v alue problem b y using finite elemen ts in time with a new time ba sis whic h we defined here so as to enforce the initial conditions on the so- lution at t he b eginning of eac h time finite elemen t. This resulted in set of sim ulta neous equations for the expansion co efficien ts. The actual algorithm employ ed here w as an LU decomp osition whic h w as v ery fast and efficien t. The relativ e effic iency of the metho d should be most apparent when implemen ted on high-p erformance computers since the algorithm is highly para llel. The metho d w as applied to sev eral matrices as a pro o f of the v alidity of the algorithm. The results of our calculations show that w e only need ab out 8 basis functions a nd 8 time steps for the matrices considered f or accuracies as great as 13 significan t digits. W e trust that this me tho d of n umerically calculating the exp o nen tia l of a matrix will b e recognized to b e a nondubious o ne! Evaluation of the matrix exp onential function using finite elements in time 7 Ac kno wlegemen ts This w ork w as supported b y the NSF C REST Cen ter for As troph ysical Science a nd T ec hnolo g y under Co o p erativ e Agreemen t HRD -06303 70. References [1] Moler C and v an Loan C 1 978 Ninete en dubious ways to c ompute the ex p onential of a m atr ix SIA M R eview 2 0 4 [2] Bellman R 1 995 Int r o duct ion to matrix analysis (SIAM) [3] Higham N J 200 8 F unctions of matric es: the ory and c omputation (SIAM) [4] W ea therford C A, Red E , and Wynn A 20 02 Solution of the time-dep endent S chr¨ odinger e quation using a b asis in time Journal of Mole cular Stru ctur e 592 47 [5] Hesthav en J S, Gottlieb S, and Gottlieb D 2007 Sp e ct r al Metho ds for Time-Dep endent Pr oblems (CAMBRIDGE UNIVERSITY P RESS) [6] William P , Sa ul T , William V and Br ian F 200 7 Numeric al r e cip es t he art of scientific c omputing (CAMBRIDGE UNIVERSITY P RESS)