smith form may be too slow

[Enrique Artal]

I have the toy example  belowto compare the direct use of smith_form and the combination of hermite_form and smith_form.

Last execution (there are random inputs) gave:

CPU times: user 2.22 s, sys: 5.93 ms, total: 2.23 s Wall time: 2.23 s CPU times: user 454 ms, sys: 990 µs, total: 455 ms Wall time: 456 ms

%%%%%%%%%%%%%%%%%%%%%

R.<t> = QQ[]
n = 3
m = 5
ds = (n, m)
M = MatrixSpace(R, n, m)
V = random_vector(R, n)
A0 = Matrix(R, n, m)
A0[0, 0] = R.random_element()
for j in range(1, n):
    A0[j, j] = A0[j -1, j - 1] * R.random_element()
def par(n):
    i1 = ZZ.random_element(0,n)
    j1 = i1
    while j1 == i1:
        j1 = ZZ.random_element(0,n)
    return (i1, j1)    
def smith_hermite(A):
    D1, U1 = A.hermite_form(transformation=True)
    D, U2, V = D1.smith_form()
    U = U2 * U1  
    return D, U, V
P = {j: identity_matrix(R, ds[j]) for j in range(2)}
s = 20
for k in range(2):
    for a in range(s):
        i, j = par(n)
        f = R.random_element()
        P[k].add_multiple_of_row(i, j ,f)
        i, j = par(n)
        P[k].swap_rows(i, j)
        i, j = par(n)
        f = R.random_element()
        P[k].add_multiple_of_column(i, j ,f)
        i, j = par(n)
        P[0].swap_columns(i, j)
A = P[0] * A0 * P[1]
%time D, U, V = A.smith_form()
%time D1, U1, V1 = smith_hermite(A)

[Dima Pasechnik]

hmm, what's Smith form of a matrix over Q ?
Q is not a PID.
>
> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/dec24250-3a00-4c3e-a2e6-e352cde6ff2bn%40googlegroups.com.

[Enrique Artal]

It is over the ring of polynomials in one variable over Q.

[Nils Bruin]

(side note: Q is generally considered a PID, but indeed it makes little sense worrying about Smith normal form over it)

Concerning the slow performance: from the documentation, R._matrix_smith_form would allow an optimized routine to kick in. That doesn't exist, though.

ZZ doesn't have it either, but there is a whole different class for that in sage/matrix/matrix_integer_dense.pyx that provides a specific-for-ZZ smith form algorithm. For QQ['t'] you end up with the generic algorithm, which is liable to be inefficient and have quirks such having large variations in runtime, depending on the specific input.

So: it makes sense it's not very optimal (yet?). It would require work to get a better implementation specific for QQ['t']. This could be worthwhile to do well. Summer of Code?

[Enrique Artal]

Sure, any field is a PID but it is not the targe of a Smith form :)

This question came from an actual computation of a colleague where the direct approach is taking days and hermite+smith is very fast, and could be useful for any one variable ring. Any way,  I agree with you that maybe a Summer of Code is a good idea. Until this, a small PR could be helpful?

[Dima Pasechnik]

On Thu, Dec 14, 2023 at 5:37 AM Nils Bruin <nbr...@sfu.ca> wrote:
>
> (side note: Q is generally considered a PID, but indeed it makes little sense worrying about Smith normal form over it)
>
> Concerning the slow performance: from the documentation, R._matrix_smith_form would allow an optimized routine to kick in. That doesn't exist, though.
> ZZ doesn't have it either, but there is a whole different class for that in sage/matrix/matrix_integer_dense.pyx that provides a specific-for-ZZ smith form algorithm. For QQ['t'] you end up with the generic algorithm, which is liable to be inefficient and have quirks such having large variations in runtime, depending on the specific input.
>
> So: it makes sense it's not very optimal (yet?). It would require work to get a better implementation specific for QQ['t']. This could be worthwhile to do well. Summer of Code?

Can't Sage just wrap PARI/gp matsnf (and mathnf) for QQ[t]? They work
for this case too, e.g.

sage: R.<x> = QQ[]
sage: M=matrix(R, [[1+x^2, 2-x+x^3],[42-x, x^3]])
sage: M.__pari__().matsnf(0)
[x^5 + x^4 - 41*x^3 - x^2 + 44*x - 84, 1]
sage: M.__pari__().mathnf(0)
[x^5 + x^4 - 41*x^3 - x^2 + 44*x - 84, 1/74088*x^4 + 43/74088*x^3 +
1765/74088*x^2 + 41/74088*x + 883/37044; 0, 1]

