Large matrices, Cython and Sage

[kroo...@gmail.com]

Hi,

Is it possible to work with large matrices (say 26000 x 15000) in sage?

I have a C function that returns a large matrix.

I would like to convert this matrix into an object that sage recognizes, so that I could run matrix-related functions on it.

The dimensions of the matrix are 26620 over 15180.

Although this is a relatively large matrix, the C code that generates it takes only a few seconds.

However, I can't seem to be able to transform it into a sage matrix object.

I have the following code inside a cython function

cdef char* resultMatrix = get_c_matrix();

answer = [[resultMatrix[(i*15180)+j] for j in range(15180)] for i in range(26620)];

The second line that tries to convert the matrix into a python object kills sage.

Even a simple statement as

A = [[0 for j in range(15180)] for i in range(26620)]

kills Sage.

Is there another way to work with my C-generated matrix in sage?

Thanks,

Avi

[Harald Schilly]

Well, if one entry is 64bit, you are maybe at 4gb memory? Converting
that to Sage, might add some overhead for each element. What's inside
the matrix, is it sparse?

-- h

> --
> You received this message because you are subscribed to the Google Groups
> "sage-cloud" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-cloud+...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sage-cloud/5f2d1d65-b9e9-4e11-91b0-971b914b2693%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Dima Pasechnik]

Is it a numeric matrix, or some kind of "special" one, e.g. 0-1 ? 

IMHO working with matrices of this size "directly" is pretty much hopeless.

Perhaps I am stating obvious to you:

usually one would resort to working with the matrix as with a "black box", that can e.g.

compute matrix-vector products.

Just in case,

Dima

Thanks,Avi

[Dan Christensen]

On Feb 23, 2016, kroo...@gmail.com wrote:

> Is it possible to work with large matrices (say 26000 x 15000) in sage?

An approach that treats each entry as a separate python object is going
to have too much overhead. I recommend using numpy. It can take the
existing C matrix and use it directly, without copying any data.

Dan

[kroo...@gmail.com]

The matrix is not sparse.

Its elements are from a small prime finite field (each entry is a byte in it's C representation).

בתאריך יום רביעי, 24 בפברואר 2016 בשעה 00:30:08 UTC+2, מאת Harald Schilly:

[kroo...@gmail.com]

Is there an easy way to compute the rank of a numpy matrix over a Galois field?

[kroo...@gmail.com]

The matrix entries are byte sized, each represents an element from a small prime field.

I want to check if the matrix's null-space contains the null-space of another (smaller matrix).

I planned to use sage's ability to compute rank over finite fields (compute the rank A and check whether the rank of the block matrix [A over B] is larger).

בתאריך יום רביעי, 24 בפברואר 2016 בשעה 01:12:12 UTC+2, מאת Dima Pasechnik:

[Dan Christensen]

On Feb 24, 2016, kroo...@gmail.com wrote:

> Is there an easy way to compute the rank of a numpy matrix over a Galois
> field?

Not that I know of. I mistakenly assumed you were working with floating
point data.

Dan

PS: Please don't copy me on replies, as indicated by the headers:

Mail-Followup-To: sage-...@googlegroups.com
Mail-Copies-To: never

[Simon King]

Hi!

On 2016-02-24, kroo...@gmail.com <kroo...@gmail.com> wrote:
> The matrix entries are byte sized, each represents an element from a small
> prime field.
> I want to check if the matrix's null-space contains the null-space of
> another (smaller matrix).
> I planned to use sage's ability to compute rank over finite fields (compute
> the rank A and check whether the rank of the block matrix [A over B] is
> larger).

Recently, support for MeatAxe matrices was added to Sage. If your
version of Sage is new, then you could install an optional MeatAxe
package by
sage -i meataxe
on the commandline. Then, recompile the parts of sage depending on the
new module, by
sage -b
If you then start Sage, MeatAxe would be used for matrices over finite
non-prime fields of order <255. Internally, a compressed representation
is used (i.e. if the field is very small then several elements are
stored in one byte).

