search

gap-system / gap

Main development repository for GAP - Groups, Algorithms, Programming, a System for Computational Discrete Algebra
https://www.gap-system.org
GNU General Public License v2.0
804 stars 163 forks source link

`IsAssociativeElement`: Documentation wrong? Bug? #2776

Open markuspf opened 6 years ago

markuspf commented 6 years ago

Documentation

31.15-1 IsAssociativeElement

  ‣ IsAssociativeElement( obj )  Category
  ‣ IsAssociativeElementCollection( obj )  Category
  ‣ IsAssociativeElementCollColl( obj )  Category

  An  element  obj  in the category IsAssociativeElement knows that the
  multiplication of any elements in the family of obj is associative. For
  example, all permutations lie in this category, as well as those ordinary
  matrices (see IsOrdinaryMatrix (24.2-2)) whose entries are also in
  IsAssociativeElement.

GAP master

gap> M := [[Z(4)]];                                                   
[ [ Z(2^2) ] ]
gap> IsOrdinaryMatrix(M);
true
gap> IsAssociativeElementCollection(M);
false
gap> IsAssociativeElementCollection( [ M ] );
false

and also

gap> M := Matrix(GF(4), [[Z(4)]]);
[ [ Z(2^2) ] ]
gap> IsOrdinaryMatrix(M);
true
gap> IsAssociativeElementCollection(M);
false
fingolfin commented 6 years ago

IsAssociativeElementCollection(M) = false; seems correct to me. Moreover, note that IsAssociativeElement(M)=true holds.

I agree that IsAssociativeElementCollection( [ M ] ) = false is surprising. However, I do not strictly see an error in the documentation you quote. Could you specify which sentence in it you consider wrong?

But let's look at the technical reasons for this, and what we might be able to do about it. You give two examples, but in master, both M are really the same, even that returned by Matrix (and there are reasons for that, but I'll get back to that below).

The problem is how we compute types of lists (and hence whether they are in a given filter or not). Clearly it is desirable to have this fast. For proper objects, this is no problem, as the type is stored. But for lists, the type has to be recomputed each time we ask for it.

So to make it fast, for homogenous lists, we essentially only take the families of the elements in to account. Consider this snippet:

gap> v:=[Z(4)];;  # a row vector
gap> M:=[v];; # = [[Z(4)]],  a "rectangular table" of ring elements
gap> L:=[M];;
gap>
gap> FamilyObj(v) = CollectionsFamily(FamilyObj(Z(4)));
true
gap> FamilyObj(M) = CollectionsFamily(FamilyObj(v));
true
gap> FamilyObj(L) = CollectionsFamily(FamilyObj(M));
true

Now, the type of the list of course inherits all filters implied by its family. But it has more: namely any filters that come from implications. For the example at hand, the following declarations and implications are relevant:

DeclareSynonym( "IsAdditiveElementWithInverseTable",
    IsAdditiveElementWithInverseCollColl   and IsTable );

DeclareSynonym( "IsRingElementTable",
        IsAdditiveElementWithInverseTable
    and IsMultiplicativeElementTable );

DeclareSynonym( "IsMatrix", IsRingElementTable );

InstallTrueMethod( IsOrdinaryMatrix, IsMatrix and IsInternalRep );

InstallTrueMethod( IsAssociativeElement, IsOrdinaryMatrix and IsAssociativeElementCollColl );

The family object of M, derived recursively as described above, implies the filter IsAssociativeElementCollColl for M. But it does not imply IsOrdinaryMatrix; instead, that comes ultimately from the fact that the kernel keeps track of lists which are actually (rectangular) "tables". So together this ensures that IsAssociativeElement(M)=true.

But crucially, the IsTable bit is not part of the family of M -- and it can't be, because of the way those families are derived: of course there are lists-of-lists-of-FFEs which are not rectangular tables, so the their family cannot imply that.

In order to get the type right here, GAP would have to iterate over all list elements, compute their types, then "intersect" the filters/flags they imply, and use that. That's a lot more expensive for a homogenous list than to just consider the type of a single element in the list.

