diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/CMakeLists.txt b/fastlib/trunk/contrib/dongryel/thesis_research/core/math/CMakeLists.txt
index 95697688dc..2ec51303b0 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/CMakeLists.txt
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/math/CMakeLists.txt
@@ -6,7 +6,6 @@ set(SOURCES
   global.cc
   math_lib.h
   math_lib_impl.h
-  linear_algebra.h
   range.h
 )
 
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
deleted file mode 100644
index 4b645c2d12..0000000000
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/math/linear_algebra.h
+++ /dev/null
@@ -1,225 +0,0 @@
-/** @file linear_algebra.h
- *
- *  @author Dongryeol Lee (dongryel@cc.gatech.edu)
- */
-
-#ifndef CORE_MATH_LINEAR_ALGEBRA_H
-#define CORE_MATH_LINEAR_ALGEBRA_H
-
-#include <armadillo>
-
-namespace core {
-namespace table {
-class DensePoint;
-class DenseMatrix;
-};
-};
-
-namespace core {
-namespace math {
-
-template<typename MatrixType>
-static void MatrixTripleProduct(
-  const MatrixType &left, const MatrixType &mid,
-  core::table::DenseMatrix *product) {
-
-  product->Init(left.n_rows(), left.n_cols());
-
-  // Use armadillo matrices to compute the triple product. This makes
-  // a light copy, so there is very little performance lost.
-  arma::mat left_alias(left.ptr(), left.n_rows(), left.n_cols());
-  arma::mat mid_alias(mid.ptr(), mid.n_rows(), mid.n_cols());
-  arma::mat product_alias(
-    product->ptr(), product->n_rows(), product->n_cols(), false);
-  product_alias = left_alias * mid_alias * arma::trans(left_alias);
-}
-
-template<typename MatrixType>
-static void MatrixTripleProduct(
-  const MatrixType &left, const MatrixType &mid, const MatrixType &right,
-  core::table::DenseMatrix *product) {
-
-  product->Init(left.n_rows(), right.n_cols());
-
-  // Use armadillo matrices to compute the triple product. This makes
-  // a light copy, so there is very little performance lost.
-  arma::mat left_alias(left.ptr(), left.n_rows(), left.n_cols());
-  arma::mat mid_alias(mid.ptr(), mid.n_rows(), mid.n_cols());
-  arma::mat right_alias(right.ptr(), right.n_rows(), right.n_cols());
-  arma::mat product_alias(
-    product->ptr(), product->n_rows(), product->n_cols(), false);
-  product_alias = left_alias * mid_alias * right_alias;
-}
-
-template<typename VectorType>
-static double Dot(const VectorType &a, const VectorType &b) {
-  arma::vec a_mat(a.ptr(), a.length());
-  arma::vec b_mat(b.ptr(), b.length());
-  return arma::dot(a_mat, b_mat);
-}
-
-template<typename VectorType>
-static double LengthEuclidean(const VectorType &a) {
-  return sqrt(core::math::Dot(a, a));
-}
-
-template<typename VectorType>
-static void CopyValues(
-  const VectorType &vec_in, VectorType *vec_out) {
-
-  for(unsigned int i = 0; i < vec_in.n_elem; i++) {
-    (*vec_out)[i] = vec_in[i];
-  }
-}
-
-template<typename VectorType>
-static void SubFrom(
-  const VectorType &vec_in, VectorType *vec_out) {
-
-  for(unsigned int i = 0; i < vec_in.n_elem; i++) {
-    (*vec_out)[i] -= vec_in[i];
-  }
-}
-
-template<typename VectorType>
-static void SubOverwrite(
-  const VectorType &sub, const VectorType &sub_from, VectorType *vec_out) {
-  for(int i = 0; i < sub.length(); i++) {
-    (*vec_out)[i] = sub_from[i] - sub[i];
-  }
-}
-
-template<typename VectorType>
-static void SubInit(
-  const VectorType &sub, const VectorType &sub_from, VectorType *vec_out) {
-  vec_out->Init(sub.length());
-  SubOverwrite(sub, sub_from, vec_out);
-}
-
-template<typename T>
-class AddExpertTrait {
-  public:
-    static void Compute(double scale, const T &vec_scaled, T *vec_add_to);
-};
-
-template<>
-class AddExpertTrait<core::table::DenseMatrix> {
-  public:
-    static void Compute(
-      double scale, const core::table::DenseMatrix &mat_scaled,
-      core::table::DenseMatrix *mat_add_to) {
-
-      arma::mat mat_scaled_alias(
-        mat_scaled.ptr(), mat_scaled.n_rows(), mat_scaled.n_cols());
-      arma::mat mat_add_to_alias(
-        mat_add_to->ptr(), mat_add_to->n_rows(), mat_add_to->n_cols(), false);
-      mat_add_to_alias = mat_add_to_alias + scale * mat_scaled_alias;
-    }
-};
-
-template<>
-class AddExpertTrait<core::table::DensePoint> {
-  public:
-    static void Compute(
-      double scale, const core::table::DensePoint &vec_scaled,
-      core::table::DensePoint *vec_add_to) {
-      arma::vec vec_scaled_alias(
-        vec_scaled.ptr(), vec_scaled.length());
-      arma::vec vec_add_to_alias(
-        vec_add_to->ptr(), vec_add_to->length(), false);
-      vec_add_to_alias = vec_add_to_alias + scale * vec_scaled_alias;
-    }
-};
-
-template<typename T>
-static void AddExpert(
-  double scale, const T &vec_scaled, T *vec_add_to) {
-  core::math::AddExpertTrait<T>::Compute(scale, vec_scaled, vec_add_to);
-}
-
-template<typename VectorType>
-static void AddTo(
-  const VectorType &vec_in, VectorType *vec_out) {
-
-  arma::vec vec_in_alias(vec_in.ptr(), vec_in.length());
-  arma::vec vec_out_alias(vec_out->ptr(), vec_in.length(), false);
-  vec_out_alias = vec_out_alias + vec_in_alias;
-}
-
-/** @brief Computes $c = c + \alpha * a b^T$.
- */
-template<typename VectorType, typename MatrixType>
-static void MulExpert(
-  double alpha, const VectorType &a, const VectorType &b, MatrixType *c) {
-
-  arma::vec a_alias(a.ptr(), a.length());
-  arma::vec b_alias(b.ptr(), b.length());
-  arma::mat c_alias(c->ptr(), c->n_rows(), c->n_cols(), false);
-  c_alias = c_alias + alpha * a_alias * arma::trans(b_alias);
-}
-
-template<typename MatrixType, typename VectorType>
-static void MulInit(
-  const MatrixType &a, const VectorType &b, VectorType *c) {
-
-  arma::mat a_alias(a.ptr(), a.n_rows(), a.n_cols());
-  arma::vec b_alias(b.ptr(), b.length());
-  c->Init(a.n_rows());
-  arma::vec c_alias(c->ptr(), c.length(), false);
-  c_alias = a_alias * b_alias;
-}
-
-template<typename VectorType>
-class ScaleTrait {
-  public:
-    static void Compute(double scale, VectorType *vec);
-};
-
-template<>
-class ScaleTrait<core::table::DensePoint> {
-  public:
-    static void Compute(double scale, core::table::DensePoint *vec) {
-      for(int i = 0; i < vec->length(); i++) {
-        (*vec)[i] *= scale;
-      }
-    }
-};
-
-template<>
-class ScaleTrait<core::table::DenseMatrix> {
-  public:
-    static void Compute(double scale, core::table::DenseMatrix *vec) {
-      for(int j = 0; j < vec->n_cols(); j++) {
-        for(int i = 0; i < vec->n_rows(); i++) {
-          vec->set(i, j, vec->get(i, j) * scale);
-        }
-      }
-    }
-};
-
-template<typename VectorType>
-static void Scale(double scale, VectorType *vec) {
-  ScaleTrait<VectorType>::Compute(scale, vec);
-}
-
-template<typename VectorType>
-static void ScaleOverwrite(
-  double scale, const VectorType &vec_in, VectorType *vec_out) {
-  arma::vec vec_in_alias(vec_in.ptr(), vec_in.length());
-  arma::vec vec_out_alias(vec_out->ptr(), vec_in.length(), false);
-  for(unsigned int i = 0; i < vec_in_alias.n_elem; i++) {
-    vec_out_alias[i] = vec_in_alias[i] * scale;
-  }
-}
-
-template<typename VectorType>
-static void ScaleInit(
-  double scale, const VectorType &vec_in, VectorType *vec_out) {
-
-  vec_out->set_size(vec_in.n_elem);
-  ScaleOverwrite(scale, vec_in, vec_out);
-}
-};
-};
-
-#endif
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/monte_carlo/mean_variance_pair_matrix.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/monte_carlo/mean_variance_pair_matrix.h
index de58dadcd3..8744eb74f8 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/monte_carlo/mean_variance_pair_matrix.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/monte_carlo/mean_variance_pair_matrix.h
@@ -6,6 +6,7 @@
 #ifndef CORE_MONTE_CARLO_MEAN_VARIANCE_PAIR_MATRIX_H
 #define CORE_MONTE_CARLO_MEAN_VARIANCE_PAIR_MATRIX_H
 
