diff --git a/Cargo.toml b/Cargo.toml
index 46d60f0..43e7fe4 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -16,5 +16,6 @@ statrs = "0.13"
 
 [dev-dependencies]
 approx = "0.4"
+blis-src = "0.2"
 intel-mkl-src = "0.5"
 time = "0.2"
diff --git a/src/kernel.rs b/src/kernel.rs
index 854f74a..501bfde 100644
--- a/src/kernel.rs
+++ b/src/kernel.rs
@@ -10,8 +10,35 @@ pub use exponential::Exponential;
 pub use matern32::Matern32;
 pub use matern52::Matern52;
 
+#[inline]
+pub(crate) fn transition(t0: f64, t1: f64, feedback: Array2<f64>) -> Array2<f64> {
+    let a = feedback * (t1 - t0);
+
+    if a.shape() == [1, 1] {
+        array![[a[(0, 0)].exp()]]
+    } else {
+        let mut b = Array2::<f64>::zeros(a.dim());
+
+        crate::expm::expm(&a, &mut b);
+
+        b
+    }
+}
+
+pub(crate) fn distance(ts1: &[f64], ts2: &[f64]) -> Array2<f64> {
+    let mut r = Array2::zeros((ts1.len(), ts2.len()));
+
+    for (i, v1) in ts1.iter().enumerate() {
+        for (j, v2) in ts2.iter().enumerate() {
+            r[(i, j)] = (v1 - v2).abs();
+        }
+    }
+
+    r
+}
+
 pub trait Kernel {
-    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> ArrayD<f64>;
+    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2<f64>;
     fn k_diag(&self, ts: &[f64]) -> Array1<f64>;
     fn order(&self) -> usize;
     fn state_mean(&self, t: f64) -> Array1<f64>;
@@ -28,14 +55,7 @@ pub trait Kernel {
     }
 
     fn transition(&self, t0: f64, t1: f64) -> Array2<f64> {
-        let f = self.feedback();
-
-        let a = f * (t1 - t0);
-        let mut b = Array2::<f64>::zeros(a.dim());
-
-        crate::expm::expm(&a, &mut b);
-
-        b
+        transition(t0, t1, self.feedback())
     }
 
     fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2<f64> {
@@ -60,7 +80,7 @@ pub trait Kernel {
 }
 
 impl Kernel for Vec<Box<dyn Kernel>> {
-    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> ArrayD<f64> {
+    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> Array2<f64> {
         unimplemented!();
     }
 
diff --git a/src/kernel/constant.rs b/src/kernel/constant.rs
index 09cba3f..f2f6b61 100644
--- a/src/kernel/constant.rs
+++ b/src/kernel/constant.rs
@@ -1,5 +1,4 @@
 use ndarray::prelude::*;
-use ndarray::IxDyn;
 
 use super::Kernel;
 
@@ -14,11 +13,11 @@ impl Constant {
 }
 
 impl Kernel for Constant {
-    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> ArrayD<f64> {
+    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2<f64> {
         let n = ts1.len();
         let m = ts2.map_or(n, |ts| ts.len());
 
-        ArrayD::ones(IxDyn(&[n, m])) * self.var
+        Array2::ones((n, m)) * self.var
     }
 
     fn k_diag(&self, ts: &[f64]) -> Array1<f64> {
@@ -60,6 +59,9 @@ impl Kernel for Constant {
 
 #[cfg(test)]
 mod tests {
+    extern crate blis_src;
+    extern crate intel_mkl_src;
+
     use approx::assert_abs_diff_eq;
     use rand::{distributions::Standard, thread_rng, Rng};
 
@@ -73,7 +75,7 @@ mod tests {
 
         assert_abs_diff_eq!(
             kernel.k_mat(&ts, None),
-            array![[2.5, 2.5, 2.5], [2.5, 2.5, 2.5], [2.5, 2.5, 2.5]].into_dyn()
+            array![[2.5, 2.5, 2.5], [2.5, 2.5, 2.5], [2.5, 2.5, 2.5]]
         );
     }
 
@@ -124,4 +126,18 @@ mod tests {
 
         assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts));
     }
+
+    #[test]
+    fn test_ssm_matrices() {
+        let kernel = Constant::new(2.5);
+
+        let deltas = [0.01, 1.0, 10.0];
+
+        for delta in &deltas {
+            assert_abs_diff_eq!(
+                crate::kernel::transition(0.0, *delta, kernel.feedback()),
+                kernel.transition(0.0, *delta)
+            );
+        }
+    }
 }
diff --git a/src/kernel/exponential.rs b/src/kernel/exponential.rs
index 83ecd9a..dc23857 100644
--- a/src/kernel/exponential.rs
+++ b/src/kernel/exponential.rs
@@ -14,8 +14,13 @@ impl Exponential {
 }
 
 impl Kernel for Exponential {
-    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> ArrayD<f64> {
-        unimplemented!();
+    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2<f64> {
+        let ts2 = ts2.unwrap_or(ts1);
+
+        let mut r = super::distance(ts1, ts2) / self.l_scale;
+        r.mapv_inplace(|v| (-v).exp());
+
+        r * self.var
     }
 
     fn k_diag(&self, ts: &[f64]) -> Array1<f64> {
@@ -42,6 +47,10 @@ impl Kernel for Exponential {
         (-1.0 / self.l_scale) * array![[1.0]]
     }
 
+    fn noise_effect(&self) -> Array2<f64> {
+        array![[1.0]]
+    }
+
     fn transition(&self, t0: f64, t1: f64) -> Array2<f64> {
         (-(t1 - t0) / self.l_scale).exp() * array![[1.0]]
     }
@@ -50,3 +59,92 @@ impl Kernel for Exponential {
         self.var * (1.0 - (-2.0 * (t1 - t0) / self.l_scale).exp()) * array![[1.0]]
     }
 }
+
+#[cfg(test)]
+mod tests {
+    extern crate blis_src;
+    extern crate intel_mkl_src;
+
+    use approx::assert_abs_diff_eq;
+    use rand::{distributions::Standard, thread_rng, Rng};
+
+    use super::*;
+
+    #[test]
+    fn test_kernel_matrix() {
+        let kernel = Exponential::new(1.1, 2.2);
+
+        let ts = [1.26, 1.46, 2.67];
+
+        assert_abs_diff_eq!(
+            kernel.k_mat(&ts, None),
+            array![
+                [1.1, 1.0044107879104887, 0.5794946136290716],
+                [1.0044107879104887, 1.1, 0.6346447914185355],
+                [0.5794946136290716, 0.6346447914185355, 1.1]
+            ]
+        );
+    }
+
+    #[test]
+    fn test_kernel_diag() {
+        let kernel = Exponential::new(1.1, 2.2);
+
+        let ts: Vec<_> = thread_rng()
+            .sample_iter::<f64, _>(Standard)
+            .take(10)
+            .map(|x| x * 10.0)
+            .collect();
+
+        assert_eq!(kernel.k_mat(&ts, None).diag(), kernel.k_diag(&ts));
+    }
+
+    #[test]
+    fn test_kernel_order() {
+        let kernel = Exponential::new(1.1, 2.2);
+
+        let m = kernel.order();
+
+        assert_eq!(kernel.state_mean(0.0).shape(), &[m]);
+        assert_eq!(kernel.state_cov(0.0).shape(), &[m, m]);
+        assert_eq!(kernel.measurement_vector().shape(), &[m]);
+        assert_eq!(kernel.feedback().shape(), &[m, m]);
+        assert_eq!(kernel.noise_effect().shape()[0], m);
+        assert_eq!(kernel.transition(0.0, 1.0).shape(), &[m, m]);
+        assert_eq!(kernel.noise_cov(0.0, 1.0).shape(), &[m, m]);
+    }
+
+    #[test]
+    fn test_ssm_variance() {
+        let kernel = Exponential::new(1.1, 2.2);
+
+        let ts: Vec<_> = thread_rng()
+            .sample_iter::<f64, _>(Standard)
+            .take(10)
+            .map(|x| x * 10.0)
+            .collect();
+
+        let h = kernel.measurement_vector();
+
+        let vars = ts
+            .iter()
+            .map(|t| h.dot(&kernel.state_cov(*t)).dot(&h))
+            .collect::<Vec<_>>();
+
+        assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts));
+    }
+
+    #[test]
+    fn test_ssm_matrices() {
+        let kernel = Exponential::new(1.1, 2.2);
+
+        let deltas = [0.01, 1.0, 10.0];
+
+        for delta in &deltas {
+            assert_abs_diff_eq!(
+                crate::kernel::transition(0.0, *delta, kernel.feedback()),
+                kernel.transition(0.0, *delta)
+            );
+        }
+    }
+}
diff --git a/src/kernel/matern32.rs b/src/kernel/matern32.rs
index 8b8037d..fc2002c 100644
--- a/src/kernel/matern32.rs
+++ b/src/kernel/matern32.rs
@@ -20,7 +20,7 @@ impl Matern32 {
 }
 
 impl Kernel for Matern32 {
-    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> ArrayD<f64> {
+    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> Array2<f64> {
         unimplemented!();
     }
 
diff --git a/src/kernel/matern52.rs b/src/kernel/matern52.rs
index cfe3c43..b4663d4 100644
--- a/src/kernel/matern52.rs
+++ b/src/kernel/matern52.rs
@@ -19,7 +19,7 @@ impl Matern52 {
 }
 
 impl Kernel for Matern52 {
-    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> ArrayD<f64> {
+    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> Array2<f64> {
         unimplemented!();
     }