Ticket #10848: trac_10848-hermitian-matrices-v3.patch
File trac_10848-hermitian-matrices-v3.patch, 7.3 KB (added by , 10 years ago) |
---|
-
sage/matrix/matrix0.pyx
# HG changeset patch # User Rob Beezer <beezer@ups.edu> # Date 1298574394 28800 # Node ID 46f71a629c3d7f0c5ca7465c7091422ae37f7c9a # Parent 8438b7c20d79c02a2ece3e1c3f7224a772ff8f07 10848: add is_hermitian for matrices diff -r 8438b7c20d79 -r 46f71a629c3d 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 conjugate-transpose. 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 double-precision floating-point 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 non-trivial. :: 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 non-real 2946 diagonal entry. :: 2947 2948 sage: A = matrix(QQbar, [[ 2, 2-I, 1+4*I], 2949 ... [ 2+I, 3+I, 2-6*I], 2950 ... [1-4*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 8438b7c20d79 -r 46f71a629c3d sage/matrix/matrix_double_dense.pyx
a b 1327 1327 return False 1328 1328 b = True 1329 1329 for i from 0 < i < self._nrows: 1330 for j from 0 <= j < i: 1330 for j from 0 <= j < i: 1331 1331 if math.fabs(self.get_unsafe(i,j) - self.get_unsafe(j,i)) > tol: 1332 1332 b = False 1333 1333 break 1334 1334 self.cache(key, b) 1335 1335 return b 1336 1336 1337 def is_hermitian(self, tol = 1e-12): 1338 r""" 1339 Return ``True`` if the matrix is equal to its conjugate-transpose. 1340 1341 INPUT: 1342 1343 - ``tol`` - default: ``1e-12`` - 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 round-off 1358 and other imprecisions. 1359 1360 The result is cached, on a per-tolerance 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.0e-17) 1388 True 1389 sage: B.is_hermitian(tol=3.0e-17) 1390 False 1391 1392 A matrix that is nearly Hermitian, but for a non-real 1393 diagonal entry. :: 1394 1395 sage: A = matrix(CDF, [[ 2, 2-I, 1+4*I], 1396 ... [ 2+I, 3+I, 2-6*I], 1397 ... [1-4*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 1337 1430 1338 1431 def cholesky(self): 1339 1432 r"""