I don't know, however, if conversion of your matrices into the MeatAxe
matrices will turn out to be a bottle neck.

Best regards,
Simon

[kcrisman]

Recently, support for MeatAxe matrices was added to Sage. If your
version of Sage is new, then you could install an optional MeatAxe
package by
   sage -i meataxe
on the commandline.

Since this message is on sage-cloud, presumably the version is the newest William wants to install ;-)

William, is meataxe already installed standard in the SMC version of Sage?

[William Stein]

Simon -- what's the status of meataxe? It's an experimental package
for some reason...

--

Best Regards,
William Stein

CEO, SageMath, Inc.

Professor of Mathematics at a university

[Dima Pasechnik]

William,

while you're at it, you asked to promote csdp to optional:

http://trac.sagemath.org/ticket/19977

[Simon King]

Hi William,

On 2016-02-24, William Stein <wst...@sagemath.com> wrote:
> On Wed, Feb 24, 2016 at 9:44 AM, kcrisman <kcri...@gmail.com> wrote:
>
> Simon -- what's the status of meataxe? It's an experimental package
> for some reason...

Isn't it the case that new packages should first be experimental,
unless there is a compelling reason to immediately make them optional?

Best regards,
Simon

[Simon King]

PS:

On 2016-02-26, Simon King <simon...@uni-koeln.de> wrote:
> Isn't it the case that new packages should first be experimental,
> unless there is a compelling reason to immediately make them optional?

Some more background information:

Upstream is the "Aachen C-MeatAxe". I could be mistaken, but I *think* I
was told that it is used in GAP (internally --- if a user creates a
matrix, MeatAxe is *not* used).

The MeatAxe project is not very actively developed, but it isn't dead.
Upstream told me that they prepare a new version (after almost 10
years).

Moreover, I have used MeatAxe inside of my group cohomology spkg, that
used to be "optional" but was broken by the release manager.

That's what speaks in favour of making the package "optional".

However, I did patch the upstream sources rather non-trivially:
- One patch prevents MeatAxe from writing multiplication tables into the
current directory.
- One patch implements error propagation (upstream would just crash on
any kind of error).
- One patch improves the performance of the computation of echelon
forms, in a rather obvious way.
- One patch implements asymptotically fast matrix multiplication, which
beats MeatAxe's schoolbook implementation even for relatively small
matrices. Thus, I use it as default algorithm in my Cython wrapper.
- One patch makes it so that MeatAxe uses asymptotically fast
multiplication internally, in some places.

I did send my patches to upstream, but there was no reaction at all (not
even acknowleding receipt). That's what speaks in favour of making it
"experimental".

Concerning test coverage:
- MeatAxe has a test suite, that can be executed during installation of
the package (of course only upon request).
- In my Cython wrapper, I have doc tests with random matrices, and compare
the results of computing echelon form and matrix products obtained by
MeatAxe with or without asymptotically fast algorithms respectively
obtained by Sage's current implementation of matrices.

MeatAxe is A *LOT* faster in matrix arithmetics and echelon computation
than the current matrix implementation in Sage, in the case of finite
non-prime fields of odd characteristic and order < 255. For these base
fields, my MeatAxe wrapper can be seen as a drop-in replacement of the
current default implementation in Sage.

So, that's another reason for making it "experimental".

I guess I should post this to sage-devel?

Best regards,
Simon

[Dima Pasechnik]

Hi Simon,

On Friday, February 26, 2016 at 8:27:52 AM UTC, Simon King wrote:

PS:

On 2016-02-26, Simon King <simon...@uni-koeln.de> wrote:
> Isn't it the case that new packages should first be experimental,
> unless there is a compelling reason to immediately make them optional?

Some more background information:

Upstream is the "Aachen C-MeatAxe". I could be mistaken, but I *think* I
was told that it is used in GAP (internally --- if a user creates a
matrix, MeatAxe is *not* used).

IMHO presently there is no meataxe C code in GAP, so it's another implementation

used, by D.Holt and S.Rees. 

 

It would be interesting to compare two implementations (meataxe and GAP) for speed.

 