Now, truth be told, in the above example, we actually are iterating over all elements -- because those lists (and lists of lists, and so on) all are mutable, so their types are not fully stored, and thus have to be recomputed every. single. time. we interact with these lists (which is a nightmare of its own, and we discussed this plenty before). But while we iterate over everything, we try to do so quite quickly; computing the full type of everything and intersecting flags etc. etc. would make it a lot slower. So this is not a viable route to go.

In theory, we could add a special case to the kernel here: for homogenous lists, inspect the type of one elements, and check for certain filters; if set, add the corresponding CategoryCollections(FILTER) (or so) to the list type we are computing.

But that's still a rather nasty hack (and I am not even sure it actually would work the way I think it might). I have doubts whether it'd be really worth it.

Which is of course one of the many reasons why we want "real" matrix and vector objects.

Which gets us back to my original remark: Why does Matrix in master still return a MatrixObj.

The reasons for this really are: no fully function MatrixObj and VectorObj types are available yet; indeed there are still lots of details to be fleshed out; existing implementations are slow or incomplete; and also, anything we put out "officially" risks becoming a hard depdenency we can never get rid of again.

Moreover, in long discussions on whether ...

  1. IsMatrixObj should imply IsMatrix;
  2. IsMatrix should IsMatrixObj;
  3. neither should imply the other;

it was finally determined that option 2 was most beneficial, and this has now been implemented in the library. But to be honest, I can't fully get the argument back together, gotta look at notes when I have had slept (and perhaps @ThomasBreuer can chime in). Part of it certainly was to make a gradual transition of existing code possible. However, now that I am writing this, I am once again wondering about it... ah well...

Anyway, from a pragmatic point of view, what we need is a highly optimized and complete "generic" MatrixObj implementation that works over arbitrary base domains, to replace uses of lists-of-lists (and the same for vectors, of course). Then, we need matrix/vector implementations for various specialized domains. For small finite fields, we of course have our compressed matrices; but there are various bottlenecks in them that should be addressed (and for compatibility reasons, we might instead want to add a new type, and leave the existing compressed matrix code for legacy support only... I could explain more of why I think that might be best, but I am already digressing too far). In addition to this, the API needs to be fleshed out more, but the best way to do that is to actually write code that uses MatrixObj, to see what APIs are really needed (to avoid inventing ones that turn out useless in practice). That in turn requires close feedback from users, including the semigroups authors! It's a feedback loop, and that's why I pushed hard to get our MatrixObj work at least in 4.10, incomplete and bare as it is...

Anyway, regarding performance: We could make IsPlistMatrixRep the default representation for Matrix, but that type has not (yet?) been optimized in any way and is painfully slow, even slower than lists-of-lists, due to object overhead (but this could be improved). Part of the problem is that it was written in a time when it was thought that all matrices should be row-major (we decided to be flexible about this), and when it was deemed sensible to require MatrixObj to support row access (i.e., M[1] returns the first row as a vector obj) -- the latter prevents many optimizations, both for memory and performance; e.g. the user of kernel functions like PROD_LIST_SCL_DEFAULT. And for compressed matrices, it causes extra overhead: instead of storing a 10x10 matrix over Z(2) in 100 bits, plus overhead, we store both rows are stored as separate, full GAP objects, for a whopping 632 bytes in total:

gap> id:=ImmutableMatrix(GF(2),IdentityMat(10,GF(2)));
<an immutable 10x10 matrix over GF2>
gap> MemoryUsage(id);
632
gap> MemoryUsage(id[1]);
48

If this was stored in a single GAP object, it'd need about 56 bytes of storage (24 for GAP meta data, 8 for a pointer to a type object, 4+4 for row/column length, 13 for the actual data, 3 for padding); perhaps some more in the end if we need to store some more meta data. So a factor of 10 in storage, and also many advantage in performance.

This is part of the reasons why I think we need to overhaul our compressed matrix type, or implement yet another one: it's better to store the matrix more compactly (also for memory locality), and use proxy objects to fake support for id[1] to get the first row (of course new code can use the right APIs to avoid proxy objects).

To give an idea about the performance of IsPlistMatrixRep:

gap> M1:=IdentityMatrix(GF(4),10);
< mutable compressed matrix 10x10 over GF(4) >
gap> M2:=IdentityMat(10)*Z(2);
[ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ]
gap> M3:=IdentityMatrix(IsPlistMatrixRep, GF(4),10);
<10x10-matrix over GF(2^2)>

