Re: Re: Questions regarding MatrixExp, and its usage

*To*: mathgroup at smc.vnet.net*Subject*: [mg63390] Re: [mg63355] Re: [mg63335] Questions regarding MatrixExp, and its usage*From*: "Michael Chang" <michael_chang86 at hotmail.com>*Date*: Wed, 28 Dec 2005 05:24:13 -0500 (EST)*Sender*: owner-wri-mathgroup at wolfram.com

Hi, >From: Pratik Desai <pdesai1 at umbc.edu> To: mathgroup at smc.vnet.net >CC: "mathgroup at smc.vnet.net" <mathgroup at smc.vnet.net>, Andrzej >Kozlowski <akoz at mimuw.edu.pl> >Subject: [mg63390] Re: [mg63355] Re: [mg63335] Questions regarding MatrixExp, and its >usage >Date: Tue, 27 Dec 2005 12:54:37 -0500 > >Michael Chang wrote: > >>Hi, >> >>>From: Andrzej Kozlowski <akoz at mimuw.edu.pl> To: mathgroup at smc.vnet.net >>>Subject: [mg63390] Re: [mg63355] Re: [mg63335] Questions regarding MatrixExp, and >>>its usage >>>Date: Tue, 27 Dec 2005 15:52:53 +0900 >>> >>> >>>On 27 Dec 2005, at 09:58, Andrzej Kozlowski wrote: >>> >>>> >>>>On 27 Dec 2005, at 09:42, Andrzej Kozlowski wrote: >>>> >>>>> >>>>>On 27 Dec 2005, at 08:19, Andrzej Kozlowski wrote: >>>>> >>>>>> Now you should know that in general this is not going to hold in >>>>>>all of complex plane (but will hold in most). >>>>> >>>>> >>>>> >>>>>I wrote the above quite thoughtlessly: obviously there no sense in >>>>>which the equality "holds in most of the complex plane". Clearly in >>>>>every sense it holds just as often as it does not. Sorry about that; I >>>>>replied too quickly. >>>>>Andrzej Kozlowski >>>>> >>>>> >>>> >>>> >>>>One more correction is needed: in a certain obvious sense the equation >>>>does not hold "more often" than it holds since Exp is a surjective >>>>mapping of the compelx plane to itself which covers it infinitely many >>>>times (the fibre is Z - the integers). >>>> >>>>Because of the holidays I am now constanlty in a rush and can't find >>>>enough free time even to write a proper reply! >>>> >>>>Andrzej Kozlowski >>> >>> >>> >>>Actually, even the above is not strictly correct: Exp is a surjective >>>mapping form the complex plane to the complex plane minus the point 0. >>>Now that I have a little bit of time I can try to analyse the entire >>>problem more carefully. (In fact, I have not taught complex analysis for >>>over 15 years and I have become a little bit rusty. So when I first saw >>>this post I thought the problem lied in the branch discontinuity of Log, >>>which is why I wrote the relations was true in "most of the complex >>>plane". Of course I was completely wrong in this respect). >>> >>>Let's again define the function >>> >>> >>>f[x_, y_] := E^(x*y) - E^(y*Log[E^x]) >>> >>>We want to investigate where in the complex plane this is 0. This is by >>>definition the same as >>> >>> >>>f[x, y] >>> >>> >>>E^(x*y) - (E^x)^y >>> >>> >>>First, this is going to be zero for any real x and and an arbitrary y: >>> >>> >>>ComplexExpand[f[x,y],{y}] >>> >>>0 >>> >>>Secondly, suppose we have any pair of complex numbers a,b where f[a,b] >>>==0. That is: >>> >>>a /: f[a, b] = 0; >>> >>>Then we have >>> >>> >>> >>>ExpandAll[FullSimplify[ >>> f[a + 2*Pi*I, b]]] >>> >>> >>>E^(a*b + 2*I*b*Pi) - (E^a)^b >>> >>> >>>This will be zero if an only if b is an integer: >>> >>> >>>Simplify[%,Element[b,Integers]] >>> >>>0 >>> >>>So for every pair (a,b) for which the identity holds and b is not an >>>integer we can generate uncountably many pairs for which it does not >>>hold by simply adding 2*Pi*I to a. For example: >>> >>> >>>f[2,3/4] >>> >>>0 >>> >>> >>>FullSimplify[f[2 + 2*Pi*I, >>> 3/4]] >>> >>> >>>(-1 - I)*E^(3/2) >>> >>>On the other hand, we can get pairs of complex numbers for which the >>>identity holds provided the imaginary part of the first complex number >>>is not large: >>> >>> >>>f[1,2+3I] >>> >>> >>>0 >>> >>>However, for complex numbers with large imaginary part: >>> >>> >>>Simplify[f[1 + 12*I, 2 + 3*I]] >>> >>>(-E^(-34 + 27*I))* >>> (-1 + E^(12*Pi)) >>> >>>it is easy in this way to give a complete description of the pairs (a,b) >>>for which f is 0, but I will skip it and turn to matrices. >>> >>>In this case, while I am not 100% sure, I tend to believe the situation >>>to be quite analogous. We are interested in the equation >>> >>>MatrixExp[B*p]==MatrixPower[MatrixExp[B],p] >> >> >>Many thanks to Pratik, Daniel, and Andrzej for their very insightful and >>expert feedback! :) >> >>>I believe this will hold for real matrices B and (probably) all complex >>>p but will not hold in general. In fact I believe most what I wrote >>>above can be generalised to this case, although the statements and >>>proofs would be more complicated. >> >> >>Hmm ... actually, from the sample example listed below, I don't believe >>that it will hold *in general* for real B *and* real p: >> >>In[1]: params={theta->Pi^Pi,p->Sqrt[2]}; >>In[2]: B=theta {{Cot[theta],Csc[theta]},{-Csc[theta],-Cot[theta]}}; >>In[3]: test1=Simplify[MatrixExp[B p]/.params]; >>In[4]: test2=Simplify[MatrixPower[MatrixExp[B],p]/.params]; >>In[5]: Simplify[test1 == test2] >>Out[5]: False >> >>Daniel has suggested that for (square matrix) B and (scalar) p both being >>real-valued, this only will hold if B is positive definite (although I >>suspect that this also may hold with B being positive semi-definite too). >>Using the above example: Hmm ... I'm beginning to get 'cold feet', and am not so sure about the 'semi' definite claim anymore ... perhaps positive definite is 'as good as it gets' ... I don't know ... >>In[6]: BLim = Limit[B,theta->0]; >>In[7]: Eigenvalues[BLim] >>Out[7]: {0, 0} >>In[8]: MatrixPower[MatrixExp[BLim],p]==MatrixExp[BLim,p] >>Out[8]: True > >This is the special case where your matrix is Nilpotent refer to >http://mathworld.wolfram.com/MatrixExponential.html Yes, the above matrix is nilpotent; however, here's what still confuses me. For any square matrix B, one can form the JordanDecomposition, J, of B, where s.J.Inverse[s]==B. Hence, for any scalar p, B p == s.p J.Inverse[s], and MatrixExp[B p] = s.MatrixExp[p J].Inverse[s]. This all makes sense -- but what is MatrixPower[MatrixExp[B], p]? One can form MatrixExp[B] == s.MatrixExp[J].Inverse[s]; When p is a (positive) integer, say, MatrixExp[B] multiplied p times is (s.MatrixExp[J].Inverse[s]).(s.MatrixExp[J].Inverse[s]). .... .(s.MatrixExp[J].Inverse[s]) and so this is equal to s.(MatrixExp[J] multiplied p times).Inverse[s]. But what expression results for p a *NON*-integral number? How does one 'handle' the s and Inverse[s] terms, for instance? >The other thing that might be of interest is the diagonalizability of the >matrices >One can easily try this >Clear[aa, theta, p] >theta = Pi^Pi >p = Rationalize[1.414] >aa = theta {{Cot[theta], Csc[theta]}, {-Csc[theta], -Cot[theta]}}; >{mata, matb} = JordanDecomposition[aa] >test1 = mata.MatrixExp[matb, p].Inverse[mata] // N // Chop >test3 = MatrixExp[aa, p] >test1 - test3 // Chop // MatrixForm > >>{{0, 0}, {0, 0}} Mathematica's documentation mentions that MatrixExp uses Putzer's method or JordanDecomposition. >The advantage of using the Jordan Decomposition is that matb will be >diagonal for full Rank case. Hence matrix exponential is essentially the >exponential of each element of the diagonal > >In fact this may be a way to find Powers of Matrix with non integer n with >out controversies and they hold >test4 = MatrixPower[aa, Sqrt[2]] // N >test5 = mata.MatrixPower[matb, Sqrt[2]].Inverse[mata] // N >test4 - test5 // Chop > >>{{0, 0}, {0, 0}} > >> >>(By the way ... does anyone know *exactly* what the second argument for >>MatrixExp does? I've emailed Wolfram, since they only document MatrixExp >>with one argument, but I've seen their own documentation examples using >>*two* arguments; empirically, thus far, it seems that: >> >>MatrixExp[B,p]==MatrixExp[B p] >> >>with p being a scalar. But I digress ...) >> >>Anyways, for B and p real, I can 'sorta' see this point from what I >>(trivially) understand of the Spectral Mapping Theorem, since, as Andrzej >>has pointed out, >> >>Exp[a]^b !=Exp[a b] >> >>in general, with a complex, and b real; hence, any (strictly >>complex-valued) eigenvalue of >> >>MatrixPower[MatrixExp[B],p] >> >>will in general *not* be equal to >> >>MatrixExp[B p] >> >>Does this seem reasonable? >> >>Overall, too, I guess that I'm still kinda perplexed by what >>MatrixPower[B,pi] *means*? Somehow, I can feel 'comfortable' with >> >>MatrixExp[B pi] >> >>but not with >> >>MatrixPower[MatrixExp[B],pi] >> >>since I tend to think of MatrixExp[B pi] as Limit[MatrixExp[B t],t->Pi] >>(and can even revert back to an infinite power series matrix sum for an >>additional 'ease of understanding'), but, unlike general 'x^y' for >>scalars, can't quite grasp what the MatrixPower[*,*] equivalent signifies >>...) :( >> >>>Let's just illustrate this in the case of a 2 by 2 random matrix. >>> >>> >>>B=Array[Random[Integer,{1,6}]&,{2,2}] >>> >>> >>>{{6,1},{5,1}} >>> >>>Let's take some complex p, e.g. 1+I >>> >>>In[65]:= >>>N[MatrixExp[B*(1 + I)]]==N[MatrixPower[MatrixExp[B],1+I]] >>> >>>Out[65]= >>>True >>> >>>To produce a case where the relationship does not hold just imitate the >>>procedure for complex numbers given above. First we add to 2Pi * times >>>the identity matrix to B: >>> >>>Z = B + 2 Pi*I IdentityMatrix[2]; >>> >>>For p take any non-integer number, real or complex: >>> >>> >>>N[MatrixExp[Z*(1/2)]]==N[MatrixPower[MatrixExp[Z],1/2]] >>> >>> >>>False >>> >>>Andrzej Kozlowski >> >> >>Do any of my comments above make sense (or does anyone have a better >>explanation of what exactly how MatrixPower[B,p] can be interpreted with p >>not an integer)? My musings simply are from a 'layman's' perspective, and >>probably not very mathematically 'strict' ... :( >> >>Regards, >> >>Michael Anyways, again, many thanks for everyone's help and feedback! Here's wishing everyone all the best for a very Happy New Year! Bonne année! ;) Regards, Michael

**Re: Re: Questions regarding MatrixExp, and its usage**

**Re: Re: Questions regarding MatrixExp, and its usage**

**Re: Re: Questions regarding MatrixExp, and its usage**

**Re: Re: Re: Questions regarding MatrixExp, and its usage**