Ticket #10791: trac_10791-fix-gram-schmidt-v2.patch

File trac_10791-fix-gram-schmidt-v2.patch, 26.8 KB (added by rbeezer, 10 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1298770076 28800
    # Node ID 7846e3628f42a420092e5b972e3d32e1a8d8a17d
    # Parent  2c5ad01c9e18f83ec0086dc382f60d5c8041f3b4
    10791: fix and expand Gram-Schmidt vector orthogonalization
    
    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    73837383            self.cache('QR_factors', QR)
    73847384        return QR
    73857385
    7386     def gram_schmidt(self):
     7386    def _gram_schmidt_noscale(self):
    73877387        r"""
    7388         Return the matrix G whose rows are obtained from the rows of self
    7389         (=A) by applying the Gram-Schmidt orthogonalization process. Also
    7390         return the coefficients mu ij, i.e., a matrix mu such that
    7391         ``(mu + 1)*G == A``.
    7392        
     7388        Performs Gram-Schmidt orthogonalization, with no scaling to unit vectors.
     7389
     7390        INPUT:
     7391
     7392        - ``self`` - is a matrix whose columns are to be orthogonalized.
     7393          The base ring of the matrix needs to have its fraction field
     7394          implemented.
     7395
    73937396        OUTPUT:
    7394        
    7395        
    7396         - ``G`` - a matrix whose rows are orthogonal
    7397        
    7398         - ``mu`` - a matrix that gives the transformation, via
    7399           the relation ``(mu + 1)*G == self``
    7400        
    7401        
     7397
     7398        Two matrices, ``Q`` and ``R`` such that if ``A`` represents ``self``:
     7399
     7400        - ``A = Q*R``
     7401        - The columns of ``Q`` are an orthogonal set which spans the
     7402          column space of ``A``.
     7403        - The conjugate-transpose of ``Q`` times ``Q`` is a diagonal matrix.
     7404        - ``R`` is a full-rank matrix, that has all entries below the
     7405          main diagonal equal to zero.
     7406
     7407        This is basically a "reduced" QR decomposition of ``self`` with
     7408        the distinction that the orthogonal column vectors of ``Q`` have
     7409        not been scaled to unit vectors, avoiding the need to take square
     7410        roots.
     7411
     7412        EXAMPLES:
     7413
     7414        A rectangular matrix whose columns have full-rank.  Notice that the
     7415        routine computes in the fraction field, requiring the column space
     7416        check to step up to ``QQ``.  ::
     7417
     7418            sage: A = matrix(ZZ, [[-1, -3,  0, -1],
     7419            ...                   [ 1,  2, -1,  2],
     7420            ...                   [-3, -6,  4, -7]])
     7421            sage: Q,R = A._gram_schmidt_noscale()
     7422            sage: Q
     7423            [    -1 -10/11      0]
     7424            [     1  -1/11   3/10]
     7425            [    -3   3/11   1/10]
     7426            sage: R
     7427            [     1  23/11 -13/11  24/11]
     7428            [     0      1  13/10 -13/10]
     7429            [     0      0      1     -1]
     7430            sage: Q*R == A
     7431            True
     7432            sage: Q.transpose()*Q
     7433            [   11     0     0]
     7434            [    0 10/11     0]
     7435            [    0     0  1/10]
     7436            sage: A.change_ring(QQ).column_space() == Q.column_space()
     7437            True
     7438
     7439        A matrix over a subfield of the complex numbers, with four
     7440        columns but rank 3, so the orthogonal set has just 3 vectors
     7441        as well.  Orthogonality comes from the Hermitian inner product
     7442        so we need to check with the conjugate-transpose.  This
     7443        example verifies that the bug on #10791 is fixed.  ::
     7444
     7445            sage: F.<a> = QuadraticField(-5)
     7446            sage: A = matrix(F, [[    1,   a - 3,   a - 2, a + 1],
     7447            ...                  [    a, 2*a + 1, 3*a + 1,     1],
     7448            ...                  [a + 1,   a - 6, 2*a - 5,     1],
     7449            ...                  [  2*a,       a,     3*a,    -3],
     7450            ...                  [    1,   a - 1,       a, a - 1]])
     7451            sage: A.rank()
     7452            3
     7453            sage: Q, R = A._gram_schmidt_noscale()
     7454            sage: Q
     7455            [                      1         25/33*a - 38/11   641/1163*a + 453/1163]
     7456            [                      a         17/11*a + 73/33  322/1163*a + 1566/1163]
     7457            [                  a + 1        10/33*a - 173/33 -784/1163*a + 1614/1163]
     7458            [                    2*a          1/11*a + 80/33  196/1163*a - 1234/1163]
     7459            [                      1         25/33*a - 16/11  855/1163*a - 1717/1163]
     7460            sage: R
     7461            [                    1         8/33*a + 5/11        8/33*a + 16/11         2/11*a + 1/33]
     7462            [                    0                     1                     1 -107/1163*a - 78/1163]
     7463            [                    0                     0                     0                     1]
     7464            sage: Q*R == A
     7465            True
     7466            sage: Q.transpose().conjugate()*Q
     7467            [        33          0          0]
     7468            [         0    2326/33          0]
     7469            [         0          0 16532/1163]
     7470            sage: Q.column_space() == A.column_space()
     7471            True
     7472
     7473        Some trivial cases.  ::
     7474
     7475            sage: A = matrix(ZZ, 3, 0)
     7476            sage: Q, R = A._gram_schmidt_noscale()
     7477            sage: Q.parent()
     7478            Full MatrixSpace of 3 by 0 dense matrices over Rational Field
     7479            sage: R.parent()
     7480            Full MatrixSpace of 0 by 0 dense matrices over Rational Field
     7481            sage: Q*R == A
     7482            True
     7483
     7484            sage: A = matrix(ZZ, 0, 3)
     7485            sage: Q, R = A._gram_schmidt_noscale()
     7486            sage: Q.parent()
     7487            Full MatrixSpace of 0 by 0 dense matrices over Rational Field
     7488            sage: R.parent()
     7489            Full MatrixSpace of 0 by 3 dense matrices over Rational Field
     7490            sage: Q*R == A
     7491            True
     7492
     7493        TESTS:
     7494
     7495        Without a fraction field, we cannot hope to proceed. ::
     7496
     7497            sage: A = matrix(Integers(6), 2, range(4))
     7498            sage: A._gram_schmidt_noscale()
     7499            Traceback (most recent call last):
     7500            ...
     7501            TypeError: Gram-Schmidt orthogonalization requires a base ring with a fraction field, not Ring of integers modulo 6
     7502
     7503        AUTHORS:
     7504
     7505        - William Stein (2007-11-18)
     7506        - Rob Beezer (2011-02-25)
     7507        """
     7508        import sage.matrix.constructor
     7509        R = self.base_ring()
     7510        try:
     7511            F = R.fraction_field()
     7512        except TypeError:
     7513            raise TypeError("Gram-Schmidt orthogonalization requires a base ring with a fraction field, not %s" % R)
     7514        n = self.ncols()
     7515        B = self.columns()
     7516        zero = F(0)
     7517        Bstar = []
     7518        R = sage.matrix.constructor.zero_matrix(F, n)
     7519        nnz = 0 # number non-zero rows in R, or number of nonzero vectors in Bstar
     7520        for i in range(n):
     7521            ortho = B[i]
     7522            for j in range(nnz):
     7523                R[j,i] = Bstar[j].hermitian_inner_product(B[i])/Bstar[j].hermitian_inner_product(Bstar[j])
     7524                ortho = ortho - R[j,i]*Bstar[j]
     7525            if ortho.hermitian_inner_product(ortho) != zero:
     7526                Bstar.append(ortho)
     7527                R[nnz, i] = 1
     7528                nnz = nnz + 1
     7529        R = R[0:nnz]
     7530        if Bstar == []:
     7531            Q = sage.matrix.constructor.matrix(F, 0, self.nrows()).transpose()
     7532        else:
     7533            Q = sage.matrix.constructor.matrix(F, Bstar).transpose()
     7534        return Q, R
     7535
     7536    def gram_schmidt(self, orthonormal=False):
     7537        r"""
     7538        Performs Gram-Schmidt orthogonalization on the rows of the matrix,
     7539        returning a new matrix and a matrix accomplishing the transformation.
     7540
     7541        INPUT:
     7542
     7543        - ``self`` - a matrix whose rows are to be orthogonalized.
     7544        - ``orthonormal`` - default: ``False`` - if ``True`` the
     7545          returned orthogonal vectors are unit vectors.  This keyword
     7546          is ignored if the matrix is over ``RDF`` or ``CDF`` and the
     7547          results are always orthonormal.
     7548
     7549        OUTPUT:
     7550
     7551        A pair of matrices, ``G`` and ``M`` such that if ``A``
     7552        represents ``self``, where the parenthetical properties occur
     7553        when ``orthonormal = True``:
     7554
     7555        - ``A = M*G``
     7556        - The rows of ``G`` are an orthogonal (resp. orthonormal)
     7557          set of vectors.
     7558        - ``G`` times the conjugate-transpose of ``G`` is a diagonal
     7559          (resp. identity) matrix.
     7560        - The row space of ``G`` equals the row space of ``A``.
     7561        - ``M`` is a full-rank matrix with zeros above the diagonal.
     7562
     7563        For exact rings, any zero vectors produced (when the original
     7564        vectors are linearly dependent) are not output, thus the
     7565        orthonormal set is linearly independent, and thus a basis for the
     7566        row space of the original matrix.
     7567
     7568        Any notion of a Gram-Schmidt procedure requires that the base
     7569        ring of the matrix has a fraction field implemented.  In order
     7570        to arrive at an orthonormal set, it must be possible to construct
     7571        square roots of the elements of the base field.  In Sage, your
     7572        best option is the field of algebraic numbers, ``QQbar``, which
     7573        properly contain the rationals and number fields.
     7574
     7575        If you have an approximate numerical matrix, then this routine
     7576        requires that your base field be the real and complex
     7577        double-precision floating point numbers, ``RDF`` and ``CDF``.
     7578        In this case, the matrix is treated as having full rank, as no
     7579        attempt is made to recognize linear dependence with approximate
     7580        calculations.
     7581
    74027582        EXAMPLES::
    7403        
     7583
     7584        Inexact Rings, Numerical Matrices:
     7585
     7586        First, the inexact rings, ``CDF`` and ``RDF``.  ::
     7587
     7588            sage: A = matrix(CDF, [[ 0.6454 + 0.7491*I, -0.8662 + 0.1489*I,  0.7656 - 0.00344*I],
     7589            ...                    [-0.2913 + 0.8057*I,  0.8321 + 0.8170*I, -0.6744 + 0.9248*I],
     7590            ...                    [ 0.2554 + 0.3517*I, -0.4454 - 0.1715*I,  0.8325 - 0.6282*I]])
     7591            sage: G, M = A.gram_schmidt()
     7592            sage: G.round(6)
     7593            [-0.422243 - 0.490087*I  0.566698 - 0.097416*I -0.500882 + 0.002251*I]
     7594            [-0.057002 - 0.495035*I  -0.35059 - 0.625323*I  0.255514 - 0.415284*I]
     7595            [ 0.394105 - 0.421778*I -0.392266 - 0.039345*I  -0.352905 + 0.62195*I]
     7596            sage: M.round(6)
     7597            [             -1.528503                      0                      0]
     7598            [  0.459974 - 0.40061*I              -1.741233                      0]
     7599            [-0.934304 + 0.148868*I   0.54833 + 0.073202*I              -0.550725]
     7600            sage: (A - M*G).zero_at(10^-12)
     7601            [0 0 0]
     7602            [0 0 0]
     7603            [0 0 0]
     7604            sage: (G*G.conjugate().transpose()).zero_at(10^-12)
     7605            [1.0   0   0]
     7606            [  0 1.0   0]
     7607            [  0   0 1.0]
     7608
     7609        A rectangular matrix.  Note that the ``orthonormal`` keyword
     7610        is ignored in these cases.  ::
     7611
     7612            sage: A = matrix(RDF, [[-0.978325, -0.751994, 0.925305, -0.200512, 0.420458],
     7613            ...                    [-0.474877, -0.983403, 0.089836,  0.132218, 0.672965]])
     7614            sage: G, M = A.gram_schmidt(orthonormal=False)
     7615            sage: G.round(6)
     7616            [-0.607223 -0.466745  0.574315 -0.124453  0.260968]
     7617            [ 0.123203 -0.617909 -0.530578  0.289773  0.487368]
     7618            sage: M.round(6)
     7619            [1.611147      0.0]
     7620            [0.958116 0.867778]
     7621            sage: (A-M*G).zero_at(10^-12)
     7622            [0.0 0.0 0.0 0.0 0.0]
     7623            [0.0 0.0 0.0 0.0 0.0]
     7624            sage: (G*G.transpose()).zero_at(10^-12)
     7625            [1.0 0.0]
     7626            [0.0 1.0]
     7627
     7628        Even though a set of vectors may be linearly dependent, no effort
     7629        is made to decide when a zero vector is really the result of a
     7630        relation of linear dependence.  So in this regard, input matrices
     7631        are treated as being of full rank.  Try one of the base rings that
     7632        provide exact results if you need exact results.  ::
     7633
     7634            sage: entries = [[1,1,2], [2,1,3], [3,1,4]]
     7635            sage: A = matrix(QQ, entries)
     7636            sage: A.rank()
     7637            2
     7638            sage: B = matrix(RDF, entries)
     7639            sage: G, M = B.gram_schmidt()
     7640            sage: G.round(6)
     7641            [-0.408248 -0.408248 -0.816497]
     7642            [ 0.707107 -0.707107      -0.0]
     7643            [ -0.57735  -0.57735   0.57735]
     7644            sage: M.round(10)
     7645            [-2.4494897428           0.0           0.0]
     7646            [-3.6742346142  0.7071067812           0.0]
     7647            [-4.8989794856  1.4142135624           0.0]
     7648
     7649        Exact Rings, Orthonormalization:
     7650
     7651        To scale a vector to unit length requires taking
     7652        a square root, which often takes us outside the base ring.
     7653        For the integers,  and rationals the field of algebraic numbers,
     7654        ``QQbar``, is big enough to contain what we need, but the price
     7655        is that the computations are very slow, hence mostly of value
     7656        for small cases or instruction. Now we need to use the
     7657        ``orthonormal`` keyword.  ::
     7658
     7659            sage: A = matrix(QQbar, [[6, -8,  1],
     7660            ...                      [4,  1,  3],
     7661            ...                      [6,  3,  3],
     7662            ...                      [7,  1, -5],
     7663            ...                      [7, -3,  5]])
     7664            sage: G, M = A.gram_schmidt(orthonormal=True)
     7665            sage: G
     7666            [ 0.5970223141259934? -0.7960297521679913? 0.09950371902099891?]
     7667            [ 0.6063218341690895?  0.5289635311888953?  0.5937772444966257?]
     7668            [ 0.5252981913594170?  0.2941669871612735?  -0.798453250866314?]
     7669            sage: M
     7670            [ 10.04987562112089?                   0                   0]
     7671            [ 1.890570661398980?  4.735582601355131?                   0]
     7672            [ 1.492555785314984?  7.006153332071100?  1.638930357041381?]
     7673            [ 2.885607851608969?  1.804330147889395?  7.963520581008761?]
     7674            [ 7.064764050490923?  5.626248468100069? -1.197679876299471?]
     7675            sage: M*G-A
     7676            [0 0 0]
     7677            [0 0 0]
     7678            [0 0 0]
     7679            [0 0 0]
     7680            [0 0 0]
     7681            sage: G*G.transpose()
     7682            [      1 0.?e-37 0.?e-35]
     7683            [0.?e-37       1 0.?e-35]
     7684            [0.?e-35 0.?e-35       1]
     7685            sage: G.row_space() == A.row_space()
     7686            True
     7687
     7688        Starting with complex numbers with rational real and imaginary parts.
     7689        Note the use of the conjugate-transpose when checking the
     7690        orthonormality. ::
     7691
     7692            sage: A = matrix(QQbar, [[  -2,    -I - 1, 4*I + 2,       -1],
     7693            ...                      [-4*I, -2*I + 17,       0,  9*I + 1],
     7694            ...                      [   1,  -2*I - 6, -I + 11, -5*I + 1]])
     7695            sage: G, M = A.gram_schmidt(orthonormal=True)
     7696            sage: G
     7697            [                         -0.3849001794597505?  -0.1924500897298753? - 0.1924500897298753?*I   0.3849001794597505? + 0.7698003589195010?*I                          -0.1924500897298753?]
     7698            [-0.06165497274852388? - 0.1387236886841787?*I  0.8188551068163327? - 0.10018933071635130?*I  0.2003786614327026? + 0.05394810115495839?*I  0.02119389688230508? + 0.5028733714801478?*I]
     7699            [  0.3842387256410419? - 0.5694103019142261?*I  0.1416892863096208? - 0.06139779741542298?*I  0.4633778333528464? - 0.01285039016180503?*I  0.02658516588219101? - 0.5373044261814995?*I]
     7700            sage: M
     7701            [                        5.196152422706632?                                          0                                          0]
     7702            [-3.079201435678004? + 3.464101615137755?*I             19.22286447225071? + 0.?e-37*I                                          0]
     7703            [  4.426352063787131? - 8.66025403784439?*I -5.117362738127481? - 3.502773139275513?*I             7.480012456446966? + 0.?e-35*I]
     7704            sage: M*G-A
     7705            [            0.?e-37 0.?e-37 + 0.?e-37*I 0.?e-37 + 0.?e-37*I             0.?e-37]
     7706            [0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-36*I 0.?e-36 + 0.?e-36*I 0.?e-36 + 0.?e-36*I]
     7707            [0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I]
     7708            sage: G*G.conjugate().transpose()
     7709            [1.000000000000000? + 0.?e-37*I            0.?e-37 + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
     7710            [           0.?e-37 + 0.?e-37*I 1.000000000000000? + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
     7711            [           0.?e-36 + 0.?e-36*I            0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-35*I]
     7712            sage: G.row_space() == A.row_space()
     7713            True
     7714
     7715        A square matrix with small rank.  The zero vectors produced as a
     7716        result of linear dependence get eliminated, so the rows of ``G``
     7717        are a basis for the row space of ``A``.  ::
     7718
     7719            sage: A = matrix(QQbar, [[2, -6, 3, 8],
     7720            ...                      [1, -3, 2, 5],
     7721            ...                      [0,  0, 2, 4],
     7722            ...                      [2, -6, 3, 8]])
     7723            sage: A.change_ring(QQ).rank()
     7724            2
     7725            sage: G, M = A.gram_schmidt(orthonormal=True)
     7726            sage: G
     7727            [ 0.1881441736767195? -0.5644325210301583?  0.2822162605150792?  0.7525766947068779?]
     7728            [-0.2502818123591464?   0.750845437077439?  0.3688363550555841?  0.4873908977520218?]
     7729            sage: M
     7730            [10.630145812734649?                   0]
     7731            [ 6.208757731331742? 0.6718090752798139?]
     7732            [ 3.574739299857670?  2.687236301119256?]
     7733            [10.630145812734649?                   0]
     7734            sage: M*G-A
     7735            [0 0 0 0]
     7736            [0 0 0 0]
     7737            [0 0 0 0]
     7738            [0 0 0 0]
     7739            sage: G*G.transpose()
     7740            [      1 0.?e-35]
     7741            [0.?e-35       1]
     7742            sage: G.row_space() == A.row_space()
     7743            True
     7744
     7745        Exact Rings, Orthogonalization:
     7746
     7747        If we forego scaling orthogonal vectors to unit vectors, we
     7748        can apply Gram-Schmidt to a much greater variety of rings.
     7749        Use the ``orthonormal=False`` keyword (or assume it as the default).
     7750        Note that now the orthogonality check creates a diagonal matrix
     7751        whose diagonal entries are the squares of the lengths of the
     7752        vectors. ::
     7753
     7754        First, in the rationals, without involving ``QQbar``.  ::
     7755
     7756            sage: A = matrix(QQ, [[-1,  3,  2,  2],
     7757            ...                   [-1,  0, -1,  0],
     7758            ...                   [-1, -2, -3, -1],
     7759            ...                   [ 1,  1,  2,  0]])
     7760            sage: A.rank()
     7761            3
     7762            sage: G, M = A.gram_schmidt()
     7763            sage: G
     7764            [    -1      3      2      2]
     7765            [-19/18    1/6   -8/9    1/9]
     7766            [  2/35  -4/35  -2/35   9/35]
     7767            sage: M
     7768            [     1      0      0]
     7769            [ -1/18      1      0]
     7770            [-13/18  59/35      1]
     7771            [   1/3 -48/35     -2]
     7772            sage: M*G-A
     7773            [0 0 0 0]
     7774            [0 0 0 0]
     7775            [0 0 0 0]
     7776            [0 0 0 0]
     7777            sage: G*G.transpose()
     7778            [   18     0     0]
     7779            [    0 35/18     0]
     7780            [    0     0  3/35]
     7781            sage: G.row_space() == A.row_space()
     7782            True
     7783
     7784        A complex subfield of the complex numbers.  ::
     7785
     7786            sage: C.<z> = CyclotomicField(5)
     7787            sage: A = matrix(C, [[              -z^3 - 2*z,             -z^3 - 1, 2*z^3 - 2*z^2 + 2*z,             1],
     7788            ...                  [         z^3 - 2*z^2 + 1, -z^3 + 2*z^2 - z - 1,                  -1,       z^2 + z],
     7789            ...                  [-1/2*z^3 - 2*z^2 + z + 1,         -z^3 + z - 2,    -2*z^3 + 1/2*z^2, 2*z^2 - z + 2]])
     7790            sage: G, M = A.gram_schmidt(orthonormal=False)
     7791            sage: G
     7792            [                                                      -z^3 - 2*z                                                         -z^3 - 1                                              2*z^3 - 2*z^2 + 2*z                                                                1]
     7793            [                   155/139*z^3 - 161/139*z^2 + 31/139*z + 13/139                 -175/139*z^3 + 180/139*z^2 - 125/139*z - 142/139                     230/139*z^3 + 124/139*z^2 + 6/139*z + 19/139                      -14/139*z^3 + 92/139*z^2 - 6/139*z - 95/139]
     7794            [-10359/19841*z^3 - 36739/39682*z^2 + 24961/39682*z - 11879/39682  -28209/39682*z^3 - 3671/19841*z^2 + 51549/39682*z - 38613/39682    -42769/39682*z^3 - 615/39682*z^2 - 1252/19841*z - 14392/19841   4895/19841*z^3 + 57885/39682*z^2 - 46094/19841*z + 65747/39682]
     7795            sage: M
     7796            [                                                           1                                                            0                                                            0]
     7797            [                14/139*z^3 + 47/139*z^2 + 145/139*z + 95/139                                                            1                                                            0]
     7798            [              -7/278*z^3 + 199/278*z^2 + 183/139*z + 175/278 -3785/39682*z^3 + 3346/19841*z^2 - 3990/19841*z + 2039/19841                                                            1]
     7799            sage: M*G - A
     7800            [0 0 0 0]
     7801            [0 0 0 0]
     7802            [0 0 0 0]
     7803            sage: G*G.conjugate().transpose()
     7804            [                               15*z^3 + 15*z^2 + 28                                                   0                                                   0]
     7805            [                                                  0                463/139*z^3 + 463/139*z^2 + 1971/139                                                   0]
     7806            [                                                  0                                                   0 230983/19841*z^3 + 230983/19841*z^2 + 1003433/39682]
     7807            sage: G.row_space() == A.row_space()
     7808            True
     7809
     7810        A slightly edited legacy example.  ::
     7811
    74047812            sage: A = matrix(ZZ, 3, [-1, 2, 5, -11, 1, 1, 1, -1, -3]); A
    74057813            [ -1   2   5]
    74067814            [-11   1   1]
     
    74117819            [  -52/5    -1/5      -2]
    74127820            [  2/187  36/187 -14/187]
    74137821            sage: mu
    7414             [     0      0      0]
    7415             [   3/5      0      0]
    7416             [  -3/5 -7/187      0]
     7822            [     1      0      0]
     7823            [   3/5      1      0]
     7824            [  -3/5 -7/187      1]
    74177825            sage: G.row(0) * G.row(1)
    74187826            0
    74197827            sage: G.row(0) * G.row(2)
    74207828            0
    74217829            sage: G.row(1) * G.row(2)
    74227830            0
    7423        
    7424         The relation between mu and A is as follows::
    7425        
    7426             sage: (mu + 1)*G == A
    7427             True
    7428         """       
    7429         from sage.modules.misc import gram_schmidt
    7430         from constructor import matrix
    7431         Bstar, mu = gram_schmidt(self.rows())
    7432         return matrix(Bstar), mu
     7831
     7832        The relation between mu and A is as follows.  ::
     7833
     7834            sage: mu*G == A
     7835            True
     7836        """
     7837        import sage.rings.real_double
     7838        import sage.rings.complex_double
     7839        R = self.base_ring()
     7840        if R in [sage.rings.real_double.RDF, sage.rings.complex_double.CDF]:
     7841            Q, R = self.transpose().QR()
     7842            m = R.nrows(); n = R.ncols()
     7843            if m > n:
     7844                Q = Q[0:m, 0:n]
     7845                R = R[0:n, 0:n]
     7846        elif R.is_exact():
     7847            if orthonormal:
     7848                Q, R = self.transpose().QR(full=False)
     7849            else:
     7850                Q, R = self.transpose()._gram_schmidt_noscale()
     7851        else:
     7852            raise NotImplementedError("Gram-Scmidt orthogonalization not implemented for matrices over inexact rings, except for RDF and CDF")
     7853        return Q.transpose(), R.transpose()
    74337854
    74347855    def jordan_form(self, base_ring=None, sparse=False, subdivide=True, transformation=False):
    74357856        r"""
  • sage/matrix/matrix_integer_dense.pyx

    diff --git a/sage/matrix/matrix_integer_dense.pyx b/sage/matrix/matrix_integer_dense.pyx
    a b  
    28572857            raise TypeError("eta must be >= 0.5")
    28582858
    28592859        # this is pretty slow
    2860         G, mu = self.gram_schmidt()
    2861  
     2860        import sage.modules.misc
     2861        G, mu = sage.modules.misc.gram_schmidt(self.rows())
    28622862        #For any $i>j$, we have $|mu_{i, j}| <= \eta$
    28632863        for e in mu.list():
    28642864            if e.abs() > eta:
    28652865                return False
    28662866
    28672867        #For any $i<d$, we have $\delta |b_i^*|^2 <= |b_{i+1}^* + mu_{i+1, i} b_i^* |^2$
    2868         norms = [G[i].norm()**2 for i in range(G.nrows())]
     2868        norms = [G[i].norm()**2 for i in range(len(G))]
    28692869        for i in xrange(1,self.nrows()):
    28702870            if norms[i] < (delta - mu[i,i-1]**2) * norms[i-1]:
    28712871                return False
  • sage/modules/misc.py

    diff --git a/sage/modules/misc.py b/sage/modules/misc.py
    a b  
    1515from sage.matrix.constructor import matrix
    1616from sage.rings.integer_ring import ZZ
    1717
     18# Function below could be replicated into
     19# sage.matrix.matrix_integer_dense.Matrix_integer_dense.is_LLL_reduced
     20# which is its only current use (2011-02-26).  Then this could
     21# be deprecated and this file removed.
     22
    1823def gram_schmidt(B):
    19     """
     24    r"""
    2025    Return the Gram-Schmidt orthogonalization of the entries in the list
    2126    B of vectors, along with the matrix mu of Gram-Schmidt coefficients.
    2227
    2328    Note that the output vectors need not have unit length. We do this
    2429    to avoid having to extract square roots.
    2530
     31    .. note::
     32
     33        Use of this function is discouraged.  It fails on linearly
     34        dependent input and its output format is not as natural as it
     35        could be.  Instead, see :meth:`sage.matrix.matrix2.Matrix2.gram_schmidt`
     36        which is safer and more general-purpose.
     37
    2638    EXAMPLES:
     39
    2740        sage: B = [vector([1,2,1/5]), vector([1,2,3]), vector([-1,0,0])]
    2841        sage: from sage.modules.misc import gram_schmidt
    2942        sage: G, mu = gram_schmidt(B)
     
    4558        sage: a = matrix([[],[],[],[]])
    4659        sage: a.gram_schmidt()
    4760         ([], [])
     61
     62    Linearly dependent input leads to a zero dot product in a denominator.
     63    This shows that Trac #10791 is fixed. ::
     64
     65        sage: from sage.modules.misc import gram_schmidt
     66        sage: V = [vector(ZZ,[1,1]), vector(ZZ,[2,2]), vector(ZZ,[1,2])]
     67        sage: gram_schmidt(V)
     68        Traceback (most recent call last):
     69        ...
     70        ValueError: linearly dependent input for module version of Gram-Schmidt
    4871    """
     72    import sage.modules.free_module_element
    4973    if len(B) == 0 or len(B[0]) == 0:
    5074        return B, matrix(ZZ,0,0,[])
    5175    n = len(B)
    5276    Bstar = [B[0]]
    5377    K = B[0].base_ring().fraction_field()
     78    zero = sage.modules.free_module_element.vector(K, len(B[0]))
     79    if Bstar[0] == zero:
     80        raise ValueError("linearly dependent input for module version of Gram-Schmidt")
    5481    mu = matrix(K, n, n)
    5582    for i in range(1,n):
    5683        for j in range(i):
    5784            mu[i,j] = B[i].dot_product(Bstar[j]) / (Bstar[j].dot_product(Bstar[j]))
    5885        Bstar.append(B[i] - sum(mu[i,j]*Bstar[j] for j in range(i)))
     86        if Bstar[i] == zero:
     87            raise ValueError("linearly dependent input for module version of Gram-Schmidt")
    5988    return Bstar, mu
    60        
     89