diff --git a/pymc_experimental/gp/__init__.py b/pymc_experimental/gp/__init__.py
index c08daef8..c8804dd4 100644
--- a/pymc_experimental/gp/__init__.py
+++ b/pymc_experimental/gp/__init__.py
@@ -13,8 +13,4 @@
 #   limitations under the License.
 
 
-from pymc_experimental.gp.latent_approx import (
-    HSGP,
-    KarhunenLoeveExpansion,
-    ProjectedProcess,
-)
+from pymc_experimental.gp.latent_approx import KarhunenLoeveExpansion, ProjectedProcess
diff --git a/pymc_experimental/gp/latent_approx.py b/pymc_experimental/gp/latent_approx.py
index 899d4545..abb9fc02 100644
--- a/pymc_experimental/gp/latent_approx.py
+++ b/pymc_experimental/gp/latent_approx.py
@@ -75,108 +75,6 @@ def conditional(self, name, Xnew, jitter=1e-6, **kwargs):
         return pm.MvNormal(name, mu=mu, chol=chol)
 
 
-class HSGP(pm.gp.Latent):
-    ## inputs: M, c
-
-    def __init__(
-        self, n_basis, c=3 / 2, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)
-    ):
-        ## TODO: specify either c or L
-        self.M = n_basis
-        self.c = c
-        super().__init__(mean_func=mean_func, cov_func=cov_func)
-
-    def _validate_cov_func(self, cov_func):
-        ## TODO: actually validate it.  Right now this fails unless cov func is exactly
-        # in the form eta**2 * pm.gp.cov.Matern12(...) and will error otherwise.
-        cov, scaling_factor = cov_func.factor_list
-        return scaling_factor, cov.ls, cov.spectral_density
-
-    def prior(self, name, X, **kwargs):
-        f, Phi, L, spd, beta, Xmu, Xsd = self._build_prior(name, X, **kwargs)
-        self.X, self.f = X, f
-        self.Phi, self.L, self.spd, self.beta = Phi, L, spd, beta
-        self.Xmu, self.Xsd = Xmu, Xsd
-        return f
-
-    def _generate_basis(self, X, L):
-        indices = pt.arange(1, self.M + 1)
-        m1 = (np.pi / (2.0 * L)) * pt.tile(L + X, self.M)
-        m2 = pt.diag(indices)
-        Phi = pt.sin(m1 @ m2) / pt.sqrt(L)
-        omega = (np.pi * indices) / (2.0 * L)
-        return Phi, omega
-
-    def _build_prior(self, name, X, **kwargs):
-        n_obs = np.shape(X)[0]
-
-        # standardize input scale
-        X = pt.as_tensor_variable(X)
-        Xmu = pt.mean(X, axis=0)
-        Xsd = pt.std(X, axis=0)
-        Xz = (X - Xmu) / Xsd
-
-        # define L using Xz and c
-        La = pt.abs(pt.min(Xz))  # .eval()?
-        Lb = pt.max(Xz)
-        L = self.c * pt.max([La, Lb])
-
-        # make basis and omega, spectral density
-        Phi, omega = self._generate_basis(Xz, L)
-        scale, ls, spectral_density = self._validate_cov_func(self.cov_func)
-        spd = scale * spectral_density(omega, ls / Xsd).flatten()
-
-        beta = pm.Normal(f"{name}_coeffs_", size=self.M)
-        f = pm.Deterministic(name, self.mean_func(X) + pt.dot(Phi * pt.sqrt(spd), beta))
-        return f, Phi, L, spd, beta, Xmu, Xsd
-
-    def _build_conditional(self, Xnew, Xmu, Xsd, L, beta):
-        Xnewz = (Xnew - Xmu) / Xsd
-        Phi, omega = self._generate_basis(Xnewz, L)
-        scale, ls, spectral_density = self._validate_cov_func(self.cov_func)
-        spd = scale * spectral_density(omega, ls / Xsd).flatten()
-        return self.mean_func(Xnew) + pt.dot(Phi * pt.sqrt(spd), beta)
-
-    def conditional(self, name, Xnew):
-        # warn about extrapolation
-        fnew = self._build_conditional(Xnew, self.Xmu, self.Xsd, self.L, self.beta)
-        return pm.Deterministic(name, fnew)
-
-
-class ExpQuad(pm.gp.cov.ExpQuad):
-    @staticmethod
-    def spectral_density(omega, ls):
-        # univariate spectral denisty, implement multi
-        return pt.sqrt(2 * np.pi) * ls * pt.exp(-0.5 * ls**2 * omega**2)
-
-
-class Matern52(pm.gp.cov.Matern52):
-    @staticmethod
-    def spectral_density(omega, ls):
-        # univariate spectral denisty, implement multi
-        # https://arxiv.org/pdf/1611.06740.pdf
-        lam = pt.sqrt(5) * (1.0 / ls)
-        return (16.0 / 3.0) * lam**5 * (1.0 / (lam**2 + omega**2) ** 3)
-
-
-class Matern32(pm.gp.cov.Matern32):
-    @staticmethod
-    def spectral_density(omega, ls):
-        # univariate spectral denisty, implement multi
-        # https://arxiv.org/pdf/1611.06740.pdf
-        lam = np.sqrt(3.0) * (1.0 / ls)
-        return 4.0 * lam**3 * (1.0 / pt.square(lam**2 + omega**2))
-
-
-class Matern12(pm.gp.cov.Matern12):
-    @staticmethod
-    def spectral_density(omega, ls):
-        # univariate spectral denisty, implement multi
-        # https://arxiv.org/pdf/1611.06740.pdf
-        lam = 1.0 / ls
-        return 2.0 * lam * (1.0 / (lam**2 + omega**2))
-
-
 class KarhunenLoeveExpansion(pm.gp.Latent):
     def __init__(
         self,