diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index c0f1d646340..1a92317115f 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -221,8 +221,8 @@ cdef class Matrix(Matrix1): return matrix([a.subs(*args, **kwds) for a in self.list()], nrows=self._nrows, ncols=self._ncols, sparse=False) - def solve_left(self, B, check=True): - """ + def solve_left(self, B, check=True, *, extend=True): + r""" Try to find a solution `X` to the equation `X A = B`. If ``self`` is a matrix `A`, then this function returns a @@ -235,16 +235,17 @@ cdef class Matrix(Matrix1): field, this method computes a least-squares solution if the system is not square. - .. NOTE:: - - In Sage one can also write ``B / A`` for - ``A.solve_left(B)``, that is, Sage implements "the - MATLAB/Octave slash operator". - INPUT: - ``B`` -- a matrix or vector + - ``extend`` -- boolean (default: ``True``); when set to ``True``, + some solvers will return solutions over a larger ring than the + base ring of the inputs (a typical case are rational solutions + for integer linear systems). When set to ``False``, a solution + over the base ring is returned, with a :class:`ValueError` + being raised if none exists. + - ``check`` -- boolean (default: ``True``); verify the answer if the system is non-square or rank-deficient, and if its entries lie in an exact ring. Meaningless over most inexact @@ -446,16 +447,16 @@ cdef class Matrix(Matrix1): """ if isinstance(B, Vector): try: - return self.transpose().solve_right(B, check=check) + return self.transpose().solve_right(B, check=check, extend=extend) except ValueError as e: raise e.__class__(str(e).replace('row', 'column')) else: try: - return self.transpose().solve_right(B.transpose(), check=check).transpose() + return self.transpose().solve_right(B.transpose(), check=check, extend=extend).transpose() except ValueError as e: raise e.__class__(str(e).replace('row', 'column')) - def solve_right(self, B, check=True): + def solve_right(self, B, check=True, *, extend=True): r""" Try to find a solution `X` to the equation `A X = B`. @@ -469,16 +470,17 @@ cdef class Matrix(Matrix1): field, this method computes a least-squares solution if the system is not square. - .. NOTE:: - - DEPRECATED. In Sage one can also write ``A \ B`` for - ``A.solve_right(B)``, that is, Sage implements "the - MATLAB/Octave backslash operator". - INPUT: - ``B`` -- a matrix or vector + - ``extend`` -- boolean (default: ``True``); when set to ``True``, + some solvers will return solutions over a larger ring than the + base ring of the inputs (a typical case are rational solutions + for integer linear systems). When set to ``False``, a solution + over the base ring is returned, with a :class:`ValueError` + being raised if none exists. + - ``check`` -- boolean (default: ``True``); verify the answer if the system is non-square or rank-deficient, and if its entries lie in an exact ring. Meaningless over most inexact @@ -897,7 +899,7 @@ cdef class Matrix(Matrix1): L = B.base_ring() # first coerce both elements to parent over same base ring P = K if L is K else coercion_model.common_parent(K, L) - if P not in _Fields and P.is_integral_domain(): + if P not in _Fields and P.is_integral_domain() and extend: # the non-integral-domain case is handled separatedly below P = P.fraction_field() if L is not P: @@ -944,6 +946,13 @@ cdef class Matrix(Matrix1): C = B.column() if b_is_vec else B + if not extend: + try: + X = self._solve_right_hermite_form(C) + except NotImplementedError: + X = self._solve_right_smith_form(C) + return X.column(0) if b_is_vec else X + if not self.is_square(): X = self._solve_right_general(C, check=check) else: @@ -952,11 +961,7 @@ cdef class Matrix(Matrix1): except NotFullRankError: X = self._solve_right_general(C, check=check) - if b_is_vec: - # Convert back to a vector - return X.column(0) - else: - return X + return X.column(0) if b_is_vec else X def _solve_right_nonsingular_square(self, B, check_rank=True): r""" @@ -1071,6 +1076,138 @@ cdef class Matrix(Matrix1): raise ValueError("matrix equation has no solutions") return X + def _solve_right_smith_form(self, B): + r""" + Solve a matrix equation over a PID using the Smith normal form. + + EXAMPLES:: + + sage: A = matrix(ZZ, 5, 7, [(1 + x^4) % 55 for x in range(5*7)]); A + [ 1 2 17 27 37 21 32] + [37 27 17 46 12 2 17] + [27 26 32 32 37 27 6] + [ 2 12 2 17 16 37 32] + [32 37 16 17 2 12 2] + + sage: y = vector(ZZ, [1, 2, 3, 4, 5]) + sage: x, = A._solve_right_smith_form(y.column()).columns() + sage: x.parent() + Ambient free module of rank 7 over the principal ideal domain Integer Ring + sage: A * x == y + True + sage: x + (10, 1530831087980480, -2969971929450215, -178745029498097, 2320752168397186, -806846536262381, -520939892126393) + + sage: z = vector(ZZ, [-4, -1, 1, 5, 14, 31, 4]) + sage: x, = A.transpose()._solve_right_smith_form(z.column()).columns() + sage: x.parent() + Ambient free module of rank 5 over the principal ideal domain Integer Ring + sage: x * A == z + True + sage: x + (-1, 0, 1, 1, -1) + + sage: X = A._solve_right_smith_form(identity_matrix(ZZ,5)) + sage: X.parent() + Full MatrixSpace of 7 by 5 dense matrices over Integer Ring + sage: X + [ -1 0 2 0 1] + [-156182670342972 -2199494166584 310625144000132 1293916 151907461896112] + [ 303010665531453 4267248023008 -602645168043247 -2510332 -294716315371323] + [ 18236418267658 256820398262 -36269645268572 -151082 -17737230430447] + [-236774176922867 -3334450741532 470910201757143 1961587 230292926737068] + [ 82318322106118 1159275026338 -163719448527234 -681977 -80065012022313] + [ 53148766104440 748485096017 -105705345467375 -440318 -51693918051894] + """ + S,U,V = self.smith_form() + + n,m = self.dimensions() + r = B.ncols() + + X_ = [] + for d, v in zip(S.diagonal(), (U*B).rows()): + if d: + X_.append(v / d) + elif v: + raise ValueError("matrix equation has no solutions") + else: + X_.append([0] * r) + + X_ += [[0] * r] * (m - n) # arbitrary + + from sage.matrix.constructor import matrix + try: + X_ = matrix(self.base_ring(), m, r, X_) + except TypeError: + raise ValueError("matrix equation has no solutions") + + return V * X_ + + def _solve_right_hermite_form(self, B): + r""" + Solve a matrix equation over a PID using the Hermite normal form. + + EXAMPLES:: + + sage: A = matrix(ZZ, 5, 7, [(1 + x^4) % 55 for x in range(5*7)]); A + [ 1 2 17 27 37 21 32] + [37 27 17 46 12 2 17] + [27 26 32 32 37 27 6] + [ 2 12 2 17 16 37 32] + [32 37 16 17 2 12 2] + + sage: y = vector(ZZ, [1, 2, 3, 4, 5]) + sage: x, = A._solve_right_hermite_form(y.column()).columns() + sage: x.parent() + Ambient free module of rank 7 over the principal ideal domain Integer Ring + sage: A * x == y + True + sage: x + (7903368738919, 1664769092987, -13561109876830, -5130794714802, 11219929764616, -2728570619520, 0) + + sage: z = vector(ZZ, [-4, -1, 1, 5, 14, 31, 4]) + sage: x, = A.transpose()._solve_right_hermite_form(z.column()).columns() + sage: x.parent() + Ambient free module of rank 5 over the principal ideal domain Integer Ring + sage: x * A == z + True + sage: x + (-1, 0, 1, 1, -1) + + sage: X = A._solve_right_hermite_form(identity_matrix(ZZ,5)) + sage: X.parent() + Full MatrixSpace of 7 by 5 dense matrices over Integer Ring + sage: X + [ 682282983347 1029536563053 550188901703 877862915448 -1147487] + [ 143716389918 216862037808 115892034032 184913433473 -241707] + [-1170705152437 -1766545243552 -944049606627 -1506293815517 1968932] + [ -442931873812 -668365722378 -357177603912 -569900577297 744938] + [ 968595469303 1461570161933 781069571508 1246248350502 -1629017] + [ -235552378240 -355438713600 -189948023680 -303074680960 396160] + [ 0 0 0 0 0] + """ + H,U = self.transpose().hermite_form(transformation=True) + H = H.transpose() + U = U.transpose() +# assert self*U == H + + n,m = self.dimensions() + r = B.ncols() + + from sage.matrix.constructor import matrix + X_ = matrix(self.base_ring(), m, r) + for i in range(min(n,m)): + v = B[i,:] + v -= H[i,:i] * X_[:i] + d = H[i][i] + try: + X_[i] = v / d + except (ZeroDivisionError, TypeError) as e: + raise ValueError("matrix equation has no solution") +# assert H*X_ == B + + return U * X_ + def prod_of_row_sums(self, cols): r""" Calculate the product of all row sums of a submatrix of `A`