Re: polyvalm problems for defective matrix input

octave-maintainers

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: polyvalm problems for defective matrix input

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.zipAdLA_polyvalmh.zip


-----
Rolf Fabian
<r dot fabian at jacobs-university dot de>

I had some trouble with converting to Matlab. However, the results Iobtained are below (some errors still persist).

I've attached the edited versions. Please check them. If you cancorrect 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 00 0 00 9.0000 + 1.0000i 9.0000 + 1.0000i0 0 00 9.0000 + 1.0000i 9.0000 + 1.0000i0 0 00 0 03.0000 0 00 0 00 9.0000 + 1.0000i 027.0000 + 3.0000i 0 00 0 15.0000 + 2.0000i27.0000 + 3.0000i 0 00 0 15.0000 + 1.0000i0 25.0000 + 3.0000i 25.0000 + 3.0000i0 0 00 0 00 9.0000 + 1.0000i 027.0000 + 3.0000i 0 00 0 15.0000 + 1.0000i27.0000 + 3.0000i 0 00 0 12.0000 + 2.0000i0 25.0000 + 3.0000i 25.0000 + 3.0000i0 0 00 0 00 0 00 0 00 25.0000 + 3.0000i 00 0 00 25.0000 + 3.0000i 032.0000 + 4.0000i 0 00 0 19.0000 + 2.0000i


 Columns 7 through 12

0 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 9.0000 + 1.0000i0 0 015.0000 + 1.0000i 0 0 15.0000 +1.0000i 12.0000 + 2.0000i 015.0000 + 2.0000i 0 0 12.0000 +2.0000i 15.0000 + 1.0000i 00 9.0000 + 1.0000i 00 0 9.0000 + 1.0000i0 0 9.0000 + 1.0000i0 0 012.0000 + 2.0000i 0 0 15.0000 +2.0000i 15.0000 + 1.0000i 015.0000 + 1.0000i 0 0 15.0000 +1.0000i 15.0000 + 2.0000i 00 9.0000 + 1.0000i 00 0 9.0000 + 1.0000i0 0 00 0 00 0 25.0000 + 3.0000i0 0 00 0 25.0000 + 3.0000i0 0 019.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 00 00 0 -10.0000i 00 00 0 0 -10.0000i0 00 0 0 0-10.0000i 00 0 00 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 arandom 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]

Re: polyvalm problems for defective matrix input, Ben Abbott <=
- Re: polyvalm problems for defective matrix input, Rolf Fabian, 2008/02/04

Prev by Date: Re: Polyfit with scaling
Next by Date: Re: Polyfit with scaling
Previous by thread: Re: Polyfit with scaling
Next by thread: Re: polyvalm problems for defective matrix input
Index(es):
- Date
- Thread