Opened 4 months ago

Last modified 10 days ago

#34597 needs_info defect

bug in cyclotomic matrix multiplication

Reported by: dimpase Owned by:
Priority: critical Milestone:
Component: linear algebra Keywords:
Cc: was Merged in:
Authors: Dima Pasechnik Reviewers:
Report Upstream: N/A Work issues:
Branch: u/dimpase/matrix/cyclot_matrix_mult_fix (Commits, GitHub, GitLab) Commit: 20f69aab69ebb4abc96d5f249a2cbc887a1c5fd2
Dependencies: Stopgaps:

GitHub link to the corresponding issue


As reported on and on

K.<z> = CyclotomicField(16)
OK = K.ring_of_integers()
L = [[-575*z^7 - 439*z^6 - 237*z^5 + 237*z^3 + 439*z^2 + 575*z + 623, 0],
    [0,     -114*z^7 - 88*z^6 - 48*z^5 + 48*z^3 + 88*z^2 + 114*z + 123]]
U = [[-1926*z^7 - 1474*z^6 - 798*z^5 + 798*z^3 + 1474*z^2 + 1926*z + 2085, 0],
    [0,   -1014*z^7 - 777*z^6 - 421*z^5 + 421*z^3 + 777*z^2 + 1014*z + 1097]]
L, U = matrix(K,L), matrix(K,U)
LU = matrix( [ [L[i].inner_product(U.transpose()[j]) for j in range(2)] for i in range(2)] )

does not report 0, but 8388593*z^7 - 8388593*z - 8388593

I checked using another cyclotomics implementation, in GAP, that basic arithmetic over K, (i.e. computation of LU) is correct.

Change History (8)

comment:1 Changed 4 months ago by dimpase

Sage's matrix multiplication is clearly wrong here:

sage: libgap(L)*libgap(U)-libgap(LU)
[ [ 0, 0 ], [ 0, 0 ] ]
sage: libgap(L)*libgap(U)-libgap(L*U)
[ [ 8388593+8388593*E(16)-8388593*E(16)^7, 0 ], [ 0, 0 ] ]

comment:2 Changed 4 months ago by dimpase

note that 8388593 is a prime, so there is some over-optimisation going on, I guess.

comment:3 Changed 4 months ago by dimpase

Cc: was added

The implementation is in src/sage/matrix/matrix_cyclo_dense.pyx

cdef _matrix_times_matrix_(self, baseMatrix right):
        Return the product of two cyclotomic dense matrices.

            self, right -- cyclotomic dense matrices with compatible
                           parents (same base ring, and compatible
                           dimensions for matrix multiplication).

            cyclotomic dense matrix

            Use a multimodular algorithm that involves multiplying the
            two matrices modulo split primes.

so this explains that appearance of a largish prime in the error is not a fluke.

comment:4 Changed 4 months ago by dimpase

the prime that pops up here is quite special:

sage: from sage.matrix.matrix_modn_dense_double import MAX_MODULUS as MAX_MODULUS_modn_dense_double
....: from sage.arith.multi_modular import MAX_MODULUS as MAX_MODULUS_multi_modular
....: MAX_MODULUS = min(MAX_MODULUS_modn_dense_double, MAX_MODULUS_multi_modular)
sage: previous_prime(MAX_MODULUS)

Could it be one actually need "more primes" in this case, whatever this exactly means?

comment:5 Changed 4 months ago by dimpase

If I take half MAX_MODULUS (which is 223 on my machine)

  • src/sage/matrix/matrix_cyclo_dense.pyx

    a b cdef class Matrix_cyclo_dense(Matrix_dense): 
    652652        bound = 1 + 2 * A.height() * B.height() * self._ncols
    654654        n = self._base_ring._n()
    655         p = previous_prime(MAX_MODULUS)
     655        p = previous_prime(2**22)
    656656        prod = 1

then this particular error goes away, while all the tests in src/sage/matrix/matrix_cyclo_dense.pyx still pass.

comment:6 Changed 4 months ago by dimpase

Authors: Dima Pasechnik
Branch: u/dimpase/matrix/cyclot_matrix_mult_fix
Commit: f4f1bf23a3bc843e62e2137a2b82c312f5180665
Status: newneeds_info

here is a fix, but why it works needs to be understood.

New commits:

f4f1bf2roughly halve maximum prime to start with, cf #34597

comment:7 Changed 4 months ago by git

Commit: f4f1bf23a3bc843e62e2137a2b82c312f518066520f69aab69ebb4abc96d5f249a2cbc887a1c5fd2

Branch pushed to git repo; I updated commit sha1. New commits:

20f69aamoved test in the right place, reduce MAX_MODULUS globally

comment:8 Changed 10 days ago by mkoeppe

Milestone: sage-9.8
Note: See TracTickets for help on using tickets.