Ticket #7852: trac_7852-solve-linear-systems-CDF.patch

File trac_7852-solve-linear-systems-CDF.patch, 14.7 KB (added by rbeezer, 10 years ago)
  • sage/matrix/matrix_complex_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1298586612 28800
    # Node ID d5c6230ee8c31e3bdc8dff95666430c537f988c4
    # Parent  3c67675a1df921b9a9fa0ef007fc54c1a1d1e8f6
    7852: solve_left, solve_right over RDF/CDF
    
    diff -r 3c67675a1df9 -r d5c6230ee8c3 sage/matrix/matrix_complex_double_dense.pyx
    a b  
    7979    ::
    8080   
    8181        sage: b = vector(CDF,[5,6])
    82         sage: m.solve_left(b)
     82        sage: m.solve_right(b)
    8383        (2.66666666667 + 0.666666666667*I, -0.333333333333 - 1.16666666667*I)
    8484   
    8585    See the commands qr, lu, and svd for QR, LU, and singular value
  • sage/matrix/matrix_double_dense.pyx

    diff -r 3c67675a1df9 -r d5c6230ee8c3 sage/matrix/matrix_double_dense.pyx
    a b  
    954954        M._matrix_numpy = scipy.linalg.lu_solve((lu, self._P_M), b)
    955955        return M
    956956
    957     def solve_left(self,vec):
     957    def solve_right(self, b):
     958        r"""
     959        Solve the vector equation ``A*x = b`` for a nonsingular ``A``.
     960
     961        INPUT:
     962
     963        - ``self`` - a square matrix that is nonsigular (of full rank).
     964        - ``b`` - a vector of the correct size.  Elements of the vector
     965          must coerce into the base ring of the coefficient matrix.  In
     966          particular, if ``b`` has entries from ``CDF`` then ``self``
     967          must have ``CDF`` as its base ring.
     968
     969        OUTPUT:
     970
     971        The unique solution ``x`` to the matrix equation ``A*x = b``,
     972        as a vector over the same base ring as ``self``.
     973
     974        ALGORITHM:
     975
     976        Uses the ``solve()`` routine from the SciPy ``scipy.linalg`` module.
     977
     978        EXAMPLES:
     979
     980        Over the reals. ::
     981
     982            sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A
     983            [ 1.0  2.0  5.0]
     984            [ 7.6  2.3  1.0]
     985            [ 1.0  2.0 -1.0]
     986            sage: b = vector(RDF,[1,2,3])
     987            sage: x = A.solve_right(b); x
     988            (-0.113695090439, 1.39018087855, -0.333333333333)
     989            sage: x.parent()
     990            Vector space of dimension 3 over Real Double Field
     991            sage: A*x
     992            (1.0, 2.0, 3.0)
     993
     994        Over the complex numbers.  ::
     995
     996            sage: A = matrix(CDF, [[      0, -1 + 2*I,  1 - 3*I,        I],
     997            ...                    [2 + 4*I, -2 + 3*I, -1 + 2*I,   -1 - I],
     998            ...                    [  2 + I,    1 - I,       -1,        5],
     999            ...                    [    3*I,   -1 - I,   -1 + I,   -3 + I]])
     1000            sage: b = vector(CDF, [2 -3*I, 3, -2 + 3*I, 8])
     1001            sage: x = A.solve_right(b); x
     1002            (1.96841637... - 1.07606761...*I, -0.614323843... + 1.68416370...*I, 0.0733985765... + 1.73487544...*I, -1.6018683... + 0.524021352...*I)
     1003            sage: x.parent()
     1004            Vector space of dimension 4 over Complex Double Field
     1005            sage: abs(A*x - b) < 1e-14
     1006            True
     1007
     1008        The vector of constants, ``b``, can be given in a
     1009        variety of forms, so long as it coerces to a vector
     1010        over the same base ring as the coefficient matrix.  ::
     1011
     1012            sage: A=matrix(CDF, 5, [1/(i+j+1) for i in range(5) for j in range(5)])
     1013            sage: A.solve_right([1]*5)
     1014            (5.0..., -120.0..., 630.0..., -1120.0..., 630.0...)
     1015
     1016        TESTS:
     1017
     1018        A degenerate case. ::
     1019
     1020            sage: A = matrix(RDF, 0, 0, [])
     1021            sage: A.solve_right(vector(RDF,[]))
     1022            ()
     1023
     1024        The coefficent matrix must be square. ::
     1025
     1026            sage: A = matrix(RDF, 2, 3, range(6))
     1027            sage: b = vector(RDF, [1,2,3])
     1028            sage: A.solve_right(b)
     1029            Traceback (most recent call last):
     1030            ...
     1031            ValueError: coefficient matrix of a system over RDF/CDF must be square, not 2 x 3
     1032
     1033        The coefficient matrix must be nonsingular.  ::
     1034
     1035            sage: A = matrix(RDF, 5, range(25))
     1036            sage: b = vector(RDF, [1,2,3,4,5])
     1037            sage: A.solve_right(b)
     1038            Traceback (most recent call last):
     1039            ...
     1040            LinAlgError: singular matrix
     1041
     1042        The vector of constants needs the correct degree.  ::
     1043
     1044            sage: A = matrix(RDF, 5, range(25))
     1045            sage: b = vector(RDF, [1,2,3,4])
     1046            sage: A.solve_right(b)
     1047            Traceback (most recent call last):
     1048            ...
     1049            TypeError: vector of constants over Real Double Field incompatible with matrix over Real Double Field
     1050
     1051        The vector of constants needs to be compatible with
     1052        the base ring of the coefficient matrix.  ::
     1053
     1054            sage: F.<a> = FiniteField(27)
     1055            sage: b = vector(F, [a,a,a,a,a])
     1056            sage: A.solve_right(b)
     1057            Traceback (most recent call last):
     1058            ...
     1059            TypeError: vector of constants over Finite Field in a of size 3^3 incompatible with matrix over Real Double Field
     1060
     1061        With a coefficient matrix over ``RDF``, a vector of constants
     1062        over ``CDF`` can be accomodated by converting the base ring
     1063        of the coefficient matrix.  ::
     1064
     1065            sage: A = matrix(RDF, 2, range(4))
     1066            sage: b = vector(CDF, [1+I,2])
     1067            sage: A.solve_right(b)
     1068            Traceback (most recent call last):
     1069            ...
     1070            TypeError: vector of constants over Complex Double Field incompatible with matrix over Real Double Field
     1071
     1072            sage: B = A.change_ring(CDF)
     1073            sage: B.solve_right(b)
     1074            (-0.5 - 1.5*I, 1.0 + 1.0*I)
    9581075        """
    959         Solve the equation A*x = b, where
    960        
     1076        if not self.is_square():
     1077            raise ValueError("coefficient matrix of a system over RDF/CDF must be square, not %s x %s " % (self.nrows(), self.ncols()))
     1078        M = self._column_ambient_module()
     1079        try:
     1080            vec = M(b)
     1081        except TypeError:
     1082            raise TypeError("vector of constants over %s incompatible with matrix over %s" % (b.base_ring(), self.base_ring()))
     1083        if vec.degree() != self.ncols():
     1084            raise ValueError("vector of constants in linear system over RDF/CDF must have degree equal to the number of columns for the coefficient matrix, not %s" % vec.degree() )
     1085
     1086        if self._ncols == 0:
     1087            return M.zero_vector()
     1088
     1089        global scipy
     1090        if scipy is None:
     1091            import scipy
     1092        import scipy.linalg
     1093        # may raise a LinAlgError for a singular matrix
     1094        return M(scipy.linalg.solve(self._matrix_numpy, vec.numpy()))
     1095
     1096    def solve_left(self, b):
     1097        r"""
     1098        Solve the vector equation ``x*A = b`` for a nonsingular ``A``.
     1099
     1100        INPUT:
     1101
     1102        - ``self`` - a square matrix that is nonsigular (of full rank).
     1103        - ``b`` - a vector of the correct size.  Elements of the vector
     1104          must coerce into the base ring of the coefficient matrix.  In
     1105          particular, if ``b`` has entries from ``CDF`` then ``self``
     1106          must have ``CDF`` as its base ring.
     1107
     1108        OUTPUT:
     1109
     1110        The unique solution ``x`` to the matrix equation ``x*A = b``,
     1111        as a vector over the same base ring as ``self``.
     1112
     1113        ALGORITHM:
     1114
     1115        Uses the ``solve()`` routine from the SciPy ``scipy.linalg`` module,
     1116        after taking the tranpose of the coefficient matrix.
     1117
    9611118        EXAMPLES:
     1119
     1120        Over the reals. ::
     1121
    9621122            sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A
    9631123            [ 1.0  2.0  5.0]
    9641124            [ 7.6  2.3  1.0]
    9651125            [ 1.0  2.0 -1.0]
    9661126            sage: b = vector(RDF,[1,2,3])
    9671127            sage: x = A.solve_left(b); x
    968             (-0.113695090439, 1.39018087855, -0.333333333333)
    969             sage: A*x
     1128            (0.666666666..., 0.0, 0.333333333...)
     1129            sage: x.parent()
     1130            Vector space of dimension 3 over Real Double Field
     1131            sage: x*A
    9701132            (1.0, 2.0, 3.0)
    9711133
     1134        Over the complex numbers.  ::
     1135
     1136            sage: A = matrix(CDF, [[      0, -1 + 2*I,  1 - 3*I,        I],
     1137            ...                    [2 + 4*I, -2 + 3*I, -1 + 2*I,   -1 - I],
     1138            ...                    [  2 + I,    1 - I,       -1,        5],
     1139            ...                    [    3*I,   -1 - I,   -1 + I,   -3 + I]])
     1140            sage: b = vector(CDF, [2 -3*I, 3, -2 + 3*I, 8])
     1141            sage: x = A.solve_left(b); x
     1142            (-1.55765124... - 0.644483985...*I, 0.183274021... + 0.286476868...*I, 0.270818505... + 0.246619217...*I, -1.69003558... - 0.828113879...*I)
     1143            sage: x.parent()
     1144            Vector space of dimension 4 over Complex Double Field
     1145            sage: abs(x*A - b) < 1e-14
     1146            True
     1147
     1148        The vector of constants, ``b``, can be given in a
     1149        variety of forms, so long as it coerces to a vector
     1150        over the same base ring as the coefficient matrix.  ::
     1151
     1152            sage: A=matrix(CDF, 5, [1/(i+j+1) for i in range(5) for j in range(5)])
     1153            sage: A.solve_left([1]*5)
     1154            (5.0..., -120.0..., 630.0..., -1120.0..., 630.0...)
     1155
    9721156        TESTS:
    973         We test two degenerate cases:
    974             sage: A = matrix(RDF, 0, 3, [])
     1157
     1158        A degenerate case. ::
     1159
     1160            sage: A = matrix(RDF, 0, 0, [])
    9751161            sage: A.solve_left(vector(RDF,[]))
    9761162            ()
    977             sage: A = matrix(RDF, 3, 0, [])
    978             sage: A.solve_left(vector(RDF,3, [1,2,3]))
    979             (0.0, 0.0, 0.0)
     1163
     1164        The coefficent matrix must be square. ::
     1165
     1166            sage: A = matrix(RDF, 2, 3, range(6))
     1167            sage: b = vector(RDF, [1,2,3])
     1168            sage: A.solve_left(b)
     1169            Traceback (most recent call last):
     1170            ...
     1171            ValueError: coefficient matrix of a system over RDF/CDF must be square, not 2 x 3
     1172
     1173        The coefficient matrix must be nonsingular.  ::
     1174
     1175            sage: A = matrix(RDF, 5, range(25))
     1176            sage: b = vector(RDF, [1,2,3,4,5])
     1177            sage: A.solve_left(b)
     1178            Traceback (most recent call last):
     1179            ...
     1180            LinAlgError: singular matrix
     1181
     1182        The vector of constants needs the correct degree.  ::
     1183
     1184            sage: A = matrix(RDF, 5, range(25))
     1185            sage: b = vector(RDF, [1,2,3,4])
     1186            sage: A.solve_left(b)
     1187            Traceback (most recent call last):
     1188            ...
     1189            TypeError: vector of constants over Real Double Field incompatible with matrix over Real Double Field
     1190
     1191        The vector of constants needs to be compatible with
     1192        the base ring of the coefficient matrix.  ::
     1193
     1194            sage: F.<a> = FiniteField(27)
     1195            sage: b = vector(F, [a,a,a,a,a])
     1196            sage: A.solve_left(b)
     1197            Traceback (most recent call last):
     1198            ...
     1199            TypeError: vector of constants over Finite Field in a of size 3^3 incompatible with matrix over Real Double Field
     1200
     1201        With a coefficient matrix over ``RDF``, a vector of constants
     1202        over ``CDF`` can be accomodated by converting the base ring
     1203        of the coefficient matrix.  ::
     1204
     1205            sage: A = matrix(RDF, 2, range(4))
     1206            sage: b = vector(CDF, [1+I,2])
     1207            sage: A.solve_left(b)
     1208            Traceback (most recent call last):
     1209            ...
     1210            TypeError: vector of constants over Complex Double Field incompatible with matrix over Real Double Field
     1211
     1212            sage: B = A.change_ring(CDF)
     1213            sage: B.solve_left(b)
     1214            (0.5 - 1.5*I, 0.5 + 0.5*I)
    9801215        """
    981         M = self._column_ambient_module()
    982         if self._nrows == 0 or self._ncols == 0:
     1216        if not self.is_square():
     1217            raise ValueError("coefficient matrix of a system over RDF/CDF must be square, not %s x %s " % (self.nrows(), self.ncols()))
     1218        M = self._row_ambient_module()
     1219        try:
     1220            vec = M(b)
     1221        except TypeError:
     1222            raise TypeError("vector of constants over %s incompatible with matrix over %s" % (b.base_ring(), self.base_ring()))
     1223        if vec.degree() != self.nrows():
     1224            raise ValueError("vector of constants in linear system over RDF/CDF must have degree equal to the number of rows for the coefficient matrix, not %s" % vec.degree() )
     1225
     1226        if self._nrows == 0:
    9831227            return M.zero_vector()
     1228
    9841229        global scipy
    9851230        if scipy is None:
    9861231            import scipy
    9871232        import scipy.linalg
    988         return M(scipy.linalg.solve(self._matrix_numpy,vec.numpy()))
    989        
     1233        # may raise a LinAlgError for a singular matrix
     1234        # call "right solve" routine with the transpose
     1235        return M(scipy.linalg.solve(self._matrix_numpy.T, vec.numpy()))
     1236
    9901237    def determinant(self):
    9911238        """
    9921239        Return the determinant of self.
  • sage/matrix/matrix_real_double_dense.pyx

    diff -r 3c67675a1df9 -r d5c6230ee8c3 sage/matrix/matrix_real_double_dense.pyx
    a b  
    7979    ::
    8080   
    8181        sage: b = vector(RDF,[5,6])
    82         sage: m.solve_left(b)
     82        sage: m.solve_right(b)
    8383        (-4.0, 4.5)
    8484   
    8585    See the commands qr, lu, and svd for QR, LU, and singular value
  • sage/rings/polynomial/polynomial_element.pyx

    diff -r 3c67675a1df9 -r d5c6230ee8c3 sage/rings/polynomial/polynomial_element.pyx
    a b  
    10021002            sage: t.inverse_mod(S.ideal((t + 1)^7)) == f
    10031003            True
    10041004       
    1005         It also works over in-exact rings, but note that due to rounding
    1006         error the product is only guaranteed to be within epsilon of the
    1007         constant polynomial 1.
    1008        
    1009         ::
     1005        This also works over inexact rings, but note that due to rounding
     1006        error the product may not always exactly equal the constant
     1007        polynomial 1 and have extra terms with coefficients close to zero. ::
    10101008       
    10111009            sage: R.<x> = RDF[]
    10121010            sage: f = inverse_mod(x^2 + 1, x^5 + x + 1); f
    10131011            0.4*x^4 - 0.2*x^3 - 0.4*x^2 + 0.2*x + 0.8
    10141012            sage: f * (x^2 + 1) % (x^5 + x + 1)
    1015             2.22044604925e-16*x^2 + 1.11022302463e-16*x + 1.0
     1013            1.0
    10161014            sage: f = inverse_mod(x^3 - x + 1, x - 2); f
    10171015            0.142857142857
    10181016            sage: f * (x^3 - x + 1) % (x - 2)
    10191017            1.0
    1020        
     1018            sage: g = 5*x^3+x-7; m = x^4-12*x+13; f = inverse_mod(g, m); f
     1019            -0.0319636125...*x^3 - 0.0383269759...*x^2 - 0.0463050900...*x + 0.346479687...
     1020            sage: f*g % m
     1021            -8.881784...e-16*x^3 + 8.881784...e-16*x^2 - 8.881784...e-16*x + 1.0
     1022
    10211023        ALGORITHM: Solve the system as + mt = 1, returning s as the inverse
    10221024        of a mod m.
    10231025