diff --git a/.github/workflows/flare.yml b/.github/workflows/flare.yml
index 9429fa7de..42a7ba6de 100644
--- a/.github/workflows/flare.yml
+++ b/.github/workflows/flare.yml
@@ -27,9 +27,9 @@ jobs:
       CXX: g++-9
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
@@ -127,7 +127,7 @@ jobs:
         sudo apt-get install python3-sphinx python3-sphinx-rtd-theme python3-breathe python3-nbsphinx
 
     - name: Run Doxygen
-      uses: mattnotmitt/doxygen-action@v1.1.0
+      uses: mattnotmitt/doxygen-action@v1.9.8
       with:
         doxyfile-path: "./Doxyfile"
         working-directory: "./docs"
@@ -141,7 +141,7 @@ jobs:
         make html
 
     - name: Publish the docs
-      uses: peaceiris/actions-gh-pages@v3
+      uses: peaceiris/actions-gh-pages@v4
       with:
         github_token: ${{ secrets.GITHUB_TOKEN }}
         publish_dir: ./docs/build/html
diff --git a/lammps_plugins/compute_flare_std_atom.cpp b/lammps_plugins/compute_flare_std_atom.cpp
index 16f01c091..d63cf09c7 100644
--- a/lammps_plugins/compute_flare_std_atom.cpp
+++ b/lammps_plugins/compute_flare_std_atom.cpp
@@ -180,10 +180,17 @@ void ComputeFlareStdAtom::compute_peratom() {
       continue;
 
     if (use_map) {
-      int power = 2;
-      compute_energy_and_u(B2_vals, B2_norm_squared, single_bond_vals, power,
-              n_species, n_max, l_max, beta_matrices[itype - 1], u, &variance, normalized);
-      variance /= sig2;
+      Eigen::VectorXd Q_desc;
+      double K_self;
+      if (normalized) {
+        K_self = 1.0;
+        double B2_norm = pow(B2_norm_squared, 0.5);
+        Q_desc = beta_matrices[itype - 1].transpose() * B2_vals / B2_norm;
+      } else {
+        K_self = B2_norm_squared; // only power 1 is supported
+        Q_desc = beta_matrices[itype - 1].transpose() * B2_vals;
+      }
+      variance = K_self - Q_desc.dot(Q_desc);
     } else {
       Eigen::VectorXd kernel_vec = Eigen::VectorXd::Zero(n_clusters);
       double K_self;
@@ -462,18 +469,14 @@ void ComputeFlareStdAtom::read_file(char *filename) {
   int beta_count = 0;
   double beta_val;
   for (int k = 0; k < n_species; k++) {
-//    for (int l = 0; l < n_species; l++) {
-
-      beta_matrix = Eigen::MatrixXd::Zero(n_descriptors, n_descriptors);
-      for (int i = 0; i < n_descriptors; i++) {
-        for (int j = 0; j < n_descriptors; j++) {
-          //beta_matrix(k * n_descriptors + i, l * n_descriptors + j) = beta[beta_count];
-          beta_matrix(i, j) = beta[beta_count];
-          beta_count++;
-        }
+    beta_matrix = Eigen::MatrixXd::Zero(n_descriptors, n_descriptors);
+    for (int i = 0; i < n_descriptors; i++) {
+      for (int j = 0; j < n_descriptors; j++) {
+        beta_matrix(i, j) = beta[beta_count];
+        beta_count++;
       }
-      beta_matrices.push_back(beta_matrix);
-//    }
+    }
+    beta_matrices.push_back(beta_matrix);
   }
 
 }
diff --git a/requirements.txt b/requirements.txt
index 2d1bf81ba..fd60d6b2a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -15,6 +15,6 @@ wandb
 docutils==0.17.1
 alabaster==0.7.12
 babel
-pygments==2.11.2
+pygments==2.15.0
 nptyping
 nbsphinx
diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp
index edd2255a0..16e671020 100644
--- a/src/ace_binding.cpp
+++ b/src/ace_binding.cpp
@@ -110,7 +110,17 @@ PYBIND11_MODULE(_C_flare, m) {
   py::class_<B1, Descriptor>(m, "B1")
       .def(py::init<const std::string &, const std::string &,
                     const std::vector<double> &, const std::vector<double> &,
-                    const std::vector<int> &>());
+                    const std::vector<int> &>())
+      .def(py::init<const std::string &, const std::string &,
+                    const std::vector<double> &, const std::vector<double> &,
+                    const std::vector<int> &,
+                    const Eigen::MatrixXd &>())
+      .def_readonly("radial_basis", &B1::radial_basis)
+      .def_readonly("cutoff_function", &B1::cutoff_function)
+      .def_readonly("radial_hyps", &B1::radial_hyps)
+      .def_readonly("cutoff_hyps", &B1::cutoff_hyps)
+      .def_readonly("cutoffs", &B1::cutoffs)
+      .def_readonly("descriptor_settings", &B1::descriptor_settings);
 
   py::class_<B2, Descriptor>(m, "B2")
       .def(py::init<const std::string &, const std::string &,
@@ -140,7 +150,17 @@ PYBIND11_MODULE(_C_flare, m) {
   py::class_<B3, Descriptor>(m, "B3")
       .def(py::init<const std::string &, const std::string &,
                     const std::vector<double> &, const std::vector<double> &,
-                    const std::vector<int> &>());
+                    const std::vector<int> &>())
+      .def(py::init<const std::string &, const std::string &,
+                    const std::vector<double> &, const std::vector<double> &,
+                    const std::vector<int> &,
+                    const Eigen::MatrixXd &>())
+      .def_readonly("radial_basis", &B3::radial_basis)
+      .def_readonly("cutoff_function", &B3::cutoff_function)
+      .def_readonly("radial_hyps", &B3::radial_hyps)
+      .def_readonly("cutoff_hyps", &B3::cutoff_hyps)
+      .def_readonly("cutoffs", &B3::cutoffs)
+      .def_readonly("descriptor_settings", &B3::descriptor_settings);
 
   // Kernel functions
   py::class_<Kernel>(m, "Kernel");
diff --git a/src/flare_pp/descriptors/b1.cpp b/src/flare_pp/descriptors/b1.cpp
index 7568d5ed4..d3ee61ede 100644
--- a/src/flare_pp/descriptors/b1.cpp
+++ b/src/flare_pp/descriptors/b1.cpp
@@ -31,9 +31,26 @@ B1 ::B1(const std::string &radial_basis, const std::string &cutoff_function,
   cutoffs = Eigen::MatrixXd::Constant(n_species, n_species, cutoff_val);
 }
 
-void B1 ::write_to_file(std::ofstream &coeff_file, int coeff_size) {
-  coeff_file << "B1" << "\n";
+B1 ::B1(const std::string &radial_basis, const std::string &cutoff_function,
+        const std::vector<double> &radial_hyps,
+        const std::vector<double> &cutoff_hyps,
+        const std::vector<int> &descriptor_settings,
+        const Eigen::MatrixXd &cutoffs) {
+
+  this->radial_basis = radial_basis;
+  this->cutoff_function = cutoff_function;
+  this->radial_hyps = radial_hyps;
+  this->cutoff_hyps = cutoff_hyps;
+  this->descriptor_settings = descriptor_settings;
+
+  set_radial_basis(radial_basis, this->radial_pointer);
+  set_cutoff(cutoff_function, this->cutoff_pointer);
+
+  // Assign cutoff matrix.
+  this->cutoffs = cutoffs;
+}
 
+void B1 ::write_to_file(std::ofstream &coeff_file, int coeff_size) {
   // Report radial basis set.
   coeff_file << radial_basis << "\n";
 
diff --git a/src/flare_pp/descriptors/b1.h b/src/flare_pp/descriptors/b1.h
index df9ef42f6..69f3f31c2 100644
--- a/src/flare_pp/descriptors/b1.h
+++ b/src/flare_pp/descriptors/b1.h
@@ -29,6 +29,12 @@ class B1 : public Descriptor {
      const std::vector<double> &cutoff_hyps,
      const std::vector<int> &descriptor_settings);
 
+  B1(const std::string &radial_basis, const std::string &cutoff_function,
+     const std::vector<double> &radial_hyps,
+     const std::vector<double> &cutoff_hyps,
+     const std::vector<int> &descriptor_settings,
+     const Eigen::MatrixXd &cutoffs);
+
   DescriptorValues compute_struc(Structure &structure);
 
   void write_to_file(std::ofstream &coeff_file, int coeff_size);
diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp
index d2325823e..d89225e51 100644
--- a/src/flare_pp/descriptors/b3.cpp
+++ b/src/flare_pp/descriptors/b3.cpp
@@ -5,6 +5,8 @@
 #include "structure.h"
 #include "wigner3j.h"
 #include "y_grad.h"
+#include <fstream> // File operations
+#include <iomanip> // setprecision
 #include <iostream>
 
 B3 ::B3() {}
@@ -24,8 +26,57 @@ B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function,
 
   set_radial_basis(radial_basis, this->radial_pointer);
   set_cutoff(cutoff_function, this->cutoff_pointer);
+
+  // Create cutoff matrix.
+  int n_species = descriptor_settings[0];
+  double cutoff_val = radial_hyps[1];
+  cutoffs = Eigen::MatrixXd::Constant(n_species, n_species, cutoff_val);
+}
+
+B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function,
+        const std::vector<double> &radial_hyps,
+        const std::vector<double> &cutoff_hyps,
+        const std::vector<int> &descriptor_settings,
+        const Eigen::MatrixXd &cutoffs) {
+
+  this->radial_basis = radial_basis;
+  this->cutoff_function = cutoff_function;
+  this->radial_hyps = radial_hyps;
+  this->cutoff_hyps = cutoff_hyps;
+  this->descriptor_settings = descriptor_settings;
+
+  set_radial_basis(radial_basis, this->radial_pointer);
+  set_cutoff(cutoff_function, this->cutoff_pointer);
+
+  // Assign cutoff matrix.
+  this->cutoffs = cutoffs;
+}
+
+void B3 ::write_to_file(std::ofstream &coeff_file, int coeff_size) {
+  // Report radial basis set.
+  coeff_file << radial_basis << "\n";
+
+  // Record number of species, nmax, lmax, and the cutoff.
+  int n_species = descriptor_settings[0];
+  int n_max = descriptor_settings[1];
+  int l_max = descriptor_settings[2];
+  double cutoff = radial_hyps[1];
+
+  coeff_file << n_species << " " << n_max << " " << l_max << " ";
+  coeff_file << coeff_size << "\n";
+  coeff_file << cutoff_function << "\n";
+
+  // Report cutoffs to 2 decimal places.
+  coeff_file << std::fixed << std::setprecision(2);
+  for (int i = 0; i < n_species; i ++){
+    for (int j = 0; j < n_species; j ++){
+      coeff_file << cutoffs(i, j) << " ";
+    }
+  }
+  coeff_file << "\n";
 }
 
+
 DescriptorValues B3 ::compute_struc(Structure &structure) {
 
   // Initialize descriptor values.
@@ -41,10 +92,10 @@ DescriptorValues B3 ::compute_struc(Structure &structure) {
   int N = descriptor_settings[1];
   int lmax = descriptor_settings[2];
 
-  complex_single_bond(single_bond_vals, force_dervs, neighbor_coords,
+  complex_single_bond_multiple_cutoffs(single_bond_vals, force_dervs, neighbor_coords,
                       unique_neighbor_count, cumulative_neighbor_count,
                       descriptor_indices, radial_pointer, cutoff_pointer, nos,
-                      N, lmax, radial_hyps, cutoff_hyps, structure);
+                      N, lmax, radial_hyps, cutoff_hyps, structure, cutoffs);
 
   // Compute descriptor values.
   Eigen::MatrixXd B3_vals, B3_force_dervs;
@@ -237,7 +288,7 @@ void compute_B3(Eigen::MatrixXd &B3_vals, Eigen::MatrixXd &B3_force_dervs,
   }
 }
 
-void complex_single_bond(
+void complex_single_bond_multiple_cutoffs(
     Eigen::MatrixXcd &single_bond_vals, Eigen::MatrixXcd &force_dervs,
     Eigen::MatrixXd &neighbor_coordinates, Eigen::VectorXi &neighbor_count,
     Eigen::VectorXi &cumulative_neighbor_count,
@@ -249,14 +300,12 @@ void complex_single_bond(
                        std::vector<double>)>
         cutoff_function,
     int nos, int N, int lmax, const std::vector<double> &radial_hyps,
-    const std::vector<double> &cutoff_hyps, const Structure &structure) {
+    const std::vector<double> &cutoff_hyps, const Structure &structure,
+    const Eigen::MatrixXd &cutoffs) {
 
   int n_atoms = structure.noa;
   int n_neighbors = structure.n_neighbors;
 
-  // TODO: Make rcut an attribute of the descriptor calculator.
-  double rcut = radial_hyps[1];
-
   // Count atoms inside the descriptor cutoff.
   neighbor_count = Eigen::VectorXi::Zero(n_atoms);
   Eigen::VectorXi store_neighbors = Eigen::VectorXi::Zero(n_neighbors);
@@ -264,9 +313,12 @@ void complex_single_bond(
   for (int i = 0; i < n_atoms; i++) {
     int i_neighbors = structure.neighbor_count(i);
     int rel_index = structure.cumulative_neighbor_count(i);
+    int central_species = structure.species[i];
     for (int j = 0; j < i_neighbors; j++) {
       int current_count = neighbor_count(i);
       int neigh_index = rel_index + j;
+      int neighbor_species = structure.neighbor_species(neigh_index);
+      double rcut = cutoffs(central_species, neighbor_species);
       double r = structure.relative_positions(neigh_index, 0);
       // Check that atom is within descriptor cutoff.
       if (r <= rcut) {
@@ -312,6 +364,10 @@ void complex_single_bond(
     int i_neighbors = structure.neighbor_count(i);
     int rel_index = structure.cumulative_neighbor_count(i);
     int neighbor_index = cumulative_neighbor_count(i);
+    int central_species = structure.species[i];
+
+    // Initialize radial hyperparameters.
+    std::vector<double> new_radial_hyps = radial_hyps;
 
     // Initialize radial and spherical harmonic vectors.
     std::vector<double> g = std::vector<double>(N, 0);
@@ -326,6 +382,8 @@ void complex_single_bond(
     int s, neigh_index, descriptor_counter, unique_ind;
     for (int j = 0; j < i_neighbors; j++) {
       neigh_index = rel_index + j;
+      int neighbor_species = structure.neighbor_species(neigh_index);
+      double rcut = cutoffs(central_species, neighbor_species);
       r = structure.relative_positions(neigh_index, 0);
       if (r > rcut)
         continue; // Skip if outside cutoff.
@@ -334,6 +392,9 @@ void complex_single_bond(
       z = structure.relative_positions(neigh_index, 3);
       s = structure.neighbor_species(neigh_index);
 
+      // Reset the endpoint of the radial basis set.
+      new_radial_hyps[1] = rcut;
+
       // Store neighbor coordinates.
       neighbor_coordinates(neighbor_index, 0) = x;
       neighbor_coordinates(neighbor_index, 1) = y;
@@ -341,7 +402,7 @@ void complex_single_bond(
 
       // Compute radial basis values and spherical harmonics.
       calculate_radial(g, gx, gy, gz, radial_function, cutoff_function, x, y, z,
-                       r, rcut, N, radial_hyps, cutoff_hyps);
+                       r, rcut, N, new_radial_hyps, cutoff_hyps);
       get_complex_Y(h, hx, hy, hz, x, y, z, lmax);
 
       // Store the products and their derivatives.
diff --git a/src/flare_pp/descriptors/b3.h b/src/flare_pp/descriptors/b3.h
index b6da7ee8f..32a144e83 100644
--- a/src/flare_pp/descriptors/b3.h
+++ b/src/flare_pp/descriptors/b3.h
@@ -22,6 +22,8 @@ class B3 : public Descriptor {
 
   std::string descriptor_name = "B3";
 
+  Eigen::MatrixXd cutoffs;
+
   B3();
 
   B3(const std::string &radial_basis, const std::string &cutoff_function,
@@ -29,8 +31,16 @@ class B3 : public Descriptor {
      const std::vector<double> &cutoff_hyps,
      const std::vector<int> &descriptor_settings);
 
+  B3(const std::string &radial_basis, const std::string &cutoff_function,
+     const std::vector<double> &radial_hyps,
+     const std::vector<double> &cutoff_hyps,
+     const std::vector<int> &descriptor_settings,
+     const Eigen::MatrixXd &cutoffs);
+
   DescriptorValues compute_struc(Structure &structure);
 
+  void write_to_file(std::ofstream &coeff_file, int coeff_size);
+
   nlohmann::json return_json();
 };
 
@@ -43,7 +53,7 @@ void compute_B3(Eigen::MatrixXd &B3_vals, Eigen::MatrixXd &B3_force_dervs,
                 const Eigen::VectorXi &descriptor_indices, int nos, int N,
                 int lmax, const Eigen::VectorXd &wigner3j_coeffs);
 
-void complex_single_bond(
+void complex_single_bond_multiple_cutoffs(
     Eigen::MatrixXcd &single_bond_vals, Eigen::MatrixXcd &force_dervs,
     Eigen::MatrixXd &neighbor_coordinates, Eigen::VectorXi &neighbor_count,
     Eigen::VectorXi &cumulative_neighbor_count,
@@ -55,6 +65,7 @@ void complex_single_bond(
                        std::vector<double>)>
         cutoff_function,
     int nos, int N, int lmax, const std::vector<double> &radial_hyps,
-    const std::vector<double> &cutoff_hyps, const Structure &structure);
+    const std::vector<double> &cutoff_hyps, const Structure &structure,
+    const Eigen::MatrixXd &cutoffs);
 
 #endif
diff --git a/src/flare_pp/kernels/dot_product.cpp b/src/flare_pp/kernels/dot_product.cpp
index f61482ef3..040b20783 100644
--- a/src/flare_pp/kernels/dot_product.cpp
+++ b/src/flare_pp/kernels/dot_product.cpp
@@ -786,8 +786,8 @@ Eigen::MatrixXd DotProduct ::compute_varmap_coefficients(
   int beta_size = p_size * (p_size + 1) / 2;
   int n_species = gp_model.sparse_descriptors[kernel_index].n_types;
   int n_sparse = gp_model.sparse_descriptors[kernel_index].n_clusters;
-  //mapping_coeffs = Eigen::MatrixXd::Zero(n_species * n_species, p_size * p_size); // can be reduced by symmetry
   mapping_coeffs = Eigen::MatrixXd::Zero(n_species, p_size * p_size); // can be reduced by symmetry
+  double jitter_root = pow(gp_model.Kuu_jitter, 0.5);
 
   // Get alpha index.
   
@@ -799,58 +799,17 @@ Eigen::MatrixXd DotProduct ::compute_varmap_coefficients(
   // Loop over types.
   for (int s = 0; s < n_species; s++){
     int n_types = gp_model.sparse_descriptors[kernel_index].n_clusters_by_type[s];
-    int c_types =
-      gp_model.sparse_descriptors[kernel_index].cumulative_type_count[s];
-    int K_ind = alpha_ind + c_types;
-
-    // Loop over clusters within each type.
-    for (int i = 0; i < n_types; i++){
-      Eigen::VectorXd pi_current =
-        gp_model.sparse_descriptors[kernel_index].descriptors[s].row(i);
-      double pi_norm =
-        gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](i);
-
-      // Skip empty environments.
-      if (pi_norm < empty_thresh)
-        continue;
-
-      // TODO: include symmetry of i & j
-      // Loop over clusters within each type.
-      for (int j = 0; j < n_types; j++){
-        Eigen::VectorXd pj_current =
-          gp_model.sparse_descriptors[kernel_index].descriptors[s].row(j);
-        double pj_norm =
-          gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](j);
-
-        // Skip empty environments.
-        if (pj_norm < empty_thresh)
-          continue;
-
-        double Kuu_inv_ij = gp_model.Kuu_inverse(K_ind + i, K_ind + j);
-        double Kuu_inv_ij_normed = Kuu_inv_ij; // / pi_norm / pj_norm;
-//        double Sigma_ij = gp_model.Sigma(K_ind + i, K_ind + j);
-//        double Sigma_ij_normed = Sigma_ij / pi_norm / pj_norm;
-        int beta_count = 0;
-
-        // First loop over descriptor values.
-        for (int k = 0; k < p_size; k++) {
-          double p_ik = pi_current(k);
-    
-          // Second loop over descriptor values.
-          for (int l = 0; l < p_size; l++){
-            double p_jl = pj_current(l);
-    
-            // Update beta vector.
-            double beta_val = sig2 * sig2 * p_ik * p_jl * (- Kuu_inv_ij_normed); // + Sigma_ij_normed); // To match the compute_cluster_uncertainty function
-            mapping_coeffs(s, beta_count) += beta_val;
-
-            if (k == l && i == 0 && j == 0) {
-              mapping_coeffs(s, beta_count) += sig2; // the self kernel term
-            }
-    
-            beta_count++;
-          }
-        }
+    Eigen::MatrixXd pi = gp_model.sparse_descriptors[kernel_index].descriptors[s];
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(p_size + n_types, n_types);
+    A.block(0, 0, p_size, n_types) = pi.transpose();
+    A.block(p_size, 0, n_types, n_types) = jitter_root * Eigen::MatrixXd::Identity(n_types, n_types) / sigma;
+    Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
+    Eigen::MatrixXd Q = qr.householderQ(); // (p_size + n_types, n_types)
+    int beta_count = 0;
+    for (int i = 0; i < p_size; i++) {
+      for (int j = 0; j < p_size; j++) {
+        if (j < n_types) mapping_coeffs(s, beta_count) += Q(i, j);
+        beta_count++;
       }
     }
   }
diff --git a/src/flare_pp/kernels/normalized_dot_product.cpp b/src/flare_pp/kernels/normalized_dot_product.cpp
index caecef892..587e93171 100644
--- a/src/flare_pp/kernels/normalized_dot_product.cpp
+++ b/src/flare_pp/kernels/normalized_dot_product.cpp
@@ -784,8 +784,8 @@ Eigen::MatrixXd NormalizedDotProduct ::compute_varmap_coefficients(
   int beta_size = p_size * (p_size + 1) / 2;
   int n_species = gp_model.sparse_descriptors[kernel_index].n_types;
   int n_sparse = gp_model.sparse_descriptors[kernel_index].n_clusters;
-  //mapping_coeffs = Eigen::MatrixXd::Zero(n_species * n_species, p_size * p_size); // can be reduced by symmetry
   mapping_coeffs = Eigen::MatrixXd::Zero(n_species, p_size * p_size); // can be reduced by symmetry
+  double jitter_root = pow(gp_model.Kuu_jitter, 0.5);
 
   // Get alpha index.
   
@@ -797,58 +797,28 @@ Eigen::MatrixXd NormalizedDotProduct ::compute_varmap_coefficients(
   // Loop over types.
   for (int s = 0; s < n_species; s++){
     int n_types = gp_model.sparse_descriptors[kernel_index].n_clusters_by_type[s];
-    int c_types =
-      gp_model.sparse_descriptors[kernel_index].cumulative_type_count[s];
-    int K_ind = alpha_ind + c_types;
-
+    Eigen::MatrixXd pi = gp_model.sparse_descriptors[kernel_index].descriptors[s];
     // Loop over clusters within each type.
     for (int i = 0; i < n_types; i++){
-      Eigen::VectorXd pi_current =
-        gp_model.sparse_descriptors[kernel_index].descriptors[s].row(i);
       double pi_norm =
         gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](i);
 
       // Skip empty environments.
       if (pi_norm < empty_thresh)
         continue;
+      pi.row(i) /= pi_norm; // normalize descriptor by norm
+    }
 
-      // TODO: include symmetry of i & j
-      // Loop over clusters within each type.
-      for (int j = 0; j < n_types; j++){
-        Eigen::VectorXd pj_current =
-          gp_model.sparse_descriptors[kernel_index].descriptors[s].row(j);
-        double pj_norm =
-          gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](j);
-
-        // Skip empty environments.
-        if (pj_norm < empty_thresh)
-          continue;
-
-        double Kuu_inv_ij = gp_model.Kuu_inverse(K_ind + i, K_ind + j);
-        double Kuu_inv_ij_normed = Kuu_inv_ij / pi_norm / pj_norm;
-//        double Sigma_ij = gp_model.Sigma(K_ind + i, K_ind + j);
-//        double Sigma_ij_normed = Sigma_ij / pi_norm / pj_norm;
-        int beta_count = 0;
-
-        // First loop over descriptor values.
-        for (int k = 0; k < p_size; k++) {
-          double p_ik = pi_current(k);
-    
-          // Second loop over descriptor values.
-          for (int l = 0; l < p_size; l++){
-            double p_jl = pj_current(l);
-    
-            // Update beta vector.
-            double beta_val = sig2 * sig2 * p_ik * p_jl * (- Kuu_inv_ij_normed); // + Sigma_ij_normed); // To match the compute_cluster_uncertainty function
-            mapping_coeffs(s, beta_count) += beta_val;
-
-            if (k == l && i == 0 && j == 0) {
-              mapping_coeffs(s, beta_count) += sig2; // the self kernel term
-            }
-    
-            beta_count++;
-          }
-        }
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(p_size + n_types, n_types);
+    A.block(0, 0, p_size, n_types) = pi.transpose();
+    A.block(p_size, 0, n_types, n_types) = jitter_root * Eigen::MatrixXd::Identity(n_types, n_types) / sigma;
+    Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
+    Eigen::MatrixXd Q = qr.householderQ(); // (p_size + n_types, n_types)
+    int beta_count = 0;
+    for (int i = 0; i < p_size; i++) {
+      for (int j = 0; j < p_size; j++) {
+        if (j < n_types) mapping_coeffs(s, beta_count) += Q(i, j);
+        beta_count++;
       }
     }
   }
diff --git a/tests/get_sgp.py b/tests/get_sgp.py
index d42d1a147..f1f869ef1 100644
--- a/tests/get_sgp.py
+++ b/tests/get_sgp.py
@@ -8,21 +8,16 @@
 from ase.build import make_supercell
 
 # Define kernel.
-sigma = 2.0
+sigma = np.random.rand() + 1.0
 power = 1.0
 dotprod_kernel = DotProduct(sigma, power)
 normdotprod_kernel = NormalizedDotProduct(sigma, power)
 
 # Define remaining parameters for the SGP wrapper.
-sigma_e = 0.3
-sigma_f = 0.2
-sigma_s = 0.1
-species_map = {6: 0, 8: 1}
-single_atom_energies = {0: -5, 1: -6}
 variance_type = "local"
 max_iterations = 20
 opt_method = "L-BFGS-B"
-bounds = [(None, None), (sigma_e, None), (None, None), (None, None)]
+bounds = [(None, None), (None, None), (None, None), (None, None)]
 
 
 def get_random_atoms(a=2.0, sc_size=2, numbers=[6, 8], set_seed: int = None):
@@ -89,6 +84,18 @@ def get_empty_sgp(n_types=2, power=2, multiple_cutoff=False, kernel_type="Normal
         cutoff_matrix,
     )
 
+    sigma_e = np.random.rand() 
+    sigma_f = np.random.rand() 
+    sigma_s = np.random.rand() 
+    print("Hyps", kernel.sigma, sigma_e, sigma_f, sigma_s)
+
+    if n_types == 1:
+        species_map = {6: 0}
+        single_atom_energies = {0: np.random.rand()}
+    elif n_types == 2:
+        species_map = {6: 0, 8: 1}
+        single_atom_energies = {0: np.random.rand(), 1: np.random.rand()}
+
     empty_sgp = SGP_Wrapper(
         [kernel],
         [b2_calc],
@@ -123,16 +130,18 @@ def get_updated_sgp(n_types=2, power=2, multiple_cutoff=False, kernel_type="Norm
     energy = training_structure.get_potential_energy()
     stress = training_structure.get_stress()
 
+    size = max(1, np.random.randint(len(training_structure)))
+    custom_range = np.random.choice(len(training_structure), size=size, replace=False).tolist()
     sgp.update_db(
         training_structure,
         forces,
-        custom_range=(1, 2, 3, 4, 5),
+        custom_range=custom_range,
         energy=energy,
         stress=stress,
         mode="specific",
-        rel_e_noise=0.1,
-        rel_f_noise=0.2,
-        rel_s_noise=0.1,
+        rel_e_noise=np.random.rand(),
+        rel_f_noise=np.random.rand(),
+        rel_s_noise=np.random.rand(),
     )
 
     # add an isolated atom to the training data
diff --git a/tests/test_lammps.py b/tests/test_lammps.py
index eaa1e6e40..f7ead128e 100644
--- a/tests/test_lammps.py
+++ b/tests/test_lammps.py
@@ -28,7 +28,8 @@
 @pytest.mark.parametrize("struc", struc_list)
 @pytest.mark.parametrize("multicut", [False, True])
 @pytest.mark.parametrize("n_cpus", n_cpus_list)
-def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus):
+@pytest.mark.parametrize("kernel_type", ["NormalizedDotProduct", "DotProduct"])
+def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus, kernel_type):
     """Test the flare pair style."""
 
     if n_species > n_types:
@@ -100,11 +101,11 @@ def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus):
             energy_lmp += sgp_model.gp_model.single_atom_energies[coded_spec]
 
     thresh = 1e-6
+    r_thresh = 1e-3
     print(energy, energy_lmp)
-    assert np.allclose(energy, energy_lmp, atol=thresh)
-    # print(forces, forces_lmp)
-    assert np.allclose(forces, forces_lmp, atol=thresh)
-    assert np.allclose(stress, stress_lmp, atol=thresh)
+    assert np.allclose(energy, energy_lmp, atol=thresh, rtol=r_thresh)
+    assert np.allclose(forces, forces_lmp, atol=thresh, rtol=r_thresh)
+    assert np.allclose(stress, stress_lmp, atol=thresh, rtol=r_thresh)
 
     # Remove files.
     os.remove(potential_name)
@@ -119,7 +120,6 @@ def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus):
 from ase.calculators.lammpsrun import LAMMPS
 from flare.md.lammps import LAMMPS_MOD, LAMMPS_MD, get_kinetic_stress
 
-
 @pytest.mark.skipif(
     not os.environ.get("lmp", False),
     reason=(