Inverse of a matrix

jlobillo · Mon Mar 22, 2010 12:40 pm

I think I've found a bug. When I execute the following code

Code:

                     SINGULAR                             /
 A Computer Algebra System for Polynomial Computations   /   version 3-1-0
                                                       0<
     by: G.-M. Greuel, G. Pfister, H. Schoenemann        \   Mar 2009
FB Mathematik der Universitaet, D-67653 Kaiserslautern    \
> LIB "matrix.lib";
// ** loaded /usr/share/Singular/LIB/matrix.lib (1.48,2009/04/10)
// ** loaded /usr/share/Singular/LIB/nctools.lib (1.54,2009/05/08)
// ** loaded /usr/share/Singular/LIB/poly.lib (1.53,2009/04/15)
// ** loaded /usr/share/Singular/LIB/general.lib (1.62,2009/04/15)
// ** loaded /usr/share/Singular/LIB/random.lib (1.20,2009/04/15)
// ** loaded /usr/share/Singular/LIB/ring.lib (1.34,2009/04/15)
// ** loaded /usr/share/Singular/LIB/primdec.lib (1.147,2009/04/15)
// ** loaded /usr/share/Singular/LIB/absfact.lib (1.7,2008/07/16)
// ** loaded /usr/share/Singular/LIB/triang.lib (1.14,2009/04/14)
// ** loaded /usr/share/Singular/LIB/elim.lib (1.34,2009/05/05)
// ** loaded /usr/share/Singular/LIB/inout.lib (1.34,2009/04/15)
> ring cuerpobase=(0,a),(dummy),dp; minpoly=a4-10a2+1;
> def sqrt6 = 1/2a2-5/2;
> def sqrt2 = 1/2a3-9/2a;
> def sqrt3 = a - sqrt2;
> vector v1=[0,0,1];
> vector v2=[0,2/3*sqrt2,-1/3];
> vector v3=[-1/3*sqrt6,-1/3*sqrt2,-1/3];
> vector v4=[1/3*sqrt6,-1/3*sqrt2,-1/3];
> matrix basesaux[3][3]=matrix(v1),matrix(v2),matrix(v3);
> matrix b1bc=transpose(basesaux);
> power(b1bc,-1);
_[1,1]=1
_[1,2]=0
_[1,3]=0
_[2,1]=0
_[2,2]=1
_[2,3]=0
_[3,1]=0
_[3,2]=0
_[3,3]=1
>

the result can't be correct since the matrix b1bc is not the unit matrix. Am I doing anything wrong?

Best,

Javier

Re: Inverse of a matrix

gorzel · Mon Mar 22, 2010 4:31 pm

Take a look, whether power from matrix.lib
is only limited to positive powers.

Use inverse from linalg.lib instead.

C. Gorzel

Re: Inverse of a matrix

seelisch · Tue Mar 23, 2010 11:10 am

Yes, 'inverse' of linalg.lib works; see below.
Indeed, 'power' of matrix.lib expects an exponent >= 0. (It is then also clear why the unit matrix is returned, because the loop inside 'power' terminates immediately with the power computed so far which is just the unit matrix.)

$ ./Singular
SINGULAR / Development
A Computer Algebra System for Polynomial Computations / version 3-1-1
0<
by: G.-M. Greuel, G. Pfister, H. Schoenemann \ Feb 2010
FB Mathematik der Universitaet, D-67653 Kaiserslautern \
// ** executing /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/.singularrc
> LIB "linalg.lib";
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/linalg.lib (12258,2009-11-09)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/general.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/poly.lib (12443,2010-01-19)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/inout.lib (12541,2010-02-09)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/elim.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/ring.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/primdec.lib (12617,2010-03-08)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/absfact.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/triang.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/random.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/matrix.lib (12231,2009-11-02)
// ** loaded /cygdrive/c/Work/SINGULAR-Branches/SINGULAR-for-fans/Singular/LIB/nctools.lib (12231,2009-11-02)
> ring cuerpobase=(0,a),(dummy),dp; minpoly=a4-10a2+1;
> def sqrt6 = 1/2a2-5/2;
> def sqrt2 = 1/2a3-9/2a;
> def sqrt3 = a - sqrt2;
> vector v1=[0,0,1];
> vector v2=[0,2/3*sqrt2,-1/3];
> vector v3=[-1/3*sqrt6,-1/3*sqrt2,-1/3];
> vector v4=[1/3*sqrt6,-1/3*sqrt2,-1/3];
> matrix basesaux[3][3]=matrix(v1),matrix(v2),matrix(v3);
> matrix b1bc=transpose(basesaux);
> inverse(b1bc);
_[1,1]=(-1/8*a^2+5/8)
_[1,2]=(1/8*a^3-9/8*a)
_[1,3]=1
_[2,1]=(-1/8*a^2+5/8)
_[2,2]=(3/8*a^3-27/8*a)
_[2,3]=0
_[3,1]=(-1/4*a^2+5/4)
_[3,2]=0
_[3,3]=0
> inverse(b1bc)*b1bc; // just checking...
_[1,1]=1
_[1,2]=0
_[1,3]=0
_[2,1]=0
_[2,2]=1
_[2,3]=0
_[3,1]=0
_[3,2]=0
_[3,3]=1
>

Re: Inverse of a matrix

jlobillo · Wed Mar 24, 2010 4:56 pm

Your answers work perfectly, however the documentation on power says that power works with integer exponents, so the manual must be updated.

And I think I've computed inverses using power(A,-1) in a different setting. If I found the code I'll post it here.

Thanks