gap> for i in [1..1000] do x:=M1*M1; od; time;
2
gap> for i in [1..1000] do x:=M2*M2; od; time;
34
gap> for i in [1..1000] do x:=M3*M3; od; time;
129

gap> for i in [1..1000] do x:=M1^5+M1; od; time;
8
gap> for i in [1..1000] do x:=M2^5+M2; od; time;
115
gap> for i in [1..1000] do x:=M3^5+M3; od; time;
346

and over the integers:

gap> I1:=IdentityMat(10);
[ [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 ],
  [ 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 ], [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ] ]
gap> I2:=IdentityMatrix(IsPlistMatrixRep,Integers,10);
<10x10-matrix over Integers>

gap> for i in [1..1000] do x:=I1*I1; od; time;
8
gap> for i in [1..1000] do x:=I2*I2; od; time;
42

gap> for i in [1..1000] do x:=I1^5+I1; od; time;
21
gap> for i in [1..1000] do x:=I2^5+I2; od; time;
126

This could be greatly improved if those IsPlistMatrixRep stored their contents simply as lists of lists, instead of lists-of-VectorObj... That is certainly doable, just nobody worked on it yet.

fingolfin commented 6 years ago

PS: I just spent a few minutes on optimizing IsPlistMatrixRep as indicated, and I got extremely encouraging results (note that I used 10 times as many iterations compared to the previous numbers). 

gap> for i in [1..10000] do x:=I1*I1; od; time;
58
gap> for i in [1..10000] do x:=I2*I2; od; time;
65

gap> for i in [1..10000] do x:=I1^5+I1; od; time;
219
gap> for i in [1..10000] do x:=I2^5+I2; od; time;
245

I'll try to make a PR out of this, and perhaps we can switch Matrix to using this, instead of lists-of-lists, but probably not for 4.10, as I expected many bugs and missing features in IsPlistMatrixRep. So this does not help with issue #2737 in GAP 4.10, but perhaps in GAP 4.11

UPDATE: used slightly different timings; those from my first run were looking overly good by chance; these numbers fluctuate quite a bit, gotta dig out my benchmarking code for better numbers when fine tuning this.

fingolfin commented 6 years ago

Oh, and of course for highest performance (beating that of plists) in the generic case, one could implement a T_MATRIX kernel type which is fixed-dimension (nxm) matrix. Just to ensure that no code gets slower by using MatrixObj...

markuspf commented 6 years ago

Disregard my comment about IsAssociativeElement, I shouldn't be opening issues when too tired. I was under the impression that IsAssociativeElement does not hold for the matrix M above. I should also start logging my GAP sessions so I can look back at what I did (I might have been playing with filters in semigroups and mucked something up).

I do think that there is some sense in a matrix being an IsAssociativeElementCollection, after all a matrix is a collection of associative elements (if the multiplication of entries is associative). Whether this is a useful filter for applications, I don't know. I'll read the rest of your post now and get back then.

ThomasBreuer commented 6 years ago

Since I was mentioned in some of the comments, just a few remarks.

stevelinton commented 6 years ago

There are two underlying problems here: (1) a matrix with mutable rows cannot remember anything which might be changed if one of the rows is mutated (in particular, in this case, if its length changes) (2) for some lists, even immutable ones, such as enumerators of FP groups, determining the length maybe arbitrarily expensive, so working out if a list of such things is rectangular would be a problem. I suspect what this really means is that IsRectangularTable should not be a Category.

fingolfin commented 6 years ago

@ThomasBreuer I agree that most of what I wrote above got quite off tangent; gotta make sure to add a link to it "somewhere" where MatrixObj design is discussed, and/or copy it there (but where is such a place? dev/TODO_matobj.txt? A wiki page somewhere?)

@stevelinton the second issue doesn't seem relevant here at all; the second is correct, but IMHO also not quite at the core of the issue -- after all, even an immutable plist cannot store a type, as T_PLIST simply has no provision to do so. Anyway, let's not fribble over finer points -- I think at the core, the only possible conclusion is that using a plist-of-plist is quite limited, and thus one really needs MatrixObj instances for all kinds of reasons, and any attempts to retro fit plists-of-plists with their full capabilities will always be limited.