diff --git a/experimental/SI/si.jl b/experimental/SI/si.jl
index dd40a2c2..da2c27b5 100644
--- a/experimental/SI/si.jl
+++ b/experimental/SI/si.jl
@@ -87,3 +87,31 @@ flag, gb_2 = Groebner.groebner_apply!(graph, G_zp)
         Groebner.groebner_apply!(graph, G_zp)
     end
 end
+
+R, x = PolynomialRing(Nemo.GF(2^31 - 1), ["x$i" for i in 1:15], ordering=:degrevlex)
+
+f = [a^rand(1:3) * b^rand(1:3) + c^rand(1:2) for a in x for b in x for c in x];
+
+graph, gb = Groebner.groebner_learn(f);
+
+@benchmark flag, gb_2 = Groebner.groebner_apply!($graph, $f)
+@benchmark Groebner.groebner($f)
+
+@myprof begin
+    for i in 1:100
+        Groebner.groebner_apply!(graph, f)
+    end
+end
+
+c = Groebner.rootn(9, ground=Nemo.GF(2^31 - 1), ordering=:degrevlex)
+graph, gb_1 = Groebner.groebner_learn(c);
+flag, gb_2 = Groebner.groebner_apply!(graph, c);
+
+@myprof begin
+    for i in 1:100
+        Groebner.groebner_apply!(graph, c)
+    end
+end
+
+c = Groebner.rootn(9, ground=Nemo.GF(2^31 - 1), ordering=:degrevlex)
+
diff --git a/src/f4/f4.jl b/src/f4/f4.jl
index 3a6c142d..2b415557 100644
--- a/src/f4/f4.jl
+++ b/src/f4/f4.jl
@@ -737,7 +737,7 @@ function basis_well_formed(key, ring, basis, hashtable)
         else
             length(basis.coeffs[i]) == length(basis.monoms[i]) && continue
             if key in (:input_f4_apply!, :output_f4_apply!)
