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/determinants-over-cyclotomic-fields-are-broken/ and on https://groups.google.com/d/msgid/sage-devel/44f0ff3c-71db-4555-a494-80b2ae222c22n%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 over-optimisation 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 223 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: | sage-9.8 |
---|
Sage's matrix multiplication is clearly wrong here: