diff --git a/src/sage/matroids/lean_matrix.pxd b/src/sage/matroids/lean_matrix.pxd
index 1cce7e863ec..a112647a0bb 100644
--- a/src/sage/matroids/lean_matrix.pxd
+++ b/src/sage/matroids/lean_matrix.pxd
@@ -54,6 +54,7 @@ cdef class BinaryMatrix(LeanMatrix):
     cdef inline list row_union(self, object L)   # Not a Sage matrix operation
 
     cdef LeanMatrix matrix_from_rows_and_columns(self, rows, columns)
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns)
 
     cdef list _character(self, bitset_t x)
     cdef BinaryMatrix _distinguish_by(self, BinaryMatrix P)
@@ -77,7 +78,7 @@ cdef class TernaryMatrix(LeanMatrix):
     cdef inline long row_inner_product(self, long i, long j)   # Not a Sage matrix operation
     cdef void row_subs(self, long x, long y)   # Not a Sage matrix operation
     cdef void _row_negate(self, long x)
-
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns)
 
 cdef class QuaternaryMatrix(LeanMatrix):
     cdef bitset_t *_M0    # _M0[i] = 1-support of row i
@@ -91,7 +92,7 @@ cdef class QuaternaryMatrix(LeanMatrix):
     cdef inline long row_len(self, long i) except -1   # Not a Sage matrix operation
     cdef inline row_inner_product(self, long i, long j)   # Not a Sage matrix operation
     cdef inline int _row_div(self, long x, object s) except -1
-
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns)
     cdef void conjugate(self)   # Not a Sage matrix operation
 
 
diff --git a/src/sage/matroids/lean_matrix.pyx b/src/sage/matroids/lean_matrix.pyx
index 55f310f5f60..df3b3bc7c30 100644
--- a/src/sage/matroids/lean_matrix.pyx
+++ b/src/sage/matroids/lean_matrix.pyx
@@ -1158,6 +1158,63 @@ cdef class BinaryMatrix(LeanMatrix):
                     bitset_add(A._M[r], c)
         return A
 
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns):
+        """
+        Return a submatrix indexed by indicated rows and columns, as well as
+        the column order of the resulting submatrix.
+        """
+        cdef BinaryMatrix A = BinaryMatrix(len(rows), len(columns))
+        cdef long r, c, lc, lg
+        cdef mp_bitcnt_t *cols
+        # deal with trivial case
+        lc = len(columns)
+        if lc == 0:
+            return A, []
+        # write [c for c in columns if c<lc] as bitset `mask` and
+        # write [c for c in columns if c>=lc] as array `cols`
+        cdef bitset_t mask
+        bitset_init(mask, lc)
+        bitset_clear(mask)
+        cols = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        g = 0
+        for c in columns:
+            if c<lc:
+                bitset_add(mask, c)
+            else:
+                cols[g] = c
+                g = g+1
+        # write [ c for c in range(lc) if c not in columns] as array `gaps`
+        cdef mp_bitcnt_t *gaps
+        gaps = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        bitset_complement(mask, mask)
+        g = 0
+        c = bitset_first(mask)
+        while c>=0:
+            gaps[g] = c
+            g = g+1
+            c =  bitset_next(mask, c+1)
+        lg = g
+        bitset_complement(mask, mask)
+        # copy relevant part of this matrix into A
+        cdef bitset_t row, row2
+        for r in xrange(len(rows)):
+            row = self._M[rows[r]]
+            row2 = A._M[r]
+            bitset_intersection(row2, row, mask) # yes, this is safe
+            for g in xrange(lg):
+                if bitset_in(row, cols[g]):
+                    bitset_add(row2, gaps[g])
+        # record order of the columns in list `order`
+        cdef list order = range(lc)
+        g = 0
+        for g in xrange(lg):
+            order[gaps[g]] = cols[g]
+        # free up the two arrays and the bitset
+        sage_free(gaps)
+        sage_free(cols)
+        bitset_free(mask)
+        return A, order
+
     cdef list _character(self, bitset_t x):   # Not a Sage matrix operation
         """
         Return the vector of intersection lengths of the rows with ``x``.
@@ -1775,6 +1832,71 @@ cdef class TernaryMatrix(LeanMatrix):
         M.resize(self._nrows)
         return M
 
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns):
+        """
+        Return a submatrix indexed by indicated rows and columns, as well as
+        the column order of the resulting submatrix.
+        """
+        cdef TernaryMatrix A = TernaryMatrix(len(rows), len(columns))
+        cdef long r, c, lc, lg
+        cdef mp_bitcnt_t *cols
+        # deal with trivial case
+        lc = len(columns)
+        if lc == 0:
+            return A, []
+        # write [c for c in columns if c<lc] as bitset `mask` and
+        # write [c for c in columns if c>=lc] as array `cols`
+        cdef bitset_t mask
+        bitset_init(mask, lc)
+        bitset_clear(mask)
+        cols = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        g = 0
+        for c in columns:
+            if c<lc:
+                bitset_add(mask, c)
+            else:
+                cols[g] = c
+                g = g+1
+        # write [ c for c in range(lc) if c not in columns] as array `gaps`
+        cdef mp_bitcnt_t *gaps
+        gaps = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        bitset_complement(mask, mask)
+        g = 0
+        c = bitset_first(mask)
+        while c>=0:
+            gaps[g] = c
+            g = g+1
+            c =  bitset_next(mask, c+1)
+        lg = g
+        bitset_complement(mask, mask)
+        # copy relevant part of this matrix into A
+        cdef bitset_t row0, row1, row0_2, row1_2
+        cdef mp_bitcnt_t p, q
+        for r in xrange(len(rows)):
+            row0 = self._M0[rows[r]]
+            row1 = self._M1[rows[r]]
+            row0_2 = A._M0[r]
+            row1_2 = A._M1[r]
+            bitset_intersection(row0_2, row0, mask) # yes, this is safe
+            bitset_intersection(row1_2, row1, mask) # yes, this is safe
+            for g in xrange(lg):
+                p = cols[g]
+                if bitset_in(row0, p):
+                    q = gaps[g]
+                    bitset_add(row0_2, q)
+                    if bitset_in(row1, p):
+                        bitset_add(row1_2, q)
+        # record order of the columns in list `order`
+        cdef list order = range(lc)
+        g = 0
+        for g in xrange(lg):
+            order[gaps[g]] = cols[g]
+        # free up the two arrays and the bitset
+        sage_free(gaps)
+        sage_free(cols)
+        bitset_free(mask)
+        return A, order
+
     def __richcmp__(left, right, op):
         """
         Compare two matrices.
