diff --git a/lib/OptimalBranchingMIS/src/OptimalBranchingMIS.jl b/lib/OptimalBranchingMIS/src/OptimalBranchingMIS.jl
index d534d32..0802a59 100644
--- a/lib/OptimalBranchingMIS/src/OptimalBranchingMIS.jl
+++ b/lib/OptimalBranchingMIS/src/OptimalBranchingMIS.jl
@@ -10,7 +10,7 @@ using SparseArrays
 using KaHyPar
 
 export MISProblem
-export MISReducer, XiaoReducer, TensorNetworkReducer
+export MISReducer, XiaoReducer, TensorNetworkReducer, SubsolverReducer, MatryoshkaReducer
 export MinBoundarySelector, MinBoundaryHighDegreeSelector,KaHyParSelector
 export TensorNetworkSolver
 export NumOfVertices, D3Measure
diff --git a/lib/OptimalBranchingMIS/src/reducer.jl b/lib/OptimalBranchingMIS/src/reducer.jl
index 22f170b..da97b12 100644
--- a/lib/OptimalBranchingMIS/src/reducer.jl
+++ b/lib/OptimalBranchingMIS/src/reducer.jl
@@ -49,6 +49,29 @@ Base.show(io::IO, reducer::TensorNetworkReducer) = print(io,
         └── recheck - $(reducer.recheck)
     """)
 
+"""
+    struct SubsolverReducer{TR, TS} <: AbstractReducer
+
+After the size of the problem is smaller than the threshold, use a subsolver to find reduction rules.
+
+# Fields
+- `reducer::TR = XiaoReducer()`: The reducer used to reduce the graph.
+- `subsolver::Symbol = :xiao`: subsolvers include `:mis2`, `:xiao` and `:ip`.
+- `threshold::Int = 100`: The threshold for using the subsolver.
+"""
+@kwdef struct SubsolverReducer <: AbstractReducer
+    reducer::AbstractReducer = XiaoReducer()
+    subsolver::Symbol = :xiao # :mis2, :xiao or :ip
+    threshold::Int = 100 # the threshold for using the subsolver
+end
+
+@kwdef struct MatryoshkaReducer <: AbstractReducer
+    reducer::AbstractReducer
+    subconfig::BranchingStrategy
+    subreducer::Union{MatryoshkaReducer, SubsolverReducer}
+    threshold::Int
+end
+
 """
     reduce_problem(::Type{R}, p::MISProblem, ::MISReducer) where R
 
@@ -71,7 +94,7 @@ The function checks the number of vertices in the graph:
 - If there are two vertices, it returns an empty instance and a count based on the presence of an edge between them.
 - For graphs with more than two vertices, it calculates the degrees of the vertices and identifies the vertex with the minimum degree to determine which vertices to remove.
 """
-function OptimalBranchingCore.reduce_problem(::Type{R}, p::MISProblem, reducer::Union{MISReducer, XiaoReducer, TensorNetworkReducer}) where R
+function OptimalBranchingCore.reduce_problem(::Type{R}, p::MISProblem, reducer::Union{MISReducer, XiaoReducer, TensorNetworkReducer, SubsolverReducer, MatryoshkaReducer}) where R
     g_new, r, _ = reduce_graph(p.g, reducer)
     if (nv(g_new) == nv(p.g)) && iszero(r)
         return p, R(0)
@@ -228,4 +251,30 @@ function best_intersect(p::MISProblem, tbl::BranchingTable, measure::AbstractMea
         end
         return best_cl
     end
+end
+
+function reduce_graph(g::SimpleGraph{Int}, reducer::SubsolverReducer)
+    if nv(g) <= reducer.threshold
+        # use the subsolver the directly solve the problem
+        if reducer.subsolver == :mis2
+            res = mis2(EliminateGraph(g))
+        elseif reducer.subsolver == :xiao
+            res = counting_xiao2013(g).size
+        elseif reducer.subsolver == :ip
+            res = ip_mis(g)
+        else
+            error("Subsolver $(reducer.subsolver) not supported, please choose from :mis2, :xiao or :ip")
+        end
+        return SimpleGraph(0), res, Int[]
+    else
+        return reduce_graph(g, reducer.reducer)
+    end
+end
+
+function reduce_graph(g::SimpleGraph{Int}, reducer::MatryoshkaReducer)
+    if nv(g) > reducer.threshold
+        return reduce_graph(g, reducer.reducer)
+    else
+        return SimpleGraph(0), mis_size(g, branching_strategy = reducer.subconfig, reducer = reducer.subreducer), Int[]
+    end
 end
\ No newline at end of file
diff --git a/lib/OptimalBranchingMIS/test/branch.jl b/lib/OptimalBranchingMIS/test/branch.jl
index 8c7c3c9..8bbaf85 100644
--- a/lib/OptimalBranchingMIS/test/branch.jl
+++ b/lib/OptimalBranchingMIS/test/branch.jl
@@ -11,7 +11,7 @@ using Test, Random
     mis_xiao = counting_xiao2013(g).size
     p = MISProblem(g)
 
-    for set_cover_solver in [IPSolver(max_itr = 10, verbose = false), LPSolver(max_itr = 10, verbose = false)], measure in [D3Measure(), NumOfVertices()], reducer in [NoReducer(), MISReducer(), XiaoReducer(), TensorNetworkReducer()], prune_by_env in [true, false], selector in [MinBoundarySelector(2), MinBoundaryHighDegreeSelector(2, 6, 0), MinBoundaryHighDegreeSelector(2, 6, 1)]
+    for set_cover_solver in [IPSolver(max_itr = 10, verbose = false), LPSolver(max_itr = 10, verbose = false)], measure in [D3Measure(), NumOfVertices()], reducer in [NoReducer(), MISReducer(), XiaoReducer(), TensorNetworkReducer(), SubsolverReducer()], prune_by_env in [true, false], selector in [MinBoundarySelector(2), MinBoundaryHighDegreeSelector(2, 6, 0), MinBoundaryHighDegreeSelector(2, 6, 1)]
         @info "set_cover_solver = $set_cover_solver, measure = $measure, reducer = $reducer, prune_by_env = $prune_by_env, selector = $selector"
         branching_strategy = BranchingStrategy(; set_cover_solver, table_solver=TensorNetworkSolver(; prune_by_env), selector=selector, measure)
         res = branch_and_reduce(p, branching_strategy, reducer, MaxSize)