-                @log level = 1 "Unlucky but probably not fatal cancellation at index $(i)" length(
+                @log level = 10^3 "Unlucky but probably not fatal cancellation at index $(i)" length(
                     basis.monoms[i]
                 ) length(basis.coeffs[i])
             else
diff --git a/src/f4/graph.jl b/src/f4/graph.jl
index a362373d..3a9899e6 100644
--- a/src/f4/graph.jl
+++ b/src/f4/graph.jl
@@ -6,9 +6,6 @@ mutable struct ComputationGraphF4{C, M, Ord}
     gb_basis::Basis{C}
     hashtable::MonomialHashtable{M, Ord}
     input_permutation::Vector{Int}
-    # The number of columns,
-    # The number of upper rows,
-    # The number of lower rows
     matrix_infos::Vector{NamedTuple{(:nup, :nlow, :ncols), Tuple{Int64, Int64, Int64}}}
     matrix_nonzeroed_rows::Vector{Vector{Int}}
     matrix_upper_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
@@ -16,15 +13,8 @@ mutable struct ComputationGraphF4{C, M, Ord}
     output_nonredundant_indices::Vector{Int}
     nonredundant_indices_before_reduce::Vector{Int}
     output_sort_indices::Vector{Int}
-    # matrix_autored_lower_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
-    # matrix_autored_upper_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
-    # F4 iteration number --> index in the basis
-    # matrix_tobereduced_rows::Vector{Vector{Int}}
-    # matrix_tobereduced_mult::Vector{Vector{MonomIdx}}
-    # # F4 iteration number --> index in the basis
-    # matrix_reducers_rows::Vector{Vector{Int}}
-    # matrix_reducers_mult::Vector{Vector{MonomIdx}}
-    # matrix_columns::Vector{Vector{Int}}
+    # hashtables_symbolic::Vector{MonomialHashtable{M, Ord}}
+    matrix_sorted_columns::Vector{Vector{Int}}
 end
 
 function initialize_computation_graph_f4(
@@ -47,7 +37,8 @@ function initialize_computation_graph_f4(
         Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
         Vector{Int}(),
         Vector{Int}(),
-        Vector{Int}()
+        Vector{Int}(),
+        Vector{Vector{Int}}()
         # Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
         # Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
     )
diff --git a/src/f4/learn-apply.jl b/src/f4/learn-apply.jl
index a70b28d7..ec33c501 100644
--- a/src/f4/learn-apply.jl
+++ b/src/f4/learn-apply.jl
@@ -1,4 +1,52 @@
 
+# NOTE: in some problems, a lot of time is spent in sorting the matrix columns,
+# hence we cache them as well
+
+function column_to_monom_mapping!(graph, matrix, symbol_ht)
+    # monoms from symbolic table represent one column in the matrix
+    hdata = symbol_ht.hashdata
+    load = symbol_ht.load
+
+    # number of pivotal cols
+    k = 0
+    @inbounds for i in (symbol_ht.offset):load
+        if hdata[i].idx == 2
+            k += 1
+        end
+    end
+
+    col2hash = matrix.col2hash
+
+    matrix.nleft = k  # CHECK!
+    # -1 as long as hashtable load is always 1 more than actual
+    matrix.nright = load - matrix.nleft - 1
+
+    # store the other direction of mapping,
+    # hash -> column
+    @inbounds for k in 1:length(col2hash)
+        hdata[col2hash[k]].idx = k
+    end
+
+    @inbounds for k in 1:(matrix.nup)
+        row = matrix.uprows[k]
+        for j in 1:length(row)
+            row[j] = hdata[row[j]].idx
+        end
+    end
+
+    @inbounds for k in 1:(matrix.nlow)
+        row = matrix.lowrows[k]
+        for j in 1:length(row)
+            row[j] = hdata[row[j]].idx
+        end
+    end
+
+    matrix.ncols = matrix.nleft + matrix.nright
+
+    @assert matrix.nleft + matrix.nright == symbol_ht.load - 1 == matrix.ncols
+    @assert matrix.nlow + matrix.nup == matrix.nrows
+end
+
 function reduction_apply!(
     graph::ComputationGraphF4,
     ring,
@@ -6,9 +54,16 @@ function reduction_apply!(
     matrix,
     ht,
     rng,
-    symbol_ht
+    symbol_ht,
+    iter
 )
-    column_to_monom_mapping!(matrix, symbol_ht)
+    if length(graph.matrix_sorted_columns) < iter
+        column_to_monom_mapping!(matrix, symbol_ht)
+        push!(graph.matrix_sorted_columns, matrix.col2hash)
+    else
+        matrix.col2hash = graph.matrix_sorted_columns[iter]
+        column_to_monom_mapping!(graph, matrix, symbol_ht)
+    end
 
     sort_matrix_upper_rows_decreasing!(matrix) # for pivots,  AB part
     sort_matrix_lower_rows_increasing!(matrix) # for reduced, CD part
@@ -101,7 +156,8 @@ function reducegb_f4_apply!(
     basis::Basis,
     matrix::MacaulayMatrix,
     hashtable::MonomialHashtable{M},
-    symbol_ht::MonomialHashtable{M}
+    symbol_ht::MonomialHashtable{M},
+    iter
 ) where {M}
     @log level = -5 "Entering apply autoreduction" basis
 
@@ -168,7 +224,14 @@ function reducegb_f4_apply!(
     matrix.nup = nup
     matrix.size = matrix.nrows
 
-    column_to_monom_mapping!(matrix, symbol_ht)
+    @log level = -2 length(graph.matrix_sorted_columns) iter
+    if length(graph.matrix_sorted_columns) < iter
+        column_to_monom_mapping!(matrix, symbol_ht)
+        push!(graph.matrix_sorted_columns, matrix.col2hash)
+    else
+        matrix.col2hash = graph.matrix_sorted_columns[iter]
+        column_to_monom_mapping!(graph, matrix, symbol_ht)
+    end
     matrix.ncols = matrix.nleft + matrix.nright
 
     sort_matrix_upper_rows_decreasing!(matrix)
@@ -258,8 +321,16 @@ function f4_apply!(graph, ring, basis::Basis{C}, params) where {C <: Coeff}
         symbolic_preprocessing!(graph, iters, basis, matrix, hashtable, symbol_ht)
         @log level = -5 "After symbolic preprocessing:" matrix
 
-        flag =
-            reduction_apply!(graph, ring, basis, matrix, hashtable, params.rng, symbol_ht)
+        flag = reduction_apply!(
+            graph,
+            ring,
+            basis,
+            matrix,
+            hashtable,
+            params.rng,
+            symbol_ht,
+            iters
+        )
         if !flag
             # Unlucky cancellation of basis coefficients happened
             return false
@@ -281,7 +352,15 @@ function f4_apply!(graph, ring, basis::Basis{C}, params) where {C <: Coeff}
     if params.reduced
         @log level = -6 "Autoreducing the final basis.."
         symbol_ht = initialize_secondary_hashtable(hashtable)
-        flag = reducegb_f4_apply!(graph, ring, basis, matrix, hashtable, symbol_ht)
+        flag = reducegb_f4_apply!(
+            graph,
+            ring,
+            basis,
+            matrix,
+            hashtable,
+            symbol_ht,
+            iters_total + 1
+        )
         if !flag
             return false
         end
diff --git a/src/f4/matrix.jl b/src/f4/matrix.jl
index e8132afc..b895ee75 100644
--- a/src/f4/matrix.jl
+++ b/src/f4/matrix.jl
@@ -415,6 +415,7 @@ function reduce_dense_row_by_known_pivots_sparse!(
     # where k - number of structural nonzeros in new reduced row, k > 0
     @assert k > 0
     j = 1
+    # TODO: Unroll by a factor of 4
     @inbounds for i in np:ncols # starting from new pivot
         if densecoeffs[i] != uzero
             newrow[j] = i