Repository: bambinos/kulprit (language: python)

Issue:
Add L1 search
## Description

This PR will add to ``kulprit`` the L1 search functionality developed to afford computationally cheaper projection search paths as is suggested by Piironen et al. (2020). We further remove the analytic solution features from the repo since the added weight and complexity of the code, and their future maintenance outweighs the benefit derived from them.

## Checklist

- [x] Tests covering the new functionality have been added
- [x] Documentation has been updated OR the change is too minor to be documented
- [ ] Changes are listed in the `CHANGELOG.md` OR changes are insignificant



Discussion / hints:

Fix patch:
diff --git a/poetry.lock b/poetry.lock
index 2f1b4b2..e5e2832 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -773,6 +773,14 @@ MarkupSafe = ">=2.0"
 [package.extras]
 i18n = ["Babel (>=2.7)"]
 
+[[package]]
+name = "joblib"
+version = "1.1.0"
+description = "Lightweight pipelining with Python functions"
+category = "main"
+optional = false
+python-versions = ">=3.6"
+
 [[package]]
 name = "jsonschema"
 version = "4.6.0"
@@ -1787,6 +1795,26 @@ packaging = "*"
 [package.extras]
 test = ["pytest (>=6,!=7.0.0,!=7.0.1)", "pytest-cov (>=3.0.0)", "pytest-qt"]
 
+[[package]]
+name = "scikit-learn"
+version = "1.0.2"
+description = "A set of python modules for machine learning and data mining"
+category = "main"
+optional = false
+python-versions = ">=3.7"
+
+[package.dependencies]
+joblib = ">=0.11"
+numpy = ">=1.14.6"
+scipy = ">=1.1.0"
+threadpoolctl = ">=2.0.0"
+
+[package.extras]
+benchmark = ["matplotlib (>=2.2.3)", "pandas (>=0.25.0)", "memory-profiler (>=0.57.0)"]
+docs = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "seaborn (>=0.9.0)", "memory-profiler (>=0.57.0)", "sphinx (>=4.0.1)", "sphinx-gallery (>=0.7.0)", "numpydoc (>=1.0.0)", "Pillow (>=7.1.2)", "sphinx-prompt (>=1.3.0)", "sphinxext-opengraph (>=0.4.2)"]
+examples = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "seaborn (>=0.9.0)"]
+tests = ["matplotlib (>=2.2.3)", "scikit-image (>=0.14.5)", "pandas (>=0.25.0)", "pytest (>=5.0.1)", "pytest-cov (>=2.9.0)", "flake8 (>=3.8.2)", "black (>=21.6b0)", "mypy (>=0.770)", "pyamg (>=4.0.0)"]
+
 [[package]]
 name = "scipy"
 version = "1.7.3"
@@ -1875,6 +1903,14 @@ tornado = ">=6.1.0"
 [package.extras]
 test = ["pre-commit", "pytest-timeout", "pytest (>=6.0)"]
 
+[[package]]
+name = "threadpoolctl"
+version = "3.1.0"
+description = "threadpoolctl"
+category = "main"
+optional = false
+python-versions = ">=3.6"
+
 [[package]]
 name = "tinycss2"
 version = "1.1.1"
