diff --git a/docs/source/guide/lre-1-intro.md b/docs/source/guide/lre-1-intro.md
index 6604837306..8e195ad1ef 100644
--- a/docs/source/guide/lre-1-intro.md
+++ b/docs/source/guide/lre-1-intro.md
@@ -78,6 +78,7 @@ mitigated = execute_with_lre(
     fold_multiplier=fold_multiplier,
 )
 
+
 print(f"Error with mitigation (LRE): {abs(ideal - mitigated):.{3}}")
 ```
 
@@ -92,7 +93,7 @@ In this section we demonstrate the use of {func}`.multivariate_layer_scaling` an
 We start by creating a number of noise-scaled circuits which we will pass to the executor.
 
 ```{code-cell} ipython3
-from mitiq.lre import multivariate_layer_scaling
+from mitiq.lre.multivariate_scaling import multivariate_layer_scaling
 
 
 noise_scaled_circuits = multivariate_layer_scaling(circuit, degree, fold_multiplier)
@@ -125,7 +126,7 @@ The penultimate step here is to fetch the coefficients we'll use to combine the
 The astute reader will note that we haven't defined or used a `degree` or `fold_multiplier` parameter, and this is where they are both needed.
 
 ```{code-cell} ipython3
-from mitiq.lre import multivariate_richardson_coefficients
+from mitiq.lre.inference import multivariate_richardson_coefficients
 
 
 coefficients = multivariate_richardson_coefficients(
@@ -144,9 +145,7 @@ mitigated = sum(
     exp_val * coeff
     for exp_val, coeff in zip(noise_scaled_exp_values, coefficients)
 )
-print(
-    f"Error with mitigation (LRE): {abs(ideal - mitigated):.{3}}"
-)
+print(f"Error with mitigation (LRE): {abs(ideal - mitigated):.{3}}")
 ```
 
 As you can see we again see a nice improvement in the accuracy using a two stage application of LRE.
diff --git a/mitiq/lre/__init__.py b/mitiq/lre/__init__.py
index 353568cdd9..1060e50c45 100644
--- a/mitiq/lre/__init__.py
+++ b/mitiq/lre/__init__.py
@@ -5,11 +5,8 @@
 
 """Methods for scaling noise in circuits by layers and using multivariate extrapolation."""
 
-from mitiq.lre.multivariate_scaling.layerwise_folding import multivariate_layer_scaling
+from mitiq.lre.lre import execute_with_lre, mitigate_executor, lre_decorator
 
-from mitiq.lre.inference.multivariate_richardson import (
-    multivariate_richardson_coefficients,
-    sample_matrix,
-)
+from mitiq.lre import multivariate_scaling
 
-from mitiq.lre.lre import execute_with_lre, mitigate_executor, lre_decorator
\ No newline at end of file
+from mitiq.lre import inference
diff --git a/mitiq/lre/inference/__init__.py b/mitiq/lre/inference/__init__.py
new file mode 100644
index 0000000000..3411c080fd
--- /dev/null
+++ b/mitiq/lre/inference/__init__.py
@@ -0,0 +1,3 @@
+from mitiq.lre.inference.multivariate_richardson import (
+    multivariate_richardson_coefficients,
+)
diff --git a/mitiq/lre/lre.py b/mitiq/lre/lre.py
index c2fd7b52af..e292e15d19 100644
--- a/mitiq/lre/lre.py
+++ b/mitiq/lre/lre.py
@@ -12,10 +12,12 @@
 from cirq import Circuit
 
 from mitiq import QPROGRAM
-from mitiq.lre import (
-    multivariate_layer_scaling,
+from mitiq.lre.inference import (
     multivariate_richardson_coefficients,
 )
+from mitiq.lre.multivariate_scaling import (
+    multivariate_layer_scaling,
+)
 from mitiq.zne.scaling import fold_gates_at_random
 
 
diff --git a/mitiq/lre/multivariate_scaling/__init__.py b/mitiq/lre/multivariate_scaling/__init__.py
new file mode 100644
index 0000000000..816dfaa208
--- /dev/null
+++ b/mitiq/lre/multivariate_scaling/__init__.py
@@ -0,0 +1,6 @@
+
+"""Methods for scaling noise in circuits by layers as required for LRE."""
+
+from mitiq.lre.multivariate_scaling.layerwise_folding import (
+    multivariate_layer_scaling,
+)