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

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

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1298770076 28800
    # Node ID 0b859949124e0d1c7752dcaea7ce7b540e6405f0
    # Parent  210ad87f9cd6a1cb713da83b8ad35dcd217b4cb7
    10791: fix and expand Gram-Schmidt vector orthogonalization
    
    diff -r 210ad87f9cd6 -r 0b859949124e sage/matrix/matrix2.pyx
    a b  
    61746174            self.cache('QR_factors', QR)
    61756175        return QR
    61766176
    6177     def gram_schmidt(self):
     6177    def _gram_schmidt_noscale(self):
    61786178        r"""
    6179         Return the matrix G whose rows are obtained from the rows of self
    6180         (=A) by applying the Gram-Schmidt orthogonalization process. Also
    6181         return the coefficients mu ij, i.e., a matrix mu such that
    6182         ``(mu + 1)*G == A``.
    6183        
     6179        Performs Gram-Schmidt orthogonalization, with no scaling to unit vectors.
     6180
     6181        INPUT:
     6182
     6183        - ``self`` - is a matrix whose columns are to be orthogonalized.
     6184          The base ring of the matrix needs to have its fraction field
     6185          implemented.
     6186
    61846187        OUTPUT:
    6185        
    6186        
    6187         - ``G`` - a matrix whose rows are orthogonal
    6188        
    6189         - ``mu`` - a matrix that gives the transformation, via
    6190           the relation ``(mu + 1)*G == self``
    6191        
    6192        
     6188
     6189        Two matrices, ``Q`` and ``R`` such that if ``A`` represents ``self``:
     6190
     6191        - ``A = Q*R``
     6192        - The columns of ``Q`` are an orthogonal set which spans the
     6193          column space of ``A``.
     6194        - The conjugate-transpose of ``Q`` times ``Q`` is a diagonal matrix.
     6195        - ``R`` is a full-rank matrix, that has all entries below the
     6196          main diagonal equal to zero.
     6197
     6198        This is basically a "reduced" QR decomposition of ``self`` with
     6199        the distinction that the orthogonal column vectors of ``Q`` have
     6200        not been scaled to unit vectors, avoiding the need to take square
     6201        roots.
     6202
     6203        EXAMPLES:
     6204
     6205        A rectangular matrix whose columns have full-rank.  Notice that the
     6206        routine computes in the fraction field, requiring the column space
     6207        check to step up to ``QQ``.  ::
     6208
     6209            sage: A = matrix(ZZ, [[-1, -3,  0, -1],
     6210            ...                   [ 1,  2, -1,  2],
     6211            ...                   [-3, -6,  4, -7]])
     6212            sage: Q,R = A._gram_schmidt_noscale()
     6213            sage: Q
     6214            [    -1 -10/11      0]
     6215            [     1  -1/11   3/10]
     6216            [    -3   3/11   1/10]
     6217            sage: R
     6218            [     1  23/11 -13/11  24/11]
     6219            [     0      1  13/10 -13/10]
     6220            [     0      0      1     -1]
     6221            sage: Q*R == A
     6222            True
     6223            sage: Q.transpose()*Q
     6224            [   11     0     0]
     6225            [    0 10/11     0]
     6226            [    0     0  1/10]
     6227            sage: A.change_ring(QQ).column_space() == Q.column_space()
     6228            True
     6229
     6230        A matrix over a subfield of the complex numbers, with four
     6231        columns but rank 3, so the orthogonal set has just 3 vectors
     6232        as well.  Orthogonality comes from the Hermitian inner product
     6233        so we need to check with the conjugate-transpose.  This
     6234        example verifies that the bug on #10791 is fixed.  ::
     6235
     6236            sage: F.<a> = QuadraticField(-5)
     6237            sage: A = matrix(F, [[    1,   a - 3,   a - 2, a + 1],
     6238            ...                  [    a, 2*a + 1, 3*a + 1,     1],
     6239            ...                  [a + 1,   a - 6, 2*a - 5,     1],
     6240            ...                  [  2*a,       a,     3*a,    -3],
     6241            ...                  [    1,   a - 1,       a, a - 1]])
     6242            sage: A.rank()
     6243            3
     6244            sage: Q, R = A._gram_schmidt_noscale()
     6245            sage: Q
     6246            [                      1         25/33*a - 38/11   641/1163*a + 453/1163]
     6247            [                      a         17/11*a + 73/33  322/1163*a + 1566/1163]
     6248            [                  a + 1        10/33*a - 173/33 -784/1163*a + 1614/1163]
     6249            [                    2*a          1/11*a + 80/33  196/1163*a - 1234/1163]
     6250            [                      1         25/33*a - 16/11  855/1163*a - 1717/1163]
     6251            sage: R
     6252            [                    1         8/33*a + 5/11        8/33*a + 16/11         2/11*a + 1/33]
     6253            [                    0                     1                     1 -107/1163*a - 78/1163]
     6254            [                    0                     0                     0                     1]
     6255            sage: Q*R == A
     6256            True
     6257            sage: Q.transpose().conjugate()*Q
     6258            [        33          0          0]
     6259            [         0    2326/33          0]
     6260            [         0          0 16532/1163]
     6261            sage: Q.column_space() == A.column_space()
     6262            True
     6263
     6264        Some trivial cases.  ::
     6265
     6266            sage: A = matrix(ZZ, 3, 0)
     6267            sage: Q, R = A._gram_schmidt_noscale()
     6268            sage: Q.parent()
     6269            Full MatrixSpace of 3 by 0 dense matrices over Rational Field
     6270            sage: R.parent()
     6271            Full MatrixSpace of 0 by 0 dense matrices over Rational Field
     6272            sage: Q*R == A
     6273            True
     6274
     6275            sage: A = matrix(ZZ, 0, 3)
     6276            sage: Q, R = A._gram_schmidt_noscale()
     6277            sage: Q.parent()
     6278            Full MatrixSpace of 0 by 0 dense matrices over Rational Field
     6279            sage: R.parent()
     6280            Full MatrixSpace of 0 by 3 dense matrices over Rational Field
     6281            sage: Q*R == A
     6282            True
     6283
     6284        TESTS:
     6285
     6286        Without a fraction field, we cannot hope to proceed. ::
     6287
     6288            sage: A = matrix(Integers(6), 2, range(4))
     6289            sage: A._gram_schmidt_noscale()
     6290            Traceback (most recent call last):
     6291            ...
     6292            TypeError: Gram-Schmidt orthogonalization requires a base ring with a fraction field, not Ring of integers modulo 6
     6293
     6294        AUTHORS:
     6295
     6296        - William Stein (2007-11-18)
     6297        - Rob Beezer (2011-02-25)
     6298        """
     6299        import sage.matrix.constructor
     6300        R = self.base_ring()
     6301        try:
     6302            F = R.fraction_field()
     6303        except:
     6304            raise TypeError("Gram-Schmidt orthogonalization requires a base ring with a fraction field, not %s" % R)
     6305        n = self.ncols()
     6306        B = self.columns()
     6307        zero = F(0)
     6308        Bstar = []
     6309        R = sage.matrix.constructor.zero_matrix(F, n)
     6310        nnz = 0 # number non-zero rows in R, or number of nonzero vectors in Bstar
     6311        for i in range(n):
     6312            ortho = B[i]
     6313            for j in range(nnz):
     6314                R[j,i] = Bstar[j].hermitian_inner_product(B[i])/Bstar[j].hermitian_inner_product(Bstar[j])
     6315                ortho = ortho - R[j,i]*Bstar[j]
     6316            if ortho.hermitian_inner_product(ortho) != zero:
     6317                Bstar.append(ortho)
     6318                R[nnz, i] = 1
     6319                nnz = nnz + 1
     6320        R = R[0:nnz]
     6321        if Bstar == []:
     6322            Q = sage.matrix.constructor.matrix(F, 0, self.nrows()).transpose()
     6323        else:
     6324            Q = sage.matrix.constructor.matrix(F, Bstar).transpose()
     6325        return Q, R
     6326
     6327    def gram_schmidt(self, orthonormal=False):
     6328        r"""
     6329        Performs Gram-Schmidt orthogonalization on the rows of the matrix,
     6330        returning a new matrix and a matrix accomplishing the transformation.
     6331
     6332        INPUT:
     6333
     6334        - ``self`` - a matrix whose rows are to be orthogonalized.
     6335        - ``orthonormal`` - default: ``False`` - if ``True`` the
     6336          returned orthogonal vectors are unit vectors.  This keyword
     6337          is ignored if the matrix is over ``RDF`` or ``CDF`` and the
     6338          results are always orthonormal.
     6339
     6340        OUTPUT:
     6341
     6342        A pair of matrices, ``G`` and ``M`` such that if ``A``
     6343        represents ``self``, where the parenthetical properties occur
     6344        when ``orthonormal = True``:
     6345
     6346        - ``A = M*G``
     6347        - The rows of ``G`` are an orthogonal (resp. orthonormal)
     6348          set of vectors.
     6349        - ``G`` times the conjugate-transpose of ``G`` is a diagonal
     6350          (resp. identity) matrix.
     6351        - The row space of ``G`` equals the row space of ``A``.
     6352        - ``M`` is a full-rank matrix with zeros above the diagonal.
     6353
     6354        For exact rings, any zero vectors produced (when the original
     6355        vectors are linearly dependent) are not output, thus the
     6356        orthonormal set is linearly independent, and thus a basis for the
     6357        row space of the original matrix.
     6358
     6359        Any notion of a Gram-Schmidt procedure requires that the base
     6360        ring of the matrix has a fraction field implemented.  In order
     6361        to arrive at an orthonormal set, it must be possible to construct
     6362        square roots of the elements of the base field.  In Sage, your
     6363        best option is the field of algebraic numbers, ``QQbar``, which
     6364        properly contain the rationals and number fields.
     6365
     6366        If you have an approximate numerical matrix, then this routine
     6367        requires that your base field be the real and complex
     6368        double-precision floating point numbers, ``RDF`` and ``CDF``.
     6369        In this case, the matrix is treated as having full rank, as no
     6370        attempt is made to recognize linear dependence with approximate
     6371        calculations.
     6372
    61936373        EXAMPLES::
    6194        
     6374
     6375        Inexact Rings, Numerical Matrices:
     6376
     6377        First, the inexact rings, ``CDF`` and ``RDF``.  ::
     6378
     6379            sage: A = matrix(CDF, [[ 0.6454 + 0.7491*I, -0.8662 + 0.1489*I,  0.7656 - 0.00344*I],
     6380            ...                    [-0.2913 + 0.8057*I,  0.8321 + 0.8170*I, -0.6744 + 0.9248*I],
     6381            ...                    [ 0.2554 + 0.3517*I, -0.4454 - 0.1715*I,  0.8325 - 0.6282*I]])
     6382            sage: G, M = A.gram_schmidt()
     6383            sage: G.round(6)
     6384            [-0.422243 - 0.490087*I  0.566698 - 0.097416*I -0.500882 + 0.002251*I]
     6385            [-0.057002 - 0.495035*I  -0.35059 - 0.625323*I  0.255514 - 0.415284*I]
     6386            [ 0.394105 - 0.421778*I -0.392266 - 0.039345*I  -0.352905 + 0.62195*I]
     6387            sage: M.round(6)
     6388            [             -1.528503                      0                      0]
     6389            [  0.459974 - 0.40061*I              -1.741233                      0]
     6390            [-0.934304 + 0.148868*I   0.54833 + 0.073202*I              -0.550725]
     6391            sage: (A - M*G).zero_at(10^-12)
     6392            [0 0 0]
     6393            [0 0 0]
     6394            [0 0 0]
     6395            sage: (G*G.conjugate().transpose()).zero_at(10^-12)
     6396            [1.0   0   0]
     6397            [  0 1.0   0]
     6398            [  0   0 1.0]
     6399
     6400        A rectangular matrix.  Note that the ``orthonormal`` keyword
     6401        is ignored in these cases.  ::
     6402
     6403            sage: A = matrix(RDF, [[-0.978325, -0.751994, 0.925305, -0.200512, 0.420458],
     6404            ...                    [-0.474877, -0.983403, 0.089836,  0.132218, 0.672965]])
     6405            sage: G, M = A.gram_schmidt(orthonormal=False)
     6406            sage: G.round(6)
     6407            [-0.607223 -0.466745  0.574315 -0.124453  0.260968]
     6408            [ 0.123203 -0.617909 -0.530578  0.289773  0.487368]
     6409            sage: M.round(6)
     6410            [1.611147      0.0]
     6411            [0.958116 0.867778]
     6412            sage: (A-M*G).zero_at(10^-12)
     6413            [0.0 0.0 0.0 0.0 0.0]
     6414            [0.0 0.0 0.0 0.0 0.0]
     6415            sage: (G*G.transpose()).zero_at(10^-12)
     6416            [1.0 0.0]
     6417            [0.0 1.0]
     6418
     6419        Even though a set of vectors may be linearly dependent, no effort
     6420        is made to decide when a zero vector is really the result of a
     6421        relation of linear dependence.  So in this regard, input matrices
     6422        are treated as being of full rank.  Try one of the base rings that
     6423        provide exact results if you need exact results.  ::
     6424
     6425            sage: entries = [[1,1,2], [2,1,3], [3,1,4]]
     6426            sage: A = matrix(QQ, entries)
     6427            sage: A.rank()
     6428            2
     6429            sage: B = matrix(RDF, entries)
     6430            sage: G, M = B.gram_schmidt()
     6431            sage: G.round(6)
     6432            [-0.408248 -0.408248 -0.816497]
     6433            [ 0.707107 -0.707107      -0.0]
     6434            [ -0.57735  -0.57735   0.57735]
     6435            sage: M.round(10)
     6436            [-2.4494897428           0.0           0.0]
     6437            [-3.6742346142  0.7071067812           0.0]
     6438            [-4.8989794856  1.4142135624           0.0]
     6439
     6440        Exact Rings, Orthonormalization:
     6441
     6442        To scale a vector to unit length requires taking
     6443        a square root, which often takes us outside the base ring.
     6444        For the integers,  and rationals the field of algebraic numbers,
     6445        ``QQbar``, is big enough to contain what we need, but the price
     6446        is that the computations are very slow, hence mostly of value
     6447        for small cases or instruction. Now we need to use the
     6448        ``orthonormal`` keyword.  ::
     6449
     6450            sage: A = matrix(QQbar, [[6, -8,  1],
     6451            ...                      [4,  1,  3],
     6452            ...                      [6,  3,  3],
     6453            ...                      [7,  1, -5],
     6454            ...                      [7, -3,  5]])
     6455            sage: G, M = A.gram_schmidt(orthonormal=True)
     6456            sage: G
     6457            [ 0.5970223141259934? -0.7960297521679913? 0.09950371902099891?]
     6458            [ 0.6063218341690895?  0.5289635311888953?  0.5937772444966257?]
     6459            [ 0.5252981913594170?  0.2941669871612735?  -0.798453250866314?]
     6460            sage: M
     6461            [ 10.04987562112089?                   0                   0]
     6462            [ 1.890570661398980?  4.735582601355131?                   0]
     6463            [ 1.492555785314984?  7.006153332071100?  1.638930357041381?]
     6464            [ 2.885607851608969?  1.804330147889395?  7.963520581008761?]
     6465            [ 7.064764050490923?  5.626248468100069? -1.197679876299471?]
     6466            sage: M*G-A
     6467            [0 0 0]
     6468            [0 0 0]
     6469            [0 0 0]
     6470            [0 0 0]
     6471            [0 0 0]
     6472            sage: G*G.transpose()
     6473            [      1 0.?e-37 0.?e-35]
     6474            [0.?e-37       1 0.?e-35]
     6475            [0.?e-35 0.?e-35       1]
     6476            sage: G.row_space() == A.row_space()
     6477            True
     6478
     6479        Starting with complex numbers with rational real and imaginary parts.
     6480        Note the use of the conjugate-transpose when checking the
     6481        orthonormality. ::
     6482
     6483            sage: A = matrix(QQbar, [[  -2,    -I - 1, 4*I + 2,       -1],
     6484            ...                      [-4*I, -2*I + 17,       0,  9*I + 1],
     6485            ...                      [   1,  -2*I - 6, -I + 11, -5*I + 1]])
     6486            sage: G, M = A.gram_schmidt(orthonormal=True)
     6487            sage: G
     6488            [                         -0.3849001794597505?  -0.1924500897298753? - 0.1924500897298753?*I   0.3849001794597505? + 0.7698003589195010?*I                          -0.1924500897298753?]
     6489            [-0.06165497274852388? - 0.1387236886841787?*I  0.8188551068163327? - 0.10018933071635130?*I  0.2003786614327026? + 0.05394810115495839?*I  0.02119389688230508? + 0.5028733714801478?*I]
     6490            [  0.3842387256410419? - 0.5694103019142261?*I  0.1416892863096208? - 0.06139779741542298?*I  0.4633778333528464? - 0.01285039016180503?*I  0.02658516588219101? - 0.5373044261814995?*I]
     6491            sage: M
     6492            [                        5.196152422706632?                                          0                                          0]
     6493            [-3.079201435678004? + 3.464101615137755?*I             19.22286447225071? + 0.?e-37*I                                          0]
     6494            [  4.426352063787131? - 8.66025403784439?*I -5.117362738127481? - 3.502773139275513?*I             7.480012456446966? + 0.?e-35*I]
     6495            sage: M*G-A
     6496            [            0.?e-37 0.?e-37 + 0.?e-37*I 0.?e-37 + 0.?e-37*I             0.?e-37]
     6497            [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]
     6498            [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]
     6499            sage: G*G.conjugate().transpose()
     6500            [1.000000000000000? + 0.?e-37*I            0.?e-37 + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
     6501            [           0.?e-37 + 0.?e-37*I 1.000000000000000? + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
     6502            [           0.?e-36 + 0.?e-36*I            0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-35*I]
     6503            sage: G.row_space() == A.row_space()
     6504            True
     6505
     6506        A square matrix with small rank.  The zero vectors produced as a
     6507        result of linear dependence get eliminated, so the rows of ``G``
     6508        are a basis for the row space of ``A``.  ::
     6509
     6510            sage: A = matrix(QQbar, [[2, -6, 3, 8],
     6511            ...                      [1, -3, 2, 5],
     6512            ...                      [0,  0, 2, 4],
     6513            ...                      [2, -6, 3, 8]])
     6514            sage: A.change_ring(QQ).rank()
     6515            2
     6516            sage: G, M = A.gram_schmidt(orthonormal=True)
     6517            sage: G
     6518            [ 0.1881441736767195? -0.5644325210301583?  0.2822162605150792?  0.7525766947068779?]
     6519            [-0.2502818123591464?   0.750845437077439?  0.3688363550555841?  0.4873908977520218?]
     6520            sage: M
     6521            [10.630145812734649?                   0]
     6522            [ 6.208757731331742? 0.6718090752798139?]
     6523            [ 3.574739299857670?  2.687236301119256?]
     6524            [10.630145812734649?                   0]
     6525            sage: M*G-A
     6526            [0 0 0 0]
     6527            [0 0 0 0]
     6528            [0 0 0 0]
     6529            [0 0 0 0]
     6530            sage: G*G.transpose()
     6531            [      1 0.?e-35]
     6532            [0.?e-35       1]
     6533            sage: G.row_space() == A.row_space()
     6534            True
     6535
     6536        Exact Rings, Orthogonalization:
     6537
     6538        If we forego scaling orthogonal vectors to unit vectors, we
     6539        can apply Gram-Schmidt to a much greater variety of rings.
     6540        Use the ``orthonormal=False`` keyword (or assume it as the default).
     6541        Note that now the orthogonality check creates a diagonal matrix
     6542        whose diagonal entries are the squares of the lengths of the
     6543        vectors. ::
     6544
     6545        First, in the rationals, without involving ``QQbar``.  ::
     6546
     6547            sage: A = matrix(QQ, [[-1,  3,  2,  2],
     6548            ...                   [-1,  0, -1,  0],
     6549            ...                   [-1, -2, -3, -1],
     6550            ...                   [ 1,  1,  2,  0]])
     6551            sage: A.rank()
     6552            3
     6553            sage: G, M = A.gram_schmidt()
     6554            sage: G
     6555            [    -1      3      2      2]
     6556            [-19/18    1/6   -8/9    1/9]
     6557            [  2/35  -4/35  -2/35   9/35]
     6558            sage: M
     6559            [     1      0      0]
     6560            [ -1/18      1      0]
     6561            [-13/18  59/35      1]
     6562            [   1/3 -48/35     -2]
     6563            sage: M*G-A
     6564            [0 0 0 0]
     6565            [0 0 0 0]
     6566            [0 0 0 0]
     6567            [0 0 0 0]
     6568            sage: G*G.transpose()
     6569            [   18     0     0]
     6570            [    0 35/18     0]
     6571            [    0     0  3/35]
     6572            sage: G.row_space() == A.row_space()
     6573            True
     6574
     6575        A complex subfield of the complex numbers.  ::
     6576
     6577            sage: C.<z> = CyclotomicField(5)
     6578            sage: A = matrix(C, [[              -z^3 - 2*z,             -z^3 - 1, 2*z^3 - 2*z^2 + 2*z,             1],
     6579            ...                  [         z^3 - 2*z^2 + 1, -z^3 + 2*z^2 - z - 1,                  -1,       z^2 + z],
     6580            ...                  [-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]])
     6581            sage: G, M = A.gram_schmidt(orthonormal=False)
     6582            sage: G
     6583            [                                                      -z^3 - 2*z                                                         -z^3 - 1                                              2*z^3 - 2*z^2 + 2*z                                                                1]
     6584            [                   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]
     6585            [-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]
     6586            sage: M
     6587            [                                                           1                                                            0                                                            0]
     6588            [                14/139*z^3 + 47/139*z^2 + 145/139*z + 95/139                                                            1                                                            0]
     6589            [              -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]
     6590            sage: M*G - A
     6591            [0 0 0 0]
     6592            [0 0 0 0]
     6593            [0 0 0 0]
     6594            sage: G*G.conjugate().transpose()
     6595            [                               15*z^3 + 15*z^2 + 28                                                   0                                                   0]
     6596            [                                                  0                463/139*z^3 + 463/139*z^2 + 1971/139                                                   0]
     6597            [                                                  0                                                   0 230983/19841*z^3 + 230983/19841*z^2 + 1003433/39682]
     6598            sage: G.row_space() == A.row_space()
     6599            True
     6600
     6601        A slightly edited legacy example.  ::
     6602
    61956603            sage: A = matrix(ZZ, 3, [-1, 2, 5, -11, 1, 1, 1, -1, -3]); A
    61966604            [ -1   2   5]
    61976605            [-11   1   1]
     
    62026610            [  -52/5    -1/5      -2]
    62036611            [  2/187  36/187 -14/187]
    62046612            sage: mu
    6205             [     0      0      0]
    6206             [   3/5      0      0]
    6207             [  -3/5 -7/187      0]
     6613            [     1      0      0]
     6614            [   3/5      1      0]
     6615            [  -3/5 -7/187      1]
    62086616            sage: G.row(0) * G.row(1)
    62096617            0
    62106618            sage: G.row(0) * G.row(2)
    62116619            0
    62126620            sage: G.row(1) * G.row(2)
    62136621            0
    6214        
    6215         The relation between mu and A is as follows::
    6216        
    6217             sage: (mu + 1)*G == A
     6622
     6623        The relation between mu and A is as follows::
     6624
     6625            sage: mu*G == A
    62186626            True
    6219         """       
    6220         from sage.modules.misc import gram_schmidt
    6221         from constructor import matrix
    6222         Bstar, mu = gram_schmidt(self.rows())
    6223         return matrix(Bstar), mu
     6627        """
     6628        import sage.rings.real_double
     6629        import sage.rings.complex_double
     6630        R = self.base_ring()
     6631        if R in [sage.rings.real_double.RDF, sage.rings.complex_double.CDF]:
     6632            Q, R = self.transpose().QR()
     6633            m = R.nrows(); n = R.ncols()
     6634            if m > n:
     6635                Q = Q[0:m, 0:n]
     6636                R = R[0:n, 0:n]
     6637        elif R.is_exact():
     6638            if orthonormal:
     6639                Q, R = self.transpose().QR(full=False)
     6640            else:
     6641                Q, R = self.transpose()._gram_schmidt_noscale()
     6642        else:
     6643            raise NotImplementedError("Gram-Scmidt orthogonalization not implemented for matrices over inexact rings, except for RDF and CDF")
     6644        return Q.transpose(), R.transpose()
    62246645
    62256646    def jordan_form(self, base_ring=None, sparse=False, subdivide=True, transformation=False):
    62266647        r"""
  • sage/matrix/matrix_integer_dense.pyx

    diff -r 210ad87f9cd6 -r 0b859949124e sage/matrix/matrix_integer_dense.pyx
    a b  
    30183018            raise TypeError("eta must be >= 0.5")
    30193019
    30203020        # this is pretty slow
    3021         G, mu = self.gram_schmidt()
    3022  
     3021        import sage.modules.misc
     3022        G, mu = sage.modules.misc.gram_schmidt(self.rows())
    30233023        #For any $i>j$, we have $|mu_{i, j}| <= \eta$
    30243024        for e in mu.list():
    30253025            if e.abs() > eta:
    30263026                return False
    30273027
    30283028        #For any $i<d$, we have $\delta |b_i^*|^2 <= |b_{i+1}^* + mu_{i+1, i} b_i^* |^2$
    3029         norms = [G[i].norm()**2 for i in range(G.nrows())]
     3029        norms = [G[i].norm()**2 for i in range(len(G))]
    30303030        for i in xrange(1,self.nrows()):
    30313031            if norms[i] < (delta - mu[i,i-1]**2) * norms[i-1]:
    30323032                return False
  • sage/modules/misc.py

    diff -r 210ad87f9cd6 -r 0b859949124e 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:
    2739        sage: B = [vector([1,2,1/5]), vector([1,2,3]), vector([-1,0,0])]
    2840        sage: from sage.modules.misc import gram_schmidt
     
    4557        sage: a = matrix([[],[],[],[]])
    4658        sage: a.gram_schmidt()
    4759         ([], [])
     60
     61    Linearly dependent input leads to a zero dot product in a denominator.
     62    This shows that Trac #10791 is fixed. ::
     63
     64        sage: from sage.modules.misc import gram_schmidt
     65        sage: V = [vector(ZZ,[1,1]), vector(ZZ,[2,2]), vector(ZZ,[1,2])]
     66        sage: gram_schmidt(V)
    4867    """
     68    import sage.modules.free_module_element
    4969    if len(B) == 0 or len(B[0]) == 0:
    5070        return B, matrix(ZZ,0,0,[])
    5171    n = len(B)
    5272    Bstar = [B[0]]
    5373    K = B[0].base_ring().fraction_field()
     74    zero = sage.modules.free_module_element.vector(K, len(B[0]))
     75    if Bstar[0] == zero:
     76        raise ValueError("linearly dependent input for module version of Gram-Schmidt")
    5477    mu = matrix(K, n, n)
    5578    for i in range(1,n):
    5679        for j in range(i):
    5780            mu[i,j] = B[i].dot_product(Bstar[j]) / (Bstar[j].dot_product(Bstar[j]))
    5881        Bstar.append(B[i] - sum(mu[i,j]*Bstar[j] for j in range(i)))
     82        if Bstar[i] == zero:
     83            raise ValueError("linearly dependent input for module version of Gram-Schmidt")
    5984    return Bstar, mu
    6085