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: 
Description
As reported on https://ask.sagemath.org/question/64194/determinantsovercyclotomicfieldsarebroken/ and on https://groups.google.com/d/msgid/sagedevel/44f0ff3c71db4555a49480b2ae222c22n%40googlegroups.com
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)] ) LU[0][0]+(L*U)[0][0]
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
comment:2 Changed 4 months ago by
note that 8388593
is a prime, so there is some overoptimisation going on, I guess.
comment:3 Changed 4 months ago by
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. INPUT: self, right  cyclotomic dense matrices with compatible parents (same base ring, and compatible dimensions for matrix multiplication). OUTPUT: cyclotomic dense matrix ALGORITHM: 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
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: MAX_MODULUS 8388608 sage: previous_prime(MAX_MODULUS) 8388593
Could it be one actually need "more primes" in this case, whatever this exactly means?
comment:5 Changed 4 months ago by
If I take half MAX_MODULUS
(which is 2^{23} on my machine)

src/sage/matrix/matrix_cyclo_dense.pyx
a b cdef class Matrix_cyclo_dense(Matrix_dense): 652 652 bound = 1 + 2 * A.height() * B.height() * self._ncols 653 653 654 654 n = self._base_ring._n() 655 p = previous_prime( MAX_MODULUS)655 p = previous_prime(2**22) 656 656 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
Authors:  → Dima Pasechnik 

Branch:  → u/dimpase/matrix/cyclot_matrix_mult_fix 
Commit:  → f4f1bf23a3bc843e62e2137a2b82c312f5180665 
Status:  new → needs_info 
here is a fix, but why it works needs to be understood.
New commits:
f4f1bf2  roughly halve maximum prime to start with, cf #34597

comment:7 Changed 4 months ago by
Commit:  f4f1bf23a3bc843e62e2137a2b82c312f5180665 → 20f69aab69ebb4abc96d5f249a2cbc887a1c5fd2 

Branch pushed to git repo; I updated commit sha1. New commits:
20f69aa  moved test in the right place, reduce MAX_MODULUS globally

comment:8 Changed 10 days ago by
Milestone:  sage9.8 

Sage's matrix multiplication is clearly wrong here: