#2190: [with bundle, positive review] implement a ZZ-module saturation
algorithm:
this is the key thing needed to compute kernels over ZZ, etc.
----------------------------------+---------------------------
Reporter: was | Owner: was
Type: enhancement | Status: closed
Priority: major | Milestone: sage-2.10.2
Component: linear algebra | Resolution: fixed
Keywords: | Merged in:
Authors: | Reviewers:
Report Upstream: N/A | Work issues:
Branch: | Commit:
Dependencies: | Stopgaps:
----------------------------------+---------------------------
Changes (by chapoton):
* upstream: => N/A
Old description:
> Saturation of a matrix, i.e., this command in Magma:
> {{{
> (<Mtrx> X) -> Mtrx
>
> The basis of the saturation w.r.t. Q of the space spanned by the rows
> of X.
> }}}
> is really important to a huge amount of my modular forms research.
> In Magma this command exists and is quite fast, because I specifically
> needed it for my work and Allan spent a lot of time on making it quite
> good. In Sage, it's built on top of PARI's {{{matkerint}}} command (one
> can -reduce computing saturation to computing the
> kernel over ZZ of the transpose matrix). As a result, Sage is almost
> completely and totally useless at doing any saturations when the size
> is at all interesting for research. For example, a very tiny example
> below is nearly a hundred times slower in Sage than in Magma!
>
> {{{
> sage: a = random_matrix(ZZ, 40, 42)
> sage: w = a.row_space()
> sage: time s =w.saturation()
> CPU times: user 0.79 s, sys: 0.07 s, total: 0.86 s
> Wall time: 0.89
> sage: b = magma(a)
> sage: magma.eval('time c := Saturation(%s);'%b.name())
> 'Time: 0.010'
> sage: 0.89 / 0.01
> 89.0000000000000
> }}}
>
> Here are some more examples:
>
> {{{
> sage: a = random_matrix(ZZ, 150,150)
> sage: w = a.row_space()
> sage: time s = w.saturation()
> CPU times: user 1.36 s, sys: 0.21 s, total: 1.57 s
> Wall time: 1.64
> sage: m = magma(a)
> sage: magma.eval('time c := Saturation(%s);'%m.name())
> 'Time: 0.290'
> }}}
>
> and a very typical one for my research would look like this:
> {{{
> sage: a = random_matrix(ZZ, 100,300)
> sage: w = a.row_space()
> sage: time s = w.saturation()
> CPU times: user 17653.78 s, sys: 360.28 s, total: 18014.06 s
> Wall time: 18241.30
> sage: 18014/3600.0 # FIVE HOURS!!
> 5.00388888888889
> sage: a = random_matrix(ZZ, 100,300)
> sage: m = magma(a)
> sage: magma.eval('time c := Saturation(%s);'%m.name())
> 'Time: 0.810'
> }}}
>
> I just want to pause to emphasize that:
> 1. there is no open source code in the world as far as I know that
> can do saturations quickly; only Magma does them right now, and
> 2. doing these reasonably quickly is utterly essential to my modular
> computational forms research.
> 3. The current situation is that what should be an easy everyday
> calculation for me -- e.g., saturating a 100x300 -- takes FIVE HOURS in
> Sage, but less than one second in Magma!
>
> Fortunately, there is a fast algorithm for saturation that depends mainly
> on kernels over GF(p) and HNF, and we now have very fast HNF because of
> #174.
>
> Here's some code in Magma that does p-satuation. This combined with
> other tricks should lead to a saturation algorithm. Porting this to Sage
> is definitely step 1.
>
> {{{
> function pAdicSaturation(B,p)
> if Type(B[1]) eq SeqEnum then
> V := RSpace(Rationals(),#B[1]);
> B := [ V | v : v in B];
> end if;
> V := Universe(B);
> n := Degree(V);
> for i in [1..#B] do
> B[i] *:= LCM([ Denominator(c) : c in Eltseq(B[i]) ]);
> end for;
> ZZ := Integers();
> FF := FiniteField(p);
> B := RMatrixSpace(ZZ,#B,n)!Matrix(B);
> m := Rank(B);
> B := Submatrix(HermiteForm(B),1,1,m,n);
> N := RMatrixSpace(FF,m,n)!B;
> while Rank(N) lt m do
> K := Kernel(N);
> vprintf pAdicSaturation :
> "Rank(N) + Rank(K) = %o + %o = %o\n", Rank(N), Rank(K), m;
> C := RMatrixSpace(ZZ,#Basis(K),n)!
> Matrix([ (1/p)*V!&+[ ZZ!u[i]*B[i] : i in [1..m] ] : u in Basis(K)
> ]);
> vtime pAdicSaturation, 2 :
> B := Submatrix(HermiteForm(VerticalJoin(B,C)),1,1,m,n);
> N := RMatrixSpace(FF,m,n)!B;
> end while;
> vprintf pAdicSaturation : "Rank(N) = %o \n", Rank(N), m;
> return [ B[i] : i in [1..m] ];
> end function;
> }}}
New description:
Saturation of a matrix, i.e., this command in Magma:
{{{
(<Mtrx> X) -> Mtrx
The basis of the saturation w.r.t. Q of the space spanned by the rows
of X.
}}}
is really important to a huge amount of my modular forms research.
In Magma this command exists and is quite fast, because I specifically
needed it for my work and Allan spent a lot of time on making it quite
good. In Sage, it's built on top of PARI's {{{matkerint}}} command (one
can -reduce computing saturation to computing the
kernel over ZZ of the transpose matrix). As a result, Sage is almost
completely and totally useless at doing any saturations when the size
is at all interesting for research. For example, a very tiny example
below is nearly a hundred times slower in Sage than in Magma!
{{{
sage: a = random_matrix(ZZ, 40, 42)
sage: w = a.row_space()
sage: time s =w.saturation()
CPU times: user 0.79 s, sys: 0.07 s, total: 0.86 s
Wall time: 0.89
sage: b = magma(a)
sage: magma.eval('time c := Saturation(%s);'%b.name())
'Time: 0.010'
sage: 0.89 / 0.01
89.0000000000000
}}}
Here are some more examples:
{{{
sage: a = random_matrix(ZZ, 150,150)
sage: w = a.row_space()
sage: time s = w.saturation()
CPU times: user 1.36 s, sys: 0.21 s, total: 1.57 s
Wall time: 1.64
sage: m = magma(a)
sage: magma.eval('time c := Saturation(%s);'%m.name())
'Time: 0.290'
}}}
and a very typical one for my research would look like this:
{{{
sage: a = random_matrix(ZZ, 100,300)
sage: w = a.row_space()
sage: time s = w.saturation()
CPU times: user 17653.78 s, sys: 360.28 s, total: 18014.06 s
Wall time: 18241.30
sage: 18014/3600.0 # FIVE HOURS!!
5.00388888888889
sage: a = random_matrix(ZZ, 100,300)
sage: m = magma(a)
sage: magma.eval('time c := Saturation(%s);'%m.name())
'Time: 0.810'
}}}
I just want to pause to emphasize that:
1. there is no open source code in the world as far as I know that
can do saturations quickly; only Magma does them right now, and
2. doing these reasonably quickly is utterly essential to my modular
computational forms research.
3. The current situation is that what should be an easy everyday
calculation for me -- e.g., saturating a 100x300 -- takes FIVE HOURS in
Sage, but less than one second in Magma!
Fortunately, there is a fast algorithm for saturation that depends mainly
on kernels over GF(p) and HNF, and we now have very fast HNF because of
#174.
Here's some code in Magma that does p-satuation. This combined with other
tricks should lead to a saturation algorithm. Porting this to Sage is
definitely step 1.
{{{
function pAdicSaturation(B,p)
if Type(B[1]) eq SeqEnum then
V := RSpace(Rationals(),#B[1]);
B := [ V | v : v in B];
end if;
V := Universe(B);
n := Degree(V);
for i in [1..#B] do
B[i] *:= LCM([ Denominator(c) : c in Eltseq(B[i]) ]);
end for;
ZZ := Integers();
FF := FiniteField(p);
B := RMatrixSpace(ZZ,#B,n)!Matrix(B);
m := Rank(B);
B := Submatrix(HermiteForm(B),1,1,m,n);
N := RMatrixSpace(FF,m,n)!B;
while Rank(N) lt m do
K := Kernel(N);
vprintf pAdicSaturation :
"Rank(N) + Rank(K) = %o + %o = %o\n", Rank(N), Rank(K), m;
C := RMatrixSpace(ZZ,#Basis(K),n)!
Matrix([ (1/p)*V!&+[ ZZ!u[i]*B[i] : i in [1..m] ] : u in Basis(K)
]);
vtime pAdicSaturation, 2 :
B := Submatrix(HermiteForm(VerticalJoin(B,C)),1,1,m,n);
N := RMatrixSpace(FF,m,n)!B;
end while;
vprintf pAdicSaturation : "Rank(N) = %o \n", Rank(N), m;
return [ B[i] : i in [1..m] ];
end function;
}}}
--
--
Ticket URL: <http://trac.sagemath.org/ticket/2190#comment:9>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/d/optout.
