Re: [sage-support] linbox bug?
William Stein
>
> ----------------------------------------------------------------------
> | Sage Version 4.0.1, Release Date: 2009-06-06 |
> | Type notebook() for the GUI, and license() for information. |
> ----------------------------------------------------------------------
> sage: A=matrix(GF(3),2,[0,0,1,2])
> sage: R.<x>=GF(3)[]
> sage: D={ x:0 , x+1:0 , x^2+x:0 }
> sage: for i in range(100000):
> ....: D[A._minpoly_linbox()]+=1
> ....:
> sage: D
> {x: 38266, x + 1: 29397, x^2 + x: 32337}
>
You're absolutely right! This *sucks* -- it seems like nothing we
have ever wrapped in Linbox is right at first. Hopefully the issue is
that somehow the algorithm is only supposed to be probabilistic, and
we're just misusing it in sage (quite possible).
Anyway, Clement Pernet will be at Sage Days next week, and we'll sort this out.
Many thanks for brining this to our attention!
This is now:
Jean-Guillaume Dumas
The minpoly algorithm used depends on which method you are using from
LinBox of course but,
If you use the solution "minpoly" you will get the blackbox algorithm
(just like if you specify "minpoly(pol, mat, Method::Blackbox())")
then (since sept 2008 and 1.1.6) we will end up using an extension field
to compute the minpoly (on my machine it will be GF(3^10)) and then
I e.g. got the following result for one try (the algorithm is still
probabilistic, but has a much larger success rate, roughly around 1/3^10):
> 99993 minimal Polynomials are x^2 +x, 3 minimal polynomial are x+1, 4
minimal polynomials are x
Now for a so small matrix it could be better to use a dense version,
which can be called by "minpoly(pol,mat,Method::Elimination())".
If i am correct this dense version is also probabilistic (choice of the
Krylov non-zero vector) and therefore should also pick vectors from an
extension.
This is not the case in 1.1.6.
Clément can you confirm this ? If so it should be easy to fix, the same
way we fixed Wiedemann.
For your example matrix in some of the cases, when vectors [1,1], and
[2,2] are chosen the Krylov space has rank 1, whereas for other non zero
vectors it has rank 2 and
thus the dense minbpoly will be x^2+x or x+1 ...
btw, the returned polynomial is always a factor of the true polynomial,
therefore to get a 1/3^{10k} probability of success it will be
sufficient to perform the lcm of k runs.
Best,
--
Jean-Guillaume Dumas.
____________________________________________________________________
Jean-Guill...@imag.fr Tél.: +33 476 514 866
Université Joseph Fourier, Grenoble I. Fax.: +33 476 631 263
Laboratoire Jean Kuntzmann, Mathématiques Appliquées et Informatique
51, avenue des Mathématiques. LJK/IMAG - BP53. 38041 Grenoble FRANCE
http://ljk.imag.fr/membres/Jean-Guillaume.Dumas
____________________________________________________________________
William Stein
This is from a linbox developer (see below). This will be fixed by:
(1) upgrading -- actually, we *already* use linbox-1.1.6 in sage, so ...
(2) making it so minpoly by default just raises a NotImplementedError, however
minpoly(proof=False) will call minpoly a bunch of times and return
the lcm of the
results.
It turns out that maybe linbox doesn't seem to have a proof=True
minpoly algorithm yet (they are hard to write), so our wrapping of
linbox is wrong, given that in Sage the default is proof=True
everywhere.
Yann -- if you want to work on improving the situation wrt any of the
above, please do.
William
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
William Stein
> William Stein wrote:
>> On Wed, Jun 10, 2009 at 6:03 PM, Yann<yannlaig...@gmail.com> wrote:
>>
>>> ----------------------------------------------------------------------
>>> | Sage Version 4.0.1, Release Date: 2009-06-06 |
>>> | Type notebook() for the GUI, and license() for information. |
>>> ----------------------------------------------------------------------
>>> sage: A=matrix(GF(3),2,[0,0,1,2])
>>> sage: R.<x>=GF(3)[]
>>> sage: D={ x:0 , x+1:0 , x^2+x:0 }
>>> sage: for i in range(100000):
>>> ....: D[A._minpoly_linbox()]+=1
>>> ....:
>>> sage: D
>>> {x: 38266, x + 1: 29397, x^2 + x: 32337}
>>>
>>>
>>
>> You're absolutely right! This *sucks* -- it seems like nothing we
>> have ever wrapped in Linbox is right at first. Hopefully the issue is
>> that somehow the algorithm is only supposed to be probabilistic, and
>> we're just misusing it in sage (quite possible).
>>
>> Anyway, Clement Pernet will be at Sage Days next week, and we'll sort this out.
>> Many thanks for brining this to our attention!
>>
>> This is now:
>>
>> http://trac.sagemath.org/sage_trac/ticket/6296
>>
> Well, I think this was corrected in linbox-1.1.6:
>
> The minpoly algorithm used depends on which method you are using from
> LinBox of course but,
> If you use the solution "minpoly" you will get the blackbox algorithm
> (just like if you specify "minpoly(pol, mat, Method::Blackbox())")
> then (since sept 2008 and 1.1.6) we will end up using an extension field
> to compute the minpoly (on my machine it will be GF(3^10)) and then
> I e.g. got the following result for one try (the algorithm is still
> probabilistic, but has a much larger success rate, roughly around 1/3^10):
void linbox_modn_dense_minpoly(mod_int modulus, mod_int **mp, size_t*
degree, size_t n, mod_int **matrix, int do_minpoly) {
ModInt F((double)modulus);
size_t m = n;
DenseMatrix<ModInt> A(linbox_new_modn_matrix( modulus, matrix, m, m));
GivPolynomial<ModInt::Element> m_A;
if (do_minpoly)
minpoly(m_A, A);
else
charpoly(m_A, A);
(*mp) = new mod_int[m_A.size()];
*degree = m_A.size() - 1;
for (size_t i=0; i <= *degree; i++) {
(*mp)[i] = (mod_int)m_A[i];
}
}
This is from the file interfaces/linbox-sage.C, which ships with linbox.
Many thanks for clarifying that minpoly fails with some probability,
and that we need to call it multiple times, take lcm's, and force the
user to give the option "proof=False" to use it.
Just out of curiosity, is there any provably correct minpoly in
linbox? We don't have one in Sage at all, so it would be useful so
we don't have to implement one from scratch.
William
Jean-Guillaume Dumas
> We're using linbox-1.1.6 in Sage.
>
> void linbox_modn_dense_minpoly(mod_int modulus, mod_int **mp, size_t*
> degree, size_t n, mod_int **matrix, int do_minpoly) {
>
> ModInt F((double)modulus);
>
> size_t m = n;
>
> DenseMatrix<ModInt> A(linbox_new_modn_matrix( modulus, matrix, m, m));
>
> GivPolynomial<ModInt::Element> m_A;
>
> if (do_minpoly)
> minpoly(m_A, A);
> else
> charpoly(m_A, A);
>
> (*mp) = new mod_int[m_A.size()];
> *degree = m_A.size() - 1;
> for (size_t i=0; i <= *degree; i++) {
> (*mp)[i] = (mod_int)m_A[i];
> }
>
> }
>
> This is from the file interfaces/linbox-sage.C, which ships with linbox.
>
> Many thanks for clarifying that minpoly fails with some probability,
> and that we need to call it multiple times, take lcm's, and force the
> user to give the option "proof=False" to use it.
>
elimination code is called and it does not have the extension trick yet,
that's why the failure probability is high for very small characteristic.
> Just out of curiosity, is there any provably correct minpoly in
> linbox?
though since it loops over minpolys until its degree is n)
I can sketch a simple one but with a pretty bad complexity (multiplied
by n) like this :
run several (k of them) probabilistic minpoly over a word size
extension, get the lcm Pi_1(X).
if the extension is of size p^t, correctness will be guaranteed except
for 1/p^{tk}.
If you want a proof of correctness, then just check that Pi_1(A) = 0. If
this is not the case lcm new minpolys to Pi_1 and check again.
William Stein
Ouch. (In sage the default for *everything* is proof=True.) Is there any
easier way to prove correctness of the minpoly (this is a basic exact
linear algebra question)? At the least I could factor Pi_1 first
before substituting in A.
William
Jean-Guillaume Dumas
> Ouch. (In sage the default for *everything* is proof=True.) Is there any
> easier way to prove correctness of the minpoly (this is a basic exact
> linear algebra question)?
answer this in a lower complexity.
We do not have those algorithms implemented yet.
In small fields they would in general require also working in an
extension, which will also kill the use of the BLAS.
Therefore in practice, I am afraid that for small finite fields blas
based lcm of non extended minpoly's then blas based evaluation of Pi_1
at A might be the best for a large range of dense matrix sizes.
> At the least I could factor Pi_1 first
> before substituting in A.
>
yes good idea (and storing the lcm'ed components too might prove also
useful).
Clement Pernet
I missed this discussion.
So I confirm, that we are using the Krylov based dense minpoly that
returns a probabilistic minpoly, and I'm the one to blame for not taking
care of this.
So far, I don't know of any better certificate for minpoly than checking
Pmin(A) = 0 (improved by using the factors, as William pointed out).
I'm also thinking of Danilevskii's elimination: although it could be
viewed as a Krylov method, it does not use any randomization, so could
it always give the whole minpoly ?
As for the use of extensions, the dense method imposes restrictions on
the size of the field (in order to use the BLAS). So this might be a
severe limiting factor.
I'm willing to sort this pb out at Sage Days 16 next week.
Clément
Jean-Guillaume Dumas a écrit :
Clement Pernet
> I'm also thinking of Danilevskii's elimination: although it could be
> viewed as a Krylov method, it does not use any randomization, so could
> it always give the whole minpoly ?
I checked that this does not solve our pb: Danilevskii's alg is
deterministic for charpoly but still can choose bad vectors for minpoly.
A conter-example is
[1 1]
[0 1]
for which Danilevskii's alg implicitly uses the first column vector
[1]
[0]
which only produces the factor X-1 of the minpoly
Clément