Ticket #10944: trac_10944-similar-matrices.patch

File trac_10944-similar-matrices.patch, 12.2 KB (added by rbeezer, 9 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1300250063 25200
    # Node ID da6aef91866b5aee2e187c72d1b1a754461d24aa
    # Parent  6a679959b54b7975af1841df4dd77d28b0faef28
    10944: similarity check for matrices
    
    diff -r 6a679959b54b -r da6aef91866b sage/matrix/matrix2.pyx
    a b  
    62996299            return J, transformation_matrix
    63006300        else:
    63016301            return J
    6302    
     6302
     6303    def is_similar(self, other, transformation=False):
     6304        r"""
     6305        Returns true if ``self`` and ``other`` are similar,
     6306        i.e. related by a change-of-basis matrix.
     6307
     6308        INPUT:
     6309
     6310        - ``other`` - a matrix, which should be square, and of the same size
     6311          as ``self``, where the entries of the matrix have a fraction field
     6312          equal to that of ``self``.  Inexact rings are not supported.
     6313
     6314        - ``transformation`` - default: ``False`` - if ``True``, the output
     6315          will include the change-of-basis matrix.  See below for an exact
     6316          description.
     6317
     6318        OUTPUT:
     6319
     6320        Two matrices, $A$ and $B$ are similar if there is an invertible
     6321        matrix $S$ such that $A=S^{-1}BS$.  $S$ can be interpreted as a
     6322        change-of-basis matrix if $A$ and $B$ are viewed as matrix
     6323        representations of the same linear transformation.
     6324
     6325        When ``transformation=False`` this method will return ``True`` if
     6326        such a matrix $S$ exists, otherwise it will return ``False``.  When
     6327        ``transformation=True`` the method returns a pair.  The first part
     6328        of the pair is ``True`` or ``False`` depending on if the matrices
     6329        are similar and the second part is the change-of-basis matrix, or
     6330        ``None`` should it not exist.
     6331
     6332        When the transformation matrix is requested, it will satisfy
     6333        ``self = S.inverse()*other*S``.
     6334
     6335        If the base rings for any of the matrices is the integers, the
     6336        rationals, or the field of algebraic numbers (``QQbar``), then the
     6337        matrices are converted to have ``QQbar`` as their base ring prior
     6338        to checking the equality of the base rings.
     6339
     6340        It is possible for this routine to fail over most fields, even when
     6341        the matrices are similar.  However, since the field of algebraic
     6342        numbers is algebraically closed, the routine will always produce
     6343        a result for matrices with rational entries.
     6344
     6345        EXAMPLES:
     6346
     6347        The two matrices in this example were constructed to be similar.
     6348        The computations happen in the field of algebraic numbers, but we
     6349        are able to convert the change-of-basis matrix back to the rationals
     6350        (which may not always be possible). ::
     6351
     6352            sage: A = matrix(ZZ, [[-5, 2, -11],
     6353            ...                   [-6, 7, -42],
     6354            ...                   [0, 1, -6]])
     6355            sage: B = matrix(ZZ, [[ 1, 12,  3],
     6356            ...                   [-1, -6, -1],
     6357            ...                   [ 0,  6,  1]])
     6358            sage: A.is_similar(B)
     6359            True
     6360            sage: _, T = A.is_similar(B, transformation=True)
     6361            sage: T
     6362            [  1.000000000000000? + 0.?e-32*I              0.?e-33 + 0.?e-33*I              0.?e-32 + 0.?e-32*I]
     6363            [-0.6666666666666667? + 0.?e-33*I  0.1666666666666667? + 0.?e-33*I -0.8333333333333333? + 0.?e-32*I]
     6364            [ 0.6666666666666667? + 0.?e-32*I              0.?e-33 + 0.?e-33*I -0.3333333333333334? + 0.?e-32*I]
     6365            sage: T.change_ring(QQ)
     6366            [   1    0    0]
     6367            [-2/3  1/6 -5/6]
     6368            [ 2/3    0 -1/3]
     6369            sage: A == T.inverse()*B*T
     6370            True
     6371
     6372        Other exact fields are supported.  ::
     6373
     6374            sage: F.<a> = FiniteField(7^2)
     6375            sage: A = matrix(F,[[2*a + 5, 6*a + 6,   a + 3],
     6376            ...                 [  a + 3, 2*a + 2, 4*a + 2],
     6377            ...                 [2*a + 6, 5*a + 5,     3*a]])
     6378            sage: B = matrix(F,[[5*a + 5, 6*a + 4,   a + 1],
     6379            ...                 [  a + 5, 4*a + 3, 3*a + 3],
     6380            ...                 [3*a + 5,   a + 4, 5*a + 6]])
     6381            sage: A.is_similar(B)
     6382            True
     6383            sage: B.is_similar(A)
     6384            True
     6385            sage: _, T = A.is_similar(B, transformation=True)
     6386            sage: T
     6387            [      1       0       0]
     6388            [6*a + 1 4*a + 3 4*a + 2]
     6389            [6*a + 3 3*a + 5 3*a + 6]
     6390            sage: A == T.inverse()*B*T
     6391            True
     6392
     6393        Two matrices with different sets of eigenvalues, so they
     6394        cannot possibly be similar. ::
     6395
     6396            sage: A = matrix(QQ, [[ 2,  3, -3, -6],
     6397            ...                   [ 0,  1, -2, -8],
     6398            ...                   [-3, -3,  4,  3],
     6399            ...                   [-1, -2,  2,  6]])
     6400            sage: B = matrix(QQ, [[ 1,  1,  2,  4],
     6401            ...                   [-1,  2, -3, -7],
     6402            ...                   [-2,  3, -4, -7],
     6403            ...                   [ 0, -1,  0,  0]])
     6404            sage: A.eigenvalues() == B.eigenvalues()
     6405            False
     6406            sage: A.is_similar(B, transformation=True)
     6407            (False, None)
     6408
     6409        Similarity is an equivalence relation, so this routine computes
     6410        a representative of the equivalence class for each matrix, the
     6411        Jordan form, as provided by :meth:`jordan_form`.  The matrices
     6412        below have identical eigenvalues (as evidenced by equal
     6413        characteristic polynomials), but slightly different Jordan forms,
     6414        and hence are not similar.
     6415
     6416            sage: A = matrix(QQ, [[ 19, -7, -29],
     6417            ...                   [-16, 11,  30],
     6418            ...                   [ 15, -7, -25]])
     6419            sage: B = matrix(QQ, [[-38, -63,  42],
     6420            ...                   [ 14,  25, -14],
     6421            ...                   [-14, -21,  18]])
     6422            sage: A.charpoly() == B.charpoly()
     6423            True
     6424            sage: A.jordan_form()
     6425            [-3| 0  0]
     6426            [--+-----]
     6427            [ 0| 4  1]
     6428            [ 0| 0  4]
     6429            sage: B.jordan_form()
     6430            [-3| 0| 0]
     6431            [--+--+--]
     6432            [ 0| 4| 0]
     6433            [--+--+--]
     6434            [ 0| 0| 4]
     6435            sage: A.is_similar(B)
     6436            False
     6437
     6438        Obtaining the Jordan form requires computing the eigenvalues of
     6439        the matrix, which may not lie in the field used for entries of
     6440        the matrix.  So the routine first checks the characteristic
     6441        polynomials - if they are unequal, then the matrices cannot be
     6442        similar. However, when the characteristic polynomials are equal,
     6443        we must examine the Jordan form. In this case, the method may fail,
     6444        EVEN when the matrices are similar.  This is not the case for
     6445        matrices over the integers, rationals or algebraic numbers,
     6446        since the computations are done in the algebraically closed
     6447        field of algebraic numbers.
     6448
     6449        Here is an example where the similarity is obvious, but the
     6450        routine fails to compute a result.  ::
     6451
     6452            sage: F.<a> = FiniteField(7^2)
     6453            sage: C = matrix(F,[[  a + 2, 5*a + 4],
     6454            ...                 [6*a + 6, 6*a + 4]])
     6455            sage: S = matrix(ZZ, [[0, 1],
     6456            ...                   [1, 0]])
     6457            sage: D = S.inverse()*C*S
     6458            sage: C.is_similar(D)
     6459            Traceback (most recent call last):
     6460            ...
     6461            ValueError: unable to compute Jordan canonical form for a matrix
     6462            sage: C.jordan_form()
     6463            Traceback (most recent call last):
     6464            ...
     6465            RuntimeError: Some eigenvalue does not exist in Finite Field in a of size 7^2.
     6466
     6467        Inexact rings and fields are also not supported.  ::
     6468
     6469            sage: A = matrix(CDF, 2, 2, range(4))
     6470            sage: B = copy(A)
     6471            sage: A.is_similar(B)
     6472            Traceback (most recent call last):
     6473            ...
     6474            ValueError: unable to compute Jordan canonical form for a matrix
     6475
     6476        Rectangular matrices and mismatched sizes return quickly.  ::
     6477
     6478            sage: A = matrix(3, 2, range(6))
     6479            sage: B = copy(A)
     6480            sage: A.is_similar(B)
     6481            False
     6482            sage: A = matrix(2, 2, range(4))
     6483            sage: B = matrix(3, 3, range(9))
     6484            sage: A.is_similar(B, transformation=True)
     6485            (False, None)
     6486
     6487        If the fraction fields of the entries are unequal, it is an error,
     6488        except in the case when the rationals gets promoted to the
     6489        algebraic numbers.  ::
     6490
     6491            sage: A = matrix(ZZ, 2, 2, range(4))
     6492            sage: B = matrix(GF(2), 2, 2, range(4))
     6493            sage: A.is_similar(B, transformation=True)
     6494            Traceback (most recent call last):
     6495            ...
     6496            TypeError: matrices need to have entries with identical fraction fields, not Algebraic Field and Finite Field of size 2
     6497            sage: A = matrix(ZZ, 2, 2, range(4))
     6498            sage: B = matrix(QQbar, 2, 2, range(4))
     6499            sage: A.is_similar(B)
     6500            True
     6501
     6502        Inputs are checked.  ::
     6503
     6504            sage: A = matrix(ZZ, 2, 2, range(4))
     6505            sage: A.is_similar('garbage')
     6506            Traceback (most recent call last):
     6507            ...
     6508            TypeError: similarity requires a matrix as an argument, not garbage
     6509            sage: B = copy(A)
     6510            sage: A.is_similar(B, transformation='junk')
     6511            Traceback (most recent call last):
     6512            ...
     6513            ValueError: transformation keyword must be True or False, not junk
     6514        """
     6515        import sage.matrix.matrix
     6516        import sage.rings.qqbar
     6517        if not sage.matrix.matrix.is_Matrix(other):
     6518            raise TypeError('similarity requires a matrix as an argument, not {0}'.format(other))
     6519        if transformation not in [True, False]:
     6520            raise ValueError('transformation keyword must be True or False, not {0}'.format(transformation))
     6521        # easy false situations
     6522        if not self.is_square() or not other.is_square:
     6523            if transformation:
     6524                return (False, None)
     6525            else:
     6526                return False
     6527        if self.nrows() != self.nrows():
     6528            if transformation:
     6529                return (False, None)
     6530            else:
     6531                return False
     6532        # convert to fraction fields for base rings
     6533        A = self.matrix_over_field()
     6534        B = other.matrix_over_field()
     6535        # move rationals to algebraically closed algebraic numbers
     6536        if A.base_ring() == QQ:
     6537            A = A.change_ring(sage.rings.qqbar.QQbar)
     6538        if B.base_ring() == QQ:
     6539            B = B.change_ring(sage.rings.qqbar.QQbar)
     6540        # require identical base fields
     6541        if A.base_ring() != B.base_ring():
     6542            raise TypeError('matrices need to have entries with identical fraction fields, not {0} and {1}'.format(A.base_ring(), B.base_ring()))
     6543        # unequal characteristic polynomials implies not similar
     6544        # and avoids any problems with eigenvalues not in the base field
     6545        if A.charpoly() != B.charpoly():
     6546            if transformation:
     6547                return (False, None)
     6548            else:
     6549                return False
     6550        # now more precisely compare Jordan form, and optionally get transformations
     6551        try:
     6552            if transformation:
     6553                JA, SA = A.jordan_form(transformation=True)
     6554            else:
     6555                JA = A.jordan_form(transformation=False)
     6556        except:
     6557            raise ValueError('unable to compute Jordan canonical form for a matrix')
     6558        try:
     6559            if transformation:
     6560                JB, SB = B.jordan_form(transformation=True)
     6561            else:
     6562                JB = B.jordan_form(transformation=False)
     6563        except:
     6564            raise ValueError('unable to compute Jordan canonical form for a matrix')
     6565        similar = (JA == JB)
     6566        transform = None
     6567        if similar and transformation:
     6568            transform = SB*SA.inverse()
     6569        if transformation:
     6570            return (similar, transform)
     6571        else:
     6572            return similar
     6573
    63036574    def symplectic_form(self):
    63046575        r"""
    63056576        Find a symplectic form for self if self is an anti-symmetric,