Ticket #10944: trac_10944-similar-matrices-v2.patch

File trac_10944-similar-matrices-v2.patch, 12.1 KB (added by ddrake, 9 years ago)

rebased for 4.7.1.alpha1; apply only this

  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1300250063 25200
    # Node ID 471f04c5d728d0af6bc86018f346fdc296cef3c5
    # Parent  e89e0526716feac93609eecfd07348f8559cacb0
    10944: similarity check for matrices
    
    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    71957195                return False
    71967196        return True
    71977197
     7198    def is_similar(self, other, transformation=False):
     7199        r"""
     7200        Returns true if ``self`` and ``other`` are similar,
     7201        i.e. related by a change-of-basis matrix.
     7202
     7203        INPUT:
     7204
     7205        - ``other`` - a matrix, which should be square, and of the same size
     7206          as ``self``, where the entries of the matrix have a fraction field
     7207          equal to that of ``self``.  Inexact rings are not supported.
     7208
     7209        - ``transformation`` - default: ``False`` - if ``True``, the output
     7210          will include the change-of-basis matrix.  See below for an exact
     7211          description.
     7212
     7213        OUTPUT:
     7214
     7215        Two matrices, $A$ and $B$ are similar if there is an invertible
     7216        matrix $S$ such that $A=S^{-1}BS$.  $S$ can be interpreted as a
     7217        change-of-basis matrix if $A$ and $B$ are viewed as matrix
     7218        representations of the same linear transformation.
     7219
     7220        When ``transformation=False`` this method will return ``True`` if
     7221        such a matrix $S$ exists, otherwise it will return ``False``.  When
     7222        ``transformation=True`` the method returns a pair.  The first part
     7223        of the pair is ``True`` or ``False`` depending on if the matrices
     7224        are similar and the second part is the change-of-basis matrix, or
     7225        ``None`` should it not exist.
     7226
     7227        When the transformation matrix is requested, it will satisfy
     7228        ``self = S.inverse()*other*S``.
     7229
     7230        If the base rings for any of the matrices is the integers, the
     7231        rationals, or the field of algebraic numbers (``QQbar``), then the
     7232        matrices are converted to have ``QQbar`` as their base ring prior
     7233        to checking the equality of the base rings.
     7234
     7235        It is possible for this routine to fail over most fields, even when
     7236        the matrices are similar.  However, since the field of algebraic
     7237        numbers is algebraically closed, the routine will always produce
     7238        a result for matrices with rational entries.
     7239
     7240        EXAMPLES:
     7241
     7242        The two matrices in this example were constructed to be similar.
     7243        The computations happen in the field of algebraic numbers, but we
     7244        are able to convert the change-of-basis matrix back to the rationals
     7245        (which may not always be possible). ::
     7246
     7247            sage: A = matrix(ZZ, [[-5, 2, -11],
     7248            ...                   [-6, 7, -42],
     7249            ...                   [0, 1, -6]])
     7250            sage: B = matrix(ZZ, [[ 1, 12,  3],
     7251            ...                   [-1, -6, -1],
     7252            ...                   [ 0,  6,  1]])
     7253            sage: A.is_similar(B)
     7254            True
     7255            sage: _, T = A.is_similar(B, transformation=True)
     7256            sage: T
     7257            [  1.000000000000000? + 0.?e-32*I              0.?e-33 + 0.?e-33*I              0.?e-32 + 0.?e-32*I]
     7258            [-0.6666666666666667? + 0.?e-33*I  0.1666666666666667? + 0.?e-33*I -0.8333333333333333? + 0.?e-32*I]
     7259            [ 0.6666666666666667? + 0.?e-32*I              0.?e-33 + 0.?e-33*I -0.3333333333333334? + 0.?e-32*I]
     7260            sage: T.change_ring(QQ)
     7261            [   1    0    0]
     7262            [-2/3  1/6 -5/6]
     7263            [ 2/3    0 -1/3]
     7264            sage: A == T.inverse()*B*T
     7265            True
     7266
     7267        Other exact fields are supported.  ::
     7268
     7269            sage: F.<a> = FiniteField(7^2)
     7270            sage: A = matrix(F,[[2*a + 5, 6*a + 6,   a + 3],
     7271            ...                 [  a + 3, 2*a + 2, 4*a + 2],
     7272            ...                 [2*a + 6, 5*a + 5,     3*a]])
     7273            sage: B = matrix(F,[[5*a + 5, 6*a + 4,   a + 1],
     7274            ...                 [  a + 5, 4*a + 3, 3*a + 3],
     7275            ...                 [3*a + 5,   a + 4, 5*a + 6]])
     7276            sage: A.is_similar(B)
     7277            True
     7278            sage: B.is_similar(A)
     7279            True
     7280            sage: _, T = A.is_similar(B, transformation=True)
     7281            sage: T
     7282            [      1       0       0]
     7283            [6*a + 1 4*a + 3 4*a + 2]
     7284            [6*a + 3 3*a + 5 3*a + 6]
     7285            sage: A == T.inverse()*B*T
     7286            True
     7287
     7288        Two matrices with different sets of eigenvalues, so they
     7289        cannot possibly be similar. ::
     7290
     7291            sage: A = matrix(QQ, [[ 2,  3, -3, -6],
     7292            ...                   [ 0,  1, -2, -8],
     7293            ...                   [-3, -3,  4,  3],
     7294            ...                   [-1, -2,  2,  6]])
     7295            sage: B = matrix(QQ, [[ 1,  1,  2,  4],
     7296            ...                   [-1,  2, -3, -7],
     7297            ...                   [-2,  3, -4, -7],
     7298            ...                   [ 0, -1,  0,  0]])
     7299            sage: A.eigenvalues() == B.eigenvalues()
     7300            False
     7301            sage: A.is_similar(B, transformation=True)
     7302            (False, None)
     7303
     7304        Similarity is an equivalence relation, so this routine computes
     7305        a representative of the equivalence class for each matrix, the
     7306        Jordan form, as provided by :meth:`jordan_form`.  The matrices
     7307        below have identical eigenvalues (as evidenced by equal
     7308        characteristic polynomials), but slightly different Jordan forms,
     7309        and hence are not similar.
     7310
     7311            sage: A = matrix(QQ, [[ 19, -7, -29],
     7312            ...                   [-16, 11,  30],
     7313            ...                   [ 15, -7, -25]])
     7314            sage: B = matrix(QQ, [[-38, -63,  42],
     7315            ...                   [ 14,  25, -14],
     7316            ...                   [-14, -21,  18]])
     7317            sage: A.charpoly() == B.charpoly()
     7318            True
     7319            sage: A.jordan_form()
     7320            [-3| 0  0]
     7321            [--+-----]
     7322            [ 0| 4  1]
     7323            [ 0| 0  4]
     7324            sage: B.jordan_form()
     7325            [-3| 0| 0]
     7326            [--+--+--]
     7327            [ 0| 4| 0]
     7328            [--+--+--]
     7329            [ 0| 0| 4]
     7330            sage: A.is_similar(B)
     7331            False
     7332
     7333        Obtaining the Jordan form requires computing the eigenvalues of
     7334        the matrix, which may not lie in the field used for entries of
     7335        the matrix.  So the routine first checks the characteristic
     7336        polynomials - if they are unequal, then the matrices cannot be
     7337        similar. However, when the characteristic polynomials are equal,
     7338        we must examine the Jordan form. In this case, the method may fail,
     7339        EVEN when the matrices are similar.  This is not the case for
     7340        matrices over the integers, rationals or algebraic numbers,
     7341        since the computations are done in the algebraically closed
     7342        field of algebraic numbers.
     7343
     7344        Here is an example where the similarity is obvious, but the
     7345        routine fails to compute a result.  ::
     7346
     7347            sage: F.<a> = FiniteField(7^2)
     7348            sage: C = matrix(F,[[  a + 2, 5*a + 4],
     7349            ...                 [6*a + 6, 6*a + 4]])
     7350            sage: S = matrix(ZZ, [[0, 1],
     7351            ...                   [1, 0]])
     7352            sage: D = S.inverse()*C*S
     7353            sage: C.is_similar(D)
     7354            Traceback (most recent call last):
     7355            ...
     7356            ValueError: unable to compute Jordan canonical form for a matrix
     7357            sage: C.jordan_form()
     7358            Traceback (most recent call last):
     7359            ...
     7360            RuntimeError: Some eigenvalue does not exist in Finite Field in a of size 7^2.
     7361
     7362        Inexact rings and fields are also not supported.  ::
     7363
     7364            sage: A = matrix(CDF, 2, 2, range(4))
     7365            sage: B = copy(A)
     7366            sage: A.is_similar(B)
     7367            Traceback (most recent call last):
     7368            ...
     7369            ValueError: unable to compute Jordan canonical form for a matrix
     7370
     7371        Rectangular matrices and mismatched sizes return quickly.  ::
     7372
     7373            sage: A = matrix(3, 2, range(6))
     7374            sage: B = copy(A)
     7375            sage: A.is_similar(B)
     7376            False
     7377            sage: A = matrix(2, 2, range(4))
     7378            sage: B = matrix(3, 3, range(9))
     7379            sage: A.is_similar(B, transformation=True)
     7380            (False, None)
     7381
     7382        If the fraction fields of the entries are unequal, it is an error,
     7383        except in the case when the rationals gets promoted to the
     7384        algebraic numbers.  ::
     7385
     7386            sage: A = matrix(ZZ, 2, 2, range(4))
     7387            sage: B = matrix(GF(2), 2, 2, range(4))
     7388            sage: A.is_similar(B, transformation=True)
     7389            Traceback (most recent call last):
     7390            ...
     7391            TypeError: matrices need to have entries with identical fraction fields, not Algebraic Field and Finite Field of size 2
     7392            sage: A = matrix(ZZ, 2, 2, range(4))
     7393            sage: B = matrix(QQbar, 2, 2, range(4))
     7394            sage: A.is_similar(B)
     7395            True
     7396
     7397        Inputs are checked.  ::
     7398
     7399            sage: A = matrix(ZZ, 2, 2, range(4))
     7400            sage: A.is_similar('garbage')
     7401            Traceback (most recent call last):
     7402            ...
     7403            TypeError: similarity requires a matrix as an argument, not garbage
     7404            sage: B = copy(A)
     7405            sage: A.is_similar(B, transformation='junk')
     7406            Traceback (most recent call last):
     7407            ...
     7408            ValueError: transformation keyword must be True or False, not junk
     7409        """
     7410        import sage.matrix.matrix
     7411        import sage.rings.qqbar
     7412        if not sage.matrix.matrix.is_Matrix(other):
     7413            raise TypeError('similarity requires a matrix as an argument, not {0}'.format(other))
     7414        if transformation not in [True, False]:
     7415            raise ValueError('transformation keyword must be True or False, not {0}'.format(transformation))
     7416        # easy false situations
     7417        if not self.is_square() or not other.is_square:
     7418            if transformation:
     7419                return (False, None)
     7420            else:
     7421                return False
     7422        if self.nrows() != self.nrows():
     7423            if transformation:
     7424                return (False, None)
     7425            else:
     7426                return False
     7427        # convert to fraction fields for base rings
     7428        A = self.matrix_over_field()
     7429        B = other.matrix_over_field()
     7430        # move rationals to algebraically closed algebraic numbers
     7431        if A.base_ring() == QQ:
     7432            A = A.change_ring(sage.rings.qqbar.QQbar)
     7433        if B.base_ring() == QQ:
     7434            B = B.change_ring(sage.rings.qqbar.QQbar)
     7435        # require identical base fields
     7436        if A.base_ring() != B.base_ring():
     7437            raise TypeError('matrices need to have entries with identical fraction fields, not {0} and {1}'.format(A.base_ring(), B.base_ring()))
     7438        # unequal characteristic polynomials implies not similar
     7439        # and avoids any problems with eigenvalues not in the base field
     7440        if A.charpoly() != B.charpoly():
     7441            if transformation:
     7442                return (False, None)
     7443            else:
     7444                return False
     7445        # now more precisely compare Jordan form, and optionally get transformations
     7446        try:
     7447            if transformation:
     7448                JA, SA = A.jordan_form(transformation=True)
     7449            else:
     7450                JA = A.jordan_form(transformation=False)
     7451        except:
     7452            raise ValueError('unable to compute Jordan canonical form for a matrix')
     7453        try:
     7454            if transformation:
     7455                JB, SB = B.jordan_form(transformation=True)
     7456            else:
     7457                JB = B.jordan_form(transformation=False)
     7458        except:
     7459            raise ValueError('unable to compute Jordan canonical form for a matrix')
     7460        similar = (JA == JB)
     7461        transform = None
     7462        if similar and transformation:
     7463            transform = SB*SA.inverse()
     7464        if transformation:
     7465            return (similar, transform)
     7466        else:
     7467            return similar
     7468
    71987469    def symplectic_form(self):
    71997470        r"""
    72007471        Find a symplectic form for self if self is an anti-symmetric,