@@ -2293,6 +2415,71 @@ cdef class QuaternaryMatrix(LeanMatrix):
         M.resize(self._nrows)
         return M
 
+    cdef matrix_from_rows_and_columns_reordered(self, rows, columns):
+        """
+        Return a submatrix indexed by indicated rows and columns, as well as
+        the column order of the resulting submatrix.
+        """
+        cdef QuaternaryMatrix A = QuaternaryMatrix(len(rows), len(columns), ring = self._gf4)
+        cdef long r, c, lc, lg
+        cdef mp_bitcnt_t *cols
+        # deal with trivial case
+        lc = len(columns)
+        if lc == 0:
+            return A, []
+        # write [c for c in columns if c<lc] as bitset `mask` and
+        # write [c for c in columns if c>=lc] as array `cols`
+        cdef bitset_t mask
+        bitset_init(mask, lc)
+        bitset_clear(mask)
+        cols = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        g = 0
+        for c in columns:
+            if c<lc:
+                bitset_add(mask, c)
+            else:
+                cols[g] = c
+                g = g+1
+        # write [ c for c in range(lc) if c not in columns] as array `gaps`
+        cdef mp_bitcnt_t *gaps
+        gaps = <mp_bitcnt_t*>sage_malloc(lc*sizeof(mp_bitcnt_t))
+        bitset_complement(mask, mask)
+        g = 0
+        c = bitset_first(mask)
+        while c>=0:
+            gaps[g] = c
+            g = g+1
+            c =  bitset_next(mask, c+1)
+        lg = g
+        bitset_complement(mask, mask)
+        # copy relevant part of this matrix into A
+        cdef bitset_t row0, row1, row0_2, row1_2
+        cdef mp_bitcnt_t p, q
+        for r in xrange(len(rows)):
+            row0 = self._M0[rows[r]]
+            row1 = self._M1[rows[r]]
+            row0_2 = A._M0[r]
+            row1_2 = A._M1[r]
+            bitset_intersection(row0_2, row0, mask) # yes, this is safe
+            bitset_intersection(row1_2, row1, mask)
+            for g in xrange(lg):
+                p = cols[g]
+                q = gaps[g]
+                if bitset_in(row0, p):
+                    bitset_add(row0_2, q)
+                if bitset_in(row1, p):
+                    bitset_add(row1_2, q)
+        # record order of the columns in list `order`
+        cdef list order = range(lc)
+        g = 0
+        for g in xrange(lg):
+            order[gaps[g]] = cols[g]
+        # free up the two arrays and the bitset
+        sage_free(gaps)
+        sage_free(cols)
+        bitset_free(mask)
+        return A, order
+    
     def __neg__(self):
         """
         Negate the matrix.
diff --git a/src/sage/matroids/linear_matroid.pyx b/src/sage/matroids/linear_matroid.pyx
index ce3d8ab3e1d..e14ef62413c 100644
--- a/src/sage/matroids/linear_matroid.pyx
+++ b/src/sage/matroids/linear_matroid.pyx
@@ -3470,9 +3470,11 @@ cdef class BinaryMatroid(LinearMatroid):
         self._move_current_basis(contractions, deletions)
         bas = list(self.basis() - contractions)
         R = [self._prow[self._idx[b]] for b in bas]
-        C = [c for c in range(len(self._E)) if self._E[c] not in deletions | contractions]
-        return BinaryMatroid(matrix=(<BinaryMatrix>self._A).matrix_from_rows_and_columns(R, C),
-                             groundset=[self._E[c] for c in C],
+        F = self.groundset() - (deletions | contractions)
+        C = [self._idx[f] for f in F]
+        A, C2 = (<BinaryMatrix>self._A).matrix_from_rows_and_columns_reordered(R, C)
+        return BinaryMatroid(matrix= A ,
+                             groundset=[self._E[c] for c in C2],
                              basis=bas,
                              keep_initial_representation=False)
 
@@ -4417,11 +4419,13 @@ cdef class TernaryMatroid(LinearMatroid):
         self._move_current_basis(contractions, deletions)
         bas = list(self.basis() - contractions)
         R = [self._prow[self._idx[b]] for b in bas]
-        C = [c for c in range(len(self._E)) if self._E[c] not in deletions | contractions]
-        return TernaryMatroid(matrix=(<TernaryMatrix>self._A).matrix_from_rows_and_columns(R, C),
-                              groundset=[self._E[c] for c in C],
-                              basis=bas,
-                              keep_initial_representation=False)
+        F = self.groundset() - (deletions | contractions)
+        C = [self._idx[f] for f in F]
+        A, C2 = (<TernaryMatrix>self._A).matrix_from_rows_and_columns_reordered(R, C)
+        return TernaryMatroid(matrix= A ,
+                             groundset=[self._E[c] for c in C2],
+                             basis=bas,
+                             keep_initial_representation=False)
 
     cpdef is_valid(self):
         r"""
@@ -5127,11 +5131,13 @@ cdef class QuaternaryMatroid(LinearMatroid):
         self._move_current_basis(contractions, deletions)
         bas = list(self.basis() - contractions)
         R = [self._prow[self._idx[b]] for b in bas]
-        C = [c for c in range(len(self._E)) if self._E[c] not in deletions | contractions]
-        return QuaternaryMatroid(matrix=(<QuaternaryMatrix>self._A).matrix_from_rows_and_columns(R, C),
-                                 groundset=[self._E[c] for c in C],
-                                 basis=bas,
-                                 keep_initial_representation=False)
+        F = self.groundset() - (deletions | contractions)
+        C = [self._idx[f] for f in F]
+        A, C2 = (<QuaternaryMatrix>self._A).matrix_from_rows_and_columns_reordered(R, C)
+        return QuaternaryMatroid(matrix= A ,
+                             groundset=[self._E[c] for c in C2],
+                             basis=bas,
+                             keep_initial_representation=False)
 
     cpdef is_valid(self):
         r"""