Re: Re: Questions regarding MatrixExp, and its usage

*To*: mathgroup at smc.vnet.net*Subject*: [mg63391] 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 Andrzej, >From: Andrzej Kozlowski <akoz at mimuw.edu.pl> To: mathgroup at smc.vnet.net >Subject: [mg63391] Re: [mg63355] Re: [mg63335] Questions regarding MatrixExp, and its >usage >Date: Wed, 28 Dec 2005 06:51:40 +0900 > > >On 28 Dec 2005, at 01:09, Michael Chang wrote: > >>Hi, >> >>>From: Andrzej Kozlowski <akoz at mimuw.edu.pl> To: mathgroup at smc.vnet.net >>><michael_chang86 at hotmail.com> >>>Subject: [mg63391] Re: [mg63355] Re: [mg63335] Questions regarding MatrixExp, and >>>its usage >>>Date: Tue, 27 Dec 2005 15:52:53 +0900 >>> >>>*This message was transferred with a trial version of CommuniGate (tm) >>>Pro* >>> >>>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: >> >>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 >> >>(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 >> >> > > >You are right of course. I was much too quick optimistic to claim that it >would hold for all real matrices. Without giving this much thought, I >imagined that this can be reduced just to the result holding for the >eigenvalues which is, at best, only the case for certain types of >matrices; in particular real normal ones. >I imagine that MatrixPower for arbitrary matrix and arbitrary exponent >is defined via the Jordan decomposition, by first you defining it for >Jordan blocks and then taking the suitable sum and finally applying the >similarity matrices. I have not considered carefully what happens for a >single Jordan block matrix, but I think the power matrix will have powers >of the eigenvalues on the diagonal and lower powers of the eigenvalue >given by differentiation z^p (where p is our exponent) with certain >coefficients above the diagonal (essentially terms of the "Taylor >expansion" of z^p). If I remember correctly, this is how one defines f[A] >for an arbitrary smooth function f and an arbitrary complex matrix A. >Thinking of this definition, however, I can see no reason, why it should >be "well behaved" in this case for non-diagonalizable matrices, so I >suspect real normal matrices (perhaps without zero eigenvalues) are the >best you can expect in general -although I am not going to have any time >to check this until well after the New Year. > >Wishing everyone a Happy New Year. > >Andrzej Kozlowski Once again, your mathematical expertise and feedback are greatly appreciated! For an arbitrary square B matrix, and taking a JordanDecomposition of B, such that s.J.Inverse[s]==B, I can see that MatrixPower[B,p] == (s.J.Inverse[s] multiplied p times) == s.(J multiplied p times).Inverse[s] where p is a (positive) integer, say. But, to my confusion, what does one obtain when p=Sqrt(Pi), say? For instance, do we still somehow get multiples of Inverse[s].s==IdentityMatrix[n] in the resulting expression, somehow? Anyways, my sincere thanks again for all of your help, and best wishes to everyone for a Happy New Year! Bonne année! ;) Regards, Michael