|

»

»

Octave - Maintainers

polyvalm problems for defective matrix input

View: New views

3 Messages — Rating Filter: Alert me

	polyvalm problems for defective matrix input by Rolf Fabian :: Rate this Message: Reply to Author \| View Threaded \| Show Only this Message 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 > AdLA_polyvalmh.zip Rolf Fabian <r dot fabian at jacobs-university dot de>
	Re: polyvalm problems for defective matrix input by Ben Abbott :: Rate this Message: Reply to Author \| View Threaded \| Show Only this Message 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 6 3.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 12 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 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 (11K) Download Attachment isnormal.m (4K) Download Attachment issquare.m (123 bytes) Download Attachment polyvalmh.m (7K) Download Attachment
	Re: polyvalm problems for defective matrix input by Rolf Fabian :: Rate this Message: Reply to Author \| View Threaded \| Show Only this Message Many many thanks to Ben Abbott who transfered the 'polyvalmh_demo.m' test script to MatLab and tested and ran it also under ML ( his results are labelled 'MATLAB RESULT' in the overview below. In summary all the matlab results for the examples given by Ben confirm the results returned by proposed 'polyvalmh.m' and discard all the results from current Octave V3.0.0 'polyvalm.m'. The most outstanding difference can be seen in EXA #3, where MatLab and 'polyvalmh.m' return a value of '5' for the upper right corner element, whereas Octave's polyvalm returns 1.8014e+016 for that element, which is more than only a bit different. EXA #6 shows that a non-trivial, non-nilpotent, regular but 'defective' x may also lead to subtle large differences between Octave's 'polyvam.m' result and both 'polyvalmh.m' and MatLab polyvalm. Again MatLab confirms the result of 'polyvalmh'. The MatLab results strongly support my suggestion to replace Octave's current 'polyvalm.m' by a renamed 'polyvalmh.m'. I already pointed out that you can do this PROVIDED that my name appears as author in the copyright notice. You might wonder whether I do insist so much on this point ? The answer is very simple. Actually it has been me who has written most of the code of the current Octave 'polyvalm.m' function. The original version by Tony Richardson (dated 1994) failed horribly without warning for even the most simple matrices. It was essentially a '3 line thing' defining only a unitary diagonalization algorithm to do polyvalm. Thus it was actually only applicable to the special case of hermitian (more precisely 'normal) input matrices and persistently failed for the more common general case, which requires a general (not only unitary) similarity transformation algorithm. Thus I decided on Aug 19, 1997 to send a bug report <http://velveeta.che.wisc.edu/octave/lists/archive//bug-octave.1997/msg00362.html> where I suggested a complete revision of 'polyvalm' for Octave 2.0.9 . The reply by jwe dated Aug. 28, 2997 <http://velveeta.che.wisc.edu/octave/lists/archive//bug-octave.1997/msg00388.html accepts the proposed function. If you have a look at the 1st URL you'll easily see that the current Octave V3.0.0 'polyvalm.m' corresponds nearly exactly character-by-character to my proposed function from Aug 19, 1997. Besides some changes in the help part (texinfo-related) The only essential difference is that you simply removed my autorship notice which I've written aside to the authorship notice of Tony (check my authorship notice in the URL). This treatment wasn't really polite and even though I don't know much about licenses I wonder whether such a policy can be conform to the terms of the GNU public license. I'm pretty much shure that such a policy concerning authorship on free software is discouraging other contributors also, not only me. Now that I've also solved the even more sophisticated problem of handling defective matrix input correctly the function doesn't have to do anything anymore with the '3-lines-code' of Tony. My offer : You can have the code, but I insist on getting single authorship. Rolf Fabian ======================================================= OVERVIEW TEST DEMOS ======================================================= octave-3.0.0.exe:1> [y_Octave,y_polyvalmh]=polyvalmh_demo(1) EXA #1 the 'classical', most simple nilpotent + defective matrix x [note: x^ (k > 2) = 0] x = 0 1 0 0 p = 1 2 3 warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 5e-293. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 3 0 0 3 y_polyvalmh = 3 2 0 3 MATLAB RESULT :> polyvalm(p,x) ans = 3 2 0 3 ======================================================= octave-3.0.0.exe:2> [y_Octave,y_polyvalmh]=polyvalmh_demo(2) EXA #2 similar to #1, but complex input matrix x = 0 + 0i 1 + 1i 0 + 0i 0 + 0i p = 1 2 3 warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 5e-293. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 3 0 0 3 y_polyvalmh = 3 + 0i 2 + 2i 0 + 0i 3 + 0i MATLAB RESULT :> polyvalm(p,x) ans = 3.0000 2.0000 + 2.0000i 0 3.0000 ======================================================= octave-3.0.0.exe:3> [y_Octave,y_polyvalmh]=polyvalmh_demo(3) EXA #3 simple defective regular matrix. p = 1 2 3 x = 1 1 1 0 1 1 0 0 1 warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 2.5e-032. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 6.0000e+000 0.0000e+000 1.8014e+016 0.0000e+000 6.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 6.0000e+000 y_polyvalmh = 6 4 5 0 6 4 0 0 6 MATLAB RESULT :> polyvalm(p,x) ans = 6 4 5 0 6 4 0 0 6 ======================================================= octave-3.0.0.exe:4> [y_Octave,y_polyvalmh]=polyvalmh_demo(4) EXA #4 same matrix as #3, scaled to very small norm. p = 1 2 3 scalefac = 1.0000e-060 warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 2.5e-032. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 3.0000e+000 0.0000e+000 9.0072e+015 0.0000e+000 3.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 3.0000e+000 y_polyvalmh = 3.00000 0.00000 0.00000 0.00000 3.00000 0.00000 0.00000 0.00000 3.00000 MATLAB RESULT :> polyvalm(p,x) ans = 3.0000 0.0000 0.0000 0 3.0000 0.0000 0 0 3.0000 ======================================================= octave-3.0.0.exe:5> [y_Octave,y_polyvalmh]=polyvalmh_demo(6) EXA #6 a 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 p = 2 3 4 5 6 7 8 9 10 11 {message from Octave for 'polyvalm(p,x)'} warning: matrix singular to machine precision, rcond = 2.22045e-017 warning: attempting to find minimum norm solution warning: dgelsd: rank deficient 4x4 matrix, rank = 3 warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 3.7e-017. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 3.5000e+001 3.2527e+001 5.2000e+001 -2.5000e+001 3.0000e+001 3.5355e+001 5.6000e+001 -2.8000e+001 -1.9984e-015 -3.2658e-015 5.3975e+000 -2.9996e+000 -8.8818e-016 -7.5364e-016 -2.7014e+000 1.5013e+000 y_polyvalmh = 35.00000 30.00000 220.00000 230.00000 30.00000 35.00000 290.00000 310.00000 0.00000 0.00000 35.00000 30.00000 0.00000 0.00000 30.00000 35.00000 MATLAB RESULT :> polyvalm(p,x) ans = 35 30 220 230 30 35 290 310 0 0 35 30 0 0 30 35 ======================================================= octave-3.0.0.exe:6> [y_Octave,y_polyvalmh]=polyvalmh_demo(8) EXA #6 matrix 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 p = 2 3 4 5 6 7 8 9 10 11 scalefac = 1.0000e-040 {message from Octave for 'polyvalm(p,x)'} warning: matrix singular to machine precision, rcond = 2.22045e-017 warning: attempting to find minimum norm solution warning: dgelsd: rank deficient 4x4 matrix, rank = 3 warning: polyvalmh: numeric underflow probable. Result should be checked. warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 3.7e-017. polyvalmh: 'Horner algorithm' used for 'defective' input. y_Octave = 1.1000e+001 2.4749e+000 1.2000e+001 -1.7500e+000 -4.4409e-016 8.1317e+000 4.0000e+000 1.0000e+000 -1.4655e-015 -2.8262e-015 8.9933e+000 -4.1586e-001 9.7700e-016 2.1981e-015 -8.4660e+000 3.9148e-001 y_polyvalmh = 11.00000 0.00000 0.00000 0.00000 0.00000 11.00000 0.00000 0.00000 0.00000 0.00000 11.00000 0.00000 0.00000 0.00000 0.00000 11.00000 MATLAB RESULT :> polyvalm(p,x) 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 ======================================================= Rolf Fabian <r dot fabian at jacobs-university dot de>