diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 852d749e1e6..f5d04784f32 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -12943,7 +12943,179 @@ cdef class Matrix(Matrix1): else: return subspace - def cholesky(self): + def _cholesky_extended_ff(self): + r""" + Performs the extended Cholesky decomposition of a Hermitian matrix over a finite field of square order. + + INPUT: + + - ``self`` -- a square matrix with entries from a finite field of square order + + OUTPUT: + + A square matrix over the same finite such that multiplying itself by its conjugate-transpose yields the input matrix + + ALGORITHM: + + First, we ensure the matrix is square and defined over a finite field of square order. Then we can perform the conjugate-symmetric + version of Gaussian elimination, but the resulting decomposition matrix `L` might not be lower triangular. + + This is a translation of ``BaseChangeToCanonical`` from the GAP ``forms`` package (for a Hermitian form). + + EXAMPLES: + + Here we use the extended decomposition, where the result may not be a lower triangular matrix:: + + sage: U = matrix(GF(17**2),[[0,1],[1,0]]) + sage: B = U._cholesky_extended_ff(); B + [13*z2 + 6 3*z2 + 16] + [13*z2 + 6 14*z2 + 1] + sage: U == B * B.H + True + sage: U = matrix(GF(13**2),[[1,4,7],[4,1,4],[7,4,1]]) + sage: B = U._cholesky_extended_ff(); B + [ 1 0 0] + [ 4 7*z2 + 3 0] + [ 7 6*z2 + 10 12*z2 + 6] + sage: U == B * B.H + True + sage: U = matrix(GF(7**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]]) + sage: B = U._cholesky_extended_ff(); B + [4*z2 + 2 6*z2 0] + [4*z2 + 2 z2 0] + [5*z2 + 6 z2 z2] + sage: U == B * B.H + True + + TESTS: + + If the matrix is not full rank, we compute the rank and throw an exception. + + sage: U = matrix(GF(3**2),[[1,4,7],[4,1,4],[7,4,1]]) + sage: U._cholesky_extended_ff() + Traceback (most recent call last): + ... + ValueError: matrix is not full rank + sage: U = matrix(GF(3**2),[[0,4,7],[4,1,4],[7,4,1]]) + sage: U._cholesky_extended_ff() + Traceback (most recent call last): + ... + ValueError: matrix is not full rank + """ + from sage.matrix.constructor import identity_matrix + + if not self.is_hermitian(): + raise ValueError("matrix is not Hermitian") + + F = self._base_ring + n = self.nrows() + if self.fetch("rank") not in [n, None]: + raise ValueError("matrix is not full rank") + if not (F.is_finite() and F.order().is_square()): + raise ValueError("the base ring must be a finite field of square order") + + q = F.order().sqrt(extend=False) + + def conj_square_root(u): + if u == 0: + return 0 + z = F.multiplicative_generator() + k = u.log(z) + if k % (q + 1) != 0: + raise ValueError(f"unable to factor: {u} is not in base field GF({q})") + return z ** (k//(q+1)) + + if self.nrows() == 1 and self.ncols() == 1: + return self.__class__(F, [conj_square_root(self[0][0])]) + + A = copy(self) + D = identity_matrix(F, n) + row = -1 + + # Diagonalize A + while True: + row += 1 + + # Look for a non-zero element on the main diagonal, starting from `row` + i = row + while i < n and A[i, i].is_zero(): + i += 1 + + if i == row: + # Do nothing since A[row, row] != 0 + pass + elif i < n: + # Swap to ensure A[row, row] != 0 + A.swap_rows(row, i) + A.swap_columns(row, i) + D.swap_rows(row, i) + else: + # All entries on the main diagonal are zero; look for an off-diagonal element + i = row + while i < n - 1: + k = i + 1 + while k < n and A[i, k].is_zero(): + k += 1 + if k == n: + i += 1 + else: + break + + if i == n - 1: + # All elements are zero; terminate + row -= 1 + r = row + 1 + break + + # Fetch the non-zero element and place it at A[row, row + 1] + if i != row: + A.swap_rows(row, i) + A.swap_columns(row, i) + D.swap_rows(row, i) + + A.swap_rows(row + 1, k) + A.swap_columns(row + 1, k) + D.swap_rows(row + 1, k) + + b = ~A[row + 1, row] + A.add_multiple_of_column(row, row + 1, b**q) + A.add_multiple_of_row(row, row + 1, b) + D.add_multiple_of_row(row, row + 1, b) + + # Eliminate below-diagonal entries in the current column + a = ~(-A[row, row]) + for i in range(row + 1, n): + b = A[i, row] * a + if not b.is_zero(): + A.add_multiple_of_column(i, row, b**q) + A.add_multiple_of_row(i, row, b) + D.add_multiple_of_row(i, row, b) + + if row == n - 1: + break + + # Count how many variables have been used + if row == n - 1: + if A[n - 1, n - 1]: # nonzero entry + r = n + else: + r = n - 1 + + if r < n: + self.cache('rank', r) + raise ValueError("matrix is not full rank") + + # Normalize diagonal elements to 1 + for i in range(r): + a = A[i, i] + if not a.is_one(): + # Find an element `b` such that `a = b*b^q = b^(q+1)` + b = conj_square_root(a) + D.rescale_row(i, 1 / b) + + return D.inverse() + + def cholesky(self, extended=False): r""" Return the Cholesky decomposition of a Hermitian matrix. @@ -12989,6 +13161,15 @@ cdef class Matrix(Matrix1): closure or the algebraic reals, depending on whether or not imaginary numbers are required. + Over finite fields, the Cholesky decomposition might + not exist, but when the field has square order (i.e., + `\GF{q^2}`), then we can perform the conjugate-symmetric + version of Gaussian elimination, but the resulting + decomposition matrix `L` might not be lower triangular. + + This is a translation of ``BaseChangeToCanonical`` from + the GAP ``forms`` package (for a Hermitian form). + EXAMPLES: This simple example has a result with entries that remain @@ -13159,6 +13340,30 @@ cdef class Matrix(Matrix1): ... ValueError: matrix is not positive definite + Here we use the extended decomposition, where the result + may not be a lower triangular matrix:: + + sage: U = matrix(GF(5**2),[[0,1],[1,0]]) + sage: B = U.cholesky(extended=True); B + [3*z2 4*z2] + [3*z2 z2] + sage: U == B * B.H + True + sage: U = matrix(GF(11**2),[[1,4,7],[4,1,4],[7,4,1]]) + sage: B = U.cholesky(extended=True); B + [ 1 0 0] + [ 4 9*z2 + 2 0] + [ 7 10*z2 + 1 3*z2 + 3] + sage: U == B * B.H + True + sage: U = matrix(GF(3**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]]) + sage: B = U.cholesky(extended=True); B + [2*z2 2 0] + [2*z2 1 0] + [ 0 1 z2] + sage: U == B * B.H + True + TESTS: This verifies that :issue:`11274` is resolved:: @@ -13206,6 +13411,9 @@ cdef class Matrix(Matrix1): ....: for R in (RR,CC,RDF,CDF,ZZ,QQ,AA,QQbar) ) True """ + if extended: + return self._cholesky_extended_ff() + cdef Matrix C # output matrix C = self.fetch('cholesky') if C is not None: