{"id":9569,"date":"2023-01-05T11:37:04","date_gmt":"2023-01-05T16:37:04","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=9569"},"modified":"2023-01-03T11:51:28","modified_gmt":"2023-01-03T16:51:28","slug":"singular-matrix-pencils-and-the-qz-algorithm","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2023\/01\/05\/singular-matrix-pencils-and-the-qz-algorithm\/","title":{"rendered":"Singular Matrix Pencils and the QZ Algorithm"},"content":{"rendered":"\r\n\r\n<div class=\"content\"><!--introduction--><p>This year, 2023, is the 50-th anniversary of the QZ algorithm for generalized matrix eignenvalue problems,<\/p><pre class=\"language-matlab\">Ax = &#955;Bx\r\n<\/pre><p>The algorithm computes these eigevalues without inverting either <tt>A<\/tt> or <tt>B<\/tt>.  And, the QZ-algorithm can help detect and analyze exceptional situaions known as <i>singular pencils<\/i>.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#19d0dd38-a456-4a63-bf2f-18e02689531d\">Matrix pencils<\/a><\/li><li><a href=\"#31c19b00-cfd8-493d-b712-08610ce3c5aa\">Example<\/a><\/li><li><a href=\"#0f5d9fb9-0792-4d4c-a9ac-c5fd3e1e7691\">Wilkinson example<\/a><\/li><li><a href=\"#cbe4eaf1-8a19-41e7-badd-1b8e69f74b61\">References<\/a><\/li><\/ul><\/div><h4>Matrix pencils<a name=\"19d0dd38-a456-4a63-bf2f-18e02689531d\"><\/a><\/h4><p>If <tt>A<\/tt> and <tt>B<\/tt> are two square matrices, the <i>linear matrix pencil<\/i> is the matrix-valued function<\/p><pre>  A - &#955;B<\/pre><p>A pencil is <i>regular<\/i> if there is at least one value of &#955; for which A - &#955;B if not singular.  The pencil is <i>singular<\/i> if both <tt>A<\/tt> and <tt>B<\/tt> are singular and, moreover, A - &#955;B is singular for all &#955;.  In other words,<\/p><pre>  det(A - &#955;B) = 0 for all &#955;.<\/pre><p>Singular pencils are more insiduous than migt appear at first glance.<\/p><h4>Example<a name=\"31c19b00-cfd8-493d-b712-08610ce3c5aa\"><\/a><\/h4><pre class=\"codeinput\">   A = [9 8 7; 6 5 4; 3 2 1]\r\n\r\n   B = [7 9 8; 4 6 5; 1 3 2]\r\n<\/pre><pre class=\"codeoutput\">\r\nA =\r\n\r\n     9     8     7\r\n     6     5     4\r\n     3     2     1\r\n\r\n\r\nB =\r\n\r\n     7     9     8\r\n     4     6     5\r\n     1     3     2\r\n\r\n<\/pre><pre class=\"codeinput\">   syms <span class=\"string\">s<\/span>\r\n\r\n   AB = A - s*B\r\n\r\n   d = det(AB)\r\n<\/pre><pre class=\"codeoutput\"> \r\nAB =\r\n \r\n[9 - 7*s, 8 - 9*s, 7 - 8*s]\r\n[6 - 4*s, 5 - 6*s, 4 - 5*s]\r\n[  3 - s, 2 - 3*s, 1 - 2*s]\r\n \r\n \r\nd =\r\n \r\n0\r\n \r\n<\/pre><pre class=\"codeinput\">   eig1 = eig(A,B)\r\n\r\n   eig2 = 1.\/eig(B,A)\r\n\r\n   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ\r\n<\/pre><pre class=\"codeoutput\">\r\neig1 =\r\n\r\n   -0.4071\r\n    1.0000\r\n    0.2439\r\n\r\n\r\neig2 =\r\n\r\n   -2.0000\r\n    1.0000\r\n    0.3536\r\n\r\n\r\nQAZ =\r\n\r\n   -1.0298  -13.0363    7.7455\r\n         0    5.6991   -4.6389\r\n         0         0    0.0000\r\n\r\n\r\nQBZ =\r\n\r\n    2.4396  -11.4948    9.6394\r\n         0    5.6991   -4.6389\r\n         0         0    0.0000\r\n\r\n<\/pre><pre class=\"codeinput\">   eig3 = eig(A',B')\r\n\r\n   eig4 = 1.\/eig(B',A')\r\n\r\n   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ\r\n<\/pre><pre class=\"codeoutput\">\r\neig3 =\r\n\r\n   -0.2169\r\n       Inf\r\n    1.0000\r\n\r\n\r\neig4 =\r\n\r\n   -0.0738\r\n         0\r\n    1.0000\r\n\r\n\r\nQATZ =\r\n\r\n   -0.0000  -15.0218    6.8390\r\n         0    2.6729   -2.2533\r\n         0         0    0.5922\r\n\r\n\r\nQBTZ =\r\n\r\n    0.0000  -15.2578    7.1280\r\n         0         0    1.0203\r\n         0         0    0.5922\r\n\r\n<\/pre><h4>Wilkinson example<a name=\"0f5d9fb9-0792-4d4c-a9ac-c5fd3e1e7691\"><\/a><\/h4><pre class=\"codeinput\">   clear\r\n\r\n   A = [4 3 2 5; 6 4 2 7; -1 -1 -2 -2; 5 3 2 6]\r\n\r\n   B = [2 1 3 4; 3 3 3 5; 0 0 -3 -2; 3 1 3 5]\r\n<\/pre><pre class=\"codeoutput\">\r\nA =\r\n\r\n     4     3     2     5\r\n     6     4     2     7\r\n    -1    -1    -2    -2\r\n     5     3     2     6\r\n\r\n\r\nB =\r\n\r\n     2     1     3     4\r\n     3     3     3     5\r\n     0     0    -3    -2\r\n     3     1     3     5\r\n\r\n<\/pre><pre class=\"codeinput\">   syms <span class=\"string\">s<\/span>\r\n\r\n   AB = A - s*B\r\n\r\n   d = det(AB)\r\n<\/pre><pre class=\"codeoutput\"> \r\nAB =\r\n \r\n[4 - 2*s,   3 - s, 2 - 3*s, 5 - 4*s]\r\n[6 - 3*s, 4 - 3*s, 2 - 3*s, 7 - 5*s]\r\n[     -1,      -1, 3*s - 2, 2*s - 2]\r\n[5 - 3*s,   3 - s, 2 - 3*s, 6 - 5*s]\r\n \r\n \r\nd =\r\n \r\n0\r\n \r\n<\/pre><pre class=\"codeinput\">   eig1 = eig(A,B)\r\n\r\n   eig2 = 1.\/eig(B,A)\r\n\r\n   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ\r\n<\/pre><pre class=\"codeoutput\">\r\neig1 =\r\n\r\n    1.2056\r\n    0.7055\r\n   -1.0000\r\n      -Inf\r\n\r\n\r\neig2 =\r\n\r\n    1.5097\r\n    0.6408\r\n         0\r\n   -1.0000\r\n\r\n\r\nQAZ =\r\n\r\n    0.7437    4.1769  -12.7279   -5.5000\r\n         0    0.0000    5.2328    2.1602\r\n         0         0    0.7857    0.0123\r\n         0         0         0   -0.2887\r\n\r\n\r\nQBZ =\r\n\r\n    0.5005    6.6143   -8.4853   -2.5000\r\n         0    0.0000    3.2668    2.0105\r\n         0         0    1.1525   -0.7904\r\n         0         0         0    0.2887\r\n\r\n<\/pre><pre class=\"codeinput\">   eig3 = eig(A',B')\r\n\r\n   eig4 = 1.\/eig(B',A')\r\n\r\n   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ\r\n<\/pre><pre class=\"codeoutput\">\r\neig3 =\r\n\r\n  -0.2141 + 0.2033i\r\n  -0.2141 - 0.2033i\r\n   0.7013 + 0.0000i\r\n   1.4508 + 0.0000i\r\n\r\n\r\neig4 =\r\n\r\n    0.3168\r\n    0.9823\r\n    1.2325\r\n         0\r\n\r\n\r\nQATZ =\r\n\r\n   0.1281 - 0.2434i   0.2665 + 0.0169i   0.2663 + 1.4905i   0.3721 + 3.5350i\r\n   0.0000 + 0.0000i   0.0587 + 0.1116i   5.2603 - 1.6197i  12.7878 - 4.0110i\r\n   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   4.1745 + 0.0000i\r\n   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.7572 + 0.0000i\r\n\r\n\r\nQBTZ =\r\n\r\n   0.9052 + 0.0000i   0.6130 - 0.6141i  -0.2443 + 0.8738i   1.2233 + 2.5485i\r\n   0.0000 + 0.0000i   0.4150 + 0.0000i   3.5658 - 1.2114i   8.0696 - 2.2671i\r\n   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   6.6127 + 0.0000i\r\n   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.5220 + 0.0000i\r\n\r\n<\/pre><h4>References<a name=\"cbe4eaf1-8a19-41e7-badd-1b8e69f74b61\"><\/a><\/h4><p>C. B. Moler and G. W. Stewart, \"An Algorithm for Generalized Matrix Eigenvalue Problems\", SIAM J.NUMER.ANAL. Vol.10, No.2, April 1973.  Also available at <a href=\"https:\/\/blogs.mathworks.com\/cleve\/files\/cbm_gws.pdf\">cbm_gws.pdf<\/a><\/p><p>J. H. Wilkinson, Kronecker's Canonical Form and the QZ Algorithm\", LINEAR ALGEBRA AND ITS APPPLICATIONS, Vol. 28, 1979.    Also available at Also available at <a href=\"https:\/\/blogs.mathworks.com\/cleve\/files\/wilkinson.pdf\">wilkinson.pdf<\/a><\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_a9bf4a4405b542598f6dba1c3d53fdbd() {\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='a9bf4a4405b542598f6dba1c3d53fdbd ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a9bf4a4405b542598f6dba1c3d53fdbd';\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 2023 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_a9bf4a4405b542598f6dba1c3d53fdbd()\"><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; R2023a<br><\/p><\/div><!--\r\na9bf4a4405b542598f6dba1c3d53fdbd ##### SOURCE BEGIN #####\r\n%% Singular Matrix Pencils and the QZ Algorithm\r\n% This year, 2023, is the 50-th anniversary of the QZ algorithm\r\n% for generalized matrix eignenvalue problems,\r\n% \r\n%   Ax = \u03bbBx\r\n%\r\n% The algorithm computes these eigevalues without inverting\r\n% either |A| or |B|.  And, the QZ-algorithm can help detect and analyze \r\n% exceptional situaions known as _singular pencils_.\r\n\r\n%% Matrix pencils\r\n% If |A| and |B| are two square matrices, the _linear matrix pencil_\r\n% is the matrix-valued function\r\n%\r\n%    A - \u03bbB\r\n%\r\n% A pencil is _regular_ if there is at least one value of \u03bb for which\r\n% A - \u03bbB if not singular.  The pencil is \r\n% _singular_ if both |A| and |B| are singular and, moreover, A - \u03bbB is \r\n% singular for all \u03bb.  In other words,\r\n% \r\n%    det(A - \u03bbB) = 0 for all \u03bb.\r\n\r\n%%\r\n% Singular pencils are more insiduous than migt appear at first glance.\r\n%% Example\r\n\r\n   A = [9 8 7; 6 5 4; 3 2 1]\r\n\r\n   B = [7 9 8; 4 6 5; 1 3 2]\r\n\r\n%%\r\n\r\n   syms s\r\n\r\n   AB = A - s*B\r\n\r\n   d = det(AB)\r\n\r\n%%\r\n\r\n   eig1 = eig(A,B)\r\n\r\n   eig2 = 1.\/eig(B,A)\r\n\r\n   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ\r\n\r\n\r\n%%\r\n\r\n   eig3 = eig(A',B')\r\n\r\n   eig4 = 1.\/eig(B',A')\r\n\r\n   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ\r\n\r\n%% Wilkinson example\r\n\r\n   clear\r\n\r\n   A = [4 3 2 5; 6 4 2 7; -1 -1 -2 -2; 5 3 2 6]\r\n\r\n   B = [2 1 3 4; 3 3 3 5; 0 0 -3 -2; 3 1 3 5]\r\n\r\n\r\n%%\r\n\r\n   syms s\r\n\r\n   AB = A - s*B\r\n\r\n   d = det(AB)\r\n\r\n%%\r\n\r\n   eig1 = eig(A,B)\r\n\r\n   eig2 = 1.\/eig(B,A)\r\n\r\n   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ\r\n\r\n%%\r\n\r\n   eig3 = eig(A',B')\r\n\r\n   eig4 = 1.\/eig(B',A')\r\n\r\n   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ\r\n\r\n%% References\r\n% C. B. Moler and G. W. Stewart,\r\n% \"An Algorithm for Generalized Matrix Eigenvalue Problems\",\r\n% SIAM J.NUMER.ANAL. Vol.10, No.2, April 1973.  Also available at \r\n% <https:\/\/blogs.mathworks.com\/cleve\/files\/cbm_gws.pdf cbm_gws.pdf>\r\n%\r\n% J. H. Wilkinson,\r\n% Kronecker's Canonical Form and the QZ Algorithm\",\r\n% LINEAR ALGEBRA AND ITS APPPLICATIONS, Vol. 28, 1979.    Also available at Also available at \r\n% <https:\/\/blogs.mathworks.com\/cleve\/files\/wilkinson.pdf wilkinson.pdf>\r\n\r\n##### SOURCE END ##### a9bf4a4405b542598f6dba1c3d53fdbd\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction--><p>This year, 2023, is the 50-th anniversary of the QZ algorithm for generalized matrix eignenvalue problems,... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2023\/01\/05\/singular-matrix-pencils-and-the-qz-algorithm\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/9569"}],"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=9569"}],"version-history":[{"count":3,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/9569\/revisions"}],"predecessor-version":[{"id":9575,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/9569\/revisions\/9575"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=9569"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=9569"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=9569"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}