Skip to content

Commit cefde12

Browse files
jacksonwalterstscrim
jacksonwalters
and
tscrim
committed
Update src/sage/matrix/matrix2.pyx
add hermitian_decomposition method given a Hermitian matrix U, returns a matrix A such that U = AA* use translated GAP source code for BaseChangeCanonical raise ValueError for non full rank matrices this method does not work for singular matrices. since we are already computing the rank in the row reduction, we can just check at the end whether the rank is full, and if not exit with a ValueError. this may be overcautious - it's not clear that this won't work for some singular matrices. there are cases where it certainly doesn't work, and we include one in the doctest. add example for singular case this example is currently failing and needs to be caught remove unnecessary import, singularity check if the matrix `U` is singular, the process will still result in a matrix but it will fail to have the expected property, namely `B*B.H == U`. Update src/sage/matrix/matrix2.pyx avoid import, directly call `sqrt` method, don't extend to keep result in `ZZ` fix example this was using the output from the `forms` package, when it should be from the GAP source translation. change examples there is a failing example due to the fact that the matrix is singular with q=3, and U = matrix(F,[[1,4,7],[4,1,4],[7,4,1]]). move _cholesky_extended_ff to hidden method move _cholesky_extended_ff to a separate hidden method, and just call it from inside cholesky if extended=true. add doctests and docs for the method separately. sparse for now, should be expanded upon. Update src/sage/matrix/matrix2.pyx add output for GF(3**2) case this case fails to have a lower triangular decomposition which has been checked by brute force. Update src/sage/matrix/matrix2.pyx move two examples over finite fields from tests the two examples computing the extended Cholesky decomposition over square order finite fields are in TESTS, and should be moved to EXAMPLES. Update src/sage/matrix/matrix2.pyx this is better and clarifies that the result of the Cholesky decomposition is lower triangular (and not possibly upper triangular). Update src/sage/matrix/matrix2.pyx add single whitespace before `extended` parameter Update src/sage/matrix/matrix2.pyx looks good. makes sense to point out the Cholesky decomposition might not exist, given that was our initial difficulty. 3 trailing whitespace move hermitian_decomposition into cholesky this method for decomposition a Hermitian matrix U = AA* is similar in spirit to the Cholesky decomposition, but extends it to work over finite fields of square order. add two examples add two examples of a Hermitian decomposition, one of which is not upper/lower triangular (so would be impossible with Cholesky decomposition) Update src/sage/matrix/matrix2.pyx yes, just forgot to update this initialize row to -1 Update src/sage/matrix/matrix2.pyx add word "package" to GAP ``forms`` Update src/sage/matrix/matrix2.pyx change \ast to * for consistency Trigger CI tests use self.__class__ within matrix.pyx file Trigger CI tests remove blank line Update src/sage/matrix/matrix2.pyx Update src/sage/matrix/matrix2.pyx Update src/sage/matrix/matrix2.pyx Update src/sage/matrix/matrix2.pyx Update src/sage/matrix/matrix2.pyx use Matrix handle 1x1 case Co-Authored-By: Travis Scrimshaw <[email protected]>
1 parent 9cd86e9 commit cefde12

File tree

1 file changed

+202
-1
lines changed

1 file changed

+202
-1
lines changed

‎src/sage/matrix/matrix2.pyx

+202-1
Original file line numberDiff line numberDiff line change
@@ -12943,7 +12943,172 @@ cdef class Matrix(Matrix1):
1294312943
else:
1294412944
return subspace
1294512945

