Ticket #13140: trac_13140_double_dense.patch

File trac_13140_double_dense.patch, 7.5 KB (added by jhpalmieri, 7 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User J. H. Palmieri <palmieri@math.washington.edu>
    # Date 1340730474 25200
    # Node ID 32efde7127c714da584ce5566ba47297e823b3f7
    # Parent  ef777986e8a5aa5ff04d778cf96746fbe8989aae
    Normalize doctests for QR decomposition for double dense matrices
    
    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b cdef class Matrix_double_dense(matrix_de 
    24742474            ...                    [3, 0, 8, 3, 3],
    24752475            ...                    [-1, 1, -6, -6, 5]])
    24762476            sage: Q, R = A.QR()
    2477             sage: Q.round(6).zero_at(10^-6)
    2478             [-0.458831  0.126051 -0.381212 -0.394574  -0.68744]
    2479             [-0.458831  -0.47269  0.051983  0.717294 -0.220963]
    2480             [ 0.229416 -0.661766 -0.661923 -0.180872  0.196411]
    2481             [ 0.688247 -0.189076  0.204468   0.09663 -0.662889]
    2482             [-0.229416 -0.535715  0.609939 -0.536422  0.024551]
    2483             sage: R.round(6).zero_at(10^-6)
     2477
     2478        At this point, ``Q`` is only well-defined up to the
     2479        signs of its columns, and similarly for ``R`` and its
     2480        rows, so we normalize them::
     2481
     2482            sage: Qnorm = Q._normalize_columns()
     2483            sage: Rnorm = R._normalize_rows()
     2484            sage: Qnorm.round(6).zero_at(10^-6)
     2485            [ 0.458831  0.126051  0.381212  0.394574   0.68744]
     2486            [ 0.458831  -0.47269 -0.051983 -0.717294  0.220963]
     2487            [-0.229416 -0.661766  0.661923  0.180872 -0.196411]
     2488            [-0.688247 -0.189076 -0.204468  -0.09663  0.662889]
     2489            [ 0.229416 -0.535715 -0.609939  0.536422 -0.024551]
     2490            sage: Rnorm.round(6).zero_at(10^-6)
    24842491            [ 4.358899 -0.458831 13.076697  6.194225  2.982405]
    2485             [      0.0 -1.670172 -0.598741   1.29202 -6.207997]
    2486             [      0.0       0.0 -5.444402 -5.468661  0.682716]
     2492            [      0.0  1.670172  0.598741  -1.29202  6.207997]
     2493            [      0.0       0.0  5.444402  5.468661 -0.682716]
    24872494            [      0.0       0.0       0.0  1.027626   -3.6193]
    24882495            [      0.0       0.0       0.0       0.0  0.024551]
    24892496            sage: (Q*Q.transpose()).zero_at(10^-14)
    cdef class Matrix_double_dense(matrix_de 
    25082515            ...                    [I + 7, 2*I + 1, -2*I + 7, -I + 1],
    25092516            ...                    [I + 2, 0, I + 12, -1]])
    25102517            sage: Q, R = A.QR()
    2511             sage: Q.round(6).zero_at(10^-6)
    2512             [             -0.730297  0.207057 + 0.538347*I -0.246305 + 0.076446*I   0.238162 - 0.10366*I]
    2513             [              0.091287 -0.207057 - 0.377878*I -0.378656 + 0.195222*I  0.701244 - 0.364371*I]
    2514             [  0.63901 + 0.091287*I  0.170822 + 0.667758*I  0.034115 - 0.040902*I  0.314017 - 0.082519*I]
    2515             [ 0.182574 + 0.091287*I  -0.036235 + 0.07247*I -0.863228 - 0.063228*I -0.449969 - 0.011612*I]
    2516             sage: R.round(6).zero_at(10^-6)
     2518            sage: Q._normalize_columns().round(6).zero_at(10^-6)
     2519            [              0.730297  0.207057 + 0.538347*I  0.246305 - 0.076446*I   0.238162 - 0.10366*I]
     2520            [             -0.091287 -0.207057 - 0.377878*I  0.378656 - 0.195222*I  0.701244 - 0.364371*I]
     2521            [ -0.63901 - 0.091287*I  0.170822 + 0.667758*I -0.034115 + 0.040902*I  0.314017 - 0.082519*I]
     2522            [-0.182574 - 0.091287*I  -0.036235 + 0.07247*I  0.863228 + 0.063228*I -0.449969 - 0.011612*I]
     2523            sage: R._normalize_rows().round(6).zero_at(10^-6)
    25172524            [             10.954451            -1.917029*I   5.385938 - 2.19089*I  -0.273861 - 2.19089*I]
    25182525            [                   0.0               4.829596 -0.869638 - 5.864879*I  0.993872 - 0.305409*I]
    2519             [                   0.0                    0.0             -12.001608  0.270953 - 0.442063*I]
     2526            [                   0.0                    0.0              12.001608 -0.270953 + 0.442063*I]
    25202527            [                   0.0                    0.0                    0.0               1.942964]
    25212528            sage: (Q.conjugate().transpose()*Q).zero_at(10^-15)
    25222529            [1.0 0.0 0.0 0.0]
    cdef class Matrix_double_dense(matrix_de 
    25452552            2
    25462553            sage: A = Arat.change_ring(CDF)
    25472554            sage: Q, R = A.QR()
    2548             sage: R.round(6).zero_at(10^-6)
    2549             [-5.567764   2.69408  -2.69408]
    2550             [      0.0 -3.569585  3.569585]
     2555            sage: R._normalize_rows().round(6).zero_at(10^-6)
     2556            [ 5.567764  -2.69408   2.69408]
     2557            [      0.0  3.569585 -3.569585]
    25512558            [      0.0       0.0       0.0]
    25522559            [      0.0       0.0       0.0]
    25532560            sage: (Q.conjugate_transpose()*Q).zero_at(10^-14)
    cdef class Matrix_double_dense(matrix_de 
    39913998        M = self._new()
    39923999        M._matrix_numpy = numpy.round(self._matrix_numpy, ndigits)
    39934000        return M
     4001
     4002    def _normalize_columns(self):
     4003        """
     4004        Returns a copy of the matrix where each column has been
     4005        multiplied by plus or minus 1, to guarantee that the real
     4006        part of the leading entry of each nonzero column is positive.
     4007
     4008        This is useful for modifying output from algorithms which
     4009        produce matrices which are only well-defined up to signs of
     4010        the columns, for example an algorithm which should produce an
     4011        orthogonal matrix.
     4012
     4013        OUTPUT:
     4014
     4015        A modified copy of the matrix.
     4016
     4017        EXAMPLES::
     4018
     4019            sage: a=matrix(CDF, [[1, -2+I, 0, -3*I], [2, 2, -2, 2], [-3, -3, -3, -2]])
     4020            sage: a
     4021            [         1.0 -2.0 + 1.0*I          0.0       -3.0*I]
     4022            [         2.0          2.0         -2.0          2.0]
     4023            [        -3.0         -3.0         -3.0         -2.0]
     4024            sage: a._normalize_columns()
     4025            [        1.0 2.0 - 1.0*I        -0.0      -3.0*I]
     4026            [        2.0        -2.0         2.0         2.0]
     4027            [       -3.0         3.0         3.0        -2.0]
     4028        """
     4029        M = self.__copy__()
     4030        cdef Py_ssize_t i, j
     4031        for j from 0 <= j < M.ncols():
     4032            for i from 0 <= i < M.column(j).degree():
     4033                a = M.column(j)[i].real()
     4034                if a != 0:
     4035                    if a < 0:
     4036                        M = M.with_rescaled_col(j, -1)
     4037                    break
     4038        return M
     4039
     4040    def _normalize_rows(self):
     4041        """
     4042        Returns a copy of the matrix where each row has been
     4043        multiplied by plus or minus 1, to guarantee that the real
     4044        part of the leading entry of each nonzero row is positive.
     4045
     4046        This is useful for modifying output from algorithms which
     4047        produce matrices which are only well-defined up to signs of
     4048        the rows, for example an algorithm which should produce an
     4049        upper triangular matrix.
     4050
     4051        OUTPUT:
     4052
     4053        A modified copy of the matrix.
     4054
     4055        EXAMPLES::
     4056
     4057            sage: a=matrix(CDF, [[1, 2, -3], [-2+I, 2, -3], [0, -2, -3], [-3*I, 2, -2]])
     4058            sage: a
     4059            [         1.0          2.0         -3.0]
     4060            [-2.0 + 1.0*I          2.0         -3.0]
     4061            [         0.0         -2.0         -3.0]
     4062            [      -3.0*I          2.0         -2.0]
     4063            sage: a._normalize_rows()
     4064            [        1.0         2.0        -3.0]
     4065            [2.0 - 1.0*I        -2.0         3.0]
     4066            [       -0.0         2.0         3.0]
     4067            [     -3.0*I         2.0        -2.0]           
     4068        """
     4069        return self.transpose()._normalize_columns().transpose()