@@ -2105,7 +2141,7 @@ testing = ["pytest (>=6)", "pytest-checkdocs (>=2.4)", "pytest-flake8", "pytest-
 [metadata]
 lock-version = "1.1"
 python-versions = ">=3.7.2, <3.10"
-content-hash = "aaac7e6fdecf83cc758d726fe3ae2d2882c68b415dca46ccee43bf20bdfd215b"
+content-hash = "d5e67c8df156b494dedffa532aff18a47988944b34e5e0cf0b786370807ac939"
 
 [metadata.files]
 aeppl = [
@@ -2518,6 +2554,10 @@ jinja2 = [
     {file = "Jinja2-3.1.2-py3-none-any.whl", hash = "sha256:6088930bfe239f0e6710546ab9c19c9ef35e29792895fed6e6e31a023a182a61"},
     {file = "Jinja2-3.1.2.tar.gz", hash = "sha256:31351a702a408a9e7595a8fc6150fc3f43bb6bf7e319770cbc0db9df9437e852"},
 ]
+joblib = [
+    {file = "joblib-1.1.0-py2.py3-none-any.whl", hash = "sha256:f21f109b3c7ff9d95f8387f752d0d9c34a02aa2f7060c2135f465da0e5160ff6"},
+    {file = "joblib-1.1.0.tar.gz", hash = "sha256:4158fcecd13733f8be669be0683b96ebdbbd38d23559f54dca7205aea1bf1e35"},
+]
 jsonschema = [
     {file = "jsonschema-4.6.0-py3-none-any.whl", hash = "sha256:1c92d2db1900b668201f1797887d66453ab1fbfea51df8e4b46236689c427baf"},
     {file = "jsonschema-4.6.0.tar.gz", hash = "sha256:9d6397ba4a6c0bf0300736057f649e3e12ecbc07d3e81a0dacb72de4e9801957"},
@@ -2815,7 +2855,6 @@ nest-asyncio = [
 ]
 netcdf4 = [
     {file = "netCDF4-1.5.8-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:fd501ccb28ebae6770112968c750a14feb39cb495d788aa67c28360f7c1f2324"},
-    {file = "netCDF4-1.5.8-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:18e257843fc29846909557a9502a29e37b381ee7760923f9280d3b26d844db8c"},
     {file = "netCDF4-1.5.8-cp310-cp310-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:8b4fac95819b0c17ca5fc1a4e8bb31116b6b808cceca0aa8b475bb50abab1063"},
     {file = "netCDF4-1.5.8-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e118bfccda464a381187b1f9c771bf9581c83c566faab309a8ec3f781668da4e"},
     {file = "netCDF4-1.5.8-cp310-cp310-win32.whl", hash = "sha256:bdba6ea34680a4c1b7018a4a7155f6112acd063289923c0c61918707e9f26910"},
@@ -2830,17 +2869,13 @@ netcdf4 = [
     {file = "netCDF4-1.5.8-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:29426faabdc017e89572ff1d7835dab426aae4c22ad1a12d1877b932e969b6ac"},
     {file = "netCDF4-1.5.8-cp37-cp37m-win32.whl", hash = "sha256:225d17f7a487ebdab99640294203b61e39e01c951b4e6a4f578d8251623f5f5a"},
     {file = "netCDF4-1.5.8-cp37-cp37m-win_amd64.whl", hash = "sha256:f86399073b582bccd278006ee0213548e7037395e1119f1af9f4faad38279b1e"},
-    {file = "netCDF4-1.5.8-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:0a33a953b60ee30dcb78db174231f7ab61923331af9645f84adff684e9add4e2"},
     {file = "netCDF4-1.5.8-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:f92b92f758dfc903af2a8a287fd68a531f73ffd3e5be72b5ad1eb3f083e7aaa2"},
-    {file = "netCDF4-1.5.8-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:b21af57acca0d70c5401f8f779409ab4e818c505fb81706eea8d9475e1f0bb9b"},
     {file = "netCDF4-1.5.8-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:7326afa46fd0c1b50d30db9764a1eefbcff576fcffa8e48bef403094590563b8"},
     {file = "netCDF4-1.5.8-cp38-cp38-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:916434a13ea317934cf248fb70dd5476c498f1def71041fc7e3fd23882ef2cda"},
     {file = "netCDF4-1.5.8-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:75ecf1ef2c841aace539f3326d101efda7c99f6c3283c48a444370aba48525af"},
     {file = "netCDF4-1.5.8-cp38-cp38-win32.whl", hash = "sha256:bd35c37342d9051d2b8fb12a9208856cc59201a94c78a742a198c81813cb00a8"},
     {file = "netCDF4-1.5.8-cp38-cp38-win_amd64.whl", hash = "sha256:bdd344d8de65849fa200f69941f1b15a2611b11b307161ce2fd8ff42148507e8"},
-    {file = "netCDF4-1.5.8-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:06f7364086fd3ae097e757d2493dc1fe006e9ae9636a109a1e4c914db05d7e18"},
     {file = "netCDF4-1.5.8-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:0f570b5b4cc0434ef8a2a648fdebfa017de695ea7c836b24ae7b216ede4e3345"},
-    {file = "netCDF4-1.5.8-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:db02d42f7b9c7d68cec351ea63ef3fc2a1ad5e7e74fc7b570b34ceb8c7645bf2"},
     {file = "netCDF4-1.5.8-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5c883a02c55fd1e5b61ad4f83dd7f11f90b894e14d120ba678d9c33d9e4b3a77"},
     {file = "netCDF4-1.5.8-cp39-cp39-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:d784d6cf5baa90909f385bf60ead91138f13ff7f870467e458fb3650ef71b48d"},
     {file = "netCDF4-1.5.8-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:49a44c7382e5e1da39d8bab5d8e406ad30d46fda9386e85a3e69491e6caaca93"},
@@ -3276,6 +3311,40 @@ qtpy = [
     {file = "QtPy-2.1.0-py3-none-any.whl", hash = "sha256:aee0586081f943029312becece9f63977b0a9e3788f77a6ac8cc74802bb173d6"},
     {file = "QtPy-2.1.0.tar.gz", hash = "sha256:ca8cd4217175186344299ee4c0f7e7adcf362c70852ba35b255a534077025c06"},
 ]
+scikit-learn = [
+    {file = "scikit-learn-1.0.2.tar.gz", hash = "sha256:b5870959a5484b614f26d31ca4c17524b1b0317522199dc985c3b4256e030767"},
+    {file = "scikit_learn-1.0.2-cp310-cp310-macosx_10_13_x86_64.whl", hash = "sha256:da3c84694ff693b5b3194d8752ccf935a665b8b5edc33a283122f4273ca3e687"},
+    {file = "scikit_learn-1.0.2-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:75307d9ea39236cad7eea87143155eea24d48f93f3a2f9389c817f7019f00705"},
+    {file = "scikit_learn-1.0.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f14517e174bd7332f1cca2c959e704696a5e0ba246eb8763e6c24876d8710049"},
+    {file = "scikit_learn-1.0.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d9aac97e57c196206179f674f09bc6bffcd0284e2ba95b7fe0b402ac3f986023"},
+    {file = "scikit_learn-1.0.2-cp310-cp310-win_amd64.whl", hash = "sha256:d93d4c28370aea8a7cbf6015e8a669cd5d69f856cc2aa44e7a590fb805bb5583"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-macosx_10_13_x86_64.whl", hash = "sha256:85260fb430b795d806251dd3bb05e6f48cdc777ac31f2bcf2bc8bbed3270a8f5"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:a053a6a527c87c5c4fa7bf1ab2556fa16d8345cf99b6c5a19030a4a7cd8fd2c0"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:245c9b5a67445f6f044411e16a93a554edc1efdcce94d3fc0bc6a4b9ac30b752"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:158faf30684c92a78e12da19c73feff9641a928a8024b4fa5ec11d583f3d8a87"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:08ef968f6b72033c16c479c966bf37ccd49b06ea91b765e1cc27afefe723920b"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:16455ace947d8d9e5391435c2977178d0ff03a261571e67f627c8fee0f9d431a"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-win32.whl", hash = "sha256:2f3b453e0b149898577e301d27e098dfe1a36943f7bb0ad704d1e548efc3b448"},
+    {file = "scikit_learn-1.0.2-cp37-cp37m-win_amd64.whl", hash = "sha256:46f431ec59dead665e1370314dbebc99ead05e1c0a9df42f22d6a0e00044820f"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-macosx_10_13_x86_64.whl", hash = "sha256:ff3fa8ea0e09e38677762afc6e14cad77b5e125b0ea70c9bba1992f02c93b028"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:9369b030e155f8188743eb4893ac17a27f81d28a884af460870c7c072f114243"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:7d6b2475f1c23a698b48515217eb26b45a6598c7b1840ba23b3c5acece658dbb"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:285db0352e635b9e3392b0b426bc48c3b485512d3b4ac3c7a44ec2a2ba061e66"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5cb33fe1dc6f73dc19e67b264dbb5dde2a0539b986435fdd78ed978c14654830"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:b1391d1a6e2268485a63c3073111fe3ba6ec5145fc957481cfd0652be571226d"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bc3744dabc56b50bec73624aeca02e0def06b03cb287de26836e730659c5d29c"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-win32.whl", hash = "sha256:a999c9f02ff9570c783069f1074f06fe7386ec65b84c983db5aeb8144356a355"},
+    {file = "scikit_learn-1.0.2-cp38-cp38-win_amd64.whl", hash = "sha256:7626a34eabbf370a638f32d1a3ad50526844ba58d63e3ab81ba91e2a7c6d037e"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-macosx_10_13_x86_64.whl", hash = "sha256:a90b60048f9ffdd962d2ad2fb16367a87ac34d76e02550968719eb7b5716fd10"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:7a93c1292799620df90348800d5ac06f3794c1316ca247525fa31169f6d25855"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:eabceab574f471de0b0eb3f2ecf2eee9f10b3106570481d007ed1c84ebf6d6a1"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:55f2f3a8414e14fbee03782f9fe16cca0f141d639d2b1c1a36779fa069e1db57"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:80095a1e4b93bd33261ef03b9bc86d6db649f988ea4dbcf7110d0cded8d7213d"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:fa38a1b9b38ae1fad2863eff5e0d69608567453fdfc850c992e6e47eb764e846"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ff746a69ff2ef25f62b36338c615dd15954ddc3ab8e73530237dd73235e76d62"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-win32.whl", hash = "sha256:e174242caecb11e4abf169342641778f68e1bfaba80cd18acd6bc84286b9a534"},
+    {file = "scikit_learn-1.0.2-cp39-cp39-win_amd64.whl", hash = "sha256:b54a62c6e318ddbfa7d22c383466d38d2ee770ebdb5ddb668d56a099f6eaf75f"},
+]
 scipy = [
     {file = "scipy-1.7.3-1-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:c9e04d7e9b03a8a6ac2045f7c5ef741be86727d8f49c45db45f244bdd2bcff17"},
     {file = "scipy-1.7.3-1-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:b0e0aeb061a1d7dcd2ed59ea57ee56c9b23dd60100825f98238c06ee5cc4467e"},
@@ -3335,6 +3404,10 @@ terminado = [
     {file = "terminado-0.15.0-py3-none-any.whl", hash = "sha256:0d5f126fbfdb5887b25ae7d9d07b0d716b1cc0ccaacc71c1f3c14d228e065197"},
     {file = "terminado-0.15.0.tar.gz", hash = "sha256:ab4eeedccfcc1e6134bfee86106af90852c69d602884ea3a1e8ca6d4486e9bfe"},
 ]
+threadpoolctl = [
+    {file = "threadpoolctl-3.1.0-py3-none-any.whl", hash = "sha256:8b99adda265feb6773280df41eece7b2e6561b772d21ffd52e372f999024907b"},
+    {file = "threadpoolctl-3.1.0.tar.gz", hash = "sha256:a335baacfaa4400ae1f0d8e3a58d6674d2f8828e3716bb2802c44955ad391380"},
+]
 tinycss2 = [
     {file = "tinycss2-1.1.1-py3-none-any.whl", hash = "sha256:fe794ceaadfe3cf3e686b22155d0da5780dd0e273471a51846d0a02bc204fec8"},
     {file = "tinycss2-1.1.1.tar.gz", hash = "sha256:b2e44dd8883c360c35dd0d1b5aad0b610e5156c2cb3b33434634e539ead9d8bf"},
diff --git a/pyproject.toml b/pyproject.toml
index 33a5de1..ac9b309 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -38,6 +38,7 @@ pymc = "4.0.0"
 python = ">=3.7.2, <3.10"
 scipy = "1.7.3"
 torch = "^1.11.0"
+scikit-learn = "1.0.2"
 
 [tool.poetry.dev-dependencies]
 autoflake = "*"
diff --git a/src/kulprit/families/__init__.py b/src/kulprit/families/__init__.py
index 1e431dd..17905f6 100644
--- a/src/kulprit/families/__init__.py
+++ b/src/kulprit/families/__init__.py
@@ -6,10 +6,6 @@
 class BaseFamily(ABC):
     """Base family class."""
 
-    @abstractmethod
-    def solve_analytic(self):  # pragma: no cover
-        pass
-
     @abstractmethod
     def solve_dispersion(self):  # pragma: no cover
         pass
diff --git a/src/kulprit/families/continuous.py b/src/kulprit/families/continuous.py
index 2682851..9287f25 100644
--- a/src/kulprit/families/continuous.py
+++ b/src/kulprit/families/continuous.py
@@ -13,97 +13,8 @@ def __init__(self, data: ModelData) -> None:
         # initialise family object with necessary attributes
         self.data = data
         self.has_dispersion_parameters = True
-        self.has_analytic_solution = True
         self.name = "gaussian"
 
-    def solve_analytic(self, submodel_structure: SubModelStructure) -> torch.tensor:
-        """Analytic solution to the reference model parameter projection.
-
-        Args:
-            submodel_structure (SubModelStructure): The submodel structure object
-
-        Returns
-            tuple: A tuple of the projection solution along with the final loss
-                value of the gradient descent
-        """
-
-        def _analytic_proj(theta_ast: np.float32) -> np.float32:
-            """Analytic solution to the point-wise parameter projection.
-
-            We separate this solution from the primary method to allow for
-            vectorisation of the projection across samples.
-
-            Args:
-                theta_ast (np.float32): The reference model posterior
-                    parameter samples
-
-            Returns:
-                np.float32: Analytic solution to the posterior parameter projection problem
-            """
-
-            f = X_ast @ theta_ast
-            jitter = np.identity(submodel_structure.num_terms) * 1e-5
-            theta_perp = np.linalg.inv(X_perp.T @ X_perp + jitter) @ X_perp.T @ f
-            return theta_perp
-
-        def _kld(theta_ast: np.float32, theta_perp: np.float32) -> np.float32:
-            """Compute the analytic KL divergence between two posterior samples.
-
-            We separate this solution from the primary method to allow for
-            vectorisation of the projection across samples.
-
-            Args:
-                theta_ast (np.float32): The reference model posterior
-                    parameter samples
-                theta_per (np.float32): The submodel posterior parameter samples
-
-            Returns:
-                np.float32: Analytic Kullback-Leibler divergence between the two
-                    samples
-            """
-
-            f = X_ast @ theta_ast
-            f_perp = X_ast @ theta_ast
-            kld = (f - f_perp).T @ (f - f_perp)
-            return kld
-
-        # extract submodel and reference model design matrices
-        X_perp = submodel_structure.X.numpy()
-        X_ast = self.data.structure.X.numpy()
-
-        # extract reference model posterior parameter samples
-        theta_ast = (
-            self.data.idata.posterior.stack(samples=("chain", "draw"))[
-                self.data.structure.term_names
-            ]
-            .to_array()
-            .transpose(*("samples", "variable"))
-            .values
-        )
-
-        # vectorise the analytic solution function
-        vec_analytic_proj = np.vectorize(
-            _analytic_proj,
-            signature="(n)->(m)",
-            doc="Vectorised `_analytic_proj` method",
-        )
-
-        # project the reference model posterior parameter samples
-        theta_perp = vec_analytic_proj(theta_ast)
-
-        # compute the Kullback-Leibler divergence between projection and truth
-        vec_kld = np.vectorize(
-            _kld,
-            signature="(n),(m)->()",
-            doc="Vectorised `_kld` method",
-        )
-        total_kld = torch.from_numpy(vec_kld(theta_ast, theta_perp)).float()
-        loss = torch.mean(total_kld).item()
-
-        # convert the projected parameters to a tensor and return
-        theta_perp = torch.from_numpy(theta_perp).float()
-        return theta_perp, loss
-
     def solve_dispersion(self, theta_perp: torch.tensor, X_perp: torch.tensor):
         """Analytic projection of the model dispersion parameters.
 
diff --git a/src/kulprit/families/family.py b/src/kulprit/families/family.py
index d6059b9..d9d5072 100644
--- a/src/kulprit/families/family.py
+++ b/src/kulprit/families/family.py
@@ -41,25 +41,6 @@ def factory_method(self) -> BaseFamily:
         family_class = self.family_dict[self.family_name]
         return family_class(self.data)
 
-    def solve_analytic(self, submodel_structure: SubModelStructure) -> torch.tensor:
-        """Analytic solution to the reference model parameter projection.
-
-        Args:
-            submodel_structure (SubModelStructure): The submodel structure object
-
-        Returns
-            tuple: A tuple of the projection solution along with the final loss
-                value of the gradient descent
-        """
-
-        # test whether or not the family has an analytic solution
-        if not self.family.has_analytic_solution:  # pragma: no cover
-            raise UserWarning(f"Cannot project {self.family.name} analytically.")
-
-        # compute the solution and return
-        solution = self.family.solve_analytic(submodel_structure=submodel_structure)
-        return solution
-
     def solve_dispersion(self, theta_perp: torch.tensor, X_perp: torch.tensor):
         """Analytic projection of the model dispersion parameters.
 
diff --git a/src/kulprit/projection/projector.py b/src/kulprit/projection/projector.py
index 6df7042..eb550be 100644
--- a/src/kulprit/projection/projector.py
+++ b/src/kulprit/projection/projector.py
@@ -6,15 +6,14 @@
 from kulprit.data.data import ModelData
 from kulprit.data.submodel import SubModelStructure, SubModelInferenceData
 from kulprit.families.family import Family
-from kulprit.projection.solvers.solver import Solver
-from kulprit.search.path import SearchPath
+from kulprit.projection.solver import Solver
 
 
 class Projector:
     def __init__(
         self,
         data: ModelData,
-        path: Optional[SearchPath] = None,
+        path: Optional[dict] = None,
         num_iters: Optional[int] = 200,
         learning_rate: Optional[float] = 0.01,
     ) -> None:
@@ -48,7 +47,6 @@ def __init__(
     def project(
         self,
         terms: Union[List[str], int],
-        method: Literal["analytic", "gradient"],
     ) -> ModelData:
         """Wrapper function for projection method.
 
@@ -57,9 +55,6 @@ def project(
                 the names of the parameters to include the submodel, or the
                 number of parameters to include in the submodel, **not**
                 including the intercept term
-            method (str): The projection method to employ, either "analytic" to
-                use the hard-coded solutions the optimisation problem, or
-                "gradient" to employ gradient descent methods
 
         Returns:
             kulprit.data.ModelData: Projected submodel ``ModelData`` object
@@ -74,12 +69,12 @@ def project(
                     + " the reference model."
                 )
             # perform projection
-            return self.project_names(term_names=terms, method=method)
+            return self.project_names(term_names=terms)
 
         # project a number of terms
         else:
             # test `model_size` input
-            if self.path is None or terms not in list(self.path.k_submodel.keys()):
+            if self.path is None or terms not in list(self.path.keys()):
                 raise UserWarning(
                     "In order to project onto an integer number of terms, please "
                     + "first complete a parameter search."
@@ -91,7 +86,6 @@ def project(
     def project_names(
         self,
         term_names: List[List[str]],
-        method: Literal["analytic", "gradient"],
     ) -> ModelData:
         """Primary projection method for GLM reference model.
 
@@ -103,7 +97,6 @@ def project_names(
         Args:
             term_names (List[str]): The names of parameters to project onto the
                 submodel, **not** including the intercept term
-            method (str): The method to use in projection
 
         Returns:
             kulprit.data.ModelData: Projected submodel ``ModelData`` object
@@ -117,7 +110,6 @@ def project_names(
         self.solver = Solver(
             data=self.data,
             family=self.family,
-            method=method,
             num_iters=self.num_iters,
             learning_rate=self.learning_rate,
         )
diff --git a/src/kulprit/projection/solvers/gradient.py b/src/kulprit/projection/solver.py
similarity index 69%
rename from src/kulprit/projection/solvers/gradient.py
rename to src/kulprit/projection/solver.py
index f0c0a97..06c734f 100644
--- a/src/kulprit/projection/solvers/gradient.py
+++ b/src/kulprit/projection/solver.py
@@ -1,4 +1,4 @@
-"""Gradient descent solver module."""
+"""Optimisation solver module."""
 
 import torch
 
@@ -6,12 +6,9 @@
 from kulprit.families.family import Family
 from kulprit.projection.architecture import Architecture
 from kulprit.projection.losses.kld import KullbackLeiblerLoss
-from kulprit.projection.solvers import BaseSolver
 
 
-class GradientDescentSolver(BaseSolver):
-    """Gradient descent projection optimisation solver."""
-
+class Solver:
     def __init__(
         self,
         data: ModelData,
@@ -19,6 +16,18 @@ def __init__(
         num_iters: int = 200,
         learning_rate: float = 0.01,
     ):
+        """Initialise solver object.
+
+        Args:
+            data (ModelData): The data object containing the model data of
+                the reference model
+            family (Family): The family object of the reference model
+            num_iters (int, optional): The number of iterations to run the
+                optimisation for, defaults to 200
+            learning_rate (float, optional): The learning rate for the optimisation
+                algorithm, defaults to 0.01
+        """
+
         # log reference model data and family
         self.data = data
         self.family = family
@@ -88,3 +97,22 @@ def optimise(self, X_perp, y_ast):
         theta_perp = list(solver.parameters())[0].data
         final_loss = loss.item()
         return theta_perp, final_loss
+
+    def solve_dispersion(
+        self, theta_perp: torch.tensor, X_perp: torch.tensor
+    ) -> torch.tensor:
+        """Analytic projection of the model dispersion parameters.
+
+        Args:
+            theta_perp (torch.tensor): A PyTorch tensor of the restricted
+                parameter draws
+            X_perp (np.ndarray): The design matrix of the restricted model we
+                are projecting onto
+
+        Returns:
+            torch.tensor: The restricted projections of the dispersion parameters
+        """
+
+        # compute the solution and return
+        solution = self.family.solve_dispersion(theta_perp=theta_perp, X_perp=X_perp)
+        return solution
diff --git a/src/kulprit/projection/solvers/__init__.py b/src/kulprit/projection/solvers/__init__.py
deleted file mode 100644
index 6e54888..0000000
--- a/src/kulprit/projection/solvers/__init__.py
+++ /dev/null
@@ -1,11 +0,0 @@
-"""Base optimisation solvers module."""
-
-from abc import ABC, abstractmethod
-
-
-class BaseSolver(ABC):
-    """Base solver class."""
-
-    @abstractmethod
-    def solve(self):  # pragma: no cover
-        pass
diff --git a/src/kulprit/projection/solvers/analytic.py b/src/kulprit/projection/solvers/analytic.py
deleted file mode 100644
index 3beffc0..0000000
--- a/src/kulprit/projection/solvers/analytic.py
+++ /dev/null
@@ -1,38 +0,0 @@
-"""Analytic solver module."""
-
-import torch
-
-from kulprit.data.data import ModelData
-from kulprit.families.family import Family
-from kulprit.projection.solvers import BaseSolver
-
-
-class AnalyticSolver(BaseSolver):
-    """Analytic projection optimisation solver."""
-
-    def __init__(
-        self,
-        data: ModelData,
-        family: Family,
-        num_iters: int = None,
-        learning_rate: float = None,
-    ):
-        # log reference model data and family
-        self.data = data
-        self.family = family
-
-    def solve(self, submodel_structure: torch.tensor) -> tuple:
-        """Perform projection by gradient descent.
-
-        Args:
-            submodel_structure (SubModelStructure): The structure of the submodel
-                being projected onto
-
-        Returns
-            tuple: A tuple of the projection solution along with the final loss
-                value of the gradient descent
-        """
-
-        # compute and return solution
-        solution = self.family.solve_analytic(submodel_structure=submodel_structure)
-        return solution
diff --git a/src/kulprit/projection/solvers/solver.py b/src/kulprit/projection/solvers/solver.py
deleted file mode 100644
index 98b9886..0000000
--- a/src/kulprit/projection/solvers/solver.py
+++ /dev/null
@@ -1,105 +0,0 @@
-"""Base solver module."""
-
-from typing import Optional
-from typing_extensions import Literal
-
-import torch
-
-from kulprit.data.data import ModelData
-from kulprit.data.submodel import SubModelStructure
-from kulprit.families.family import Family
-from kulprit.projection.solvers import BaseSolver
-from kulprit.projection.solvers.analytic import AnalyticSolver
-from kulprit.projection.solvers.gradient import GradientDescentSolver
-
-
-class Solver:
-    def __init__(
-        self,
-        data: ModelData,
-        family: Family,
-        method: Literal["analytic", "gradient"],
-        num_iters: Optional[int] = 200,
-        learning_rate: Optional[float] = 0.01,
-    ) -> None:
-        """Solver factory class.
-
-        Args:
-            data (ModelData): The reference model ModelData object
-            method (str): The method of projection to use in the procedure
-            num_iters (int): The number of iterations to run gradient descent
-            learning_rate (float): The learning rate of the gradient descent
-                optimiser
-        """
-
-        # log method, data, and family of the problem
-        self.data = data
-        self.method = method
-        self.family = family
-
-        # store gradient descent parameters if provided
-        self.num_iters = num_iters
-        self.learning_rate = learning_rate
-
-        # define all available KL divergence loss classes
-        self.method_dict = {
-            "analytic": AnalyticSolver,
-            "gradient": GradientDescentSolver,
-        }
-
-        # test valid solution method
-        if self.method not in self.method_dict:
-            raise UserWarning(
-                "Please either solve the projection analytically or with "
-                + "gradient descent."
-            )
-
-    def factory_method(self) -> BaseSolver:
-        """Choose the appropriate family class given the model."""
-
-        # return appropriate solver class given method
-        solver_class = self.method_dict[self.method]
-        return solver_class(
-            data=self.data,
-            family=self.family,
-            num_iters=self.num_iters,
-            learning_rate=self.learning_rate,
-        )
-
-    def solve(self, submodel_structure: SubModelStructure) -> tuple:
-        """Solve the optimisation (projection) problem central to the procedure.
-
-        Args:
-            submodel_structure (SubModelStructure): The structure of the submodel
-                being projected onto
-
-        Returns
-            tuple: A tuple of the projection solution along with the final loss
-                value of the gradient descent
-        """
-
-        # initialise family class
-        solver_class = self.factory_method()
-
-        # compute the solution and return
-        solution, loss = solver_class.solve(submodel_structure=submodel_structure)
-        return solution, loss
-
-    def solve_dispersion(
-        self, theta_perp: torch.tensor, X_perp: torch.tensor
-    ) -> torch.tensor:
-        """Analytic projection of the model dispersion parameters.
-
-        Args:
-            theta_perp (torch.tensor): A PyTorch tensor of the restricted
-                parameter draws
-            X_perp (np.ndarray): The design matrix of the restricted model we
-                are projecting onto
-
-        Returns:
-            torch.tensor: The restricted projections of the dispersion parameters
-        """
-
-        # compute the solution and return
-        solution = self.family.solve_dispersion(theta_perp=theta_perp, X_perp=X_perp)
-        return solution
diff --git a/src/kulprit/reference.py b/src/kulprit/reference.py
index 66940c5..4c40092 100644
--- a/src/kulprit/reference.py
+++ b/src/kulprit/reference.py
@@ -10,7 +10,6 @@
 
 from kulprit.data import ModelData, ModelStructure
 from kulprit.projection.projector import Projector
-from kulprit.search.path import SearchPath
 from kulprit.search.searcher import Searcher
 
 
@@ -67,13 +66,12 @@ def __init__(
             num_iters=num_iters,
             learning_rate=learning_rate,
         )
-        self.searcher = Searcher(self.data, self.projector)
+        self.searcher = Searcher(self.projector)
         self.path = None
 
     def project(
         self,
         terms: Union[List[str], int],
-        method: Literal["analytic", "gradient"] = "analytic",
     ) -> ModelData:
         """Projection the reference model onto a variable subset.
 
@@ -82,23 +80,20 @@ def project(
                 the names of the parameters to include the submodel, or the
                 number of parameters to include in the submodel, **not**
                 including the intercept term
-            method (str): The projection method to employ, either "analytic" to
-                use the hard-coded solutions the optimisation problem, or
-                "gradient" to employ gradient descent methods
 
         Returns:
             kulprit.data.ModelData: Projected submodel ``ModelData`` object
         """
 
         # project the reference model onto a subset of covariates
-        sub_model = self.projector.project(terms=terms, method=method)
+        sub_model = self.projector.project(terms=terms)
         return sub_model
 
     def search(
         self,
         max_terms: Optional[int] = None,
-        method: Literal["analytic", "gradient"] = "analytic",
-    ) -> SearchPath:
+        method: Literal["forward", "l1"] = "forward",
+    ) -> dict:
         """Model search method through parameter space.
 
         If ``max_terms`` is not provided, then the search path runs from the
@@ -108,12 +103,13 @@ def search(
         Args:
             max_terms (int): The number of parameters of the largest submodel in
                 the search path, **not** including the intercept term
-            method (str): The projection method to employ, either "analytic" to
-                use the hard-coded solutions the optimisation problem, or
-                "gradient" to employ gradient descent methods
+            method (str): The search method to employ, either "forward" to
+                employ a forward search heuristic through the space, or "l1" to
+                use the L1-regularised search path
 
         Returns:
-            kulprit.search.SearchPath: The model selection procedure search path
+            dict: The model selection procedure search path, containing the
+                submodels along the search path, keyed by their model size
         """
 
         # set default `max_terms` value
@@ -129,7 +125,6 @@ def search(
             )
 
         self.path = self.searcher.search(max_terms=max_terms, method=method)
-        self.projector.path = self.path  # feed path result through to the projector
         return self.path
 
     def loo_compare(
diff --git a/src/kulprit/search/__init__.py b/src/kulprit/search/__init__.py
index e83df38..17af863 100644
--- a/src/kulprit/search/__init__.py
+++ b/src/kulprit/search/__init__.py
@@ -1,1 +1,12 @@
 """Base search module."""
+
+
+from abc import ABC, abstractmethod
+
+
+class SearchPath(ABC):
+    """Base search module."""
+
+    @abstractmethod
+    def search(self, max_terms: int) -> None:
+        """Initialise search path class."""
diff --git a/src/kulprit/search/forward.py b/src/kulprit/search/forward.py
new file mode 100644
index 0000000..b294861
--- /dev/null
+++ b/src/kulprit/search/forward.py
@@ -0,0 +1,128 @@
+"""Forward search module."""
+
+from typing import List
+
+import pandas as pd
+from kulprit.data.data import ModelData
+from kulprit.projection.projector import Projector
+from kulprit.search import SearchPath
+
+
+class ForwardSearchPath(SearchPath):
+    def __init__(self, projector: Projector) -> None:
+        """Initialise search path class."""
+
+        # log the names of the terms in the reference model
+        self.ref_terms = projector.data.structure.term_names
+
+        # log the projector object
+        self.projector = projector
+
+        # initialise search
+        self.k_term_idx = {}
+        self.k_term_names = {}
+        self.k_submodel = {}
+        self.k_dist = {}
+
+    def __repr__(self) -> str:
+        """String representation of the search path."""
+
+        path_dict = {
+            k: [submodel.structure.term_names, submodel.dist_to_ref_model]
+            for k, submodel in self.k_submodel.items()
+        }
+        df = pd.DataFrame.from_dict(
+            path_dict, orient="index", columns=["Terms", "Distance from reference model"]
+        )
+        df.index.name = "Model size"
+        string = df.to_string()
+        return string
+
+    def add_submodel(
+        self,
+        k: int,
+        k_term_names: list,
+        k_submodel: ModelData,
+        k_dist: float,
+    ) -> None:
+        """Update search path with new submodel."""
+
+        self.k_term_names[k] = k_term_names
+        self.k_submodel[k] = k_submodel
+        self.k_dist[k] = k_dist
+
+    def get_candidates(self, k: int) -> List[List]:
+        """Method for extracting a list of all candidate submodels.
+
+        Args:
+            k (int): The number of terms in the previous submodel, from which we
+                wish to find all possible candidate submodels
+
+        Returns:
+            List[List]: A list of lists, each containing the terms of all candidate
+                submodels
+        """
+
+        prev_subset = self.k_term_names[k]
+        candidate_additions = list(set(self.ref_terms).difference(prev_subset))
+        candidates = [prev_subset + [addition] for addition in candidate_additions]
+        return candidates
+
+    def search(self, max_terms: int) -> dict:
+        """Forward search through the parameter space.
+
+        Args:
+            max_terms (int): Number of terms to perform the forward search
+                up to.
+
+        Returns:
+            dict: A dictionary of submodels, keyed by the number of terms in the
+                submodel.
+        """
+
+        # initial intercept-only subset
+        k = 0
+        k_term_names = []
+        k_submodel = self.projector.project(terms=k_term_names)
+        k_dist = k_submodel.dist_to_ref_model
+
+        # add submodel to search path
+        self.add_submodel(
+            k=k,
+            k_term_names=k_term_names,
+            k_submodel=k_submodel,
+            k_dist=k_dist,
+        )
+
+        # perform forward search through parameter space
+        while k < max_terms:
+            # get list of candidate submodels, project onto them, and compute
+            # their distances
+            k_candidates = self.get_candidates(k=k)
+            k_projections = [
+                self.projector.project(terms=candidate) for candidate in k_candidates
+            ]
+
+            # identify the best candidate by distance from reference model
+            best_submodel = min(
+                k_projections, key=lambda projection: projection.sort_index
+            )
+            best_dist = best_submodel.sort_index
+
+            # retrieve the best candidate's term names and indices
+            k_term_names = best_submodel.structure.common_terms
+
+            # increment number of parameters
+            k = len(k_term_names)
+
+            # add best candidate to search path
+            self.add_submodel(
+                k=k,
+                k_term_names=k_term_names,
+                k_submodel=best_submodel,
+                k_dist=best_dist,
+            )
+
+        # toggle indicator variable and return search path
+        self.search_completed = True
+        return self.k_submodel
diff --git a/src/kulprit/search/l1.py b/src/kulprit/search/l1.py
new file mode 100644
index 0000000..a3b04a3
--- /dev/null
+++ b/src/kulprit/search/l1.py
@@ -0,0 +1,123 @@
+"""L1 search path module."""
+
+import numpy as np
+import pandas as pd
+from sklearn.linear_model import lasso_path
+
+from kulprit.data.data import ModelData
+from kulprit.projection.projector import Projector
+from kulprit.search import SearchPath
+
+
+class L1SearchPath(SearchPath):
+    """L1 search path class."""
+
+    def __init__(self, projector: Projector) -> None:
+        """Initialise L1 search path class."""
+
+        # log the projector object
+        self.projector = projector
+
+        # log the model data object of the reference model
+        self.data = projector.data
+
+        # initialise search
+        self.k_term_names = {}
+        self.k_submodel = {}
+        self.k_dist = {}
+
+        self.search_completed = True
+
+    def __repr__(self) -> str:
+        """String representation of the search path."""
+
+        path_dict = {
+            k: [submodel.structure.term_names, submodel.dist_to_ref_model]
+            for k, submodel in self.k_submodel.items()
+        }
+        df = pd.DataFrame.from_dict(
+            path_dict, orient="index", columns=["Terms", "Distance from reference model"]
+        )
+        df.index.name = "Model size"
+        string = df.to_string()
+        return string
+
+    def first_non_zero_idx(self, arr):
+        """Find the index of the first non-zero element in each row of a matrix.
+
+        Args:
+            arr (np.ndarray): A matrix.
+
+        Returns:
+            dict: Dictionary keyed by the row number where each value is the index
+                of the first non-zero element in that row."""
+
+        # initialise dictionary of indices
+        idx_dict = {}
+
+        # loop through each row and find first non-zero element
+        for i, j in zip(*np.where(arr != 0)):
+            if i in idx_dict:
+                continue
+            idx_dict[i] = j
+
+        # identify which keys are missing and set their values to infinity
+        if len(idx_dict) < arr.shape[0]:
+            missing_keys = set(range(arr.shape[0])) - set(idx_dict.keys())
+            for key in missing_keys:
+                idx_dict[key] = np.inf
+
+        return idx_dict
+
+    def compute_path(self) -> None:
+        """Compute the L1 search path.
+
+        We compute the L1 search path for a given data object in the latent space
+        and return the coefficients for each model size. This use of the latent
+        space is a bit of a hack, but Catalina et al. (2021) show that this
+        space remains informative in terms of model selection and results in
+        faster computation. We use ``sklearn`` to compute the L1 path.
+
+        Returns:
+            np.ndarray: The coefficients for each model size
+        """
+
+        # extract reference model data and latent predictor
+        X = np.column_stack(
+            [
+                self.data.structure.design.common[term]
+                for term in self.data.structure.common_terms
+            ]
+        )
+        eta = self.data.structure.link.link(self.data.structure.y.numpy())
+
+        # compute L1 path in the latent space
+        _, coef_path, _ = lasso_path(X, eta)
+        cov_order = self.first_non_zero_idx(coef_path)
+
+        return cov_order
+
+    def search(self, max_terms: int) -> None:
+        """Perform L1 search through the parameter."""
+
+        # compute L1 path for each model size
+        coef_path = self.compute_path()
+
+        # sort the covariates according to their L1 ordering
+        cov_lasso = {
+            k: v for k, v in sorted(coef_path.items(), key=lambda item: item[1])
+        }
+        sorted_covs = [self.data.structure.common_terms[k] for k in cov_lasso]
+        sorted_covs = [["Intercept"] + sorted_covs[:i] for i in range(max_terms + 1)]
+
+        # produce submodels for each model size
+        self.k_term_names = {len(terms) - 1: terms for terms in sorted_covs}
+
+        # project the reference model on each of the submodels
+        for k, term_names in self.k_term_names.items():
+            self.k_submodel[k] = self.projector.project(term_names)
+            self.k_dist[k] = self.k_submodel[k].dist_to_ref_model
+
+        # toggle indicator variable and return search path
+        self.search_completed = True
+        return self.k_submodel
diff --git a/src/kulprit/search/path.py b/src/kulprit/search/path.py
deleted file mode 100644
index 882420e..0000000
--- a/src/kulprit/search/path.py
+++ /dev/null
@@ -1,72 +0,0 @@
-"""Search path module."""
-
-from typing import List
-
-import pandas as pd
-
-from kulprit.data.data import ModelData
-
-
-class SearchPath:
-    """Search path object for model selection procedure."""
-
-    def __init__(self, ref_terms: list) -> None:
-        """Initialise search path class."""
-
-        # set limit to number of terms in search path
-        self.ref_terms = ref_terms
-
-        # initialise search
-        self.k_term_idx = {}
-        self.k_term_names = {}
-        self.k_submodel = {}
-        self.k_dist = {}
-
-    def __getitem__(self, k: int) -> ModelData:
-        """Return the submodel in the search path with k terms."""
-
-        return self.k_submodel[k]
-
-    def __repr__(self) -> str:
-        """String representation of the search path."""
-
-        path_dict = {
-            k: [submodel.structure.term_names, submodel.dist_to_ref_model]
-            for k, submodel in self.k_submodel.items()
-        }
-        df = pd.DataFrame.from_dict(
-            path_dict, orient="index", columns=["Terms", "Distance from reference model"]
-        )
-        df.index.name = "Model size"
-        string = df.to_string()
-        return string
-
-    def add_submodel(
-        self,
-        k: int,
-        k_term_names: list,
-        k_submodel: ModelData,
-        k_dist: float,
-    ) -> None:
-        """Update search path with new submodel."""
-
-        self.k_term_names[k] = k_term_names
-        self.k_submodel[k] = k_submodel
-        self.k_dist[k] = k_dist
-
-    def get_candidates(self, k: int) -> List[List]:
-        """Method for extracting a list of all candidate submodels.
-
-        Args:
-            k (int): The number of terms in the previous submodel, from which we
-                wish to find all possible candidate submodels
-
-        Returns:
-            List[List]: A list of lists, each containing the terms of all candidate
-                submodels
-        """
-
-        prev_subset = self.k_term_names[k]
-        candidate_additions = list(set(self.ref_terms).difference(prev_subset))
-        candidates = [prev_subset + [addition] for addition in candidate_additions]
-        return candidates
diff --git a/src/kulprit/search/searcher.py b/src/kulprit/search/searcher.py
index 7763beb..753c499 100644
--- a/src/kulprit/search/searcher.py
+++ b/src/kulprit/search/searcher.py
@@ -6,24 +6,30 @@
 import pandas as pd
 import arviz as az
 
-from kulprit.data.data import ModelData
 from kulprit.projection.projector import Projector
-from kulprit.search.path import SearchPath
+from kulprit.search.forward import ForwardSearchPath
+from kulprit.search.l1 import L1SearchPath
 
 
 class Searcher:
-    def __init__(self, data: ModelData, projector: Projector) -> None:
-        """Initialise forward search class."""
+    def __init__(self, projector: Projector) -> None:
+        """Initialise forward search class.
 
-        # initialise reference model and set of reference model term names
-        self.data = data
-        self.term_names = data.structure.common_terms
+        Args:
+            projector (Projector): Projector object.
+
+        Raises:
+            UserWarning: If method is not "forward" or "l1".
+        """
 
         # initialise projector for the search procedure
         self.projector = projector
 
-        # initialise search path
-        self.path = SearchPath(ref_terms=self.term_names)
+        # define all available search heuristics
+        self.method_dict = {
+            "forward": ForwardSearchPath,
+            "l1": L1SearchPath,
+        }
 
         # indicator variable tracking whether or not a search has been run
         self.search_completed = False
@@ -34,8 +40,8 @@ def __repr__(self) -> str:
         return self.path.__repr__()
 
     def search(
-        self, max_terms: int, method: Literal["analytic", "gradient"]
-    ) -> SearchPath:
+        self, max_terms: int, method: Literal["forward", "l1"] = "forward"
+    ) -> dict:
         """Primary search method of the procedure.
 
         Performs forward search through the parameter space.
@@ -43,58 +49,31 @@ def search(
         Args:
             max_terms (int): Number of terms to perform the forward search
                 up to.
-            method (str): The projection method to employ, either "analytic" to
-                use the hard-coded solutions the optimisation problem, or
-                "gradient" to employ gradient descent methods
+            method (str): Method to use for search.
+
+        Returns:
+            dict: A dictionary of submodels, keyed by the number of terms in the
+                submodel.
         """
 
-        # initial intercept-only subset
-        k = 0
-        k_term_names = []
-        k_submodel = self.projector.project(terms=k_term_names, method=method)
-        k_dist = k_submodel.dist_to_ref_model
-
-        # add submodel to search path
-        self.path.add_submodel(
-            k=k,
-            k_term_names=k_term_names,
-            k_submodel=k_submodel,
-            k_dist=k_dist,
-        )
+        # test valid solution method
+        if method not in self.method_dict:
+            raise UserWarning("Please either select either forward search or L1 search.")
+
+        # initialise search path
+        self.path = self.method_dict[method](self.projector)
 
-        # perform forward search through parameter space
-        while k < max_terms:
-            # get list of candidate submodels, project onto them, and compute
-            # their distances
-            k_candidates = self.path.get_candidates(k=k)
-            k_projections = [
-                self.projector.project(terms=candidate, method=method)
-                for candidate in k_candidates
-            ]
-
-            # identify the best candidate by distance from reference model
-            best_submodel = min(
-                k_projections, key=lambda projection: projection.sort_index
-            )
-            best_dist = best_submodel.sort_index
-
-            # retrieve the best candidate's term names and indices
-            k_term_names = best_submodel.structure.common_terms
-
-            # increment number of parameters
-            k = len(k_term_names)
-
-            # add best candidate to search path
-            self.path.add_submodel(
-                k=k,
-                k_term_names=k_term_names,
-                k_submodel=best_submodel,
-                k_dist=best_dist,
-            )
+        # perform the search according to the chosen heuristic
+        k_submodels = self.path.search(max_terms=max_terms)
 
         # toggle indicator variable and return search path
         self.search_completed = True
-        return self.path
+
+        # feed path result through to the projector
+        self.projector.path = k_submodels
+
+        # return the final search path
+        return k_submodels
 
     def loo_compare(
         self,


Test patch:
diff --git a/tests/test_project.py b/tests/test_project.py
index 6ab44fc..f46ef1a 100644
--- a/tests/test_project.py
+++ b/tests/test_project.py
@@ -45,7 +45,7 @@ def test_architecture_forward(self, ref_model):
             ref_model.data.structure.num_obs,
         )
 
-    def test_analytic_project(self, ref_model):
+    def test_project(self, ref_model):
         """Test that the analytic projection method works."""
 
         # project the reference model to some parameter subset
@@ -55,16 +55,6 @@ def test_analytic_project(self, ref_model):
         assert sub_model.structure.num_terms == 2
         assert sub_model.structure.model_size == 1
 
-    def test_gradient_project(self, ref_model):
-        """Test that the gradient projection method works."""
-
-        # project the reference model to some parameter subset
-        sub_model = ref_model.project(terms=["x"], method="gradient")
-
-        assert sub_model.structure.X.shape == (ref_model.data.structure.num_obs, 2)
-        assert sub_model.structure.num_terms == 2
-        assert sub_model.structure.model_size == 1
-
     def test_projected_idata_dims(self, ref_model, bambi_model_idata):
         """Test that the dimensions of the projected inference data are correct."""
 
diff --git a/tests/test_search.py b/tests/test_search.py
index cdda950..ce25fc6 100644
--- a/tests/test_search.py
+++ b/tests/test_search.py
@@ -1,5 +1,8 @@
 import bambi as bmb
 import kulprit as kpt
+from kulprit.search.l1 import L1SearchPath
+
+import numpy as np
 
 import copy
 import pytest
@@ -15,7 +18,32 @@ def test_forward(self, ref_model):
 
         ref_model_copy = copy.copy(ref_model)
         ref_model_copy.search()
-        assert list(ref_model_copy.path.k_submodel.keys()) == [0, 1, 2]
+        assert list(ref_model_copy.path.keys()) == [0, 1, 2]
+
+    def test_l1(self, ref_model):
+        """Test that L1 search gives expected result."""
+
+        ref_model_copy = copy.copy(ref_model)
+        ref_model_copy.search(method="l1")
+        assert list(ref_model_copy.path.keys()) == [0, 1, 2]
+
+    def test_l1_utils(self, ref_model):
+        """Test that L1 utility methods return expected result."""
+
+        copy.copy(ref_model)
+        proj = ref_model.projector
+        searcher = L1SearchPath(proj)
+        arr = np.array(
+            [[0.0, 1.0, 2.0, 3.0], [0.0, 0.0, -1.0, 2.0], [0.0, 0.0, 0.0, 0.0]]
+        )
+        assert searcher.first_non_zero_idx(arr) == {0: 1, 1: 2, 2: np.inf}
+
+    def test_bad_search_method(self, ref_model):
+        """Test that an error is raised when an invalid search method is used."""
+
+        with pytest.raises(UserWarning):
+            ref_model_copy = copy.copy(ref_model)
+            ref_model_copy.search(method="bad_method")
 
     def test_search_too_many_terms(self, ref_model):
         """Test than an error is raise when too many terms are used."""
@@ -45,14 +73,22 @@ def test_loo_with_no_search_path(self, ref_model):
             ref_model = kpt.ReferenceModel(model=bambi_model, idata=idata)
             ref_model.loo_compare()
 
-    def test_repr(self, ref_model):
-        """Test the string representation of the search path."""
+    def test_forward_repr(self, ref_model):
+        """Test the string representation of the forward search path."""
 
         ref_model_copy = copy.copy(ref_model)
         path = ref_model_copy.search()
         assert type(ref_model_copy.searcher.__repr__()) == str
         assert type(path.__repr__()) == str
 
+    def test_l1_repr(self, ref_model):
+        """Test the string representation of the L1 search path."""
+
+        ref_model_copy = copy.copy(ref_model)
+        path = ref_model_copy.search(method="l1")
+        assert type(ref_model_copy.searcher.__repr__()) == str
+        assert type(path.__repr__()) == str
+
     def test_plotting(self, ref_model):
         """Test the LOO compare plotting method."""
 
diff --git a/tests/test_solver.py b/tests/test_solver.py
deleted file mode 100644
index a8c9dc6..0000000
--- a/tests/test_solver.py
+++ /dev/null
@@ -1,19 +0,0 @@
-from kulprit.data.data import ModelData
-from kulprit.families import BaseFamily
-from kulprit.families.family import Family
-from kulprit.projection.solvers.solver import Solver
-
-import pytest
-
-from tests import KulpritTest
-
-
-class TestSolver(KulpritTest):
-    """Test the solver object methods."""
-
-    def test_wrong_method(self, ref_model):
-        """Test that the correct error is raised when a wrong method is used."""
-
-        with pytest.raises(UserWarning):
-            family = Family(ref_model.data)
-            Solver(ref_model.data, family, method="random")