The MeatAxe project is not very actively developed, but it isn't dead.
Upstream told me that they prepare a new version (after almost 10
years).

Moreover, I have used MeatAxe inside of my group cohomology spkg, that
used to be "optional" but was broken by the release manager.

That's what speaks in favour of making the package "optional".

However, I did patch the upstream sources rather non-trivially:
- One patch prevents MeatAxe from writing multiplication tables into the
  current directory.
- One patch implements error propagation (upstream would just crash on
  any kind of error).
- One patch improves the performance of the computation of echelon
  forms, in a rather obvious way.
- One patch implements asymptotically fast matrix multiplication, which
  beats MeatAxe's schoolbook implementation even for relatively small
  matrices. Thus, I use it as default algorithm in my Cython wrapper.
- One patch makes it so that MeatAxe uses asymptotically fast
  multiplication internally, in some places.

I did send my patches to upstream, but there was no reaction at all (not
even acknowleding receipt). That's what speaks in favour of making it
"experimental".

did they (Lux and Ringe) even receive them? 

One possible way to make it upstream is to create a GAP package with this meataxe...

Dima

[Simon King]

Hi Dima,

On 2016-02-26, Dima Pasechnik <dim...@gmail.com> wrote:
> IMHO presently there is no meataxe C code in GAP, so it's another
> implementation
> used, by D.Holt and S.Rees.

Really? A possible reason of confusion might be that the non-matrix part
of MeatAxe is used in GAP, or am I mistaken in that point as well? In
any case, when David Green gave me his group cohomology code 10ish years
ago, he told me that he used MeatAxe as a matrix backend, since it was
used in GAP.

> It would be interesting to compare two implementations (meataxe and GAP)
> for speed.

No problem, as soon as you tell me how to use create a matrix in GAP
that uses a fast backend (no, M:=[[...],[...]] is not enough).

>> I did send my patches to upstream, but there was no reaction at all (not
>> even acknowleding receipt). That's what speaks in favour of making it
>> "experimental".
>>
>
> did they (Lux and Ringe) even receive them?

I suppose so. There was some mail exchange back and forth while I
created the patches.

> One possible way to make it upstream is to create a GAP package with this
> meataxe...

I only learned recently how to create a new style Sage package. I don't
like to create a GAP package now...
More seriously, creating a new-style version of the group cohomology
package (which requires a major re-write of David Green's C-sources) is
of higher priority to me.

Best regards,
Simon

[kcrisman]

>
> Simon -- what's the status of meataxe?  It's an experimental package
> for some reason...

Isn't it the case that new packages should first be experimental,
unless there is a compelling reason to immediately make them optional?

?  No, experimental just means it doesn't build everywhere, optional means it builds and works on all "supported" platforms, for some value of the word "supported".  (Or it's supposed to mean that, it may not in practice.)

[Dima Pasechnik]

On Friday, February 26, 2016 at 1:09:22 PM UTC, Simon King wrote:

Hi Dima,

On 2016-02-26, Dima Pasechnik <dim...@gmail.com> wrote:
> IMHO presently there is no meataxe C code in GAP, so it's another
> implementation
> used, by D.Holt and S.Rees.

Really? A possible reason of confusion might be that the non-matrix part
of MeatAxe is used in GAP, or am I mistaken in that point as well?

GAP manual says:

 69.12 The Smash MeatAxe

  

  The standard MeatAxe provided in the GAP library is based on the MeatAxe in the GAP 3 package Smash, originally written by Derek Holt and Sarah Rees [HR94]. It is accessible via the variable SMTX to which MTX

  (69.3-1) is assigned by default.  [...]

 

In
any case, when David Green gave me his group cohomology code 10ish years
ago, he told me that he used MeatAxe as a matrix backend, since it was
used in GAP.

there was C-Meataxe package in GAP3

http://gap-system.org/Gap3/Packages3/meataxe.html

> It would be interesting to compare two implementations (meataxe and GAP)
> for speed.

