Ticket #10848: trac_10848hermitianmatricesv4.patch
File trac_10848hermitianmatricesv4.patch, 7.0 KB (added by , 10 years ago) 


sage/matrix/matrix0.pyx
# HG changeset patch # User Rob Beezer <beezer@ups.edu> # Date 1298574394 28800 # Node ID 07a2e975404a08cddf524ef0ad6ac59217519672 # Parent aff88d266c1d2bca6b7be97d579e198ebf466863 10848: add is_hermitian for matrices diff r aff88d266c1d r 07a2e975404a sage/matrix/matrix0.pyx
a b 2896 2896 return False 2897 2897 return True 2898 2898 2899 def is_hermitian(self): 2900 r""" 2901 Returns ``True`` if the matrix is equal to its conjugatetranspose. 2902 2903 OUTPUT: 2904 2905 ``True`` if the matrix is square and equal to the transpose 2906 with every entry conjugated, and ``False`` otherwise. 2907 2908 Note that if conjugation has no effect on elements of the base 2909 ring (such as for integers), then the :meth:`is_symmetric` 2910 method is equivalent and faster. 2911 2912 This routine is for matrices over exact rings and so may not 2913 work properly for matrices over ``RR`` or ``CC``. For matrices with 2914 approximate entries, the rings of doubleprecision floatingpoint 2915 numbers, ``RDF`` and ``CDF``, are a better choice since the 2916 :meth:`sage.matrix.matrix_double_dense.Matrix_double_dense.is_hermitian` 2917 method has a tolerance parameter. This provides control over 2918 allowing for minor discrepancies between entries when checking 2919 equality. 2920 2921 EXAMPLES:: 2922 2923 sage: A = matrix(QQbar, [[ 1 + I, 1  6*I, 1  I], 2924 ... [3  I, 4*I, 2], 2925 ... [1 + I, 2  8*I, 2 + I]]) 2926 sage: A.is_hermitian() 2927 False 2928 sage: B = A*A.conjugate_transpose() 2929 sage: B.is_hermitian() 2930 True 2931 2932 Sage has several fields besides the entire complex numbers 2933 where conjugation is nontrivial. :: 2934 2935 sage: F.<b> = QuadraticField(7) 2936 sage: C = matrix(F, [[2*b  3, 7*b  6, b + 3], 2937 ... [2*b  3, 3*b + 2, 2*b], 2938 ... [ b + 1, 0, 2]]) 2939 sage: C.is_hermitian() 2940 False 2941 sage: C = C*C.conjugate_transpose() 2942 sage: C.is_hermitian() 2943 True 2944 2945 A matrix that is nearly Hermitian, but for a nonreal 2946 diagonal entry. :: 2947 2948 sage: A = matrix(QQbar, [[ 2, 2I, 1+4*I], 2949 ... [ 2+I, 3+I, 26*I], 2950 ... [14*I, 2+6*I, 5]]) 2951 sage: A.is_hermitian() 2952 False 2953 sage: A[1,1] = 132 2954 sage: A.is_hermitian() 2955 True 2956 2957 Rectangular matrices are never Hermitian. :: 2958 2959 sage: A = matrix(QQbar, 3, 4) 2960 sage: A.is_hermitian() 2961 False 2962 """ 2963 if not self.is_square(): 2964 return False 2965 2966 cdef Py_ssize_t i,j 2967 2968 for i from 0 <= i < self._nrows: 2969 for j from 0 <= j <= i: 2970 if self.get_unsafe(i,j) != self.get_unsafe(j,i).conjugate(): 2971 return False 2972 return True 2973 2899 2974 def is_skew_symmetric(self): 2900 2975 """ 2901 2976 Returns True if this is a skew symmetric matrix. 
sage/matrix/matrix_double_dense.pyx
diff r aff88d266c1d r 07a2e975404a sage/matrix/matrix_double_dense.pyx
a b 1334 1334 self.cache(key, b) 1335 1335 return b 1336 1336 1337 def is_hermitian(self, tol = 1e12): 1338 r""" 1339 Return ``True`` if the matrix is equal to its conjugatetranspose. 1340 1341 INPUT: 1342 1343  ``tol``  default: ``1e12``  the largest value of the 1344 absolute value of the difference between two matrix entries 1345 for which they will still be considered equal. 1346 1347 OUTPUT: 1348 1349 ``True`` if the matrix is square and equal to the transpose 1350 with every entry conjugated, and ``False`` otherwise. 1351 1352 Note that if conjugation has no effect on elements of the base 1353 ring (such as for integers), then the :meth:`is_symmetric` 1354 method is equivalent and faster. 1355 1356 The tolerance parameter is used to allow for numerical values 1357 to be equal if there is a slight difference due to roundoff 1358 and other imprecisions. 1359 1360 The result is cached, on a pertolerance basis. 1361 1362 EXAMPLES:: 1363 1364 sage: A = matrix(CDF, [[ 1 + I, 1  6*I, 1  I], 1365 ... [3  I, 4*I, 2], 1366 ... [1 + I, 2  8*I, 2 + I]]) 1367 sage: A.is_hermitian() 1368 False 1369 sage: B = A*A.conjugate_transpose() 1370 sage: B.is_hermitian() 1371 True 1372 1373 We get a unitary matrix from the SVD routine and use this 1374 numerical matrix to create a matrix that should be Hermitian 1375 (indeed it should be the identity matrix), but with some 1376 imprecision. We use this to illustrate that if the tolerance 1377 is set too small, then we can be too strict about the equality 1378 of entries and achieve the wrong result. :: 1379 1380 sage: A = matrix(CDF, [[ 1 + I, 1  6*I, 1  I], 1381 ... [3  I, 4*I, 2], 1382 ... [1 + I, 2  8*I, 2 + I]]) 1383 sage: U, _, _ = A.SVD() 1384 sage: B=U*U.conjugate_transpose() 1385 sage: B.is_hermitian() 1386 True 1387 sage: B.is_hermitian(tol=5.0e17) 1388 True 1389 sage: B.is_hermitian(tol=3.0e17) 1390 False 1391 1392 A matrix that is nearly Hermitian, but for a nonreal 1393 diagonal entry. :: 1394 1395 sage: A = matrix(CDF, [[ 2, 2I, 1+4*I], 1396 ... [ 2+I, 3+I, 26*I], 1397 ... [14*I, 2+6*I, 5]]) 1398 sage: A.is_hermitian() 1399 False 1400 sage: A[1,1] = 132 1401 sage: A.is_hermitian() 1402 True 1403 1404 Rectangular matrices are never Hermitian. :: 1405 1406 sage: A = matrix(CDF, 3, 4) 1407 sage: A.is_hermitian() 1408 False 1409 """ 1410 global numpy 1411 tol = float(tol) 1412 key = 'hermitian_%s'%tol 1413 b = self.fetch(key) 1414 if not b is None: 1415 return b 1416 if not self.is_square(): 1417 self.cache(key, False) 1418 return False 1419 if numpy is None: 1420 import numpy 1421 cdef Py_ssize_t i, j 1422 hermitian = True 1423 for i from 0 < i < self._nrows: 1424 for j from 0 <= j <= i: 1425 if numpy.absolute(self.get_unsafe(i,j)  self.get_unsafe(j,i).conjugate()) > tol: 1426 hermitian = False 1427 break 1428 self.cache(key, hermitian) 1429 return hermitian 1430 1337 1431 def schur(self, base_ring=None): 1338 1432 r""" 1339 1433 Returns the Schur decomposition of the matrix.