12946-
def cholesky(self):
12946+
def _cholesky_extended_ff(self):
12947+
r"""
12948+
Performs the extended Cholesky decomposition of a Hermitian matrix over a finite field of square order.
12949+
12950+
INPUT:
12951+
12952+
- ``self`` -- a square matrix with entries from a finite field of square order
12953+
12954+
OUTPUT:
12955+
12956+
A square matrix over the same finite such that multiplying itself by its conjugate-transpose yields the input matrix
12957+
12958+
ALGORITHM:
12959+
12960+
First, we ensure the matrix is square and defined over a finite field of square order. Then we can perform the conjugate-symmetric
12961+
version of Gaussian elimination, but the resulting decomposition matrix `L` might not be lower triangular.
12962+
12963+
This is a translation of ``BaseChangeToCanonical`` from the GAP ``forms`` package (for a Hermitian form).
12964+
12965+
EXAMPLES:
12966+
12967+
Here we use the extended decomposition, where the result may not be a lower triangular matrix::
12968+
12969+
sage: U = matrix(GF(17**2),[[0,1],[1,0]])
12970+
sage: B = U._cholesky_extended_ff(); B
12971+
[13*z2 + 6 3*z2 + 16]
12972+
[13*z2 + 6 14*z2 + 1]
12973+
sage: U == B * B.H
12974+
True
12975+
sage: U = matrix(GF(13**2),[[1,4,7],[4,1,4],[7,4,1]])
12976+
sage: B = U._cholesky_extended_ff(); B
12977+
[ 1 0 0]
12978+
[ 4 7*z2 + 3 0]
12979+
[ 7 6*z2 + 10 12*z2 + 6]
12980+
sage: U == B * B.H
12981+
True
12982+
sage: U = matrix(GF(7**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
12983+
sage: B = U._cholesky_extended_ff(); B
12984+
[4*z2 + 2 6*z2 0]
12985+
[4*z2 + 2 z2 0]
12986+
[5*z2 + 6 z2 z2]
12987+
sage: U == B * B.H
12988+
True
12989+
sage: U = matrix(GF(3**2),[[1,4,7],[4,1,4],[7,4,1]])
12990+
sage: B = U._cholesky_extended_ff(); B
12991+
Traceback (most recent call last)
12992+
...
12993+
ValueError: matrix is not full rank
12994+
sage: U = matrix(GF(3**2),[[0,4,7],[4,1,4],[7,4,1]])
12995+
sage: B = U._cholesky_extended_ff(); B
12996+
Traceback (most recent call last)
12997+
...
12998+
ValueError: matrix is not full rank
12999+
"""
13000+
from sage.matrix.constructor import identity_matrix
13001+
13002+
if not self.is_hermitian():
13003+
raise ValueError("matrix is not Hermitian")
13004+
13005+
F = self._base_ring
13006+
n = self.nrows()
13007+
13008+
if not (F.is_finite() and F.order().is_square()):
13009+
raise ValueError("the base ring must be a finite field of square order")
13010+
13011+
q = F.order().sqrt(extend=False)
13012+
13013+
def conj_square_root(u):
13014+
if u == 0:
13015+
return 0
13016+
z = F.multiplicative_generator()
13017+
k = u.log(z)
13018+
if k % (q + 1) != 0:
13019+
raise ValueError(f"unable to factor: {u} is not in base field GF({q})")
13020+
return z ** (k//(q+1))
13021+
13022+
if self.nrows() == 1 and self.ncols() == 1:
13023+
return self.__class__(F, [conj_square_root(self[0][0])])
13024+
13025+
A = copy(self)
13026+
D = identity_matrix(F, n)
13027+
row = -1
13028+
13029+
# Diagonalize A
13030+
while True:
13031+
row += 1
13032+
13033+
# Look for a non-zero element on the main diagonal, starting from `row`
13034+
i = row
13035+
while i < n and A[i, i].is_zero():
13036+
i += 1
13037+
13038+
if i == row:
13039+
# Do nothing since A[row, row] != 0
13040+
pass
13041+
elif i < n:
13042+
# Swap to ensure A[row, row] != 0
13043+
A.swap_rows(row, i)
13044+
A.swap_columns(row, i)
13045+
D.swap_rows(row, i)
13046+
else:
13047+
# All entries on the main diagonal are zero; look for an off-diagonal element
13048+
i = row
13049+
while i < n - 1:
13050+
k = i + 1
13051+
while k < n and A[i, k].is_zero():
13052+
k += 1
13053+
if k == n:
13054+
i += 1
13055+
else:
13056+
break
13057+
13058+
if i == n - 1:
13059+
# All elements are zero; terminate
13060+
row -= 1
13061+
r = row + 1
13062+
break
13063+
13064+
# Fetch the non-zero element and place it at A[row, row + 1]
13065+
if i != row:
13066+
A.swap_rows(row, i)
13067+
A.swap_columns(row, i)
13068+
D.swap_rows(row, i)
13069+
13070+
A.swap_rows(row + 1, k)
13071+
A.swap_columns(row + 1, k)
13072+
D.swap_rows(row + 1, k)
13073+
13074+
b = ~A[row + 1, row]
13075+
A.add_multiple_of_column(row, row + 1, b**q)
13076+
A.add_multiple_of_row(row, row + 1, b)
13077+
D.add_multiple_of_row(row, row + 1, b)
13078+
13079+
# Eliminate below-diagonal entries in the current column
13080+
a = ~(-A[row, row])
13081+
for i in range(row + 1, n):
13082+
b = A[i, row] * a
13083+
if not b.is_zero():
13084+
A.add_multiple_of_column(i, row, b**q)
13085+
A.add_multiple_of_row(i, row, b)
13086+
D.add_multiple_of_row(i, row, b)
13087+
13088+
if row == n - 1:
13089+
break
13090+
13091+
# Count how many variables have been used
13092+
if row == n - 1:
13093+
if A[n - 1, n - 1]: # nonzero entry
13094+
r = n
13095+
else:
13096+
r = n - 1
13097+
13098+
if r < n:
13099+
raise ValueError("matrix is not full rank")
13100+
13101+
# Normalize diagonal elements to 1
13102+
for i in range(r):
13103+
a = A[i, i]
13104+
if not a.is_one():
13105+
# Find an element `b` such that `a = b*b^q = b^(q+1)`
13106+
b = conj_square_root(a)
13107+
D.rescale_row(i, 1 / b)
13108+
13109+
return D.inverse()
13110+
13111+
def cholesky(self, extended=False):
1294713112
r"""
1294813113
Return the Cholesky decomposition of a Hermitian matrix.
1294913114

@@ -12989,6 +13154,15 @@ cdef class Matrix(Matrix1):
1298913154
closure or the algebraic reals, depending on whether or not
1299013155
imaginary numbers are required.
1299113156

13157+
Over finite fields, the Cholesky decomposition might
13158+
not exist, but when the field has square order (i.e.,
13159+
`\GF{q^2}`), then we can perform the conjugate-symmetric
13160+
version of Gaussian elimination, but the resulting
13161+
decomposition matrix `L` might not be lower triangular.
13162+
13163+
This is a translation of ``BaseChangeToCanonical`` from
13164+
the GAP ``forms`` package (for a Hermitian form).
13165+
1299213166
EXAMPLES:
1299313167

1299413168
This simple example has a result with entries that remain
@@ -13159,6 +13333,30 @@ cdef class Matrix(Matrix1):
1315913333
...
1316013334
ValueError: matrix is not positive definite
1316113335

13336+
Here we use the extended decomposition, where the result
13337+
may not be a lower triangular matrix::
13338+
13339+
sage: U = matrix(GF(5**2),[[0,1],[1,0]])
13340+
sage: B = U.cholesky(extended=True); B
13341+
[3*z2 4*z2]
13342+
[3*z2 z2]
13343+
sage: U == B * B.H
13344+
True
13345+
sage: U = matrix(GF(11**2),[[1,4,7],[4,1,4],[7,4,1]])
13346+
sage: B = U.cholesky(extended=True); B
13347+
[ 1 0 0]
13348+
[ 4 9*z2 + 2 0]
13349+
[ 7 10*z2 + 1 3*z2 + 3]
13350+
sage: U == B * B.H
13351+
True
13352+
sage: U = matrix(GF(3**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
13353+
sage: B = U.cholesky(extended=True); B
13354+
[2*z2 2 0]
13355+
[2*z2 1 0]
13356+
[ 0 1 z2]
13357+
sage: U == B * B.H
13358+
True
13359+
1316213360
TESTS:
1316313361

1316413362
This verifies that :issue:`11274` is resolved::
@@ -13206,6 +13404,9 @@ cdef class Matrix(Matrix1):
1320613404
....: for R in (RR,CC,RDF,CDF,ZZ,QQ,AA,QQbar) )
1320713405
True
1320813406
"""
13407+
if extended:
13408+
return self._cholesky_extended_ff()
13409+
1320913410
cdef Matrix C # output matrix
1321013411
C = self.fetch('cholesky')
1321113412
if C is not None:

0 commit comments

Comments
 (0)