polyvalm problems for defective matrix input
Ben Abbott
bpabbott at mac.com
Fri Feb 1 13:00:28 CST 2008
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);
>>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polyvalmh_demo.m
Type: application/octet-stream
Size: 6570 bytes
Desc: not available
Url : https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20080201/8c401602/attachment-0004.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: isnormal.m
Type: application/octet-stream
Size: 2758 bytes
Desc: not available
Url : https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20080201/8c401602/attachment-0005.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: issquare.m
Type: application/octet-stream
Size: 114 bytes
Desc: not available
Url : https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20080201/8c401602/attachment-0006.obj
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polyvalmh.m
Type: application/octet-stream
Size: 4666 bytes
Desc: not available
Url : https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20080201/8c401602/attachment-0007.obj
More information about the Octave-maintainers
mailing list