{"id":212,"date":"2012-07-23T13:40:49","date_gmt":"2012-07-23T18:40:49","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=212"},"modified":"2013-05-02T10:04:48","modified_gmt":"2013-05-02T15:04:48","slug":"a-balancing-act-for-the-matrix-exponential","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/23\/a-balancing-act-for-the-matrix-exponential\/","title":{"rendered":"A Balancing Act for the Matrix Exponential"},"content":{"rendered":"<!DOCTYPE html\r\n  PUBLIC \"-\/\/W3C\/\/DTD HTML 4.01 Transitional\/\/EN\">\r\n<style type=\"text\/css\">\r\n\r\nh1 { font-size:18pt; }\r\nh2.titlebg { font-size:13pt; }\r\nh3 { color:#4A4F55; padding:0px; margin:5px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:11pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\nh4 { color:#4A4F55; padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:10pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\n   \r\np { padding:0px; margin:0px 0px 20px; }\r\nimg { padding:0px; margin:0px 0px 20px; border:none; }\r\np img, pre img, tt img, li img { margin-bottom:0px; } \r\n\r\nul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }\r\nul li { padding:0px; margin:0px 0px 7px 0px; background:none; }\r\nul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }\r\nul li ol li { list-style:decimal; }\r\nol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }\r\nol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }\r\nol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }\r\nol li ol li { list-style-type:lower-alpha; }\r\nol li ul { padding-top:7px; }\r\nol li ul li { list-style:square; }\r\n\r\npre, tt, code { font-size:12px; }\r\npre { margin:0px 0px 20px; }\r\npre.error { color:red; }\r\npre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }\r\npre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }\r\n\r\n@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }\r\n\r\nspan.keyword { color:#0000FF }\r\nspan.comment { color:#228B22 }\r\nspan.string { color:#A020F0 }\r\nspan.untermstring { color:#B20000 }\r\nspan.syscmd { color:#B28C00 }\r\n\r\n.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }\r\n.footer p { margin:0px; }\r\n\r\n  <\/style><div class=\"content\"><!--introduction--><p>I have been interested in the computation of the matrix exponential, $e^A$, for a long time. A recent query from a user provides a provocative example.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#abdda55e-1fbd-460f-856a-ae8ee6c475a5\">Nineteen Dubious Ways<\/a><\/li><li><a href=\"#7401fe8a-5a7d-40df-92d7-1ae34f45adf2\">A Query from a User<\/a><\/li><li><a href=\"#12557a94-3f14-472e-a47b-b9964a056794\">Symbolic<\/a><\/li><li><a href=\"#ae3622c8-ad6c-4a19-b44c-c98a7ab1d537\">Classic MATLAB<\/a><\/li><li><a href=\"#54e46fbf-74ad-48e6-9f38-1842208fefd1\">The Three Demos<\/a><\/li><li><a href=\"#ff08226e-a285-43eb-a408-119ccabb44c0\">Scaling, Squaring, and Pade Approximations<\/a><\/li><li><a href=\"#49528121-70e4-41ab-b759-e7463d083dd5\">Balancing<\/a><\/li><li><a href=\"#de132ce4-54c5-4623-8d58-3bd6bae80aae\">Balancing the Exponential<\/a><\/li><li><a href=\"#647d2910-0031-4c32-bcaa-dd69f793787a\">Condition<\/a><\/li><li><a href=\"#6d4d7a14-97a4-481b-8eb3-103551ff3a31\">Should We Use Balancing?<\/a><\/li><li><a href=\"#5c82e9f2-c47d-4b38-8806-7e0bd33f9985\">Overscaling<\/a><\/li><li><a href=\"#36518ea3-7605-44f2-845e-23cf536a3156\">expm_new<\/a><\/li><li><a href=\"#76f13812-8b49-401d-9807-a739fa982195\">What Did I Learn?<\/a><\/li><\/ul><\/div><h4>Nineteen Dubious Ways<a name=\"abdda55e-1fbd-460f-856a-ae8ee6c475a5\"><\/a><\/h4><p>In 1978, Charlie Van Loan and I published a paper in <i>SIAM Review<\/i> entitled \"Nineteen Dubious Ways to Compute the Exponential of a Matrix\". The paper does not pick a \"best of the 19\", but cautiously suggests that the \"scaling and squaring\" algorithm might be OK. This was about the time I was tinkering with the first MATLAB and consequently every version of MATLAB has had an <tt>expm<\/tt> function, based on scaling and squaring. The <i>SIAM Review<\/i> paper proved to be very popular and in 2003 we published a followup, \"Nineteen Dubious Ways ..., Twenty-Five Years Later\". A <a title=\"http:\/\/www.cs.cornell.edu\/cv\/researchpdf\/19ways+.pdf (link no longer works)\">PDF is available<\/a> from Charlie's web site.<\/p><p>Our colleague Nick Higham reconsided the matrix exponential in 2005. Nick did a careful error analysis of scaling and squaring, improved the efficiency of the algorithm, and wrote a paper for the <i>SIAM Journal on Numerical Analysis<\/i>, \"The scaling and squaring method for the matrix exponential revisited\". A <a href=\"http:\/\/eprints.ma.man.ac.uk\/634\/01\/covered\/MIMS_ep2006_394.pdf\">PDF is available<\/a> from the University of Manchester's web site. The current version of <tt>expm<\/tt> in MATLAB is Nick's implementation of scaling and squaring.<\/p><p>A more recent review of Nick's work on the matrix exponential is provided by these <a href=\"http:\/\/www.maths.manchester.ac.uk\/~higham\/talks\/rome08_talk.pdf\">slides<\/a> for a talk he gave at a meeting in Rome in 2008.<\/p><h4>A Query from a User<a name=\"7401fe8a-5a7d-40df-92d7-1ae34f45adf2\"><\/a><\/h4><p>A few weeks ago, MathWorks Tech Support received a query from a user about the following matrix. Note that the elements of <tt>A<\/tt> range over 18 orders of magnitude.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span> <span class=\"string\">g<\/span>\r\na = 2e10;\r\nb = 4e8\/6;\r\nc = 200\/3;\r\nd = 3;\r\ne = 1e-8;\r\nA = [0 e 0; -(a+b) -d a; c 0 -c]\r\n<\/pre><pre class=\"codeoutput\">\r\nA =\r\n\r\n                         0                     1e-08                         0\r\n         -20066666666.6667                        -3               20000000000\r\n          66.6666666666667                         0         -66.6666666666667\r\n\r\n<\/pre><p>The computed matrix exponential has huge elements.<\/p><pre class=\"codeinput\">E = expm(A)\r\n<\/pre><pre class=\"codeoutput\">\r\nE =\r\n\r\n       1.7465684381715e+17         -923050477.783131     -1.73117355055901e+17\r\n     -3.07408665108297e+25      1.62463553675545e+17      3.04699053651329e+25\r\n      1.09189154376804e+17         -577057840.468934     -1.08226721572342e+17\r\n\r\n<\/pre><p>The report claimed that the right answer, obtained from a MATLAB competitor, differs from <tt>E<\/tt> by many orders of magnitude.<\/p><pre class=\"codeinput\">[  0.446849, 1.54044*10^-9, 0.462811,\r\n  -5.74307*10^6, -0.015283, -4.52654*10^6\r\n   0.447723, 1.5427*10^-9, 0.463481]\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n                  0.446849               1.54044e-09                  0.462811\r\n                  -5743070                 -0.015283                  -4526540\r\n                  0.447723                1.5427e-09                  0.463481\r\n\r\n<\/pre><h4>Symbolic<a name=\"12557a94-3f14-472e-a47b-b9964a056794\"><\/a><\/h4><p>Let's generate the symbolic representation of <tt>A<\/tt>.<\/p><pre class=\"codeinput\">a = sym(2e10);\r\nb = sym(4e8)\/6;\r\nc = sym(200)\/3;\r\nd = sym(3);\r\ne = sym(1e-8);\r\nS = [0 e 0; -(a+b) -d a; c 0 -c]\r\n<\/pre><pre class=\"codeoutput\"> \r\nS =\r\n \r\n[              0, 1\/100000000,           0]\r\n[ -60200000000\/3,          -3, 20000000000]\r\n[          200\/3,           0,      -200\/3]\r\n \r\n<\/pre><p>Now have the Symbolic Toolbox compute the matrix exponential, then convert the result to floating point. We can regard this as the \"right answer\". We see that it agrees with the user's expectations.<\/p><pre class=\"codeinput\">X = real(double(expm(S)))\r\n<\/pre><pre class=\"codeoutput\">\r\nX =\r\n\r\n         0.446849468283175      1.54044157383952e-09         0.462811453558774\r\n         -5743067.77947947       -0.0152830038686819         -4526542.71278401\r\n         0.447722977849494      1.54270484519591e-09         0.463480648837651\r\n\r\n<\/pre><h4>Classic MATLAB<a name=\"ae3622c8-ad6c-4a19-b44c-c98a7ab1d537\"><\/a><\/h4><p>I ran my old Fortran MATLAB from 1980.  Here is the output. It got the right answer.<\/p><pre class=\"codeinput\"><span class=\"comment\">% &lt;&gt;<\/span>\r\n<span class=\"comment\">% A = &lt;0 e 0; -(a+b) -d a; c 0 -c&gt;<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%  A     =<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%     0.000000000000000D+00   1.000000000000000D-08   0.000000000000000D+00<\/span>\r\n<span class=\"comment\">%    -2.006666666666667D+10  -3.000000000000000D+00   2.000000000000000D+10<\/span>\r\n<span class=\"comment\">%     6.666666666666667D+01   0.000000000000000D+00  -6.666666666666667D+01<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">% &lt;&gt;<\/span>\r\n<span class=\"comment\">% exp(A)<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%  ANS   =<\/span>\r\n<span class=\"comment\">%<\/span>\r\n<span class=\"comment\">%     4.468494682831735D-01   1.540441573839520D-09   4.628114535587735D-01<\/span>\r\n<span class=\"comment\">%    -5.743067779479621D+06  -1.528300386868247D-02  -4.526542712784168D+06<\/span>\r\n<span class=\"comment\">%     4.477229778494929D-01   1.542704845195912D-09   4.634806488376499D-01<\/span>\r\n<\/pre><h4>The Three Demos<a name=\"54e46fbf-74ad-48e6-9f38-1842208fefd1\"><\/a><\/h4><p>In addition to <tt>expm<\/tt>, MATLAB has for many years provided three demo functions that illustrate popular methods for computing $e^A$. The function <tt>expmdemo1<\/tt> is a MATLAB implementation of the scaling and squaring algorithm that was used in the builtin <tt>expm<\/tt> before Higham's improvements. The function <tt>expmdemo2<\/tt> implements the Taylor power series that is often the definition of $e^A$, but which is one of the worst of the nineteen ways because it is slow and numerically unreliable. The function <tt>expmdemo3<\/tt> uses eigenvalues and eigenvectors, which is OK only if the eigenvector matrix is well conditioned. (MATLAB also has a function <tt>expm1<\/tt>, which computes the scalar function $e^x\\!-\\!1$ without computing $e^x$. The <tt>m<\/tt> in the name is for minus, not matrix. This function has nothing to do with the matrix exponential.)<\/p><p>Let's see what the three demo functions do with our example.<\/p><pre class=\"codeinput\"><span class=\"comment\">%  Scaling and squaring<\/span>\r\nE1 = expmdemo1(A)\r\n<\/pre><pre class=\"codeoutput\">\r\nE1 =\r\n\r\n         0.446848323199335      1.54043901480671e-09         0.462810666904014\r\n         -5743177.01871262       -0.0152833835375292         -4526656.46142213\r\n         0.447721814330828      1.54270222301338e-09          0.46347984316225\r\n\r\n<\/pre><pre class=\"codeinput\"><span class=\"comment\">%  Taylor series<\/span>\r\nE2 = expmdemo2(A)\r\n<\/pre><pre class=\"codeoutput\">\r\nE2 =\r\n\r\n         -3627968682.81884         0.502451507654604         -3062655286.68657\r\n     -1.67974375988037e+19          3498209047.28622     -2.27506724048955e+19\r\n           15580992163.692          7.53393732504015          4987630142.66227\r\n\r\n<\/pre><pre class=\"codeinput\"><span class=\"comment\">%  Eigenvalues and eigenvectors<\/span>\r\nE3 = expmdemo3(A)\r\n<\/pre><pre class=\"codeoutput\">\r\nE3 =\r\n\r\n         0.446849468283181      1.54044157383954e-09         0.462811453558778\r\n         -5743067.77947891         -0.01528300386868         -4526542.71278343\r\n           0.4477229778495      1.54270484519593e-09         0.463480648837654\r\n\r\n<\/pre><p>You can see that both <tt>eigdemo1<\/tt>, the outdated scaling and squaring, and <tt>eigdemo3<\/tt>, eigenvectors, get the right answer, while <tt>eigdemo2<\/tt>, Taylor series, blows up.<\/p><h4>Scaling, Squaring, and Pade Approximations<a name=\"ff08226e-a285-43eb-a408-119ccabb44c0\"><\/a><\/h4><p>In outline, the scaling and squaring algorithm for computing $e^A$ is:<\/p><div><ul><li>Pick an integer $s$ and $\\sigma = 2^s$ so that $||A\/\\sigma|| \\approx 1$.<\/li><li>Find a Pade approximation, $P \\approx \\mbox{exp}(A\/\\sigma)$.<\/li><li>Use repeated squaring to compute $e^A \\approx P^\\sigma$.<\/li><\/ul><\/div><p>We have two implementations of scaling and squaring, the outdated one in <tt>expmdemo1<\/tt> and the current one in <tt>expm<\/tt>. It turns out that, for this matrix, the old implementation decides the scale factor should be <tt>2^37<\/tt> while the current implementation chooses <tt>2^32<\/tt>. Using the new scale factor  will save five matrix multiplications in the unscaling by repeated squaring.<\/p><p>The key to the comparison of these two implementations lies in the eigenvalues of the Pade approximants.<\/p><pre class=\"codeinput\">P = expmdemo1(A\/2^37);\r\ne = eig(P)\r\n\r\nP = expm(A\/2^32);\r\ne = eig(P)\r\n<\/pre><pre class=\"codeoutput\">\r\ne =\r\n\r\n         0.999999999539043\r\n         0.999999999954888\r\n         0.999999999999176\r\n\r\n\r\ne =\r\n\r\n         0.999999974526488\r\n          1.00000000930672\r\n         0.999999999946262\r\n\r\n<\/pre><p>In this case, the old code produces eigenvalues that are less than one. Powers of these eigenvalues, and hence powers of <tt>P<\/tt>, remain bounded. But the current code happens to produce an eigenvalue slightly larger than one.  The powers of <tt>e<\/tt> and of <tt>P<\/tt> blow up.<\/p><pre class=\"codeinput\">e.^(2^32)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n      3.05317714952674e-48\r\n      2.28895048607366e+17\r\n         0.793896973586281\r\n\r\n<\/pre><h4>Balancing<a name=\"49528121-70e4-41ab-b759-e7463d083dd5\"><\/a><\/h4><p>One cure, at least in this instance, is balancing. Balancing is a diagonal similarity transformation that tries to make the matrix closer to symmetric by making the row norms equal to the column norms.  This may improve the accuracy of computed eigenvalues, but seriously alter the eigenvectors. The MATLAB documentation has a good discussion of the effect of balancing on eigenvectors.<\/p><pre class=\"codeinput\"><span class=\"comment\">% doc balance<\/span>\r\n<\/pre><h4>Balancing the Exponential<a name=\"de132ce4-54c5-4623-8d58-3bd6bae80aae\"><\/a><\/h4><p>Balancing can have sometimes have a beneficial effect in the computation of $e^A$. For the example in this blog the elements of the diagonal similarity transform are powers of 2 that vary over a wide range.<\/p><pre class=\"codeinput\">[T,B] = balance(A);\r\nT\r\nlog2T = diag(log2(diag(T)))\r\n<\/pre><pre class=\"codeoutput\">\r\nT =\r\n\r\n       9.5367431640625e-07                         0                         0\r\n                         0                      2048                         0\r\n                         0                         0       1.9073486328125e-06\r\n\r\n\r\nlog2T =\r\n\r\n   -20     0     0\r\n     0    11     0\r\n     0     0   -19\r\n\r\n<\/pre><p>In the balanced matrix $B = T^{-1} A T$, the tiny element $a_{1,2}=10^{-8}$ has been magnified to be comparable with the other elements.<\/p><pre class=\"codeinput\">B\r\n<\/pre><pre class=\"codeoutput\">\r\nB =\r\n\r\n                         0               21.47483648                         0\r\n          -9.3442698319753                        -3          18.6264514923096\r\n          33.3333333333333                         0         -66.6666666666667\r\n\r\n<\/pre><p>Computing the $e^B$ presents no difficulties. The final result obtained by reversing the scaling is what we have come to expect for this example.<\/p><pre class=\"codeinput\">M = T*expm(B)\/T\r\n<\/pre><pre class=\"codeoutput\">\r\nM =\r\n\r\n         0.446849468283175      1.54044157383952e-09         0.462811453558774\r\n         -5743067.77947947       -0.0152830038686819           -4526542.712784\r\n         0.447722977849495      1.54270484519591e-09          0.46348064883765\r\n\r\n<\/pre><h4>Condition<a name=\"647d2910-0031-4c32-bcaa-dd69f793787a\"><\/a><\/h4><p>Nick Higham has contributed his <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/20820-the-matrix-function-toolbox\">Matrix Function Toolbox<\/a> to MATLAB Central.  The Toolbox has many useful functions, including <tt>expm_cond<\/tt>, which computes the condition number of the matrix exponential function.  Balancing improves the conditioning of this example by 16 orders of magnitude.<\/p><pre class=\"codeinput\">addpath(<span class=\"string\">'..\/..\/MFToolbox'<\/span>)\r\nexpm_cond(A)\r\nexpm_cond(B)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n      3.16119437847582e+18\r\n\r\n\r\nans =\r\n\r\n          238.744549689702\r\n\r\n<\/pre><h4>Should We Use Balancing?<a name=\"6d4d7a14-97a4-481b-8eb3-103551ff3a31\"><\/a><\/h4><p>Bob Ward, in the original 1977 paper on scaling and squaring, recommended balancing. Nick Higham includes balancing in the pseudocode algorithm in his 2005 paper, but in recent email with me he was reluctant to recommend it. I am also reluctant to draw any conclusions from this one case. Its scaling is too bizarre. Besides, there is a better solution, avoid overscaling.<\/p><h4>Overscaling<a name=\"5c82e9f2-c47d-4b38-8806-7e0bd33f9985\"><\/a><\/h4><p>In 2009 Nick Higham's Ph. D. Student, Awad H. Al-Mohy, wrote a dissertation entitled \"A New Scaling and Squaring Algorithm for the Matrix Exponential\". The dissertation described a MATLAB function <tt>expm_new<\/tt>. A PDF of the dissertation and a zip file with the code <a href=\"http:\/\/eprints.ma.man.ac.uk\/1442\">are available<\/a> from the University of Manchester's web site.<\/p><p>If $A$ is not a normal matrix, then the norm of the power, $||A^k||$, can grow much more slowly than the power of the norm, $||A||^k$. As a result, it is possible to suffer significant roundoff error in the repeated squaring. The choice of the scale factor $2^s$ involves a delicate compromise between the accuracy of the Pade approximation and the number of required squarings.<\/p><p>Here is one experiment with our example and various choices of $s$. The function <tt>padexpm<\/tt> is the order 13 Pade approximation taken from <tt>expm<\/tt>. We see that $s = 32$, which is the choice made by <tt>expm<\/tt>, is the worst possible choice.  The old choice, $s = 37$, is much better. These large values of $s$ result from the fact that this particular matrix has a large norm.  For this example, values of $s$ less than 10 are much better.<\/p><pre class=\"codeinput\">warning(<span class=\"string\">'off'<\/span>,<span class=\"string\">'MATLAB:nearlySingularMatrix'<\/span>)\r\nerr = zeros(40,1);\r\n<span class=\"keyword\">for<\/span> s = 1:40\r\n   P = padexpm(A\/2^s);\r\n   <span class=\"keyword\">for<\/span> k = 1:s\r\n      P = P*P;\r\n   <span class=\"keyword\">end<\/span>\r\n   err(s) = norm((P-X).\/X,inf);\r\n<span class=\"keyword\">end<\/span>\r\nsemilogy(err)\r\nxlabel(<span class=\"string\">'s'<\/span>)\r\nylabel(<span class=\"string\">'error'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/balancing_01.png\" alt=\"\"> <h4>expm_new<a name=\"36518ea3-7605-44f2-845e-23cf536a3156\"><\/a><\/h4><p>So how does the latest <tt>expm<\/tt> from Manchester do on this example? It chooses $s = 8$ and does a fine job.<\/p><pre class=\"codeinput\">expm_new(A)\r\n<\/pre><pre class=\"codeoutput\">\r\nans =\r\n\r\n         0.446849468283145      1.54044157383964e-09          0.46281145355881\r\n         -5743067.77947979       -0.0152830038686872         -4526542.71278561\r\n         0.447722977849464      1.54270484519604e-09         0.463480648837687\r\n\r\n<\/pre><h4>What Did I Learn?<a name=\"76f13812-8b49-401d-9807-a739fa982195\"><\/a><\/h4><p>This example is an extreme outlier, but it is instructive. The condition number of the problem is terrible. Small changes in the data might make huge changes in the result, but I haven't investigated that. The computed result might be the exact result for some matrix near the given one, but I haven't pursued that. The current version of <tt>expm<\/tt> in MATLAB computed an awful result, but it was sort of unlucky. We have seen at least half a dozen other functions, including classic MATLAB and <tt>expm<\/tt> with balancing, that get the right answer. It looks like <tt>expm_new<\/tt> should find its way into MATLAB.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_300aba04d9ef4d7f91b79bf4313bce3e() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='300aba04d9ef4d7f91b79bf4313bce3e ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 300aba04d9ef4d7f91b79bf4313bce3e';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2012 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_300aba04d9ef4d7f91b79bf4313bce3e()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.14<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; 7.14<br><\/p><\/div><!--\r\n300aba04d9ef4d7f91b79bf4313bce3e ##### SOURCE BEGIN #####\r\n%% A Balancing Act for the Matrix Exponential\r\n% I have been interested in the computation of the matrix exponential,\r\n% $e^A$, for a long time.\r\n% A recent query from a user provides a provocative example.\r\n\r\n%% Nineteen Dubious Ways\r\n% In 1978, Charlie Van Loan and I published a paper in _SIAM Review_\r\n% entitled \"Nineteen Dubious Ways to Compute the Exponential of a Matrix\".\r\n% The paper does not pick a \"best of the 19\", but cautiously suggests\r\n% that the \"scaling and squaring\" algorithm might be OK.\r\n% This was about the time I was tinkering with the first MATLAB and\r\n% consequently every version of MATLAB has had an |expm| function,\r\n% based on scaling and squaring.\r\n% The _SIAM Review_ paper proved to be very popular and in 2003 we\r\n% published a followup, \"Nineteen Dubious Ways ..., Twenty-Five Years Later\".\r\n% A <http:\/\/www.cs.cornell.edu\/cv\/researchpdf\/19ways+.pdf PDF is available>\r\n% from Charlie's web site.\r\n\r\n%%\r\n% Our colleague Nick Higham reconsided the matrix exponential in 2005.\r\n% Nick did a careful error analysis of scaling and squaring,\r\n% improved the efficiency of the algorithm, and wrote a paper\r\n% for the _SIAM Journal on Numerical Analysis_,\r\n% \"The scaling and squaring method for the matrix exponential revisited\".\r\n% A <http:\/\/eprints.ma.man.ac.uk\/634\/01\/covered\/MIMS_ep2006_394.pdf\r\n% PDF is available> from the University of Manchester's web site.\r\n% The current version of |expm| in MATLAB is Nick's implementation\r\n% of scaling and squaring.\r\n\r\n%%\r\n% A more recent review of Nick's work on the matrix exponential is\r\n% provided by these\r\n% <http:\/\/www.maths.manchester.ac.uk\/~higham\/talks\/rome08_talk.pdf\r\n% slides> for a talk he gave at a meeting in Rome in 2008.\r\n\r\n%% A Query from a User\r\n% A few weeks ago, MathWorks Tech Support received a query from a user\r\n% about the following matrix.\r\n% Note that the elements of |A| range over 18 orders of magnitude.\r\n\r\nformat long g\r\na = 2e10;\r\nb = 4e8\/6;\r\nc = 200\/3;\r\nd = 3;\r\ne = 1e-8;\r\nA = [0 e 0; -(a+b) -d a; c 0 -c]\r\n\r\n%%\r\n% The computed matrix exponential has huge elements.\r\nE = expm(A)\r\n\r\n%%\r\n% The report claimed that the right answer, obtained from a MATLAB competitor,\r\n% differs from |E| by many orders of magnitude.\r\n[  0.446849, 1.54044*10^-9, 0.462811,\r\n  -5.74307*10^6, -0.015283, -4.52654*10^6\r\n   0.447723, 1.5427*10^-9, 0.463481]\r\n\r\n%% Symbolic\r\n% Let's generate the symbolic representation of |A|.\r\na = sym(2e10);\r\nb = sym(4e8)\/6;\r\nc = sym(200)\/3;\r\nd = sym(3);\r\ne = sym(1e-8);\r\nS = [0 e 0; -(a+b) -d a; c 0 -c]\r\n\r\n%%\r\n% Now have the Symbolic Toolbox compute the matrix exponential,\r\n% then convert the result to floating point.\r\n% We can regard this as the \"right answer\".\r\n% We see that it agrees with the user's expectations.\r\nX = real(double(expm(S)))\r\n\r\n%% Classic MATLAB\r\n% I ran my old Fortran MATLAB from 1980.  Here is the output.\r\n% It got the right answer.  \r\n\r\n% <>\r\n% A = <0 e 0; -(a+b) -d a; c 0 -c>\r\n% \r\n%  A     =\r\n% \r\n%     0.000000000000000D+00   1.000000000000000D-08   0.000000000000000D+00\r\n%    -2.006666666666667D+10  -3.000000000000000D+00   2.000000000000000D+10\r\n%     6.666666666666667D+01   0.000000000000000D+00  -6.666666666666667D+01\r\n% \r\n% <>\r\n% exp(A)\r\n% \r\n%  ANS   =\r\n% \r\n%     4.468494682831735D-01   1.540441573839520D-09   4.628114535587735D-01\r\n%    -5.743067779479621D+06  -1.528300386868247D-02  -4.526542712784168D+06\r\n%     4.477229778494929D-01   1.542704845195912D-09   4.634806488376499D-01\r\n\r\n%% The Three Demos\r\n% In addition to |expm|, MATLAB has for many years provided three demo\r\n% functions that illustrate popular methods for computing $e^A$.\r\n% The function |expmdemo1| is a MATLAB implementation of the scaling and\r\n% squaring algorithm that was used in the builtin |expm| before Higham's\r\n% improvements.\r\n% The function |expmdemo2| implements the Taylor power series that is often the\r\n% definition of $e^A$, but which is one of the worst of the nineteen ways\r\n% because it is slow and numerically unreliable.\r\n% The function |expmdemo3| uses eigenvalues and eigenvectors, which is OK\r\n% only if the eigenvector matrix is well conditioned.\r\n% (MATLAB also has a function |expm1|, which computes the scalar\r\n% function $e^x\\!-\\!1$ without computing $e^x$.\r\n% The |m| in the name is for minus, not matrix.\r\n% This function has nothing to do with the matrix exponential.)\r\n\r\n%%\r\n% Let's see what the three demo functions do with our example.\r\n\r\n%  Scaling and squaring\r\nE1 = expmdemo1(A)\r\n\r\n%%\r\n\r\n%  Taylor series\r\nE2 = expmdemo2(A)\r\n\r\n%%\r\n\r\n%  Eigenvalues and eigenvectors\r\nE3 = expmdemo3(A)\r\n\r\n%%\r\n% You can see that both |eigdemo1|, the outdated scaling and squaring,\r\n% and |eigdemo3|, eigenvectors, get the right answer, while |eigdemo2|,\r\n% Taylor series, blows up.\r\n\r\n%% Scaling, Squaring, and Pade Approximations\r\n% In outline, the scaling and squaring algorithm for computing $e^A$ is:\r\n%\r\n% * Pick an integer $s$ and $\\sigma = 2^s$ so that $||A\/\\sigma|| \\approx 1$.\r\n% * Find a Pade approximation, $P \\approx \\mbox{exp}(A\/\\sigma)$.\r\n% * Use repeated squaring to compute $e^A \\approx P^\\sigma$.\r\n\r\n%%\r\n% We have two implementations of scaling and squaring,\r\n% the outdated one in |expmdemo1| and the current one in |expm|.\r\n% It turns out that, for this matrix, the old implementation decides the scale\r\n% factor should be |2^37| while the current implementation chooses |2^32|.\r\n% Using the new scale factor  will save five matrix multiplications\r\n% in the unscaling by repeated squaring.\r\n\r\n%%\r\n% The key to the comparison of these two implementations lies in\r\n% the eigenvalues of the Pade approximants.\r\n\r\nP = expmdemo1(A\/2^37);\r\ne = eig(P)\r\n\r\nP = expm(A\/2^32);\r\ne = eig(P)\r\n\r\n%%\r\n% In this case, the old code produces eigenvalues that are less than one.\r\n% Powers of these eigenvalues, and hence powers of |P|, remain bounded.\r\n% But the current code happens to produce an eigenvalue slightly larger\r\n% than one.  The powers of |e| and of |P| blow up.\r\n\r\ne.^(2^32)\r\n\r\n%% Balancing\r\n% One cure, at least in this instance, is balancing.\r\n% Balancing is a diagonal similarity transformation that tries\r\n% to make the matrix closer to symmetric by making the row norms equal\r\n% to the column norms.  This may improve the accuracy of computed\r\n% eigenvalues, but seriously alter the eigenvectors.\r\n% The MATLAB documentation has a good discussion of the effect\r\n% of balancing on eigenvectors.\r\n\r\n% doc balance\r\n\r\n%% Balancing the Exponential\r\n% Balancing can have sometimes have a beneficial effect in the\r\n% computation of $e^A$.\r\n% For the example in this blog the elements of the diagonal similarity\r\n% transform are powers of 2 that vary over a wide range.\r\n[T,B] = balance(A);\r\nT\r\nlog2T = diag(log2(diag(T)))\r\n\r\n%%\r\n% In the balanced matrix $B = T^{-1} A T$, the tiny element $a_{1,2}=10^{-8}$\r\n% has been magnified to be comparable with the other elements.\r\nB\r\n\r\n%%\r\n% Computing the $e^B$ presents no difficulties.\r\n% The final result obtained by reversing the scaling is what\r\n% we have come to expect for this example.\r\nM = T*expm(B)\/T\r\n\r\n%% Condition\r\n% Nick Higham has contributed his \r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/20820-the-matrix-function-toolbox Matrix Function Toolbox>\r\n% to MATLAB Central.  The Toolbox has many useful functions, including\r\n% |expm_cond|, which computes the condition number of the matrix exponential\r\n% function.  Balancing improves the conditioning of this example by\r\n% 16 orders of magnitude.\r\n\r\naddpath('..\/..\/MFToolbox')\r\nexpm_cond(A)\r\nexpm_cond(B)\r\n\r\n%% Should We Use Balancing?\r\n% Bob Ward, in the original 1977 paper on scaling and squaring,\r\n% recommended balancing.\r\n% Nick Higham includes balancing in the pseudocode algorithm in his\r\n% 2005 paper, but in recent email with me he was reluctant to recommend it.\r\n% I am also reluctant to draw any conclusions from this one case.\r\n% Its scaling is too bizarre.\r\n% Besides, there is a better solution, avoid overscaling.\r\n\r\n%% Overscaling\r\n% In 2009 Nick Higham's Ph. D. Student, Awad H. Al-Mohy, wrote a dissertation\r\n% entitled \"A New Scaling and Squaring Algorithm for the Matrix Exponential\".\r\n% The dissertation described a MATLAB function |expm_new|.\r\n% A PDF of the dissertation and a zip file with the code\r\n% <http:\/\/eprints.ma.man.ac.uk\/1442 are available>\r\n% from the University of Manchester's web site.\r\n\r\n%%\r\n% If $A$ is not a normal matrix, then the norm of the power, $||A^k||$,\r\n% can grow much more slowly than the power of the norm, $||A||^k$.\r\n% As a result, it is possible to suffer significant roundoff error\r\n% in the repeated squaring.\r\n% The choice of the scale factor $2^s$ involves a delicate compromise\r\n% between the accuracy of the Pade approximation and the number of\r\n% required squarings.\r\n\r\n%%\r\n% Here is one experiment with our example and various choices of $s$.\r\n% The function |padexpm| is the order 13 Pade approximation taken from |expm|.\r\n% We see that $s = 32$, which is the choice made by |expm|, is the worst\r\n% possible choice.  The old choice, $s = 37$, is much better. \r\n% These large values of $s$ result from the fact that this particular matrix\r\n% has a large norm.  For this example, values of $s$ less than 10 are much\r\n% better.\r\n\r\nwarning('off','MATLAB:nearlySingularMatrix')\r\nerr = zeros(40,1);\r\nfor s = 1:40\r\n   P = padexpm(A\/2^s);\r\n   for k = 1:s\r\n      P = P*P;\r\n   end\r\n   err(s) = norm((P-X).\/X,inf);\r\nend\r\nsemilogy(err)\r\nxlabel('s')\r\nylabel('error')\r\n\r\n%% expm_new\r\n% So how does the latest |expm| from Manchester do on this example?\r\n% It chooses $s = 8$ and does a fine job.\r\nexpm_new(A)\r\n\r\n%% What Did I Learn?\r\n% This example is an extreme outlier, but it is instructive.\r\n% The condition number of the problem is terrible.\r\n% Small changes in the data might make huge changes in the result,\r\n% but I haven't investigated that.\r\n% The computed result might be the exact result for some matrix near\r\n% the given one, but I haven't pursued that.\r\n% The current version of |expm| in MATLAB computed an awful result,\r\n% but it was sort of unlucky.\r\n% We have seen at least half a dozen other functions, including\r\n% classic MATLAB and |expm| with balancing, that get the right answer.\r\n% It looks like |expm_new| should find its way into MATLAB.\r\n\r\n\r\n##### SOURCE END ##### 300aba04d9ef4d7f91b79bf4313bce3e\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction--><p>I have been interested in the computation of the matrix exponential, $e^A$, for a long time. A recent query from a user provides a provocative example.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/23\/a-balancing-act-for-the-matrix-exponential\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[6],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/212"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=212"}],"version-history":[{"count":8,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/212\/revisions"}],"predecessor-version":[{"id":650,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/212\/revisions\/650"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=212"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=212"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=212"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}