The MatrixObj project

Max Horn

Matrices in GAP have some issues...

They usually...

  1. ... are assumed to be "lists of row vectors";
  2. ... don't "know" their base ring;
  3. ... don't "know" their dimensions;
  4. ... can incur massive method dispatch overhead.

Let's look at these a bit closer!

The GAP manual on matrices

  • "Matrices are represented in GAP by lists of row vectors [...]" 
    
    gap> mat:=[[1,2,3],[4,5,6],[7,8,9]];;
    gap> IsMatrix(mat);
    true
  • "The vectors must all have the same length, and their elements must lie in a common ring."
  • For performance reasons, this is not enforced. 
    
    gap> mat:=[[1,2],[3]];;
    gap> IsMatrix(mat);
    true
    gap> IsRectangularTable(mat);
    false
    gap> DimensionsMat(mat);
    fail

GAP matrices have no base ring

  • What's the base ring of this matrix? 
    
    gap> mat:=[[0,1],[2,3]];;
  • Affects questions like "Is this matrix invertible?" ("yes" over $\mathbb{Q}$, "no" over $\mathbb{Z}$) 
    
    gap> mat^-1;
    [ [ -3/2, 1/2 ], [ 1, 0 ] ]
  • Some operations, like Eigenvalues, at least allow you to specify a base field or ring.

... no base ring [cont.]

  • Perhaps we can ask GAP? Hmmm, it kinda assumes matrices are over fields: 
    
    gap> DefaultFieldOfMatrix(mat);
    Rationals
  • Alternatively, let 
    
    gap> R:=DefaultRing(mat); LeftActingDomain(R);
    ( Integers^[ 2, 2 ] )
    Integers
  • Base "field" is to restrictive; perhaps "ring" also?
  • To compute "the" (or rather "a") base ring, GAP may need to inspect all matrix entries

Fun with empty matrices!


        gap> id:=IdentityMat(3);;
        gap> m1:=[];;
        gap> m2:=[[]];; 
        gap> m3:=EmptyMatrix(0);
        EmptyMatrix( 0 )
        gap> m4:=[ EmptyRowVector(FamilyObj(42)) ];
        [ <empty row vector> ]
        gap> m5:=IdentityMat(0);  # NullMapMatrix 
        <null map>
        

        gap> TransposedMat(m1);
        [  ]
        gap> TransposedMat(m2);
        [  ]
        gap> TransposedMat(m3);
        [  ]
        gap> TransposedMat(m4);
        [  ]
        gap> TransposedMat(m5);
        [  ]

Fun with empty matrices!


        gap> id:=IdentityMat(3);;
        gap> m1:=[];;
        gap> m2:=[[]];; 
        gap> m3:=EmptyMatrix(0);
        EmptyMatrix( 0 )
        gap> m4:=[ EmptyRowVector(FamilyObj(42)) ];
        [ <empty row vector> ]
        gap> m5:=IdentityMat(0);  # NullMapMatrix 
        <null map>
        

        gap> DimensionsMat(m1);
        Error, no method found!
        gap> DimensionsMat(m2);
        Error, no method found!
        gap> DimensionsMat(m3);
        [ 0, 0 ]
        gap> DimensionsMat(m4);
        [ 1, 0 ]
        gap> DimensionsMat(m5);
        Error, no method found!

Fun with empty matrices!


        gap> id:=IdentityMat(3);;
        gap> m1:=[];;
        gap> m2:=[[]];; 
        gap> m3:=EmptyMatrix(0);
        EmptyMatrix( 0 )
        gap> m4:=[ EmptyRowVector(FamilyObj(42)) ];
        [ <empty row vector> ]
        gap> m5:=IdentityMat(0);  # NullMapMatrix 
        <null map>
        

        gap> id*m1;
        Error, Inner product multiplication of lists: no summands
        gap> id*m2;
        [ [  ], [  ], [  ] ]
        gap> id*m3;
        EmptyMatrix( 0 )
        gap> id*m4;
        [ <empty row vector>, <empty row vector>, <empty row vector> ]
        gap> id*m5;
        [ [  ], [  ], [  ] ]

Fun with empty matrices!


        gap> id:=IdentityMat(3);;
        gap> m1:=[];;
        gap> m2:=[[]];; 
        gap> m3:=EmptyMatrix(0);
        EmptyMatrix( 0 )
        gap> m4:=[ EmptyRowVector(FamilyObj(42)) ];
        [ <empty row vector> ]
        gap> m5:=IdentityMat(0);  # NullMapMatrix 
        <null map>
        

        gap> m1*id;
        Error, Inner product multiplication of lists: no summands
        gap> m2*id;
        Error, Inner product multiplication of lists: no summands
        gap> m3*id;
        EmptyMatrix( 0 )
        gap> m4*id;
        Error, no method found!
        gap> m5*id;
        <null map>

Method dispatch overhead

  • Quiz: why are there both First and FirstOp? 

    gap> list := List([1..1000], i -> i^2);;
    gap> for i in [1..10^4] do First(list,ReturnFalse); od; time;
    205
    gap> for i in [1..10^4] do FirstOp(list,ReturnFalse); od; time;
    209
    

    gap> list := List([1..1000], i -> [i^2]);;
    gap> for i in [1..10^4] do First(list,ReturnFalse); od; time;
    208
    gap> for i in [1..10^4] do FirstOp(list,ReturnFalse); od; time;
    499
    

    gap> list := List([1..1000], i -> [[i^2]]);;
    gap> for i in [1..10^4] do First(list,ReturnFalse); od; time;
    201
    gap> for i in [1..10^4] do FirstOp(list,ReturnFalse); od; time;
    1227

More on dispatch overhead


    gap> list := List([1..1000], i -> [[i^2]]);;
    gap> TNUM_OBJ(list);
    [ 22, "list (plain)" ]
    

    gap> KnownTruePropertiesOfObject(list);
    [ "IsFinite", "IsSmallList", "IsRectangularTable" ]
    

    gap> TNUM_OBJ(list);
    [ 26, "list (plain,dense)" ]
    

    gap> MakeImmutable(list);;
    gap> TNUM_OBJ(list);
    [ 27, "list (plain,dense,imm)" ]
    

    gap> KnownTruePropertiesOfObject(list);
    [ "IsFinite", "IsSmallList", "IsRectangularTable" ]
    gap> TNUM_OBJ(list);
    [ 49, "list (plain,rect table,imm)" ]
    

    gap> for i in [1..10^4] do First(list,ReturnFalse); od; time;
    210
    gap> for i in [1..10^4] do FirstOp(list,ReturnFalse); od; time;
    202

List-of-rows is too restrictive

There are many ways to store a matrix, e.g.

  • as a list of row vectors;
  • as a list of column vectors;
  • as a flat list with $rows \times columns$ entries;
  • as a "sparse" matrix; e.g.
    • a list of tuples $[i,j,a_{ij}]$ for non-zero entries
  • as a lazily computed matrix
  • ...

Solution:

A dedicated matrix (and vector) interface