No problem, as soon as you tell me how to use create a matrix in GAP
that uses a fast backend (no, M:=[[...],[...]] is not enough).

ImmutableMatrix (sect 24.14-1 ImmutableMatrix)

Dima

[William Stein]

Hi,

I tried "sage -i meataxe" on sage-6.10 (the one installed SMC). That
seemed to work fine. I then did "sage -b", which seemed to work. I
then made a random matrix over GF(9), and it is of type
Matrix_generic_dense, so I don't think meataxe is used. I then
reverted what I did. That this modifies code in the core library
makes me extra nervous...

William

> --
> You received this message because you are subscribed to the Google Groups
> "sage-cloud" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sage-cloud+...@googlegroups.com.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sage-cloud/95ae8eaa-e4ce-4242-aa32-638d0403319e%40googlegroups.com.

>
> For more options, visit https://groups.google.com/d/optout.

[Simon King]

Hi William,

On 2016-02-26, William Stein <wst...@sagemath.com> wrote:
> I tried "sage -i meataxe" on sage-6.10 (the one installed SMC). That
> seemed to work fine. I then did "sage -b", which seemed to work. I
> then made a random matrix over GF(9), and it is of type
> Matrix_generic_dense, so I don't think meataxe is used.

I don't know what version the MeatAxe wrapper was merged into. See
whether SAGE_ROOT/src/sage/matrix/matrix_gfpn_dense.pyx and ...pxd is
present or not. If it isn't, then the wrapper is missing and "sage -i
meataxe" has no visible effect.

> I then
> reverted what I did. That this modifies code in the core library
> makes me extra nervous...

The *wrapper* is in Sage, regardless whether you did "sage -i meataxe"
or not. It is an optional extension module. Look into
SAGE_ROOT/src/module_list.py, where it says (if you have a recent Sage
version):
OptionalExtension("sage.matrix.matrix_gfpn_dense",
sources = ['sage/matrix/matrix_gfpn_dense.pyx'],
libraries = ['mtx'],
package = 'meataxe'),

And if you further look into sage.matrix.matrix_space, you'll find that
it is *tried* to import matrix_gfpn_dense. On success, it is used for
certain base rings. On error, the generic (very slow) implementation is
used.

Hope that clearifies how it works!
Cheers,
Simon

[Simon King]

Hi Karl-Dieter,

