Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Re: Re: Questions regarding MatrixExp, and its usage
  • Next by Date: Re: Re: Questions regarding MatrixExp, and its usage
  • Previous by thread: Re: Re: Questions regarding MatrixExp, and its usage
  • Next by thread: Re: Re: Re: Questions regarding MatrixExp, and its usage