+#include <armadillo>
 #include "core/monte_carlo/mean_variance_pair.h"
 
 namespace core {
@@ -55,8 +56,8 @@ class MeanVariancePairVector {
              new core::monte_carlo::MeanVariancePair[n_elements_];
     }
 
-    void sample_means(core::table::DensePoint *point_out) const {
-      point_out->Init(n_elements_);
+    void sample_means(arma::vec *point_out) const {
+      point_out->set_size(n_elements_);
       for(int i = 0; i < n_elements_; i++) {
         (*point_out)[i] = ptr_[i].sample_mean();
       }
@@ -110,11 +111,11 @@ class MeanVariancePairMatrix {
       }
     }
 
-    void sample_means(core::table::DenseMatrix *point_out) const {
-      point_out->Init(n_rows_, n_cols_);
+    void sample_means(arma::mat *point_out) const {
+      point_out->set_size(n_rows_, n_cols_);
       for(int j = 0; j < n_cols_; j++) {
         for(int i = 0; i < n_rows_; i++) {
-          point_out->set(i, j, this->get(i, j).sample_mean());
+          point_out->at(i, j) = this->get(i, j).sample_mean();
         }
       }
     }
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h
index e613f04c21..cdb9876ad5 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region.h
@@ -14,11 +14,16 @@ namespace core {
 namespace optimization {
 template<typename FunctionType>
 class TrustRegion {
+  public:
+    enum TrustRegionSearchMethod {CAUCHY, DOGLEG, STEIHAUG};
+
   private:
     double max_radius_;
 
     FunctionType *function_;
 
+    TrustRegionSearchMethod search_method_;
+
   private:
 
     void Evaluate_(const arma::vec &iterate);
@@ -53,7 +58,8 @@ class TrustRegion {
       max_radius_ = max_radius_in;
     }
 
-    void Init(FunctionType &function_in);
+    void Init(
+      FunctionType &function_in, TrustRegionSearchMethod search_method_in);
 
     void Optimize(int num_iterations, arma::vec *iterate);
 };
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
index 753960217e..7b6533a9b1 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/optimization/trust_region_dev.h
@@ -16,6 +16,12 @@ template<typename FunctionType>
 TrustRegion<FunctionType>::TrustRegion() {
   max_radius_ = 10.0;
   function_ = NULL;
+  search_method_ = CAUCHY;
+}
+
+template<typename FunctionType>
+double TrustRegion<FunctionType>::ReductionRatio_(const arma::vec &iterate, const arma::vec &step,) {
+
 }
 
 template<typename FunctionType>
@@ -266,12 +272,16 @@ template<typename FunctionType>
 void TrustRegion<FunctionType>::ObtainStepDirection_(
   arma::vec *step_direction, double *step_direction_norm) {
 
+  switch(search_method_) {
 
+  }
 }
 
 template<typename FunctionType>
-void TrustRegion<FunctionType>::Init(FunctionType &function_in) {
+void TrustRegion<FunctionType>::Init(
+  FunctionType &function_in, TrustRegionSearchMethod search_method_in) {
   function_ = &function_in;
+  search_method_ = search_method_in;
 }
 
 template<typename FunctionType>
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
index 59dee7fe50..556c287cbd 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/dense_point.h
@@ -8,6 +8,7 @@
 #ifndef CORE_TABLE_DENSE_POINT_H
 #define CORE_TABLE_DENSE_POINT_H
 
+#include <armadillo>
 #include <vector>
 #include <boost/interprocess/offset_ptr.hpp>
 #include <boost/serialization/serialization.hpp>
@@ -175,6 +176,16 @@ class DensePoint {
       printf("\n");
     }
 };
+
+template<typename DensePointType>
+static void DensePointToArmaVec(
+  const DensePointType &point_in, arma::vec *vec_out) {
+  const_cast<arma::u32 &>(vec_out->n_rows) = point_in.length();
+  const_cast<arma::u32 &>(vec_out->n_cols) = 1;
+  const_cast<arma::u32 &>(vec_out->n_elem) = point_in.length();
+  const_cast<bool &>(vec_out->use_aux_mem) = true;
+  const_cast<double *&>(vec_out->mem) = const_cast<double *>(point_in.ptr());
+}
 };
 };
 
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
index 41d57b044f..31b9665be1 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/table/distributed_table.h
@@ -152,7 +152,9 @@ class DistributedTable: public boost::noncopyable {
     void ReplenishNodes_(std::vector<TreeType *> &top_leaf_nodes) {
 
       core::table::DensePoint tmp_point;
+      arma::vec tmp_point_alias;
       tmp_point.Init(top_leaf_nodes[0]->bound().center().length());
+      core::table::DensePointToArmaVec(tmp_point, &tmp_point_alias);
       int num_additional = table_outbox_group_comm_size_ -
                            top_leaf_nodes.size();
       int num_samples = std::max(1, core::math::RandInt(top_leaf_nodes.size()));
@@ -160,15 +162,17 @@ class DistributedTable: public boost::noncopyable {
       // Randomly add new dummy nodes with randomly chosen centroids
       // averaged.
       for(int j = 0; j < num_additional; j++) {
-        tmp_point.SetZero();
+        tmp_point_alias.zeros();
         for(int i = 0; i < num_samples; i++) {
-          tmp_point.SetZero();
-          core::math::AddTo(
+          arma::vec random_node_center;
+          core::table::DensePointToArmaVec(
             top_leaf_nodes[
               core::math::RandInt(top_leaf_nodes.size())]->bound().center(),
-            &tmp_point);
+            &random_node_center);
+          tmp_point_alias += random_node_center;
         }
-        core::math::Scale(1.0 / static_cast<double>(num_samples), &tmp_point);
+        tmp_point_alias = (1.0 / static_cast<double>(num_samples)) *
+                          tmp_point_alias;
         top_leaf_nodes.push_back(new TreeType());
         top_leaf_nodes[ top_leaf_nodes.size() - 1 ]->bound().center().Copy(
           tmp_point);
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
index 699bf0de92..7bdcee803a 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/distributed_local_kmeans.h
@@ -6,9 +6,9 @@
 #ifndef CORE_TREE_DISTRIBUTED_LOCAL_KMEANS_H
 #define CORE_TREE_DISTRIBUTED_LOCAL_KMEANS_H
 
+#include <armadillo>
 #include <boost/mpi.hpp>
 #include <boost/serialization/string.hpp>
-#include "core/math/linear_algebra.h"
 #include "core/table/dense_point.h"
 
 namespace core {
@@ -25,12 +25,15 @@ class DistributedLocalKMeans {
 
         core::table::DensePoint centroid_;
 
+        arma::vec centroid_alias_;
+
         int num_points_;
 
       public:
 
         void Init(int length_in) {
           centroid_.Init(length_in);
+          core::table::DensePointToArmaVec(centroid_, &centroid_alias_);
         }
 
         const core::table::DensePoint &centroid() const {
@@ -53,22 +56,26 @@ class DistributedLocalKMeans {
           double factor =
             static_cast<double>(num_points_) /
             static_cast<double>(num_points_ + centroid_in.num_points());
-          core::math::Scale(factor, &centroid_);
-          core::math::AddExpert(
-            1.0 - factor, centroid_in.centroid(), &centroid_);
+          arma::vec centroid_in_alias;
+          core::table::DensePointToArmaVec(
+            centroid_in.centroid(), &centroid_in_alias);
+          centroid_alias_ = factor * centroid_alias_ +
+                            (1.0 - factor) * centroid_in_alias;
           num_points_ = num_points_ + centroid_in.num_points();
         }
 
         void Add(const core::table::DensePoint &point_in) {
           double factor = static_cast<double>(num_points_) /
                           static_cast<double>(num_points_ + 1);
-          core::math::Scale(factor, &centroid_);
-          core::math::AddExpert(1.0 - factor, point_in, &centroid_);
+          arma::vec point_in_alias;
+          core::table::DensePointToArmaVec(point_in, &point_in_alias);
+          centroid_alias_ = factor * centroid_alias_;
+          centroid_alias_ += (1.0 - factor) * point_in_alias;
           num_points_++;
         }
 
         void Reset() {
-          centroid_.SetZero();
+          centroid_alias_.zeros();
           num_points_ = 0;
         }
 
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
index d1fb53c219..b0c035ca10 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/core/tree/gen_metric_tree.h
@@ -9,7 +9,6 @@
 #include <vector>
 #include "ball_bound.h"
 #include "general_spacetree.h"
-#include "core/math/linear_algebra.h"
 #include "core/metric_kernels/abstract_metric.h"
 #include "core/table/dense_matrix.h"
 #include "core/table/memory_mapped_file.h"
@@ -69,11 +68,15 @@ class GenMetricTree {
 
       int end = begin + count;
       core::table::DensePoint col_point;
+      arma::vec bound_ref;
+      core::table::DensePointToArmaVec(bounds->center(), &bound_ref);
       for(int i = begin; i < end; i++) {
         matrix.MakeColumnVector(i, &col_point);
-        core::math::AddTo(col_point, &(bounds->center()));
+        arma::vec col_point_ref;
+        core::table::DensePointToArmaVec(col_point, &col_point_ref);
+        bound_ref += col_point_ref;
       }
-      core::math::Scale(1.0 / static_cast<double>(count), &(bounds->center()));
+      bound_ref = (1.0 / static_cast<double>(count)) * bound_ref;
 
       double furthest_distance;
       FurthestColumnIndex_(
@@ -88,12 +91,17 @@ class GenMetricTree {
       TreeType *node, TreeType *left, TreeType *right) {
 
       // Compute the weighted sum of the two pivots
-      node->bound().center().CopyValues(left->bound().center());
-      core::math::Scale(left->count(), &(node->bound().center()));
-      core::math::AddExpert(
-        right->count(), right->bound().center(), & (node->bound().center()));
-      core::math::Scale(
-        1.0 / static_cast<double>(node->count()), & (node->bound().center()));
+      arma::vec bound_ref;
+      core::table::DensePointToArmaVec(node->bound().center(), &bound_ref);
+      arma::vec left_bound_ref;
+      core::table::DensePointToArmaVec(left->bound().center(), &left_bound_ref);
+      arma::vec right_bound_ref;
+      core::table::DensePointToArmaVec(
+        right->bound().center(), &right_bound_ref);
+      bound_ref = left->count() * left_bound_ref +
+                  right->count() * right_bound_ref;
+      bound_ref =
+        (1.0 / static_cast<double>(node->count())) * bound_ref;
 
       double left_max_dist, right_max_dist;
       FurthestColumnIndex_(
diff --git a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h
index 79db4233de..03839fd44b 100644
--- a/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h
+++ b/fastlib/trunk/contrib/dongryel/thesis_research/mlpack/mixed_logit_dcm/dcm_table.h
@@ -6,10 +6,10 @@
 #ifndef MLPACK_MIXED_LOGIT_DCM_DCM_TABLE_H
 #define MLPACK_MIXED_LOGIT_DCM_DCM_TABLE_H
 
+#include <armadillo>
 #include <algorithm>
 #include <vector>
 #include "core/table/table.h"
-#include "core/math/linear_algebra.h"
 #include "core/monte_carlo/mean_variance_pair.h"
 #include "core/monte_carlo/mean_variance_pair_matrix.h"
 #include "mlpack/mixed_logit_dcm/mixed_logit_dcm_distribution.h"
@@ -92,8 +92,8 @@ class DCMTable {
      *         form.
      */
     void ComputeChoiceProbabilities_(
-      int person_index, const core::table::DensePoint &parameter_vector,
-      core::table::DensePoint *choice_probabilities) {
+      int person_index, const arma::vec &parameter_vector,
+      arma::vec *choice_probabilities) {
 
       int num_discrete_choices = this->num_discrete_choices(person_index);
       choice_probabilities->Init(num_discrete_choices);
@@ -109,8 +109,12 @@ class DCMTable {
         core::table::DensePoint attribute_for_discrete_choice;
         this->get_attribute_vector(
           person_index, discrete_choice_index, &attribute_for_discrete_choice);
-        double dot_product = core::math::Dot(
-                               parameter_vector, attribute_for_discrete_choice);
+        arma::vec attribute_for_discrete_choice_alias;
+        core::table::DensePointToArmaVec(
+          attribute_for_discrete_choice, attribute_for_discrete_choice_alias);
+        double dot_product =
+          arma::dot(
+            parameter_vector, attribute_for_discrete_choice_alias);
         double unnormalized_probability = exp(dot_product);
         normalizing_sum += unnormalized_probability;
         (*choice_probabilities)[discrete_choice_index] =
@@ -152,11 +156,11 @@ class DCMTable {
      *         computation of Equation 8.14 in the paper.
      */
     void SimulatedLoglikelihoodHessian(
-      core::table::DenseMatrix *likelihood_hessian) const {
+      arma::mat *likelihood_hessian) const {
 
-      likelihood_hessian->Init(
+      likelihood_hessian->set_size(
         distribution_->num_parameters(), distribution_->num_parameters());
-      likelihood_hessian->SetZero();
+      likelihood_hessian->zeros();
 
       // For each active person,
       for(int i = 0; i < num_active_people_; i++) {
@@ -178,23 +182,23 @@ class DCMTable {
           simulated_loglikelihood_hessians_[person_index].first;
         const core::monte_carlo::MeanVariancePairVector &hessian_second_part =
           simulated_loglikelihood_hessians_[person_index].second;
-        core::table::DenseMatrix hessian_first;
-        core::table::DensePoint hessian_second;
+        arma::mat hessian_first;
+        arma::vec hessian_second;
         hessian_first_part.sample_means(&hessian_first);
         hessian_second_part.sample_means(&hessian_second);
 
         // Construct the contribution on the fly.
-        core::math::AddExpert(
-          inverse_simulated_choice_probability,
-          hessian_first, likelihood_hessian);
-        core::math::MulExpert(
-          - core::math::Sqr(inverse_simulated_choice_probability),
-          hessian_second, hessian_second, likelihood_hessian);
+        (*likelihood_hessian) += inverse_simulated_choice_probability *
+                                 hessian_first;
+        (*likelihood_hessian) +=
+          (- core::math::Sqr(inverse_simulated_choice_probability)) *
+          hessian_second * arma::trans(hessian_second);
       }
 
       // Divide by the number of people.
-      core::math::Scale(
-        1.0 / static_cast<double>(num_active_people_), likelihood_hessian);
+      (*likelihood_hessian) =
+        (1.0 / static_cast<double>(num_active_people_)) *
+        (*likelihood_hessian);
     }
 
     /** @brief Return the gradient of the current simulated log
@@ -202,10 +206,10 @@ class DCMTable {
      *         in the paper.
      */
     void SimulatedLoglikelihoodGradient(
-      core::table::DensePoint *likelihood_gradient) const {
+      arma::vec *likelihood_gradient) const {
 
-      likelihood_gradient->Init(distribution_->num_parameters());
-      likelihood_gradient->SetZero();
+      likelihood_gradient->set_size(distribution_->num_parameters());
+      likelihood_gradient->zeros();
 
       // For each active person,
       for(int i = 0; i < num_active_people_; i++) {
@@ -230,14 +234,14 @@ class DCMTable {
 
         // Add the inverse probability weighted gradient vector for
         // the current person to the total tally.
-        core::math::AddExpert(
-          inverse_simulated_choice_probability,
-          gradient_vector, likelihood_gradient);
+        (*likelihood_gradient) +=
+          inverse_simulated_choice_probability * gradient_vector;
       }
 
       // Divide by the number of people.
-      core::math::Scale(
-        1.0 / static_cast<double>(num_active_people_), likelihood_gradient);
+      (*likelihood_gradient) =
+        (1.0 / static_cast<double>(num_active_people_)) *
+        (*likelihood_gradient);
     }
 
     /** @brief Return the current simulated log likelihood score.
@@ -378,8 +382,8 @@ class DCMTable {
 
       // Simulated log-likelihood gradient update by the simulated
       // choice probabilty scaled gradient product.
-      core::math::Scale(
-        choice_probabilities[discrete_choice_index], &beta_gradient_product);
+      beta_gradient_product = choice_probabilities[discrete_choice_index] *
+                              beta_gradient_product;
       simulated_loglikelihood_gradients_[person_index].push_back(
         beta_gradient_product);