On 2016-02-26, kcrisman <kcri...@gmail.com> wrote:
> ? No, experimental just means it doesn't build everywhere, optional means
> it builds and works on all "supported" platforms, for some value of the
> word "supported". (Or it's supposed to mean that, it may not in practice.)

Really? Then I am confident that it can be promoted to "optional", and
will ask for it on sage-devel. What I use in my old group cohomology
spkg is based on a different MeatAxe version, but it did work on all
supported platforms (big and little endian and so on).

Best regards,
Simon

[William Stein]

I defined "experimental" packages at the top of
http://files.sagemath.org/spkg/experimental/ to be "These are
EXPERIMENTAL! They probably won't work at all for you! Use at your own
risk! Many of these have *never* been successfully built on any
platform! (But still, if you can figure out how to build them, I'd
like to know about it.) These also may not be available under a
GPL-compatible license."

>
> Best regards,
> Simon
>
> --
> You received this message because you are subscribed to the Google Groups "sage-cloud" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-cloud+...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-cloud/naqngm%245mn%242%40ger.gmane.org.

> For more options, visit https://groups.google.com/d/optout.

[Simon King]

Hi Dima,

On 2016-02-26, Dima Pasechnik <dim...@gmail.com> wrote:
>> No problem, as soon as you tell me how to use create a matrix in GAP
>> that uses a fast backend (no, M:=[[...],[...]] is not enough).
>>
>
> ImmutableMatrix (sect 24.14-1 ImmutableMatrix)

Sorry, it still is a problem. Let me rephrase my question:
Let M be a native Sage matrix. How to quickly convert M into a matrix in
the libgap (or gap) interface?

sage: M = random_matrix(GF(9,'x'),4000,4000)
sage: Mg = libgap.ImmutableMatrix(M.base_ring(),M)
Traceback (most recent call last):
...
AttributeError: No such attribute: ImmutableMatrix.
sage: Mg = gap.ImmutableMatrix(M.base_ring(),M)
<hitting Ctrl-c after several minutes>

So, libgap doesn't seem to know ImmutableMatrix and gap takes ages.

Best regards,
Simon

[Dima Pasechnik]

it is perfectly within your means to let it know it.

Edit 

src/sage/libs/gap/gap_functions.py

 --- a/src/sage/libs/gap/gap_functions.py

+++ b/src/sage/libs/gap/gap_functions.py

@@ -242,6 +242,7 @@ common_gap_functions = [

   'Identity',

   'Image',

   'Images',

+  'ImmutableMatrix',

   'Index',

   'InfoAlgebra',

   'InfoAttributes',

and then 'sage -br'

(and add this to your next ticket to make it permanent...)

There is also libgap's function_factory facility to do more or less the same.

Having said this, libgap does not seem to be fast enough here, either.

sage: M1 = random_matrix(GF(9,'x'),400,400)

sage: %time libgap.ImmutableMatrix(M1.base_ring(),M1)

CPU times: user 19.5 s, sys: 20.2 ms, total: 19.5 s

Wall time: 19.5 s

< immutable compressed matrix 400x400 over GF(9) >

I guess it's due to very slow conversion between Sage and GAP finite fields,

especially non-prime ones.

(if I do the same for ZZ, I get 10-fold speedup; if I do with GF(13), I get 6-fold speedup)

Not sure how easy is to fix this.

Dima

Best regards,
Simon

[Simon King]

Hi Dima,

On 2016-02-27, Dima Pasechnik <dim...@gmail.com> wrote:
> Edit
> src/sage/libs/gap/gap_functions.py
>
> --- a/src/sage/libs/gap/gap_functions.py
> +++ b/src/sage/libs/gap/gap_functions.py
> @@ -242,6 +242,7 @@ common_gap_functions = [
> 'Identity',
> 'Image',
> 'Images',
> + 'ImmutableMatrix',
> 'Index',
> 'InfoAlgebra',
> 'InfoAttributes',
>
> and then 'sage -br'

Thanks! I get this:

sage: M = random_matrix(GF(9,'x'),4000,4000)

sage: K = M.base_ring()
sage: D = dict((a,libgap(a)) for a in K)
sage: %time gL = libgap([libgap([D[a] for a in row]) for row in M])
CPU times: user 57 s, sys: 40 ms, total: 57 s
Wall time: 56.9 s
sage: gM = libgap.ImmutableMatrix(K, gL)
sage: gM
< immutable compressed matrix 4000x4000 over GF(9) >

> I guess it's due to very slow conversion between Sage and GAP finite fields,
> especially non-prime ones.

That's why above I use a dictionary, to avoid repeated conversion.

It is a bit strange that small finite fields cache their elements
but do not cache their conversions.

Anyway.

Timings for some matrix computations:
sage: %time gM2 = gM*gM
CPU times: user 2min 5s, sys: 0 ns, total: 2min 5s
Wall time: 2min 5s
sage: %time M2 = M*M
CPU times: user 25.1 s, sys: 4 ms, total: 25.1 s
Wall time: 25.1 s

and
sage: %time gM2 = gM^(-1)
CPU times: user 2min 46s, sys: 4 ms, total: 2min 46s
Wall time: 2min 46s
sage: %time M2 = M^(-1)
CPU times: user 1min, sys: 12 ms, total: 1min
Wall time: 1min

The good timing for multiplication is due to my patches: I added
asymptotically fast Winograd-Strassen matrix multiplication with
a memory efficient schedule to meataxe.

I also added some trivial improvements to the implementation of Gaussian
elimination in meataxe, but I did not bother to implement an asymptotically
fast algorithm. Still it's faster than what is in gap. And way faster than
what was in Sage before.

Alas, it is still way slower than Magma, if I recall correctly.

Best regards,
Simon