For some reason
sage: M.__pari__().matsnf(0).sage()

cannot deal with conversion back to Sage, which is probably a bug:

NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 M.__pari__().matsnf(Integer(0)).sage()

File cypari2/gen.pyx:1965, in cypari2.gen.Gen.sage()

File /mnt/opt/Sage/sage-dev/src/sage/libs/pari/convert_sage.pyx:41, in
sage.libs.pari.convert_sage.gen_to_sage()
39
40
---> 41 cpdef gen_to_sage(Gen z, locals=None) noexcept:
42 """
43 Convert a PARI gen to a Sage/Python object.

File /mnt/opt/Sage/sage-dev/src/sage/libs/pari/convert_sage.pyx:300,
in sage.libs.pari.convert_sage.gen_to_sage()
298 return K([gen_to_sage(real), gen_to_sage(imag)])
299 elif t == t_VEC or t == t_COL:
--> 300 return [gen_to_sage(x, locals) for x in z.python_list()]
301 elif t == t_VECSMALL:
302 return z.python_list_small()

File /mnt/opt/Sage/sage-dev/src/sage/libs/pari/convert_sage.pyx:326,
in sage.libs.pari.convert_sage.gen_to_sage()
324 from sage.misc.sage_eval import sage_eval
325 locals = {} if locals is None else locals
--> 326 return sage_eval(str(z), locals=locals)
327
328

File /mnt/opt/Sage/sage-dev/src/sage/misc/sage_eval.py:199, in
sage_eval(source, locals, cmds, preparse)
197 return locals['_sage_eval_returnval_']
198 else:
--> 199 return eval(source, sage.all.__dict__, locals)

File <string>:1

NameError: name 'x' is not defined

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/eae0bd95-4ba8-4e8e-a5c7-c6ec8d83518an%40googlegroups.com.

[Dima Pasechnik]

On Thu, Dec 14, 2023 at 12:52 PM Dima Pasechnik <dim...@gmail.com> wrote:
>
> On Thu, Dec 14, 2023 at 5:37 AM Nils Bruin <nbr...@sfu.ca> wrote:
> >
> > (side note: Q is generally considered a PID, but indeed it makes little sense worrying about Smith normal form over it)
> >
> > Concerning the slow performance: from the documentation, R._matrix_smith_form would allow an optimized routine to kick in. That doesn't exist, though.
> > ZZ doesn't have it either, but there is a whole different class for that in sage/matrix/matrix_integer_dense.pyx that provides a specific-for-ZZ smith form algorithm. For QQ['t'] you end up with the generic algorithm, which is liable to be inefficient and have quirks such having large variations in runtime, depending on the specific input.
> >
> > So: it makes sense it's not very optimal (yet?). It would require work to get a better implementation specific for QQ['t']. This could be worthwhile to do well. Summer of Code?
>
> Can't Sage just wrap PARI/gp matsnf (and mathnf) for QQ[t]? They work
> for this case too, e.g.
>
> sage: R.<x> = QQ[]
> sage: M=matrix(R, [[1+x^2, 2-x+x^3],[42-x, x^3]])
> sage: M.__pari__().matsnf(0)
> [x^5 + x^4 - 41*x^3 - x^2 + 44*x - 84, 1]
> sage: M.__pari__().mathnf(0)
> [x^5 + x^4 - 41*x^3 - x^2 + 44*x - 84, 1/74088*x^4 + 43/74088*x^3 +
> 1765/74088*x^2 + 41/74088*x + 883/37044; 0, 1]
>

There are few caveats here.

matsnf only takes square matrices in this case
(something that can be fixed by padding with zero rows/columns)

a bit of testing shows that on square random matrices, of size up to
7x7, testing along the lines of
the code posted here shows that, if one doesn't do conversion, then
PARI is, on average, twice, or even more as fast as smith_hermite() -
but sometimes it's a slower, in between smith_hermite(A) and
A.smith_form()

Interestingly, Sage's hermite_form(), which seems to be an optimised
Cython implementation,
using _hermite_form_euclidean in matrix2.pyx, beats PARI's mathnf in
terms of speed.
(and mathnf is usually slower than matsnf!)

Conversion back to Sage matrices doesn't work as there is no conevsion
from PARI's polynomials to Sage
polynomials implemented for gen_to_sage(), so it's not a bug per se,
just a gap in functionality.

Dima