|
From: | Ben Abbott |
Subject: | Re: polyvalm problems for defective matrix input |
Date: | Fri, 1 Feb 2008 14:00:28 -0500 |
On Jan 30, 2008, at 9:32 AM, Rolf Fabian wrote:
Some of Octave's square matrix functions like 'logm' and 'polyvalm' have the problem that they do not correctly recognize defective matrix input and consequently their output may differ essentially from expectation. Attached zip-file contains an excerpt from AdLA (Advanced Linear Algebra) package which I started to develop in 1997. There are 3 m-files included : 'polyvalmh.m' 'isnormal.m' ( auxilliary function ) 'polyvalmh_demo.m' All script functions are heavily commented and 'polyvalmh' provides two optional switches 'verb' and 'warn' which I've activated prior to upload in order to increase verbosity. If you are interested in the topic and like to see several examples where Octave's output from 'polyvalm' is significantly in error, please run something like: '[y_Octave, y_AdLA] = polyvalmh_demo ( N )' where N = 1, 2, .... 13 . If the demos convince you that 'polyvalmh.m' (incl. auxilliary function 'isnormal.m') represent a better approach to matrix polynomial evaluation than the currently used 'polyvalm.m', you might rename and implement it into Octave, as long as the copyright notice is retained. May be somebody can check the examples also against MatLab, which is not available for me. Rolf Fabian < r dot fabian at jacobs-university dot de >http://www.nabble.com/file/p15183200/AdLA_polyvalmh.zip AdLA_polyvalmh.zip----- Rolf Fabian <r dot fabian at jacobs-university dot de>
I had some trouble with converting to Matlab. However, the results I obtained are below (some errors still persist).
I've attached the edited versions. Please check them. If you can correct the remaining problems, I'll run them for you.
Ben >> polyvalmh_demo(1) EXA %1\nthe "classical", most simple nilpotent + defective matrix x, note: x^ (k>2) = 0 x = 0 1 0 0 p = 1 2 3 \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 3 2 0 3 >> polyvalmh_demo(2) EXA %2\nsimilar to %1, but complex input matrix x = 0 1.0000 + 1.0000i 0 0 p = 1 2 3 \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 3.0000 2.0000 + 2.0000i 0 3.0000 >> polyvalmh_demo(3) EXA %3\nsimple defective regular matrix. \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 6 4 5 0 6 4 0 0 6 >> polyvalmh_demo(4) EXA %4\nsame matrix as %3, scaled to very small norm. \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 3.0000 0.0000 0.0000 0 3.0000 0.0000 0 0 3.0000 >> polyvalmh_demo(5) EXA %5\nmuch more sophisticated nilpotent + defective + matrix x note the extremely large noise in Octave's core function output. x = -0.3568 0.1661 0.0556 0.1435 0.1830 -0.6527 -0.0866 0.7122 0.8908 0.7510 0.1923 -1.1334 -0.3927 -0.7526 -1.0773 -0.8711 0.3456 0.1227 0.4245 1.1486 -0.9055 0.4768 0.4519 0.6874 0.4116 x5 = 1.0e-13 * 0.0167 0.0002 -0.0059 -0.0092 -0.0103 0.0291 0.0056 -0.0058 -0.0097 -0.0132 -0.1130 0.0031 0.0454 0.0694 0.0769 0.0761 -0.0056 -0.0296 -0.0469 -0.0536 0.0427 0.0044 -0.0124 -0.0195 -0.0235 \n ??? Error using ==> power Matrix dimensions must agree. Error in ==> polyvalmh at 80 p = p .* nfro.^(po:-1:0); Error in ==> polyvalmh_demo at 81 yh = polyvalmh (p,x); % AdLA >> polyvalmh_demo(6)EXA %6\na non-trivial, non-nilpotent, regular (det(x)=1) but 'defective' x
x = 0 1 0 0 1 0 2 3 0 0 0 1 0 0 1 0 \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 35 30 220 230 30 35 290 310 0 0 35 30 0 0 30 35 >> polyvalmh_demo(7) EXA %7\nmatrix from example %6 scaled to very large norm. x = 0 1 0 0 1 0 2 3 0 0 0 1 0 0 1 0 polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. \n HANGS ON OCTAVE V3.0.0 'i686-pc-msdosmsvc' (Windows XP) ans = [] >> polyvalmh_demo(8) EXA %6\nmatrix from example %6 scaled to very small norm. x = 0 1 0 0 1 0 2 3 0 0 0 1 0 0 1 0 \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = 11.0000 0.0000 0.0000 0.0000 0.0000 11.0000 0.0000 0.0000 0 0 11.0000 0.0000 0 0 0.0000 11.0000 >> polyvalmh_demo(9) EXA %9\nnon-trivial 16x16, singular and defective x \n polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. ans = Columns 1 through 63.0000 0 0 0 0 0 0 9.0000 + 1.0000i 9.0000 + 1.0000i 0 0 0 0 9.0000 + 1.0000i 9.0000 + 1.0000i 0 0 0 0 0 0 3.0000 0 0 0 0 0 0 9.0000 + 1.0000i 0 27.0000 + 3.0000i 0 0 0 0 15.0000 + 2.0000i 27.0000 + 3.0000i 0 0 0 0 15.0000 + 1.0000i 0 25.0000 + 3.0000i 25.0000 + 3.0000i 0 0 0 0 0 0 0 9.0000 + 1.0000i 0 27.0000 + 3.0000i 0 0 0 0 15.0000 + 1.0000i 27.0000 + 3.0000i 0 0 0 0 12.0000 + 2.0000i 0 25.0000 + 3.0000i 25.0000 + 3.0000i 0 0 0 0 0 0 0 0 0 0 0 0 0 25.0000 + 3.0000i 0 0 0 0 0 25.0000 + 3.0000i 0 32.0000 + 4.0000i 0 0 0 0 19.0000 + 2.0000i
Columns 7 through 120 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9.0000 + 1.0000i 0 0 0 15.0000 + 1.0000i 0 0 15.0000 + 1.0000i 12.0000 + 2.0000i 0 15.0000 + 2.0000i 0 0 12.0000 + 2.0000i 15.0000 + 1.0000i 0 0 9.0000 + 1.0000i 0 0 0 9.0000 + 1.0000i 0 0 9.0000 + 1.0000i 0 0 0 12.0000 + 2.0000i 0 0 15.0000 + 2.0000i 15.0000 + 1.0000i 0 15.0000 + 1.0000i 0 0 15.0000 + 1.0000i 15.0000 + 2.0000i 0 0 9.0000 + 1.0000i 0 0 0 9.0000 + 1.0000i 0 0 0 0 0 0 0 0 25.0000 + 3.0000i 0 0 0 0 0 25.0000 + 3.0000i 0 0 0 19.0000 + 2.0000i 0 0 19.0000 + 2.0000i 19.0000 + 2.0000i 0
Columns 13 through 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3.0000 0 0 0 0 9.0000 + 1.0000i 9.0000 + 1.0000i 0 0 9.0000 + 1.0000i 9.0000 + 1.0000i 0 0 0 0 0 >> polyvalmh_demo(10) EXA %10\nzero input matrix polyvalmh: zero matrix x input: 'p(end).*eye(dim(x))' returned. \n ans =0 -10.0000i 0 0 0 0 0 0 -10.0000i 0 0 0 0 0 0 -10.0000i 0 0 0 0 0 0 -10.0000i 0 0 0 0 0 0 -10.0000i
>> polyvalmh_demo(11)EXA %11\nspeed test using a complex, 'normal' matrix of larger size [200x200] As a 'normal' matrix example, we use a unitary mat resulting from a random complex schur decomposition.
t_polyvalm_tictoc = 0.1535 t_polyvalm_cputime = 0.1600 ??? Error using ==> power Matrix dimensions must agree. Error in ==> polyvalmh at 80 p = p .* nfro.^(po:-1:0); Error in ==> polyvalmh_demo at 162 yo = polyvalmh (p,x); % AdLA >> polyvalmh_demo(12) EXA %12\ninput matrix contains non-finite elements (polyvalm) leads to an error! x = 1 Inf 0 NaN ans = NaN NaN NaN NaN >> polyvalmh_demo(13) EXA %13\ninput matrix contains non-finite elements (polyvalmh) leads to an error! x = 1 Inf 0 NaN ??? NaN's cannot be converted to logicals. Error in ==> polyvalmh at 55 if (nfro) Error in ==> polyvalmh_demo at 179 yh = polyvalmh (p,x); >>
polyvalmh_demo.m
Description: Binary data
isnormal.m
Description: Binary data
issquare.m
Description: Binary data
polyvalmh.m
Description: Binary data
[Prev in Thread] | Current Thread | [Next in Thread] |