diff --git a/doc/index.md b/doc/index.md
index 607a2ae214..cc71edfed4 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -123,6 +123,10 @@ Transform data from one space to another.
  * [`AMF`](user/methods/amf.md): alternating matrix factorization
  * [`LocalCoordinateCoding`](user/methods/local_coordinate_coding.md): local
    coordinate coding with dictionary learning
+ * [`LMNN`](user/methods/lmnn.md): large margin nearest neighbor (distance
+   metric learning)
+ * [`NCA`](user/methods/nca.md): neighborhood components analysis (distance
+   metric learning)
  * [`NMF`](user/methods/nmf.md): non-negative matrix factorization
  * [`PCA`](user/methods/pca.md): principal components analysis
  * [`SparseCoding`](user/methods/sparse_coding.md): sparse coding with
diff --git a/doc/sidebar.html b/doc/sidebar.html
index 568496310e..86bedf3e6d 100644
--- a/doc/sidebar.html
+++ b/doc/sidebar.html
@@ -170,6 +170,16 @@ when the sidebar is built for each page.
                 <code>LocalCoordinateCoding</code>
               </a>
             </li>
+            <li>
+              <a href="LINKROOTuser/methods/lmnn.html">
+                <code>LMNN</code>
+              </a>
+            </li>
+            <li>
+              <a href="LINKROOTuser/methods/nca.html">
+                <code>NCA</code>
+              </a>
+            </li>
             <li>
               <a href="LINKROOTuser/methods/nmf.html">
                 <code>NMF</code>
diff --git a/doc/user/bindings/cli.md b/doc/user/bindings/cli.md
index bc02653831..adfd721f22 100644
--- a/doc/user/bindings/cli.md
+++ b/doc/user/bindings/cli.md
@@ -1677,9 +1677,9 @@ $ mlpack_linear_svm --input_model_file lsvm_model.bin --test_file test.csv
 $ mlpack_lmnn [--batch_size 50] [--center] [--distance_file <string>]
         [--help] [--info <string>] --input_file <string> [--k 1] [--labels_file
         <string>] [--linear_scan] [--max_iterations 100000] [--normalize]
-        [--optimizer 'amsgrad'] [--passes 50] [--print_accuracy] [--range 1]
-        [--rank 0] [--regularization 0.5] [--seed 0] [--step_size 0.01]
-        [--tolerance 1e-07] [--verbose] [--version] [--centered_data_file
+        [--optimizer 'amsgrad'] [--passes 50] [--print_accuracy] [--rank 0]
+        [--regularization 0.5] [--seed 0] [--step_size 0.01] [--tolerance 1e-07]
+        [--update_interval 1] [--verbose] [--version] [--centered_data_file
         <string>] [--output_file <string>] [--transformed_data_file <string>]
 ```
 
@@ -1706,12 +1706,12 @@ An implementation of Large Margin Nearest Neighbors (LMNN), a distance learning
 | `--optimizer (-O)` | [`string`](#doc_string) | Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or 'lbfgs'. | `'amsgrad'` |
 | `--passes (-p)` | [`int`](#doc_int) | Maximum number of full passes over dataset for AMSGrad, BB_SGD and SGD. | `50` |
 | `--print_accuracy (-P)` | [`flag`](#doc_flag) | Print accuracies on initial and transformed dataset |  |
-| `--range (-R)` | [`int`](#doc_int) | Number of iterations after which impostors needs to be recalculated | `1` |
 | `--rank (-A)` | [`int`](#doc_int) | Rank of distance matrix to be optimized.  | `0` |
 | `--regularization (-r)` | [`double`](#doc_double) | Regularization for LMNN objective function  | `0.5` |
 | `--seed (-s)` | [`int`](#doc_int) | Random seed.  If 0, 'std::time(NULL)' is used. | `0` |
 | `--step_size (-a)` | [`double`](#doc_double) | Step size for AMSGrad, BB_SGD and SGD (alpha). | `0.01` |
 | `--tolerance (-t)` | [`double`](#doc_double) | Maximum tolerance for termination of AMSGrad, BB_SGD, SGD or L-BFGS. | `1e-07` |
+| `--update_interval (-R)` | [`int`](#doc_int) | Number of iterations after which impostors need to be recalculated. | `1` |
 | `--verbose (-v)` | [`flag`](#doc_flag) | Display informational messages and the full list of parameters and timers at the end of execution. |  |
 | `--version (-V)` | [`flag`](#doc_flag) | Display the version of mlpack.  <span class="special">Only exists in CLI binding.</span> |  |
 
@@ -1731,7 +1731,7 @@ This program implements Large Margin Nearest Neighbors, a distance learning tech
 
 To work, this algorithm needs labeled data.  It can be given as the last row of the input dataset (specified with `--input_file (-i)`), or alternatively as a separate matrix (specified with `--labels_file (-l)`).  Additionally, a starting point for optimization (specified with `--distance_file (-d)`can be given, having (r x d) dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank matrix will be optimized. Alternatively, Low-Rank distance can be learned by specifying the `--rank (-A)`parameter (A Low-Rank matrix with uniformly distributed values will be used as initial learning point). 
 
-The program also requires number of targets neighbors to work with ( specified with `--k (-k)`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `--regularization (-r)`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `--range (-R)`).
+The program also requires number of targets neighbors to work with ( specified with `--k (-k)`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `--regularization (-r)`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `--update_interval (-R)`).
 
 Output can either be the learned distance matrix (specified with `--output_file (-o)`), or the transformed dataset  (specified with `--transformed_data_file (-D)`), or both. Additionally mean-centered dataset (specified with `--centered_data_file (-c)`) can be accessed given mean-centering (specified with `--center (-C)`) is performed on the dataset. Accuracy on initial dataset and final transformed dataset can be printed by specifying the `--print_accuracy (-P)`parameter. 
 
@@ -1755,10 +1755,10 @@ $ mlpack_lmnn --input_file iris.csv --labels_file iris_labels.csv --k 3
   --optimizer bbsgd --output_file output.csv
 ```
 
-An another program call making use of range & regularization parameter with dataset having labels as last column can be made as: 
+Another program call making use of update interval & regularization parameter with dataset having labels as last column can be made as: 
 
 ```bash
-$ mlpack_lmnn --input_file letter_recognition.csv --k 5 --range 10
+$ mlpack_lmnn --input_file letter_recognition.csv --k 5 --update_interval 10
   --regularization 0.4 --output_file output.csv
 ```
 
diff --git a/doc/user/bindings/go.md b/doc/user/bindings/go.md
index c9baae764f..71878e2430 100644
--- a/doc/user/bindings/go.md
+++ b/doc/user/bindings/go.md
@@ -2069,12 +2069,12 @@ param.Normalize = false
 param.Optimizer = "amsgrad"
 param.Passes = 50
 param.PrintAccuracy = false
-param.Range = 1
 param.Rank = 0
 param.Regularization = 0.5
 param.Seed = 0
 param.StepSize = 0.01
 param.Tolerance = 1e-07
+param.UpdateInterval = 1
 param.Verbose = false
 
 centered_data, output, transformed_data := mlpack.Lmnn(input, param)
@@ -2102,12 +2102,12 @@ There are two types of input options: required options, which are passed directl
 | `Optimizer` | [`string`](#doc_string) | Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or 'lbfgs'. | `"amsgrad"` |
 | `Passes` | [`int`](#doc_int) | Maximum number of full passes over dataset for AMSGrad, BB_SGD and SGD. | `50` |
 | `PrintAccuracy` | [`bool`](#doc_bool) | Print accuracies on initial and transformed dataset | `false` |
-| `Range` | [`int`](#doc_int) | Number of iterations after which impostors needs to be recalculated | `1` |
 | `Rank` | [`int`](#doc_int) | Rank of distance matrix to be optimized.  | `0` |
 | `Regularization` | [`float64`](#doc_float64) | Regularization for LMNN objective function  | `0.5` |
 | `Seed` | [`int`](#doc_int) | Random seed.  If 0, 'std::time(NULL)' is used. | `0` |
 | `StepSize` | [`float64`](#doc_float64) | Step size for AMSGrad, BB_SGD and SGD (alpha). | `0.01` |
 | `Tolerance` | [`float64`](#doc_float64) | Maximum tolerance for termination of AMSGrad, BB_SGD, SGD or L-BFGS. | `1e-07` |
+| `UpdateInterval` | [`int`](#doc_int) | Number of iterations after which impostors need to be recalculated. | `1` |
 | `Verbose` | [`bool`](#doc_bool) | Display informational messages and the full list of parameters and timers at the end of execution. | `false` |
 
 ### Output options
@@ -2127,7 +2127,7 @@ This program implements Large Margin Nearest Neighbors, a distance learning tech
 
 To work, this algorithm needs labeled data.  It can be given as the last row of the input dataset (specified with `Input`), or alternatively as a separate matrix (specified with `Labels`).  Additionally, a starting point for optimization (specified with `Distance`can be given, having (r x d) dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank matrix will be optimized. Alternatively, Low-Rank distance can be learned by specifying the `Rank`parameter (A Low-Rank matrix with uniformly distributed values will be used as initial learning point). 
 
-The program also requires number of targets neighbors to work with ( specified with `K`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `Regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `Range`).
+The program also requires number of targets neighbors to work with ( specified with `K`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `Regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `UpdateInterval`).
 
 Output can either be the learned distance matrix (specified with `Output`), or the transformed dataset  (specified with `TransformedData`), or both. Additionally mean-centered dataset (specified with `CenteredData`) can be accessed given mean-centering (specified with `Center`) is performed on the dataset. Accuracy on initial dataset and final transformed dataset can be printed by specifying the `PrintAccuracy`parameter. 
 
@@ -2156,13 +2156,13 @@ param.Optimizer = "bbsgd"
 _, output, _ := mlpack.Lmnn(iris, param)
 ```
 
-An another program call making use of range & regularization parameter with dataset having labels as last column can be made as: 
+Another program call making use of update interval & regularization parameter with dataset having labels as last column can be made as: 
 
 ```go
 // Initialize optional parameters for Lmnn().
 param := mlpack.LmnnOptions()
 param.K = 5
-param.Range = 10
+param.UpdateInterval = 10
 param.Regularization = 0.4
 
 _, output, _ := mlpack.Lmnn(letter_recognition, param)
diff --git a/doc/user/bindings/julia.md b/doc/user/bindings/julia.md
index 5d9e699354..6b81648246 100644
--- a/doc/user/bindings/julia.md
+++ b/doc/user/bindings/julia.md
@@ -1696,9 +1696,9 @@ julia> using mlpack: lmnn
 julia> centered_data, output, transformed_data = lmnn(input;
           batch_size=50, center=false, distance=zeros(0, 0), k=1, labels=Int[],
           linear_scan=false, max_iterations=100000, normalize=false,
-          optimizer="amsgrad", passes=50, print_accuracy=false, range=1, rank=0,
+          optimizer="amsgrad", passes=50, print_accuracy=false, rank=0,
           regularization=0.5, seed=0, step_size=0.01, tolerance=1e-07,
-          verbose=false)
+          update_interval=1, verbose=false)
 ```
 
 An implementation of Large Margin Nearest Neighbors (LMNN), a distance learning technique.  Given a labeled dataset, this learns a transformation of the data that improves k-nearest-neighbor performance; this can be useful as a preprocessing step. [Detailed documentation](#lmnn_detailed-documentation).
@@ -1722,12 +1722,12 @@ An implementation of Large Margin Nearest Neighbors (LMNN), a distance learning
 | `optimizer` | [`String`](#doc_String) | Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or 'lbfgs'. | `"amsgrad"` |
 | `passes` | [`Int`](#doc_Int) | Maximum number of full passes over dataset for AMSGrad, BB_SGD and SGD. | `50` |
 | `print_accuracy` | [`Bool`](#doc_Bool) | Print accuracies on initial and transformed dataset | `false` |
-| `range` | [`Int`](#doc_Int) | Number of iterations after which impostors needs to be recalculated | `1` |
 | `rank` | [`Int`](#doc_Int) | Rank of distance matrix to be optimized.  | `0` |
 | `regularization` | [`Float64`](#doc_Float64) | Regularization for LMNN objective function  | `0.5` |
 | `seed` | [`Int`](#doc_Int) | Random seed.  If 0, 'std::time(NULL)' is used. | `0` |
 | `step_size` | [`Float64`](#doc_Float64) | Step size for AMSGrad, BB_SGD and SGD (alpha). | `0.01` |
 | `tolerance` | [`Float64`](#doc_Float64) | Maximum tolerance for termination of AMSGrad, BB_SGD, SGD or L-BFGS. | `1e-07` |
+| `update_interval` | [`Int`](#doc_Int) | Number of iterations after which impostors need to be recalculated. | `1` |
 | `verbose` | [`Bool`](#doc_Bool) | Display informational messages and the full list of parameters and timers at the end of execution. | `false` |
 
 ### Output options
@@ -1747,7 +1747,7 @@ This program implements Large Margin Nearest Neighbors, a distance learning tech
 
 To work, this algorithm needs labeled data.  It can be given as the last row of the input dataset (specified with `input`), or alternatively as a separate matrix (specified with `labels`).  Additionally, a starting point for optimization (specified with `distance`can be given, having (r x d) dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank matrix will be optimized. Alternatively, Low-Rank distance can be learned by specifying the `rank`parameter (A Low-Rank matrix with uniformly distributed values will be used as initial learning point). 
 
-The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `range`).
+The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `update_interval`).
 
 Output can either be the learned distance matrix (specified with `output`), or the transformed dataset  (specified with `transformed_data`), or both. Additionally mean-centered dataset (specified with `centered_data`) can be accessed given mean-centering (specified with `center`) is performed on the dataset. Accuracy on initial dataset and final transformed dataset can be printed by specifying the `print_accuracy`parameter. 
 
@@ -1774,13 +1774,13 @@ julia> _, output, _ = lmnn(iris; k=3, labels=iris_labels,
             optimizer="bbsgd")
 ```
 
-An another program call making use of range & regularization parameter with dataset having labels as last column can be made as: 
+Another program call making use of update interval & regularization parameter with dataset having labels as last column can be made as: 
 
 ```julia
 julia> using CSV
 julia> letter_recognition = CSV.read("letter_recognition.csv")
-julia> _, output, _ = lmnn(letter_recognition; k=5, range=10,
-            regularization=0.4)
+julia> _, output, _ = lmnn(letter_recognition; k=5,
+            regularization=0.4, update_interval=10)
 ```
 
 ### See also
diff --git a/doc/user/bindings/python.md b/doc/user/bindings/python.md
index a0a777fc72..8399b13531 100644
--- a/doc/user/bindings/python.md
+++ b/doc/user/bindings/python.md
@@ -1715,8 +1715,8 @@ Then, to use that model to predict classes for the dataset '`'test'`', storing t
         copy_all_inputs=False, distance=np.empty([0, 0]), input_=np.empty([0,
         0]), k=1, labels=np.empty([0], dtype=np.uint64), linear_scan=False,
         max_iterations=100000, normalize=False, optimizer='amsgrad', passes=50,
-        print_accuracy=False, range=1, rank=0, regularization=0.5, seed=0,
-        step_size=0.01, tolerance=1e-07, verbose=False)
+        print_accuracy=False, rank=0, regularization=0.5, seed=0,
+        step_size=0.01, tolerance=1e-07, update_interval=1, verbose=False)
 >>> centered_data = d['centered_data']
 >>> output = d['output']
 >>> transformed_data = d['transformed_data']
@@ -1744,12 +1744,12 @@ An implementation of Large Margin Nearest Neighbors (LMNN), a distance learning
 | `optimizer` | [`str`](#doc_str) | Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or 'lbfgs'. | `'amsgrad'` |
 | `passes` | [`int`](#doc_int) | Maximum number of full passes over dataset for AMSGrad, BB_SGD and SGD. | `50` |
 | `print_accuracy` | [`bool`](#doc_bool) | Print accuracies on initial and transformed dataset | `False` |
-| `range` | [`int`](#doc_int) | Number of iterations after which impostors needs to be recalculated | `1` |
 | `rank` | [`int`](#doc_int) | Rank of distance matrix to be optimized.  | `0` |
 | `regularization` | [`float`](#doc_float) | Regularization for LMNN objective function  | `0.5` |
 | `seed` | [`int`](#doc_int) | Random seed.  If 0, 'std::time(NULL)' is used. | `0` |
 | `step_size` | [`float`](#doc_float) | Step size for AMSGrad, BB_SGD and SGD (alpha). | `0.01` |
 | `tolerance` | [`float`](#doc_float) | Maximum tolerance for termination of AMSGrad, BB_SGD, SGD or L-BFGS. | `1e-07` |
+| `update_interval` | [`int`](#doc_int) | Number of iterations after which impostors need to be recalculated. | `1` |
 | `verbose` | [`bool`](#doc_bool) | Display informational messages and the full list of parameters and timers at the end of execution. | `False` |
 
 ### Output options
@@ -1769,7 +1769,7 @@ This program implements Large Margin Nearest Neighbors, a distance learning tech
 
 To work, this algorithm needs labeled data.  It can be given as the last row of the input dataset (specified with `input_`), or alternatively as a separate matrix (specified with `labels`).  Additionally, a starting point for optimization (specified with `distance`can be given, having (r x d) dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank matrix will be optimized. Alternatively, Low-Rank distance can be learned by specifying the `rank`parameter (A Low-Rank matrix with uniformly distributed values will be used as initial learning point). 
 
-The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `range`).
+The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `update_interval`).
 
 Output can either be the learned distance matrix (specified with `output`), or the transformed dataset  (specified with `transformed_data`), or both. Additionally mean-centered dataset (specified with `centered_data`) can be accessed given mean-centering (specified with `center`) is performed on the dataset. Accuracy on initial dataset and final transformed dataset can be printed by specifying the `print_accuracy`parameter. 
 
@@ -1793,10 +1793,10 @@ Example - Let's say we want to learn distance on iris dataset with number of tar
 >>> output = output['output']
 ```
 
-An another program call making use of range & regularization parameter with dataset having labels as last column can be made as: 
+Another program call making use of update interval & regularization parameter with dataset having labels as last column can be made as: 
 
 ```python
->>> output = lmnn(input_=letter_recognition, k=5, range=10,
+>>> output = lmnn(input_=letter_recognition, k=5, update_interval=10,
   regularization=0.4)
 >>> output = output['output']
 ```
diff --git a/doc/user/bindings/r.md b/doc/user/bindings/r.md
index 2aa0037959..b6a242b521 100644
--- a/doc/user/bindings/r.md
+++ b/doc/user/bindings/r.md
@@ -1689,9 +1689,9 @@ R> library(mlpack)
 R> d <- lmnn(batch_size=50, center=FALSE, distance=matrix(numeric(), 0,
         0), input=matrix(numeric(), 0, 0), k=1, labels=matrix(integer(), 0, 0),
         linear_scan=FALSE, max_iterations=100000, normalize=FALSE,
-        optimizer="amsgrad", passes=50, print_accuracy=FALSE, range=1, rank=0,
+        optimizer="amsgrad", passes=50, print_accuracy=FALSE, rank=0,
         regularization=0.5, seed=0, step_size=0.01, tolerance=1e-07,
-        verbose=getOption("mlpack.verbose", FALSE))
+        update_interval=1, verbose=getOption("mlpack.verbose", FALSE))
 R> centered_data <- d$centered_data
 R> output <- d$output
 R> transformed_data <- d$transformed_data
@@ -1718,12 +1718,12 @@ An implementation of Large Margin Nearest Neighbors (LMNN), a distance learning
 | `optimizer` | [`character`](#doc_character) | Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or 'lbfgs'. | `"amsgrad"` |
 | `passes` | [`integer`](#doc_integer) | Maximum number of full passes over dataset for AMSGrad, BB_SGD and SGD. | `50` |
 | `print_accuracy` | [`logical`](#doc_logical) | Print accuracies on initial and transformed dataset | `FALSE` |
-| `range` | [`integer`](#doc_integer) | Number of iterations after which impostors needs to be recalculated | `1` |
 | `rank` | [`integer`](#doc_integer) | Rank of distance matrix to be optimized.  | `0` |
 | `regularization` | [`numeric`](#doc_numeric) | Regularization for LMNN objective function  | `0.5` |
 | `seed` | [`integer`](#doc_integer) | Random seed.  If 0, 'std::time(NULL)' is used. | `0` |
 | `step_size` | [`numeric`](#doc_numeric) | Step size for AMSGrad, BB_SGD and SGD (alpha). | `0.01` |
 | `tolerance` | [`numeric`](#doc_numeric) | Maximum tolerance for termination of AMSGrad, BB_SGD, SGD or L-BFGS. | `1e-07` |
+| `update_interval` | [`integer`](#doc_integer) | Number of iterations after which impostors need to be recalculated. | `1` |
 | `verbose` | [`logical`](#doc_logical) | Display informational messages and the full list of parameters and timers at the end of execution. | `getOption("mlpack.verbose", FALSE)` |
 
 ### Output options
@@ -1743,7 +1743,7 @@ This program implements Large Margin Nearest Neighbors, a distance learning tech
 
 To work, this algorithm needs labeled data.  It can be given as the last row of the input dataset (specified with `input`), or alternatively as a separate matrix (specified with `labels`).  Additionally, a starting point for optimization (specified with `distance`can be given, having (r x d) dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank matrix will be optimized. Alternatively, Low-Rank distance can be learned by specifying the `rank`parameter (A Low-Rank matrix with uniformly distributed values will be used as initial learning point). 
 
-The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `range`).
+The program also requires number of targets neighbors to work with ( specified with `k`), A regularization parameter can also be passed, It acts as a trade of between the pulling and pushing terms (specified with `regularization`), In addition, this implementation of LMNN includes a parameter to decide the interval after which impostors must be re-calculated (specified with `update_interval`).
 
 Output can either be the learned distance matrix (specified with `output`), or the transformed dataset  (specified with `transformed_data`), or both. Additionally mean-centered dataset (specified with `centered_data`) can be accessed given mean-centering (specified with `center`) is performed on the dataset. Accuracy on initial dataset and final transformed dataset can be printed by specifying the `print_accuracy`parameter. 
 
@@ -1767,10 +1767,10 @@ R> output <- lmnn(input=iris, labels=iris_labels, k=3, optimizer="bbsgd")
 R> output <- output$output
 ```
 
-An another program call making use of range & regularization parameter with dataset having labels as last column can be made as: 
+Another program call making use of update interval & regularization parameter with dataset having labels as last column can be made as: 
 
 ```R
-R> output <- lmnn(input=letter_recognition, k=5, range=10,
+R> output <- lmnn(input=letter_recognition, k=5, update_interval=10,
   regularization=0.4)
 R> output <- output$output
 ```
diff --git a/doc/user/core.md b/doc/user/core.md
index 0caaea7459..c482688d39 100644
--- a/doc/user/core.md
+++ b/doc/user/core.md
@@ -862,9 +862,9 @@ including:
 
  * [`NeighborSearch`](/src/mlpack/methods/neighbor_search/neighbor_search.hpp)
  * [`RangeSearch`](/src/mlpack/methods/range_search/range_search.hpp)
- * [`LMNN`](/src/mlpack/methods/lmnn/lmnn.hpp)
+ * [`LMNN`](methods/lmnn.md)
  * [`EMST`](/src/mlpack/methods/emst/emst.hpp)
- * [`NCA`](/src/mlpack/methods/nca/nca.hpp)
+ * [`NCA`](methods/nca.md)
  * [`RANN`](/src/mlpack/methods/rann/rann.hpp)
  * [`KMeans`](/src/mlpack/methods/kmeans/kmeans.hpp)
 
@@ -1281,6 +1281,11 @@ std::cout << "Squared Mahalanobis distance on 32-bit floating point data:"
     << std::endl;
 std::cout << " - Points 3 and 5:   " << d1 << "." << std::endl;
 std::cout << " - Points 11 and 31: " << d2 << "." << std::endl;
+
+// Note that an equivalent transformation matrix can be recovered from Q with
+// an upper Cholesky decomposition (Q -> R.t() * R).
+arma::mat recoveredW = arma::chol(md.Q(), "lower");
+// A transformed dataset can be created with `(recoveredW * dataset)`.
 ```
 
 ---
diff --git a/doc/user/methods/lmnn.md b/doc/user/methods/lmnn.md
new file mode 100644
index 0000000000..322a5f683e
--- /dev/null
+++ b/doc/user/methods/lmnn.md
@@ -0,0 +1,454 @@
+## LMNN
+
+The `LMNN` class implements large margin nearest neighbor, which can be used
+as both a linear dimensionality reduction technique and a distance learning
+technique (also called metric learning).  LMNN finds a linear transformation of
+the dataset that improves `k`-nearest-neighbor classification performance.
+
+#### Simple usage example:
+
+```c++
+// Learn a distance metric that improves kNN classification performance.
+
+// All data and labels are uniform random; 10 dimensional data, 5 classes.
+// Replace with a data::Load() call or similar for a real application.
+arma::mat dataset(10, 1000, arma::fill::randu); // 1000 points.
+arma::Row<size_t> labels =
+    arma::randi<arma::Row<size_t>>(1000, arma::distr_param(0, 4));
+
+mlpack::LMNN lmnn(3 /* neighbors to consider */); // Step 1: create object.
+arma::mat distance;
+lmnn.LearnDistance(dataset, labels, distance);    // Step 2: learn distance.
+
+// `distance` can now be used as a transformation matrix for the data.
+arma::mat transformedData = distance * dataset;
+// Or, you can create a MahalanobisDistance to evaluate points in the
+// transformed dataset space.
+arma::mat q = distance.t() * distance;
+mlpack::MahalanobisDistance d(std::move(q));
+
+std::cout << "Distance between points 0 and 1:" << std::endl;
+std::cout << " - Before LMNN: "
+    << mlpack::EuclideanDistance::Evaluate(dataset.col(0), dataset.col(1))
+    << "." << std::endl;
+std::cout << " - After LMNN:  "
+    << d.Evaluate(dataset.col(0), dataset.col(1)) << "." << std::endl;
+```
+<p style="text-align: center; font-size: 85%"><a href="#simple-examples">More examples...</a></p>
+
+#### Quick links:
+
+ * [Constructors](#constructors): create `LMNN` objects.
+ * [`LearnDistance()`](#learning-distances): learn distance metrics.
+ * [Other functionality](#other-functionality) for loading and saving.
+ * [Examples](#simple-examples) of simple usage and integration with other
+   techniques.
+
+#### See also:
+
+<!-- TODO: link to kNN -->
+
+ * [mlpack distance metrics](../core.md#distances)
+ * [`NCA`](nca.md)
+ * [Metric learning on Wikipedia](https://en.wikipedia.org/wiki/Similarity_learning#Metric_learning)
+ * [Large margin nearest neighbor on Wikipedia](https://en.wikipedia.org/wiki/Large_margin_nearest_neighbor)
+ * [Distance metric learning for Large Margin Nearest Neighbor Classification (pdf)](https://proceedings.neurips.cc/paper_files/paper/2005/file/a7f592cef8b130a6967a90617db5681b-Paper.pdf)
+
+### Constructors
+
+ * `lmnn = LMNN(k, regularization=0.5, updateInterval=1)`
+   - Create an `LMNN` object considering the specified number `k` of neighbors.
+   - Optionally, specify the regularization to be applied to the LMNN cost
+     function (a `double`), and the number of iterations between recomputation
+     of neighbors (`updateInterval`, a `size_t`).
+
+---
+
+ * `lmnn = LMNN<DistanceType>(k, regularization=0.5, updateInterval=1)`
+ * `lmnn = LMNN<DistanceType>(k, regularization,     updateInterval,   distance)`
+   - Create an `LMNN` object using a custom
+     [`DistanceType`](../core.md#distances).
+   - `k` specifies the number of neighbors to consider.
+   - `regularization` specifies the regularization penalty to be applied to the
+     LMNN cost function (a `double`).
+   - `updateInterval` specifies the number of iterations between recomputation
+     of neighbors (a `size_t`).
+   - An instantiated `DistanceType` can optionally be passed with the `distance`
+     parameter.
+   - Using a custom `DistanceType` means that `LearnDistance()` will learn a
+     linear transformation for the data *in the metric space of the custom
+     `DistanceType`*.
+     * This means any learned distance may not necessarily improve
+       classification performance with the
+       [Euclidean distance](../core.md#lmetric).
+     * Instead, classification performance will be improved when the learned
+       distance is used with the given `DistanceType` only.
+   - Any mlpack `DistanceType` can be used as a drop-in replacement, or a
+     [custom `DistanceType`](../../developer/distances.md).
+     * A list of mlpack's provided distance metrics can be found
+       [here](../core.md#distances).
+   - ***Note: be sure that you understand the implications of a custom
+     `DistanceType` before using this version.***
+
+---
+
+***Notes***:
+
+ - A larger `k` will cause `LearnDistance()` to take longer to compute, but will
+   give more accurate results.  It is generally suggested to keep `k` in roughly
+   the `3` to `5` range, depending on the dataset.  Using `k = 1` can provide
+   fast convergence, but the learned distance metric may be of lower quality.
+
+ - `regularization` controls the balance between encouraging small distances for
+   points of the same class and penalizing small distances for points of
+   different classes.  When `regularization` is increased, small distances for
+   points of different classes are further penalized.
+
+ - Setting `updateInterval` greater than `1` will allow the LMNN algorithm to
+   take multiple steps without the expensive recomputation of neighbors, but
+   this means that subsequent optimization steps may not be using the true
+   nearest neighbors.
+   * If using an SGD-like algorithm (i.e. an optimizer for a
+     [differentiable separable function](https://www.ensmallen.org/docs.html#differentiable-separable-functions)),
+     this can often be set to a relatively high value (100 is not unreasonable).
+   * If using an optimizer like L-BFGS (i.e. a full-batch optimizer for
+     [differentiable functions](https://www.ensmallen.org/docs.html#differentiable-functions)),
+     this should be kept relatively low (going above 10 is not advised).
+   * It is worth cross-validating different values of the parameter to see what
+     works for your dataset.
+
+---
+
+### Learning Distances
+
+Once an `LMNN` object has been created, the `LearnDistance()` method can be used
+to learn a distance.
+
+ * `lmnn.LearnDistance(data, labels, distance,            [callbacks...])`
+ * `lmnn.LearnDistance(data, labels, distance, optimizer, [callbacks...])`
+   - Learn a distance metric on the given `data` and `labels`, filling
+     `distance` with a transformation matrix that can be used to map the data
+     into the space of the learned distance.
+   - Optionally, pass an instantiated
+     [ensmallen optimizer](https://www.ensmallen.org) and/or
+     [ensmallen callbacks](https://www.ensmallen.org/docs.html#callback-documentation)
+     to be used for the learning process.
+   - If no optimizer is passed,
+     [`ens::AMSGrad`](https://www.ensmallen.org/docs.html#amsgrad) is used.
+   - If `distance` already has size `r` x `data.n_rows` for some `r` less than
+     or equal to `data.n_rows`, it will be used as the starting point for
+     optimization.  Otherwise, the identity matrix with size `data.n_rows` x
+     `data.n_rows` will be used.
+   - When optimization is complete, `distance` will have size `r` x
+     `data.n_rows`, where `r` is less than or equal to `data.n_rows`.
+     * *Note*: If `r < data.n_rows`, then LMNN has learned a distance metric
+       that also reduces the dimensionality of the data.  See the
+       [last example](#simple-examples).
+
+To use `distance`, either:
+
+ * Compute a new transformed dataset as `distance * data`, or
+ * Use an instantiated [`MahalanobisDistance`](../core.md#mahalanobisdistance)
+   with `distance.t() * distance` as the `Q` matrix.
+
+See the [examples section](#simple-examples) for more details.
+
+#### `LearnDistance()` Parameters:
+
+| **name** | **type** | **description** |
+|----------|----------|-----------------|
+| `data` | [`arma::mat`](../matrices.md) | [Column-major](../matrices.md#representing-data-in-mlpack) training matrix. |
+| `labels` | [`arma::Row<size_t>`](../matrices.md) | Training labels, [between `0` and `numClasses - 1`](../load_save.md#normalizing-labels) (inclusive).  Should have length `data.n_cols`.  |
+| `distance` | [`arma::mat`](../matrices.md) | Output matrix to store transformation matrix representing learned distance. |
+| `optimizer` | [any ensmallen optimizer](https://www.ensmallen.org) | Instantiated ensmallen optimizer for [differentiable functions](https://www.ensmallen.org/docs.html#differentiable-functions) or [differentiable separable functions](https://www.ensmallen.org/docs.html#differentiable-separable-functions). | `ens::AMSGrad()` |
+| `callbacks...` | [any set of ensmallen callbacks](https://www.ensmallen.org/docs.html#callback-documentation) | Optional callbacks for the ensmallen optimizer, such as e.g. `ens::ProgressBar()`, `ens::Report()`, or others. | _(N/A)_ |
+
+***Note***: any matrix type can be used for `data` and `distance`, so long as
+that type implements the Armadillo API.  So, e.g., `arma::fmat` can be used.
+
+### Other Functionality
+
+ * An `LMNN` object can be serialized with
+   [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects).
+   Note that this is only meaningful if a custom `DistanceType` is being used,
+   and that custom `DistanceType` has state to be saved.
+
+ * `lmnn.K()` returns the number of neighbors used by LMNN, and `lmnn.K() = k`
+   will set the number of neighbors to use to `k`.
+
+ * `lmnn.Regularization()` returns the current regularization value of the LMNN
+   object (as a `double`), and `lmnn.Regularization() = r` can be used to set
+   the regularization value to `r`.
+
+ * `lmnn.UpdateInterval()` returns the current number of iterations between
+   neighbor recomputation (as a `size_t`), and `lmnn.UpdateInterval() = i` sets
+   the number of iterations between neighbor recomputation to `i`.
+
+ * `lmnn.Distance()` will return the `DistanceType` being used for learning.
+   Unless a custom `DistanceType` was specified in the constructor,
+   this simply returns a [`SquaredEuclideanDistance`](../core.md#lmetric)
+   object.
+
+### Simple Examples
+
+Learn a distance metric to improve classification performance on the iris
+dataset, and show improved performance when using
+[`NaiveBayesClassifier`](naive_bayes_classifier.md).
+
+```c++
+// See https://datasets.mlpack.org/satellite.test.csv.
+// (We are using the test set here just because it is a little smaller and
+// we want this example to run quickly.)
+arma::mat dataset;
+mlpack::data::Load("satellite.test.csv", dataset, true);
+// See https://datasets.mlpack.org/satellite.test.labels.csv.
+arma::Row<size_t> labels;
+mlpack::data::Load("satellite.test.labels.csv", labels, true);
+
+// Create an LMNN object using 5 nearest neighbors and learn a distance.
+arma::mat distance;
+mlpack::LMNN lmnn(5);
+lmnn.LearnDistance(dataset, labels, distance);
+
+// The distance matrix has size equal to the dimensionality of the data.
+std::cout << "Learned distance size: " << distance.n_rows << " x "
+    << distance.n_cols << "." << std::endl;
+
+// Learn a NaiveBayesClassifier model on the data and print the performance.
+mlpack::NaiveBayesClassifier nbc1(dataset, labels, 2);
+arma::Row<size_t> predictions;
+nbc1.Classify(dataset, predictions);
+std::cout << "Naive Bayes Classifier without LMNN: "
+    << arma::accu(labels == predictions) << " of " << labels.n_elem
+    << " correct." << std::endl;
+
+// Now transform the data and learn another NaiveBayesClassifier.
+arma::mat transformedDataset = distance * dataset;
+mlpack::NaiveBayesClassifier nbc2(transformedDataset, labels, 2);
+nbc2.Classify(transformedDataset, predictions);
+std::cout << "Naive Bayes Classifier with LMNN:    "
+    << arma::accu(labels == predictions) << " of " << labels.n_elem
+    << " correct." << std::endl;
+```
+
+---
+
+Learn a distance metric on the vehicle dataset, using 32-bit floating point to
+represent the data and metric.
+
+```c++
+// See https://datasets.mlpack.org/vehicle.csv.
+arma::fmat dataset;
+mlpack::data::Load("vehicle.csv", dataset, true);
+
+// The labels are contained as the last row of the dataset.
+arma::Row<size_t> labels =
+    arma::conv_to<arma::Row<size_t>>::from(dataset.row(dataset.n_rows - 1));
+dataset.shed_row(dataset.n_rows - 1);
+
+// Create an LMNN object with k=1 and learn distance on float32 data.
+// Set updateInterval to a large value (100) because we are using the default
+// AMSGrad optimizer (which will take very many small steps).
+arma::fmat distance;
+mlpack::LMNN lmnn(1, 0.5, 100);
+
+lmnn.LearnDistance(dataset, labels, distance, ens::ProgressBar());
+
+// We want to compute six quantities:
+//
+//  - Average distance to points of the same class before LMNN.
+//  - Average distance to points of the same class after LMNN, using
+//    MahalanobisDistance.
+//  - Average distance to points of the same class after LMNN, using the
+//    transformed dataset.
+//
+//  - The same three quantities above, but for points of the other class.
+//
+// LMNN should reduce the average distance to points in the same class, while
+// increasing the average distance to points in other classes.
+float distSums[6] = { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f };
+size_t sameCount = 0;
+arma::fmat q = distance.t() * distance;
+mlpack::MahalanobisDistance md(std::move(q));
+arma::fmat transformedDataset = distance * dataset;
+for (size_t i = 1; i < dataset.n_cols; ++i)
+{
+  const double d1 = mlpack::EuclideanDistance::Evaluate(
+      dataset.col(0), dataset.col(i));
+  const double d2 = md.Evaluate(dataset.col(0), dataset.col(i));
+  const double d3 = mlpack::EuclideanDistance::Evaluate(
+      transformedDataset.col(0), transformedDataset.col(i));
+
+  // Determine whether the point has the same label as point 0.
+  if (labels[i] == labels[0])
+  {
+    distSums[0] += d1;
+    distSums[1] += d2;
+    distSums[2] += d3;
+    ++sameCount;
+  }
+  else
+  {
+    distSums[3] += d1;
+    distSums[4] += d2;
+    distSums[5] += d3;
+  }
+}
+
+// Turn the results into average distances across the class.
+distSums[0] /= sameCount;
+distSums[1] /= sameCount;
+distSums[2] /= sameCount;
+distSums[3] /= (dataset.n_cols - sameCount);
+distSums[4] /= (dataset.n_cols - sameCount);
+distSums[5] /= (dataset.n_cols - sameCount);
+
+// Print the results.
+std::cout << "Average distance between point 0 and other points of the same "
+    << "class:" << std::endl;
+std::cout << " - Before LMNN:                           " << distSums[0] << "."
+    << std::endl;
+std::cout << " - After LMNN (with MahalanobisDistance): " << distSums[1] << "."
+    << std::endl;
+std::cout << " - After LMNN (with transformed dataset): " << distSums[2] << "."
+    << std::endl;
+std::cout << std::endl;
+
+std::cout << "Average distance between point 0 and points of other classes: "
+    << std::endl;
+std::cout << " - Before LMNN:                           " << distSums[3] << "."
+    << std::endl;
+std::cout << " - After LMNN (with MahalanobisDistance): " << distSums[4] << "."
+    << std::endl;
+std::cout << " - After LMNN (with transformed dataset): " << distSums[5] << "."
+    << std::endl;
+std::cout << std::endl;
+
+std::cout << "Ratio of other-class to same-class distances:" << std::endl;
+std::cout << "(We expect this to go up.)" << std::endl;
+std::cout << " - Before LMNN: " << (distSums[3] / distSums[0]) << "."
+    << std::endl;
+std::cout << " - After LMNN:  " << (distSums[5] / distSums[2]) << "."
+    << std::endl;
+```
+
+---
+
+Learn a distance metric on the iris dataset, using the L-BFGS optimizer with
+callbacks.
+
+```c++
+// See https://datasets.mlpack.org/iris.csv.
+arma::mat dataset;
+mlpack::data::Load("iris.csv", dataset, true);
+// See https://datasets.mlpack.org/iris.labels.csv.
+arma::Row<size_t> labels;
+mlpack::data::Load("iris.labels.csv", labels, true);
+
+// Learn a distance with ensmallen's L-BFGS optimizer.
+ens::L_BFGS lbfgs;
+lbfgs.NumBasis() = 5;
+lbfgs.MaxIterations() = 1000;
+
+// Use 5 neighbors for LMNN, and leave updateInterval at the default of 1,
+// because we are using L-BFGS (a full-back optimizer).
+mlpack::LMNN lmnn(5);
+
+// Use a callback that prints a final optimization report.
+arma::mat distance;
+lmnn.LearnDistance(dataset, labels, distance, lbfgs, ens::Report());
+```
+
+---
+
+Learn a distance metric on the vehicle dataset, but instead of using the
+Euclidean distance as the underlying metric, use the Manhattan distance.  This
+means that LMNN is optimizing k-NN performance under the Manhattan distance, not
+under the Euclidean distance.
+
+```c++
+// See https://datasets.mlpack.org/vehicle.csv.
+arma::mat dataset;
+mlpack::data::Load("vehicle.csv", dataset, true);
+
+// The labels are contained as the last row of the dataset.
+arma::Row<size_t> labels =
+    arma::conv_to<arma::Row<size_t>>::from(dataset.row(dataset.n_rows - 1));
+dataset.shed_row(dataset.n_rows - 1);
+
+// Create the LMNN object and optimize.  Use k=3 and Nesterov momentum SGD,
+// printing a progress bar during optimization.  Because Nesterov momentum SGD
+// is an ensmallen optimizer for differentiable separable functions, we increase
+// updateInterval to reduce the number of neighbor recomputations.  We also set
+// the regularization parameter to 1.0 to increase the penalty for nearby
+// neighbors of a different class.
+mlpack::LMNN<mlpack::ManhattanDistance> lmnn(3, 1.0, 100);
+arma::mat distance;
+ens::NesterovMomentumSGD opt(0.000001 /* step size */,
+                             32 /* batch size */,
+                             20 * dataset.n_cols /* 20 epochs */);
+lmnn.LearnDistance(dataset, labels, distance, opt, ens::ProgressBar());
+
+// Now inspect distances between points with the Euclidean distance and with the
+// inner product distance.
+arma::mat transformedDataset = distance * dataset;
+
+// Points 0 and 1 have the same label (0).  See their original distance---with
+// both the Euclidean and Manhattan distances---and their transformed distances.
+// We expect these points to get closer together, in the Manhattan distance.
+const double d1 = mlpack::ManhattanDistance::Evaluate(
+    dataset.col(0), dataset.col(1));
+const double d2 = mlpack::ManhattanDistance::Evaluate(
+    transformedDataset.col(0), transformedDataset.col(1));
+
+std::cout << "Distance between points 0 and 1 (same class):" << std::endl;
+std::cout << " - Manhattan distance:" << std::endl;
+std::cout << "   * Before LMNN: " << d1 << std::endl;
+std::cout << "   * After LMNN:  " << d2 << std::endl;
+std::cout << std::endl;
+
+// Point 3 has a different label.  We therefore expect this point to get further
+// from point 0 with the Manhattan distance, but not necessarily with the
+// Euclidean distance.
+const double d3 = mlpack::ManhattanDistance::Evaluate(
+    dataset.col(0), dataset.col(3));
+const double d4 = mlpack::ManhattanDistance::Evaluate(
+    transformedDataset.col(0), transformedDataset.col(3));
+
+std::cout << "Distance between points 0 and 3 (different class):" << std::endl;
+std::cout << " - Manhattan distance:" << std::endl;
+std::cout << "   * Before LMNN: " << d3 << std::endl;
+std::cout << "   * After LMNN:  " << d4 << std::endl;
+
+// Note that point 3 has been moved further away from point 0 than point 1.
+```
+
+---
+
+Learn a distance metric while also performing dimensionality reduction, reducing
+the dimensionality of the satellite dataset by 3 dimensions.
+
+```c++
+// See https://datasets.mlpack.org/satellite.train.csv.
+arma::mat dataset;
+mlpack::data::Load("satellite.train.csv", dataset, true);
+// See https://datasets.mlpack.org/satellite.labels.csv.
+arma::Row<size_t> labels;
+mlpack::data::Load("satellite.train.labels.csv", labels, true);
+
+// Use a random initialization for the distance transformation, with the
+// specified output dimensionality.
+arma::mat distance(dataset.n_rows - 3, dataset.n_rows, arma::fill::randu);
+mlpack::LMNN lmnn(3);
+ens::L_BFGS opt;
+opt.MaxIterations() = 10; // You may want more in a real application.
+lmnn.LearnDistance(dataset, labels, distance, opt, ens::Report());
+
+// Now transform the dataset.
+arma::mat transformedData = distance * dataset;
+
+std::cout << "Original data has size " << dataset.n_rows << " x "
+    << dataset.n_cols << "." << std::endl;
+std::cout << "Transformed data has size " << transformedData.n_rows << " x "
+    << transformedData.n_cols << "." << std::endl;
+```
diff --git a/doc/user/methods/nca.md b/doc/user/methods/nca.md
new file mode 100644
index 0000000000..4870907af8
--- /dev/null
+++ b/doc/user/methods/nca.md
@@ -0,0 +1,440 @@
+## NCA
+
+The `NCA` class implements neighborhood components analysis, which can be used
+as both a linear dimensionality reduction technique and a distance learning
+technique (also called metric learning).  Neighborhood components analysis finds
+a linear transformation of the dataset that improves `k`-nearest-neighbor
+classification performance.
+
+Note that `NCA` is a computationally intensive technique (each optimization
+iteration takes time quadratic in the data size!), and may be slow to run even
+for datasets of only moderate size.  See [`LMNN`](lmnn.md) for another distance
+learning technique that scales better to larger datasets.
+
+#### Simple usage example:
+
+```c++
+// Learn a distance metric that improves kNN classification performance.
+
+// All data and labels are uniform random; 10 dimensional data, 5 classes.
+// Replace with a data::Load() call or similar for a real application.
+arma::mat dataset(10, 1000, arma::fill::randu); // 1000 points.
+arma::Row<size_t> labels =
+    arma::randi<arma::Row<size_t>>(1000, arma::distr_param(0, 4));
+
+mlpack::NCA nca;                              // Step 1: create object.
+arma::mat distance;
+nca.LearnDistance(dataset, labels, distance); // Step 2: learn distance.
+
+// `distance` can now be used as a transformation matrix for the data.
+arma::mat transformedData = distance * dataset;
+// Or, you can create a MahalanobisDistance to evaluate points in the
+// transformed dataset space.
+arma::mat q = distance.t() * distance;
+mlpack::MahalanobisDistance d(std::move(q));
+
+std::cout << "Distance between points 0 and 1:" << std::endl;
+std::cout << " - Before NCA: "
+    << mlpack::EuclideanDistance::Evaluate(dataset.col(0), dataset.col(1))
+    << "." << std::endl;
+std::cout << " - After NCA:  "
+    << d.Evaluate(dataset.col(0), dataset.col(1)) << "." << std::endl;
+```
+<p style="text-align: center; font-size: 85%"><a href="#simple-examples">More examples...</a></p>
+
+#### Quick links:
+
+ * [Constructors](#constructors): create `NCA` objects.
+ * [`LearnDistance()`](#learning-distances): learn distance metrics.
+ * [Other functionality](#other-functionality) for loading and saving.
+ * [Examples](#simple-examples) of simple usage and integration with other
+   techniques.
+
+#### See also:
+
+<!-- TODO: link to kNN -->
+
+ * [mlpack distance metrics](../core.md#distances)
+ * [`LMNN`](lmnn.md)
+ * [Metric learning on Wikipedia](https://en.wikipedia.org/wiki/Similarity_learning#Metric_learning)
+ * [Neighborhood Components Analysis on Wikipedia](https://en.wikipedia.org/wiki/Neighbourhood_components_analysis)
+ * [Neighbourhood Components Analysis (pdf)](https://proceedings.neurips.cc/paper_files/paper/2004/file/42fe880812925e520249e808937738d2-Paper.pdf)
+
+### Constructors
+
+ * `nca = NCA()`
+   - Create an `NCA` object with default parameters.
+
+---
+
+ * `nca = NCA<DistanceType>()`
+ * `nca = NCA<DistanceType>(distance)`
+   - Create an `NCA` object using a custom
+     [`DistanceType`](../core.md#distances).
+   - An instantiated `DistanceType` can optionally be passed with the `distance`
+     parameter.
+   - Using a custom `DistanceType` means that `LearnDistance()` will learn a
+     linear transformation for the data *in the metric space of the custom
+     `DistanceType`*.
+     * This means any learned distance may not necessarily improve
+       classification performance with the
+       [Euclidean distance](../core.md#lmetric).
+     * Instead, classification performance will be improved when the learned
+       distance is used with the given `DistanceType` only.
+   - Any mlpack `DistanceType` can be used as a drop-in replacement, or a
+     [custom `DistanceType`](../../developer/distances.md).
+     * A list of mlpack's provided distance metrics can be found
+       [here](../core.md#distances).
+   - ***Note: be sure that you understand the implications of a custom
+     `DistanceType` before using this version.***
+
+---
+
+### Learning Distances
+
+Once an `NCA` object has been created, the `LearnDistance()` method can be used
+to learn a distance.
+
+ * `nca.LearnDistance(data, labels, distance,            [callbacks...])`
+ * `nca.LearnDistance(data, labels, distance, optimizer, [callbacks...])`
+   - Learn a distance metric on the given `data` and `labels`, filling
+     `distance` with a transformation matrix that can be used to map the data
+     into the space of the learned distance.
+   - Optionally, pass an instantiated
+     [ensmallen optimizer](https://www.ensmallen.org) and/or
+     [ensmallen callbacks](https://www.ensmallen.org/docs.html#callback-documentation)
+     to be used for the learning process.
+   - If `distance` already has size `r` x `data.n_rows` for some `r` less than
+     or equal to `data.n_rows`, it will be used as the starting point for
+     optimization.  Otherwise, the identity matrix with size `data.n_rows` x
+     `data.n_rows` will be used.
+   - When optimization is complete, `distance` will have size `r` x
+     `data.n_rows`, where `r` is less than or equal to `data.n_rows`.
+     * *Note*: If `r < data.n_rows`, then NCA has learned a distance metric that
+       also reduces the dimensionality of the data.  See the
+       [last example](#simple-examples).
+
+To use `distance`, either:
+
+ * Compute a new transformed dataset as `distance * data`, or
+ * Use an instantiated [`MahalanobisDistance`](../core.md#mahalanobisdistance)
+   with `distance.t() * distance` as the `Q` matrix.
+
+See the [examples section](#simple-examples) for more details.
+
+***Caveat:*** NCA operates by repeatedly computing expressions of the form
+`exp(-distance.Evaluate(data.col(i), data.col(j)))` (that is, the exponential of
+the negative distance between two points).  When distances are very large, this
+***quantity underflows to 0*** and results will not be reasonable.
+
+ - This situation can be detected, usually by a result where `distance` is equal
+   to the identity matrix.
+ - Alternately, if the [`ens::ProgressBar()`
+   callback](https://www.ensmallen.org/docs.html#progressbar) is used, a loss of
+   0 often means this situation has occurred.
+ - To mitigate the problem, consider scaling data such that the maximum pairwise
+   distance is less than 10.  See the [simple examples](#simple-examples) that
+   use the `vehicle` dataset.
+
+#### `LearnDistance()` Parameters:
+
+| **name** | **type** | **description** |
+|----------|----------|-----------------|
+| `data` | [`arma::mat`](../matrices.md) | [Column-major](../matrices.md#representing-data-in-mlpack) training matrix. |
+| `labels` | [`arma::Row<size_t>`](../matrices.md) | Training labels, [between `0` and `numClasses - 1`](../load_save.md#normalizing-labels) (inclusive).  Should have length `data.n_cols`.  |
+| `distance` | [`arma::mat`](../matrices.md) | Output matrix to store transformation matrix representing learned distance. |
+| `optimizer` | [any ensmallen optimizer](https://www.ensmallen.org) | Instantiated ensmallen optimizer for [differentiable functions](https://www.ensmallen.org/docs.html#differentiable-functions) or [differentiable separable functions](https://www.ensmallen.org/docs.html#differentiable-separable-functions). | `ens::StandardSGD()` |
+| `callbacks...` | [any set of ensmallen callbacks](https://www.ensmallen.org/docs.html#callback-documentation) | Optional callbacks for the ensmallen optimizer, such as e.g. `ens::ProgressBar()`, `ens::Report()`, or others. | _(N/A)_ |
+
+***Note***: any matrix type can be used for `data` and `distance`, so long as
+that type implements the Armadillo API.  So, e.g., `arma::fmat` can be used.
+
+### Other Functionality
+
+ * An `NCA` object can be serialized with
+   [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects).
+   Note that this is only meaningful if a custom `DistanceType` is being used,
+   and that custom `DistanceType` has state to be saved.
+
+ * `nca.Distance()` will return the `DistanceType` being used for learning.
+   Unless a custom `DistanceType` was specified in the constructor,
+   this simply returns a [`SquaredEuclideanDistance`](../core.md#lmetric)
+   object.
+
+### Simple Examples
+
+Learn a distance metric to improve classification performance on the iris
+dataset, and show improved performance when using
+[`NaiveBayesClassifier`](naive_bayes_classifier.md).
+
+```c++
+// See https://datasets.mlpack.org/iris.csv.
+arma::mat dataset;
+mlpack::data::Load("iris.csv", dataset, true);
+// See https://datasets.mlpack.org/iris.labels.csv.
+arma::Row<size_t> labels;
+mlpack::data::Load("iris.labels.csv", labels, true);
+
+// Create an NCA object and learn a distance.
+arma::mat distance;
+mlpack::NCA nca;
+nca.LearnDistance(dataset, labels, distance);
+
+// The distance matrix has size equal to the dimensionality of the data.
+std::cout << "Learned distance size: " << distance.n_rows << " x "
+    << distance.n_cols << "." << std::endl;
+
+// Learn a NaiveBayesClassifier model on the data and print the performance.
+mlpack::NaiveBayesClassifier nbc1(dataset, labels, 3);
+arma::Row<size_t> predictions;
+nbc1.Classify(dataset, predictions);
+std::cout << "Naive Bayes Classifier without NCA: "
+    << arma::accu(labels == predictions) << " of " << labels.n_elem
+    << " correct." << std::endl;
+
+// Now transform the data and learn another NaiveBayesClassifier.
+arma::mat transformedDataset = distance * dataset;
+mlpack::NaiveBayesClassifier nbc2(transformedDataset, labels, 3);
+nbc2.Classify(transformedDataset, predictions);
+std::cout << "Naive Bayes Classifier with NCA:    "
+    << arma::accu(labels == predictions) << " of " << labels.n_elem
+    << " correct." << std::endl;
+```
+
+---
+
+Learn a distance metric on the ionosphere dataset, using 32-bit floating point
+to represent the data and metric.
+
+```c++
+// See https://datasets.mlpack.org/ionosphere.csv.
+arma::fmat dataset;
+mlpack::data::Load("ionosphere.csv", dataset, true);
+
+// The labels are the last row of the dataset.
+arma::Row<size_t> labels =
+    arma::conv_to<arma::Row<size_t>>::from(dataset.row(dataset.n_rows - 1));
+dataset.shed_row(dataset.n_rows - 1);
+
+// Create an NCA object and learn distance on float32 data.
+// To keep computation time down, we use an instantiated optimizer that will
+// only perform 10 epochs of training.  (In a real application you may want to
+// train for longer!)
+arma::fmat distance;
+mlpack::NCA nca;
+
+ens::StandardSGD opt;
+opt.MaxIterations() = 10 * dataset.n_cols;
+nca.LearnDistance(dataset, labels, distance, opt, ens::ProgressBar());
+
+// We want to compute six quantities:
+//
+//  - Average distance to points of the same class before NCA.
+//  - Average distance to points of the same class after NCA, using
+//    MahalanobisDistance.
+//  - Average distance to points of the same class after NCA, using the
+//    transformed dataset.
+//
+//  - The same three quantities above, but for points of the other class.
+//
+// NCA should reduce the average distance to points in the same class, while
+// increasing the average distance to points in other classes.
+float distSums[6] = { 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f };
+size_t sameCount = 0;
+arma::fmat q = distance.t() * distance;
+mlpack::MahalanobisDistance md(std::move(q));
+arma::fmat transformedDataset = distance * dataset;
+for (size_t i = 1; i < dataset.n_cols; ++i)
+{
+  const double d1 = mlpack::EuclideanDistance::Evaluate(
+      dataset.col(0), dataset.col(i));
+  const double d2 = md.Evaluate(dataset.col(0), dataset.col(i));
+  const double d3 = mlpack::EuclideanDistance::Evaluate(
+      transformedDataset.col(0), transformedDataset.col(i));
+
+  // Determine whether the point has the same label as point 0.
+  if (labels[i] == labels[0])
+  {
+    distSums[0] += d1;
+    distSums[1] += d2;
+    distSums[2] += d3;
+    ++sameCount;
+  }
+  else
+  {
+    distSums[3] += d1;
+    distSums[4] += d2;
+    distSums[5] += d3;
+  }
+}
+
+// Turn the results into average distances across the class.
+distSums[0] /= sameCount;
+distSums[1] /= sameCount;
+distSums[2] /= sameCount;
+distSums[3] /= (dataset.n_cols - sameCount);
+distSums[4] /= (dataset.n_cols - sameCount);
+distSums[5] /= (dataset.n_cols - sameCount);
+
+// Print the results.
+std::cout << "Average distance between point 0 and other points of the same "
+    << "class:" << std::endl;
+std::cout << " - Before NCA:                           " << distSums[0] << "."
+    << std::endl;
+std::cout << " - After NCA (with MahalanobisDistance): " << distSums[1] << "."
+    << std::endl;
+std::cout << " - After NCA (with transformed dataset): " << distSums[2] << "."
+    << std::endl;
+std::cout << std::endl;
+
+std::cout << "Average distance between point 0 and points of other classes: "
+    << std::endl;
+std::cout << " - Before NCA:                           " << distSums[3] << "."
+    << std::endl;
+std::cout << " - After NCA (with MahalanobisDistance): " << distSums[4] << "."
+    << std::endl;
+std::cout << " - After NCA (with transformed dataset): " << distSums[5] << "."
+    << std::endl;
+std::cout << std::endl;
+
+std::cout << "Ratio of other-class to same-class distances:" << std::endl;
+std::cout << "(We expect this to go up.)" << std::endl;
+std::cout << " - Before NCA: " << (distSums[3] / distSums[0]) << "."
+    << std::endl;
+std::cout << " - After NCA:  " << (distSums[5] / distSums[2]) << "."
+    << std::endl;
+```
+
+---
+
+Learn a distance metric on the iris dataset, using the L-BFGS optimizer with
+callbacks.
+
+```c++
+// See https://datasets.mlpack.org/iris.csv.
+arma::mat dataset;
+mlpack::data::Load("iris.csv", dataset, true);
+// See https://datasets.mlpack.org/iris.labels.csv.
+arma::Row<size_t> labels;
+mlpack::data::Load("iris.labels.csv", labels, true);
+
+// Learn a distance with ensmallen's L-BFGS optimizer.
+ens::L_BFGS lbfgs;
+lbfgs.NumBasis() = 5;
+lbfgs.MaxIterations() = 1000;
+
+arma::mat distance;
+mlpack::NCA nca;
+
+// Use a callback that prints a final optimization report.
+nca.LearnDistance(dataset, labels, distance, lbfgs, ens::Report());
+```
+
+---
+
+<!-- TODO: actually use a kNN classifier here... once we have it implemented! -->
+
+Learn a distance metric on the vehicle dataset, but instead of using the
+Euclidean distance as the underlying metric, use the Manhattan distance.  This
+means that NCA is optimizing k-NN performance under the Manhattan distance, not
+under the Euclidean distance.
+
+```c++
+// See https://datasets.mlpack.org/vehicle.csv.
+arma::mat dataset;
+mlpack::data::Load("vehicle.csv", dataset, true);
+
+// The labels are contained as the last row of the dataset.
+arma::Row<size_t> labels =
+    arma::conv_to<arma::Row<size_t>>::from(dataset.row(dataset.n_rows - 1));
+dataset.shed_row(dataset.n_rows - 1);
+
+// Because typical distances between points in the vehicle dataset are large,
+// we will center the dataset and scale it to have points in the unit ball.
+// (That is, all points will have values in each dimension between -1 and 1.)
+// This means that the maximum pairwise distance is 2.
+dataset.each_col() -= arma::mean(dataset, 1);
+dataset /= arma::max(arma::max(arma::abs(dataset)));
+
+// Create the NCA object and optimize.  Use Nesterov momentum SGD, printing a
+// progress bar during optimization.
+mlpack::NCA<mlpack::ManhattanDistance> nca;
+arma::mat distance;
+ens::NesterovMomentumSGD opt(0.01 /* step size */,
+                             32 /* batch size */,
+                             20 * dataset.n_cols /* 20 epochs */);
+nca.LearnDistance(dataset, labels, distance, opt, ens::ProgressBar());
+
+// Now inspect distances between points with the Euclidean distance and with the
+// inner product distance.
+arma::mat transformedDataset = distance * dataset;
+
+// Points 0 and 1 have the same label (0).  See their original distance---with
+// both the Euclidean and Manhattan distances---and their transformed distances.
+// We expect these points to get closer together, in the Manhattan distance.
+const double d1 = mlpack::ManhattanDistance::Evaluate(
+    dataset.col(0), dataset.col(1));
+const double d2 = mlpack::ManhattanDistance::Evaluate(
+    transformedDataset.col(0), transformedDataset.col(1));
+
+std::cout << "Distance between points 0 and 1 (same class):" << std::endl;
+std::cout << " - Manhattan distance:" << std::endl;
+std::cout << "   * Before NCA: " << d1 << std::endl;
+std::cout << "   * After NCA:  " << d2 << std::endl;
+std::cout << std::endl;
+
+// Point 3 has a different label.  We therefore expect this point to get further
+// from point 0 with the Manhattan distance, but not necessarily with the
+// Euclidean distance.
+const double d3 = mlpack::ManhattanDistance::Evaluate(
+    dataset.col(0), dataset.col(3));
+const double d4 = mlpack::ManhattanDistance::Evaluate(
+    transformedDataset.col(0), transformedDataset.col(3));
+
+std::cout << "Distance between points 0 and 3 (different class):" << std::endl;
+std::cout << " - Manhattan distance:" << std::endl;
+std::cout << "   * Before NCA: " << d3 << std::endl;
+std::cout << "   * After NCA:  " << d4 << std::endl;
+
+// Note that point 3 has been moved further away from point 0 than point 1.
+```
+
+---
+
+Learn a distance metric while also performing dimensionality reduction, reducing
+the dimensionality of the vehicle dataset by 2 dimensions.
+
+```c++
+// See https://datasets.mlpack.org/vehicle.csv.
+arma::mat dataset;
+mlpack::data::Load("vehicle.csv", dataset, true);
+
+// The labels are contained as the last row of the dataset.
+arma::Row<size_t> labels =
+    arma::conv_to<arma::Row<size_t>>::from(dataset.row(dataset.n_rows - 1));
+dataset.shed_row(dataset.n_rows - 1);
+
+// Because typical distances between points in the vehicle dataset are large,
+// we will center the dataset and scale it to have points in the unit ball.
+// (That is, all points will have values in each dimension between -1 and 1.)
+// This means that the maximum pairwise distance is 2.
+dataset.each_col() -= arma::mean(dataset, 1);
+dataset /= arma::max(arma::max(arma::abs(dataset)));
+
+// Use a random initialization for the distance transformation, with the
+// specified output dimensionality.
+arma::mat distance(dataset.n_rows - 2, dataset.n_rows, arma::fill::randu);
+mlpack::NCA nca;
+ens::L_BFGS opt;
+opt.MaxIterations() = 10; // You may want more in a real application.
+nca.LearnDistance(dataset, labels, distance, opt);
+
+// Now transform the dataset.
+arma::mat transformedData = distance * dataset;
+
+std::cout << std::endl << std::endl;
+std::cout << "Original data has size " << dataset.n_rows << " x "
+    << dataset.n_cols << "." << std::endl;
+std::cout << "Transformed data has size " << transformedData.n_rows << " x "
+    << transformedData.n_cols << "." << std::endl;
+```
diff --git a/src/mlpack/core.hpp b/src/mlpack/core.hpp
index 73eeaaece4..7205587ab5 100644
--- a/src/mlpack/core.hpp
+++ b/src/mlpack/core.hpp
@@ -37,6 +37,7 @@
 // Now the core mlpack classes.
 #include <mlpack/core/util/arma_traits.hpp>
 #include <mlpack/core/util/ens_traits.hpp>
+#include <mlpack/core/util/first_element_is_arma.hpp>
 #include <mlpack/core/util/using.hpp>
 #include <mlpack/core/util/conv_to.hpp>
 #include <mlpack/core/util/log.hpp>
diff --git a/src/mlpack/core/util/first_element_is_arma.hpp b/src/mlpack/core/util/first_element_is_arma.hpp
new file mode 100644
index 0000000000..60adcc0fe3
--- /dev/null
+++ b/src/mlpack/core/util/first_element_is_arma.hpp
@@ -0,0 +1,44 @@
+/**
+ * @file core/util/first_element_is_arma.hpp
+ * @author Ryan Curtin
+ *
+ * Utility struct to detect whether the first element in a parameter pack is an
+ * Armadillo type.
+ */
+#ifndef MLPACK_CORE_UTIL_FIRST_ELEMENT_IS_ARMA_HPP
+#define MLPACK_CORE_UTIL_FIRST_ELEMENT_IS_ARMA_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+
+// This utility struct returns the first type of a parameter pack.
+template<typename... CallbackTypes>
+struct First
+{
+  typedef void type;
+};
+
+// This matches whenever CallbackTypes has one or more elements.
+template<typename T, typename... CallbackTypes>
+struct First<T, CallbackTypes...>
+{
+  typedef T type;
+};
+
+// This utility template struct detects whether the first element in a
+// parameter pack is an Armadillo type.  It is entirely for the deprecated
+// constructor below and can be removed when that is removed during the
+// release of mlpack 5.0.0.
+template<typename... CallbackTypes>
+struct FirstElementIsArma
+{
+  static constexpr bool value = arma::is_arma_type<
+      typename std::remove_reference<
+          typename First<CallbackTypes...>::type
+      >::type>::value;
+};
+
+}
+
+#endif
diff --git a/src/mlpack/methods/lmnn/constraints.hpp b/src/mlpack/methods/lmnn/constraints.hpp
index 082b8961ed..7704d239fe 100644
--- a/src/mlpack/methods/lmnn/constraints.hpp
+++ b/src/mlpack/methods/lmnn/constraints.hpp
@@ -27,12 +27,25 @@ namespace mlpack {
  * data point) and Triplets() (Generates sets of {dataset, target neighbors,
  * impostors} tripltets.)
  */
-template<typename DistanceType = SquaredEuclideanDistance>
+template<typename MatType = arma::mat,
+         typename LabelsType = arma::Row<size_t>,
+         typename DistanceType = SquaredEuclideanDistance>
 class Constraints
 {
  public:
   //! Convenience typedef.
-  typedef NeighborSearch<NearestNeighborSort, DistanceType> KNN;
+  typedef NeighborSearch<NearestNeighborSort, DistanceType, MatType> KNN;
+
+  // Convenience typedef for element type of data.
+  typedef typename MatType::elem_type ElemType;
+  // Convenience typedef for column vector of data.
+  typedef typename GetColType<MatType>::type VecType;
+  // Convenience typedef for cube of data.
+  typedef typename GetCubeType<MatType>::type CubeType;
+  // Convenience typedef for dense matrix of indices.
+  typedef typename GetUDenseMatType<MatType>::type UMatType;
+  // Convenience typedef for dense vector of indices.
+  typedef typename GetColType<UMatType>::type UVecType;
 
   /**
    * Constructor for creating a Constraints instance.
@@ -41,8 +54,8 @@ class Constraints
    * @param labels Input dataset labels.
    * @param k Number of target neighbors, impostors & triplets.
    */
-  Constraints(const arma::mat& dataset,
-              const arma::Row<size_t>& labels,
+  Constraints(const MatType& dataset,
+              const LabelsType& labels,
               const size_t k);
 
   /**
@@ -54,10 +67,10 @@ class Constraints
    * @param labels Input dataset labels.
    * @param norms Input dataset norms.
    */
-  void TargetNeighbors(arma::Mat<size_t>& outputMatrix,
-                       const arma::mat& dataset,
-                       const arma::Row<size_t>& labels,
-                       const arma::vec& norms);
+  void TargetNeighbors(UMatType& outputMatrix,
+                       const MatType& dataset,
+                       const LabelsType& labels,
+                       const VecType& norms);
 
   /**
    * Calculates k similar labeled nearest neighbors for a batch of dataset and
@@ -70,10 +83,10 @@ class Constraints
    * @param begin Index of the initial point of dataset.
    * @param batchSize Number of data points to use.
    */
-  void TargetNeighbors(arma::Mat<size_t>& outputMatrix,
-                       const arma::mat& dataset,
-                       const arma::Row<size_t>& labels,
-                       const arma::vec& norms,
+  void TargetNeighbors(UMatType& outputMatrix,
+                       const MatType& dataset,
+                       const LabelsType& labels,
+                       const VecType& norms,
                        const size_t begin,
                        const size_t batchSize);
 
@@ -86,10 +99,10 @@ class Constraints
    * @param labels Input dataset labels.
    * @param norms Input dataset norms.
    */
-  void Impostors(arma::Mat<size_t>& outputMatrix,
-                 const arma::mat& dataset,
-                 const arma::Row<size_t>& labels,
-                 const arma::vec& norms);
+  void Impostors(UMatType& outputMatrix,
+                 const MatType& dataset,
+                 const LabelsType& labels,
+                 const VecType& norms);
 
   /**
    * Calculates k differently labeled nearest neighbors & distances to
@@ -101,11 +114,11 @@ class Constraints
    * @param labels Input dataset labels.
    * @param norms Input dataset norms.
    */
-  void Impostors(arma::Mat<size_t>& outputNeighbors,
-                 arma::mat& outputDistance,
-                 const arma::mat& dataset,
-                 const arma::Row<size_t>& labels,
-                 const arma::vec& norms);
+  void Impostors(UMatType& outputNeighbors,
+                 MatType& outputDistance,
+                 const MatType& dataset,
+                 const LabelsType& labels,
+                 const VecType& norms);
 
   /**
    * Calculates k differently labeled nearest neighbors for a batch of dataset
@@ -118,10 +131,10 @@ class Constraints
    * @param begin Index of the initial point of dataset.
    * @param batchSize Number of data points to use.
    */
-  void Impostors(arma::Mat<size_t>& outputMatrix,
-                 const arma::mat& dataset,
-                 const arma::Row<size_t>& labels,
-                 const arma::vec& norms,
+  void Impostors(UMatType& outputMatrix,
+                 const MatType& dataset,
+                 const LabelsType& labels,
+                 const VecType& norms,
                  const size_t begin,
                  const size_t batchSize);
 
@@ -137,11 +150,11 @@ class Constraints
    * @param begin Index of the initial point of dataset.
    * @param batchSize Number of data points to use.
    */
-  void Impostors(arma::Mat<size_t>& outputNeighbors,
-                 arma::mat& outputDistance,
-                 const arma::mat& dataset,
-                 const arma::Row<size_t>& labels,
-                 const arma::vec& norms,
+  void Impostors(UMatType& outputNeighbors,
+                 MatType& outputDistance,
+                 const MatType& dataset,
+                 const LabelsType& labels,
+                 const VecType& norms,
                  const size_t begin,
                  const size_t batchSize);
 
@@ -158,12 +171,12 @@ class Constraints
    * @param points Indices of data points to calculate impostors on.
    * @param numPoints Number of points to actually calculate impostors on.
    */
-  void Impostors(arma::Mat<size_t>& outputNeighbors,
-                 arma::mat& outputDistance,
-                 const arma::mat& dataset,
-                 const arma::Row<size_t>& labels,
-                 const arma::vec& norms,
-                 const arma::uvec& points,
+  void Impostors(UMatType& outputNeighbors,
+                 MatType& outputDistance,
+                 const MatType& dataset,
+                 const LabelsType& labels,
+                 const VecType& norms,
+                 const UVecType& points,
                  const size_t numPoints);
 
   /**
@@ -175,10 +188,10 @@ class Constraints
    * @param labels Input dataset labels.
    * @param norms Input dataset norms.
    */
-  void Triplets(arma::Mat<size_t>& outputMatrix,
-                const arma::mat& dataset,
-                const arma::Row<size_t>& labels,
-                const arma::vec& norms);
+  void Triplets(UMatType& outputMatrix,
+                const MatType& dataset,
+                const LabelsType& labels,
+                const VecType& norms);
 
   //! Get the number of target neighbors (k).
   const size_t& K() const { return k; }
@@ -195,13 +208,13 @@ class Constraints
   size_t k;
 
   //! Store unique labels.
-  arma::Row<size_t> uniqueLabels;
+  LabelsType uniqueLabels;
 
   //! Store indices of data points having similar label.
-  std::vector<arma::uvec> indexSame;
+  std::vector<UVecType> indexSame;
 
   //! Store indices of data points having different label.
-  std::vector<arma::uvec> indexDiff;
+  std::vector<UVecType> indexDiff;
 
   //! False if nothing has ever been precalculated.
   bool precalculated;
@@ -210,15 +223,15 @@ class Constraints
   * Precalculate the unique labels, and indices of similar
   * and different datapoints on the basis of labels.
   */
-  inline void Precalculate(const arma::Row<size_t>& labels);
+  inline void Precalculate(const LabelsType& labels);
 
   /**
   * Re-order neighbors on the basis of increasing norm in case
   * of ties among distances.
   */
-  inline void ReorderResults(const arma::mat& distances,
-                             arma::Mat<size_t>& neighbors,
-                             const arma::vec& norms);
+  inline void ReorderResults(const MatType& distances,
+                             UMatType& neighbors,
+                             const VecType& norms);
 };
 
 } // namespace mlpack
diff --git a/src/mlpack/methods/lmnn/constraints_impl.hpp b/src/mlpack/methods/lmnn/constraints_impl.hpp
index 6f227fd970..27ed28417a 100644
--- a/src/mlpack/methods/lmnn/constraints_impl.hpp
+++ b/src/mlpack/methods/lmnn/constraints_impl.hpp
@@ -17,10 +17,10 @@
 
 namespace mlpack {
 
-template<typename DistanceType>
-Constraints<DistanceType>::Constraints(
-    const arma::mat& /* dataset */,
-    const arma::Row<size_t>& labels,
+template<typename MatType, typename LabelsType, typename DistanceType>
+Constraints<MatType, LabelsType, DistanceType>::Constraints(
+    const MatType& /* dataset */,
+    const LabelsType& labels,
     const size_t k) :
     k(k),
     precalculated(false)
@@ -36,11 +36,11 @@ Constraints<DistanceType>::Constraints(
   }
 }
 
-template<typename DistanceType>
-inline void Constraints<DistanceType>::ReorderResults(
-    const arma::mat& distances,
-    arma::Mat<size_t>& neighbors,
-    const arma::vec& norms)
+template<typename MatType, typename LabelsType, typename DistanceType>
+inline void Constraints<MatType, LabelsType, DistanceType>::ReorderResults(
+    const MatType& distances,
+    UMatType& neighbors,
+    const VecType& norms)
 {
   // Shortcut...
   if (neighbors.n_rows == 1)
@@ -64,24 +64,21 @@ inline void Constraints<DistanceType>::ReorderResults(
       if (start != end)
       {
         // We must sort these elements by norm.
-        arma::Col<size_t> newNeighbors =
-            neighbors.col(i).subvec(start, end - 1);
-        arma::uvec indices = ConvTo<arma::uvec>::From(newNeighbors);
-
-        arma::uvec order = arma::sort_index(norms.elem(indices));
-        neighbors.col(i).subvec(start, end - 1) =
-            newNeighbors.elem(order);
+        UVecType indices = neighbors.col(i).subvec(start, end - 1);
+        UVecType order = arma::sort_index(norms.elem(indices));
+        neighbors.col(i).subvec(start, end - 1) = indices.elem(order);
       }
     }
   }
 }
 
 // Calculates k similar labeled nearest neighbors.
-template<typename DistanceType>
-void Constraints<DistanceType>::TargetNeighbors(arma::Mat<size_t>& outputMatrix,
-                                                const arma::mat& dataset,
-                                                const arma::Row<size_t>& labels,
-                                                const arma::vec& norms)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::TargetNeighbors(
+    UMatType& outputMatrix,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
@@ -89,8 +86,8 @@ void Constraints<DistanceType>::TargetNeighbors(arma::Mat<size_t>& outputMatrix,
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -114,28 +111,29 @@ void Constraints<DistanceType>::TargetNeighbors(arma::Mat<size_t>& outputMatrix,
 
 // Calculates k similar labeled nearest neighbors  on a
 // batch of data points.
-template<typename DistanceType>
-void Constraints<DistanceType>::TargetNeighbors(arma::Mat<size_t>& outputMatrix,
-                                                const arma::mat& dataset,
-                                                const arma::Row<size_t>& labels,
-                                                const arma::vec& norms,
-                                                const size_t begin,
-                                                const size_t batchSize)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::TargetNeighbors(
+    UMatType& outputMatrix,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms,
+    const size_t begin,
+    const size_t batchSize)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
 
-  arma::mat subDataset = dataset.cols(begin, begin + batchSize - 1);
-  arma::Row<size_t> sublabels = labels.cols(begin, begin + batchSize - 1);
+  MatType subDataset = dataset.cols(begin, begin + batchSize - 1);
+  LabelsType sublabels = labels.cols(begin, begin + batchSize - 1);
 
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   // Vectors to store indices.
-  arma::uvec subIndexSame;
+  UVecType subIndexSame;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -161,11 +159,12 @@ void Constraints<DistanceType>::TargetNeighbors(arma::Mat<size_t>& outputMatrix,
 }
 
 // Calculates k differently labeled nearest neighbors.
-template<typename DistanceType>
-void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputMatrix,
-                                          const arma::mat& dataset,
-                                          const arma::Row<size_t>& labels,
-                                          const arma::vec& norms)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Impostors(
+    UMatType& outputMatrix,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
@@ -173,8 +172,8 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputMatrix,
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -198,12 +197,13 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputMatrix,
 
 // Calculates k differently labeled nearest neighbors. The function
 // writes back calculated neighbors & distances to passed matrices.
-template<typename DistanceType>
-void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
-                                          arma::mat& outputDistance,
-                                          const arma::mat& dataset,
-                                          const arma::Row<size_t>& labels,
-                                          const arma::vec& norms)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Impostors(
+    UMatType& outputNeighbors,
+    MatType& outputDistance,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
@@ -211,8 +211,8 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -237,28 +237,29 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
 
 // Calculates k differently labeled nearest neighbors on a
 // batch of data points.
-template<typename DistanceType>
-void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputMatrix,
-                                          const arma::mat& dataset,
-                                          const arma::Row<size_t>& labels,
-                                          const arma::vec& norms,
-                                          const size_t begin,
-                                          const size_t batchSize)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Impostors(
+    UMatType& outputMatrix,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms,
+    const size_t begin,
+    const size_t batchSize)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
 
-  arma::mat subDataset = dataset.cols(begin, begin + batchSize - 1);
-  arma::Row<size_t> sublabels = labels.cols(begin, begin + batchSize - 1);
+  MatType subDataset = dataset.cols(begin, begin + batchSize - 1);
+  LabelsType sublabels = labels.cols(begin, begin + batchSize - 1);
 
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   // Vectors to store indices.
-  arma::uvec subIndexSame;
+  UVecType subIndexSame;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -285,29 +286,30 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputMatrix,
 
 // Calculates k differently labeled nearest neighbors & distances on a
 // batch of data points.
-template<typename DistanceType>
-void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
-                                          arma::mat& outputDistance,
-                                          const arma::mat& dataset,
-                                          const arma::Row<size_t>& labels,
-                                          const arma::vec& norms,
-                                          const size_t begin,
-                                          const size_t batchSize)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Impostors(
+    UMatType& outputNeighbors,
+    MatType& outputDistance,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms,
+    const size_t begin,
+    const size_t batchSize)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
 
-  arma::mat subDataset = dataset.cols(begin, begin + batchSize - 1);
-  arma::Row<size_t> sublabels = labels.cols(begin, begin + batchSize - 1);
+  MatType subDataset = dataset.cols(begin, begin + batchSize - 1);
+  LabelsType sublabels = labels.cols(begin, begin + batchSize - 1);
 
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   // Vectors to store indices.
-  arma::uvec subIndexSame;
+  UVecType subIndexSame;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -335,14 +337,15 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
 
 // Calculates k differently labeled nearest neighbors & distances over some
 // data points.
-template<typename DistanceType>
-void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
-                                          arma::mat& outputDistance,
-                                          const arma::mat& dataset,
-                                          const arma::Row<size_t>& labels,
-                                          const arma::vec& norms,
-                                          const arma::uvec& points,
-                                          const size_t numPoints)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Impostors(
+    UMatType& outputNeighbors,
+    MatType& outputDistance,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms,
+    const UVecType& points,
+    const size_t numPoints)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
@@ -350,11 +353,11 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
   // KNN instance.
   KNN knn;
 
-  arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  UMatType neighbors;
+  MatType distances;
 
   // Vectors to store indices.
-  arma::uvec subIndexSame;
+  UVecType subIndexSame;
 
   for (size_t i = 0; i < uniqueLabels.n_cols; ++i)
   {
@@ -384,31 +387,35 @@ void Constraints<DistanceType>::Impostors(arma::Mat<size_t>& outputNeighbors,
 
 // Generates {data point, target neighbors, impostors} triplets using
 // TargetNeighbors() and Impostors().
-template<typename DistanceType>
-void Constraints<DistanceType>::Triplets(arma::Mat<size_t>& outputMatrix,
-                                         const arma::mat& dataset,
-                                         const arma::Row<size_t>& labels,
-                                         const arma::vec& norms)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void Constraints<MatType, LabelsType, DistanceType>::Triplets(
+    UMatType& outputMatrix,
+    const MatType& dataset,
+    const LabelsType& labels,
+    const VecType& norms)
 {
   // Perform pre-calculation. If neccesary.
   Precalculate(labels);
 
   size_t N = dataset.n_cols;
 
-  arma::Mat<size_t> impostors(k, dataset.n_cols);
+  UMatType impostors(k, dataset.n_cols);
   Impostors(impostors, dataset, labels, norms);
 
-  arma::Mat<size_t> targetNeighbors(k, dataset.n_cols);;
+  UMatType targetNeighbors(k, dataset.n_cols);;
   TargetNeighbors(targetNeighbors, dataset, labels, norms);
 
-  outputMatrix = arma::Mat<size_t>(3, k * k * N , arma::fill::zeros);
+  outputMatrix = UMatType(3, k * k * N , arma::fill::zeros);
 
-  for (size_t i = 0, r = 0; i < N; ++i)
+  #pragma omp parallel for collapse(3)
+  for (size_t i = 0; i < N; ++i)
   {
     for (size_t j = 0; j < k; ++j)
     {
-      for (size_t l = 0; l < k; l++, r++)
+      for (size_t l = 0; l < k; l++)
       {
+        const size_t r = i * (k * k) + j * k + l;
+
         // Generate triplets.
         outputMatrix(0, r) = i;
         outputMatrix(1, r) = targetNeighbors(j, i);
@@ -418,9 +425,9 @@ void Constraints<DistanceType>::Triplets(arma::Mat<size_t>& outputMatrix,
   }
 }
 
-template<typename DistanceType>
-inline void Constraints<DistanceType>::Precalculate(
-    const arma::Row<size_t>& labels)
+template<typename MatType, typename LabelsType, typename DistanceType>
+inline void Constraints<MatType, LabelsType, DistanceType>::Precalculate(
+    const LabelsType& labels)
 {
   // Make sure the calculation is necessary.
   if (precalculated)
@@ -431,6 +438,7 @@ inline void Constraints<DistanceType>::Precalculate(
   indexSame.resize(uniqueLabels.n_elem);
   indexDiff.resize(uniqueLabels.n_elem);
 
+  #pragma omp parallel for
   for (size_t i = 0; i < uniqueLabels.n_elem; ++i)
   {
     // Store same and diff indices.
diff --git a/src/mlpack/methods/lmnn/lmnn.hpp b/src/mlpack/methods/lmnn/lmnn.hpp
index 92646df11a..8dadc5b31e 100644
--- a/src/mlpack/methods/lmnn/lmnn.hpp
+++ b/src/mlpack/methods/lmnn/lmnn.hpp
@@ -49,7 +49,7 @@ namespace mlpack {
  * @tparam OptimizerType Optimizer to use for developing distance.
  */
 template<typename DistanceType = SquaredEuclideanDistance,
-         typename OptimizerType = ens::AMSGrad>
+         typename DeprecatedOptimizerType = ens::AMSGrad>
 class LMNN
 {
  public:
@@ -63,11 +63,27 @@ class LMNN
    * @param k Number of targets to consider.
    * @param distance Type of distance metric used for computation.
    */
+  [[deprecated("Will be removed in mlpack 5.0.0.  Pass the dataset directly to "
+               "LearnDistance() instead.")]]
   LMNN(const arma::mat& dataset,
        const arma::Row<size_t>& labels,
        const size_t k,
        const DistanceType distance = DistanceType());
 
+  /**
+   * Construct the LMNN object, optionally with an instantiated distance metric.
+   *
+   * @param k Number of target neighbors to consider.
+   * @param regularization Penalty to apply to objective function.
+   * @param updateInterval Number of iterations between each recomputation of
+   *     true neighbors and impostors.
+   * @param distance Instantiated distance metric for computation.
+   */
+  LMNN(const size_t k,
+       const double regularization = 0.5,
+       const size_t updateInterval = 1,
+       DistanceType distance = DistanceType());
+
 
   /**
    * Perform Large Margin Nearest Neighbors metric learning. The output
@@ -80,25 +96,99 @@ class LMNN
    * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
    *      See https://www.ensmallen.org/docs.html#callback-documentation.
    */
-  template<typename... CallbackTypes>
+  template<typename... CallbackTypes,
+           typename = typename std::enable_if<IsEnsCallbackTypes<
+               CallbackTypes...
+           >::value>::type,
+           typename = typename std::enable_if<
+               !FirstElementIsArma<CallbackTypes...>::value
+           >::type>
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the version that takes a "
+               "dataset as a parameter.")]]
   void LearnDistance(arma::mat& outputMatrix, CallbackTypes&&... callbacks);
 
+  /**
+   * Perform Large Margin Nearest Neighbors metric learning. The output
+   * distance matrix is written into the passed reference. If the
+   * LearnDistance() is called with an outputMatrix with correct dimensions,
+   * then that matrix will be used as the starting point for optimization.
+   *
+   * @param dataset Dataset to learn distance metric on.
+   * @param labels Labels for dataset.
+   * @param outputMatrix Covariance matrix of Mahalanobis distance.
+   * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
+   *      See https://www.ensmallen.org/docs.html#callback-documentation.
+   */
+  template<typename MatType,
+           typename LabelsType,
+           typename... CallbackTypes,
+           typename = typename std::enable_if<!IsEnsOptimizer<
+               typename First<CallbackTypes...>::type,
+               LMNNFunction<MatType, LabelsType, DistanceType>,
+               MatType
+           >::value>::type,
+           typename = typename std::enable_if<IsEnsCallbackTypes<
+               CallbackTypes...
+           >::value>::type>
+  void LearnDistance(const MatType& dataset,
+                     const LabelsType& labels,
+                     MatType& outputMatrix,
+                     CallbackTypes&&... callbacks) const;
+
+  /**
+   * Perform Large Margin Nearest Neighbors metric learning. The output
+   * distance matrix is written into the passed reference. If the
+   * LearnDistance() is called with an outputMatrix with correct dimensions,
+   * then that matrix will be used as the starting point for optimization.
+   *
+   * @param dataset Dataset to learn distance metric on.
+   * @param labels Labels for dataset.
+   * @param optimizer Instantiated ensmallen optimizer to use for LMNN.
+   * @param outputMatrix Covariance matrix of Mahalanobis distance.
+   * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
+   *      See https://www.ensmallen.org/docs.html#callback-documentation.
+   */
+  template<typename MatType,
+           typename LabelsType,
+           typename OptimizerType,
+           typename... CallbackTypes,
+           typename = typename std::enable_if<IsEnsOptimizer<
+               OptimizerType,
+               LMNNFunction<MatType, LabelsType, DistanceType>,
+               MatType
+           >::value>::type>
+  void LearnDistance(const MatType& dataset,
+                     const LabelsType& labels,
+                     MatType& outputMatrix,
+                     OptimizerType& optimizer,
+                     CallbackTypes&&... callbacks) const;
 
   //! Get the dataset reference.
-  const arma::mat& Dataset() const { return dataset; }
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  const arma::mat& Dataset() const { return *dataset; }
 
   //! Get the labels reference.
-  const arma::Row<size_t>& Labels() const { return labels; }
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  const arma::Row<size_t>& Labels() const { return *labels; }
 
   //! Access the regularization value.
   const double& Regularization() const { return regularization; }
   //! Modify the regularization value.
   double& Regularization() { return regularization; }
 
-  //! Access the range value.
-  const size_t& Range() const { return range; }
-  //! Modify the range value.
-  size_t& Range() { return range; }
+  //! Access the iteration update interval value.
+  const size_t& UpdateInterval() const { return updateInterval; }
+  //! Modify the iteration update interval value.
+  size_t& UpdateInterval() { return updateInterval; }
+
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use UpdateInterval() "
+               "instead.")]]
+  const size_t& Range() const { return updateInterval; }
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use UpdateInterval() "
+               "instead.")]]
+  size_t& Range() { return updateInterval; }
 
   //! Access the value of k.
   const size_t& K() const { return k; }
@@ -106,15 +196,23 @@ class LMNN
   size_t K() { return k; }
 
   //! Get the optimizer.
-  const OptimizerType& Optimizer() const { return optimizer; }
-  OptimizerType& Optimizer() { return optimizer; }
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  const DeprecatedOptimizerType& Optimizer() const { return optimizer; }
+  //! Modify the optimizer.
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  DeprecatedOptimizerType& Optimizer() { return optimizer; }
+
+  // Serialize the LMNN object.
+  template<typename Archive>
+  void serialize(Archive& ar, const unsigned int /* version */);
 
  private:
-  //! Dataset reference.
-  const arma::mat& dataset;
-
-  //! Labels reference.
-  const arma::Row<size_t>& labels;
+  //! Dataset pointer (will be removed in mlpack 5.0.0).
+  const arma::mat* dataset;
+  //! Labels pointer (will be removed in mlpack 5.0.0).
+  const arma::Row<size_t>* labels;
 
   //! Number of target points.
   size_t k;
@@ -122,14 +220,14 @@ class LMNN
   //! Regularization value.
   double regularization;
 
-  //! Range after which impostors need to be recalculated.
-  size_t range;
+  //! Number of iterations after which impostors need to be recalculated.
+  size_t updateInterval;
 
   //! Distance to be used.
   DistanceType distance;
 
-  //! The optimizer to use.
-  OptimizerType optimizer;
+  //! The optimizer to use (will be removed in mlpack 5.0.0).
+  DeprecatedOptimizerType optimizer;
 }; // class LMNN
 
 } // namespace mlpack
diff --git a/src/mlpack/methods/lmnn/lmnn_function.hpp b/src/mlpack/methods/lmnn/lmnn_function.hpp
index f35fcf8cd0..e64b8ae0d0 100644
--- a/src/mlpack/methods/lmnn/lmnn_function.hpp
+++ b/src/mlpack/methods/lmnn/lmnn_function.hpp
@@ -41,9 +41,22 @@ namespace mlpack {
  * operate on one point in the dataset.  This is useful for optimizers like
  * stochastic gradient descent (see ens::SGD).
  */
-template<typename DistanceType = SquaredEuclideanDistance>
+template<typename MatType = arma::mat,
+         typename LabelsType = arma::Row<size_t>,
+         typename DistanceType = SquaredEuclideanDistance>
 class LMNNFunction
 {
+  // Convenience typedef for element type of data.
+  typedef typename MatType::elem_type ElemType;
+  // Convenience typedef for column vector of data.
+  typedef typename GetColType<MatType>::type VecType;
+  // Convenience typedef for cube of data.
+  typedef typename GetCubeType<MatType>::type CubeType;
+  // Convenience typedef for dense matrix of indices.
+  typedef typename GetUDenseMatType<MatType>::type UMatType;
+  // Convenience typedef for dense vector of indices.
+  typedef typename GetColType<UMatType>::type UVecType;
+
  public:
   /**
    * Constructor for LMNNFunction class.
@@ -52,14 +65,14 @@ class LMNNFunction
    * @param labels Input dataset labels.
    * @param k Number of target neighbors to be used.
    * @param regularization Regularization value.
-   * @param range Range after which impostors need to be recalculated.
+   * @param updateInterval Number of iterations before impostors are recomputed.
    * @param distance Type of distance metric used for computation.
    */
-  LMNNFunction(const arma::mat& dataset,
-               const arma::Row<size_t>& labels,
+  LMNNFunction(const MatType& dataset,
+               const LabelsType& labels,
                size_t k,
                double regularization,
-               size_t range,
+               size_t updateInterval,
                DistanceType distance = DistanceType());
 
 
@@ -69,13 +82,13 @@ class LMNNFunction
   void Shuffle();
 
   /**
-   * Evaluate the LMNN function for the given transformation matrix.  This is the
-   * non-separable implementation, where the objective function is not
+   * Evaluate the LMNN function for the given transformation matrix.  This is
+   * the non-separable implementation, where the objective function is not
    * decomposed into the sum of several objective functions.
    *
    * @param transformation Transformation matrix of Mahalanobis distance.
    */
-  double Evaluate(const arma::mat& transformation);
+  ElemType Evaluate(const MatType& transformation);
 
   /**
    * Evaluate the LMNN objective function for the given transformation matrix on
@@ -89,9 +102,9 @@ class LMNNFunction
    * @param begin Index of the initial point to use for objective function.
    * @param batchSize Number of points to use for objective function.
    */
-  double Evaluate(const arma::mat& transformation,
-                  const size_t begin,
-                  const size_t batchSize = 1);
+  ElemType Evaluate(const MatType& transformation,
+                    const size_t begin,
+                    const size_t batchSize = 1);
 
   /**
    * Evaluate the gradient of the LMNN function for the given transformation
@@ -103,7 +116,7 @@ class LMNNFunction
    * @param gradient Matrix to store the calculated gradient in.
    */
   template<typename GradType>
-  void Gradient(const arma::mat& transformation, GradType& gradient);
+  void Gradient(const MatType& transformation, GradType& gradient);
 
   /**
    * Evaluate the gradient of the LMNN function for the given transformation
@@ -121,7 +134,7 @@ class LMNNFunction
    * @param batchSize Number of points to use for objective function.
    */
   template<typename GradType>
-  void Gradient(const arma::mat& transformation,
+  void Gradient(const MatType& transformation,
                 const size_t begin,
                 GradType& gradient,
                 const size_t batchSize = 1);
@@ -137,8 +150,8 @@ class LMNNFunction
    * @param gradient Matrix to store the calculated gradient in.
    */
   template<typename GradType>
-  double EvaluateWithGradient(const arma::mat& transformation,
-                              GradType& gradient);
+  ElemType EvaluateWithGradient(const MatType& transformation,
+                                GradType& gradient);
 
   /**
    * Evaluate the LMNN objective function together with gradient for the given
@@ -156,13 +169,13 @@ class LMNNFunction
    * @param batchSize Number of points to use for objective function.
    */
   template<typename GradType>
-  double EvaluateWithGradient(const arma::mat& transformation,
-                              const size_t begin,
-                              GradType& gradient,
-                              const size_t batchSize = 1);
+  ElemType EvaluateWithGradient(const MatType& transformation,
+                                const size_t begin,
+                                GradType& gradient,
+                                const size_t batchSize = 1);
 
   //! Return the initial point for the optimization.
-  const arma::mat& GetInitialPoint() const { return initialPoint; }
+  const MatType& GetInitialPoint() const { return initialPoint; }
 
   /**
    * Get the number of functions the objective function can be decomposed into.
@@ -171,7 +184,7 @@ class LMNNFunction
   size_t NumFunctions() const { return dataset.n_cols; }
 
   //! Return the dataset passed into the constructor.
-  const arma::mat& Dataset() const { return dataset; }
+  const MatType& Dataset() const { return dataset; }
 
   //! Access the regularization value.
   const double& Regularization() const { return regularization; }
@@ -183,26 +196,26 @@ class LMNNFunction
   //! Modify the value of k.
   size_t& K() { return k; }
 
-  //! Access the value of range.
-  const size_t& Range() const { return range; }
-  //! Modify the value of k.
-  size_t& Range() { return range; }
+  //! Access the number of iterations between impostor recomputation.
+  const size_t& UpdateInterval() const { return updateInterval; }
+  //! Modify the number of iterations between impostor recomputation..
+  size_t& UpdateInterval() { return updateInterval; }
 
  private:
   //! data.  This will be an alias until Shuffle() is called.
-  arma::mat dataset;
+  MatType dataset;
   //! labels.  This will be an alias until Shuffle() is called.
-  arma::Row<size_t> labels;
+  LabelsType labels;
   //! Initial parameter point.
-  arma::mat initialPoint;
+  MatType initialPoint;
   //! Store transformed dataset.
-  arma::mat transformedDataset;
+  MatType transformedDataset;
   //! Store target neighbors of data points.
-  arma::Mat<size_t> targetNeighbors;
+  UMatType targetNeighbors;
   //! Initial impostors.
-  arma::Mat<size_t> impostors;
+  UMatType impostors;
   //! Cache distance. Used to avoid repetive calculation.
-  arma::mat distanceMat;
+  MatType distanceMat;
   //! Number of target neighbors.
   size_t k;
   //! The instantiated distance metric.
@@ -211,28 +224,28 @@ class LMNNFunction
   double regularization;
   //! Keep iterations count.
   size_t iteration;
-  //! Range after which impostors need to be recalculated.
-  size_t range;
+  //! Number of iterations before impostors need to be recalculated.
+  size_t updateInterval;
   //! Constraints Object.
-  Constraints<DistanceType> constraint;
+  Constraints<MatType, LabelsType, DistanceType> constraint;
   //! Holds pre-calculated cij.
-  arma::mat pCij;
+  MatType pCij;
   //! Holds the norm of each data point.
-  arma::vec norm;
+  VecType norm;
   //! Hold previous eval values for each datapoint.
-  arma::cube evalOld;
+  CubeType evalOld;
   //! Hold previous maximum norm of impostor.
-  arma::mat maxImpNorm;
+  MatType maxImpNorm;
   //! Holds previous transformation matrix. Used for L-BFGS like optimizer.
-  arma::mat transformationOld;
+  MatType transformationOld;
   //! Holds previous transformation matrices.
-  std::vector<arma::mat> oldTransformationMatrices;
+  std::vector<MatType> oldTransformationMatrices;
   //! Holds number of points which are using each transformation matrix.
   std::vector<size_t> oldTransformationCounts;
   //! Holds points to transformation matrix mapping.
-  arma::vec lastTransformationIndices;
+  VecType lastTransformationIndices;
   //! Used for storing points to re-calculate impostors for.
-  arma::uvec points;
+  UVecType points;
   //! Flag for controlling use of bounds over impostors.
   bool impBounds;
   /**
@@ -242,12 +255,12 @@ class LMNNFunction
   */
   inline void Precalculate();
   //! Update cache transformation matrices.
-  inline void UpdateCache(const arma::mat& transformation,
+  inline void UpdateCache(const MatType& transformation,
                           const size_t begin,
                           const size_t batchSize);
   //! Calculate norm of change in transformation.
-  inline void TransDiff(std::map<size_t, double>& transformationDiffs,
-                        const arma::mat& transformation,
+  inline void TransDiff(std::unordered_map<size_t, ElemType>& transDiffs,
+                        const MatType& transformation,
                         const size_t begin,
                         const size_t batchSize);
 };
diff --git a/src/mlpack/methods/lmnn/lmnn_function_impl.hpp b/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
index 01aa77afea..58dc54011c 100644
--- a/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
+++ b/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
@@ -18,18 +18,19 @@
 
 namespace mlpack {
 
-template<typename DistanceType>
-LMNNFunction<DistanceType>::LMNNFunction(const arma::mat& datasetIn,
-                                         const arma::Row<size_t>& labelsIn,
-                                         size_t k,
-                                         double regularization,
-                                         size_t range,
-                                         DistanceType distance) :
+template<typename MatType, typename LabelsType, typename DistanceType>
+LMNNFunction<MatType, LabelsType, DistanceType>::LMNNFunction(
+    const MatType& datasetIn,
+    const LabelsType& labelsIn,
+    size_t k,
+    double regularization,
+    size_t updateInterval,
+    DistanceType distance) :
     k(k),
     distance(distance),
     regularization(regularization),
     iteration(0),
-    range(range),
+    updateInterval(updateInterval),
     constraint(datasetIn, labelsIn, k),
     points(datasetIn.n_cols),
     impBounds(false)
@@ -60,7 +61,7 @@ LMNNFunction<DistanceType>::LMNNFunction(const arma::mat& datasetIn,
   lastTransformationIndices.zeros();
 
   // Reserve the first element of cache.
-  arma::mat emptyMat;
+  MatType emptyMat;
   oldTransformationMatrices.push_back(emptyMat);
   oldTransformationCounts.push_back(dataset.n_cols);
 
@@ -92,18 +93,18 @@ LMNNFunction<DistanceType>::LMNNFunction(const arma::mat& datasetIn,
 }
 
 //! Shuffle the dataset.
-template<typename DistanceType>
-void LMNNFunction<DistanceType>::Shuffle()
+template<typename MatType, typename LabelsType, typename DistanceType>
+void LMNNFunction<MatType, LabelsType, DistanceType>::Shuffle()
 {
-  arma::mat newDataset = dataset;
-  arma::Mat<size_t> newLabels = labels;
-  arma::cube newEvalOld = evalOld;
-  arma::vec newlastTransformationIndices = lastTransformationIndices;
-  arma::mat newMaxImpNorm = maxImpNorm;
-  arma::vec newNorm = norm;
+  MatType newDataset = dataset;
+  LabelsType newLabels = labels;
+  CubeType newEvalOld = evalOld;
+  VecType newlastTransformationIndices = lastTransformationIndices;
+  MatType newMaxImpNorm = maxImpNorm;
+  VecType newNorm = norm;
 
   // Generate ordering.
-  arma::uvec ordering = arma::shuffle(arma::linspace<arma::uvec>(0,
+  UVecType ordering = arma::shuffle(arma::linspace<UVecType>(0,
       dataset.n_cols - 1, dataset.n_cols));
 
   ClearAlias(dataset);
@@ -126,9 +127,9 @@ void LMNNFunction<DistanceType>::Shuffle()
 }
 
 // Update cache transformation matrices.
-template<typename DistanceType>
-inline void LMNNFunction<DistanceType>::UpdateCache(
-    const arma::mat& transformation,
+template<typename MatType, typename LabelsType, typename DistanceType>
+inline void LMNNFunction<MatType, LabelsType, DistanceType>::UpdateCache(
+    const MatType& transformation,
     const size_t begin,
     const size_t batchSize)
 {
@@ -162,31 +163,13 @@ inline void LMNNFunction<DistanceType>::UpdateCache(
   }
 
   oldTransformationCounts[index] += batchSize;
-
-  #ifdef DEBUG
-    size_t total = 0;
-    for (size_t i = 1; i < oldTransformationCounts.size(); ++i)
-    {
-      std::ostringstream oss;
-      oss << "transformation counts for matrix " << i
-          << " invalid (" << oldTransformationCounts[i] << ")!";
-      Log::Assert(oldTransformationCounts[i] <= dataset.n_cols, oss.str());
-      total += oldTransformationCounts[i];
-    }
-
-    std::ostringstream oss;
-    oss << "total count for transformation matrices invalid (" << total
-        << ", " << "should be " << dataset.n_cols << "!";
-    if (begin + batchSize == dataset.n_cols)
-      Log::Assert(total == dataset.n_cols, oss.str());
-  #endif
 }
 
 // Calculate norm of change in transformation.
-template<typename DistanceType>
-inline void LMNNFunction<DistanceType>::TransDiff(
-    std::map<size_t, double>& transformationDiffs,
-    const arma::mat& transformation,
+template<typename MatType, typename LabelsType, typename DistanceType>
+inline void LMNNFunction<MatType, LabelsType, DistanceType>::TransDiff(
+    std::unordered_map<size_t, ElemType>& transformationDiffs,
+    const MatType& transformation,
     const size_t begin,
     const size_t batchSize)
 {
@@ -209,22 +192,24 @@ inline void LMNNFunction<DistanceType>::TransDiff(
 }
 
 //! Evaluate cost over whole dataset.
-template<typename DistanceType>
-double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
+template<typename MatType, typename LabelsType, typename DistanceType>
+typename MatType::elem_type
+LMNNFunction<MatType, LabelsType, DistanceType>::Evaluate(
+    const MatType& transformation)
 {
-  double cost = 0;
+  ElemType cost = 0;
 
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
-  double transformationDiff = 0;
+  ElemType transformationDiff = 0;
   if (!transformationOld.is_empty())
   {
     // Calculate norm of change in transformation.
     transformationDiff = arma::norm(transformation - transformationOld);
   }
 
-  if (!transformationOld.is_empty() && iteration++ % range == 0)
+  if (!transformationOld.is_empty() && iteration++ % updateInterval == 0)
   {
     if (impBounds)
     {
@@ -251,7 +236,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
           norm);
     }
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
     // Re-calculate impostors on transformed dataset.
     constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -263,7 +248,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate cost due to distance between target neighbors & data point.
-      double eval = distance.Evaluate(transformedDataset.col(i),
+      ElemType eval = distance.Evaluate(transformedDataset.col(i),
                           transformedDataset.col(targetNeighbors(j, i)));
       cost += (1 - regularization) * eval;
     }
@@ -276,7 +261,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (!transformationOld.is_empty() && evalOld(l, j, i) < -1)
@@ -292,7 +277,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -338,21 +323,23 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation)
 }
 
 //! Calculate cost over batches.
-template<typename DistanceType>
-double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
-                                            const size_t begin,
-                                            const size_t batchSize)
+template<typename MatType, typename LabelsType, typename DistanceType>
+typename MatType::elem_type
+LMNNFunction<MatType, LabelsType, DistanceType>::Evaluate(
+    const MatType& transformation,
+    const size_t begin,
+    const size_t batchSize)
 {
-  double cost = 0;
+  ElemType cost = 0;
 
   // Calculate norm of change in transformation.
-  std::map<size_t, double> transformationDiffs;
+  std::unordered_map<size_t, ElemType> transformationDiffs;
   TransDiff(transformationDiffs, transformation, begin, batchSize);
 
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
-  if (impBounds && iteration++ % range == 0)
+  if (impBounds && iteration++ % updateInterval == 0)
   {
     // Track number of data points to use for impostors calculatiom.
     size_t numPoints = 0;
@@ -378,7 +365,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
     constraint.Impostors(impostors, distanceMat,
         transformedDataset, labels, norm, points, numPoints);
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
     // Re-calculate impostors on transformed dataset.
     constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -390,7 +377,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate cost due to distance between target neighbors & data point.
-      double eval = distance.Evaluate(transformedDataset.col(i),
+      ElemType eval = distance.Evaluate(transformedDataset.col(i),
                           transformedDataset.col(targetNeighbors(j, i)));
       cost += (1 - regularization) * eval;
     }
@@ -403,7 +390,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (lastTransformationIndices(i) && evalOld(l, j, i) < -1)
@@ -419,7 +406,7 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -467,16 +454,16 @@ double LMNNFunction<DistanceType>::Evaluate(const arma::mat& transformation,
 }
 
 //! Compute gradient over whole dataset.
-template<typename DistanceType>
+template<typename MatType, typename LabelsType, typename DistanceType>
 template<typename GradType>
-void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
-                                          GradType& gradient)
+void LMNNFunction<MatType, LabelsType, DistanceType>::Gradient(
+    const MatType& transformation, GradType& gradient)
 {
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
-  double transformationDiff = 0;
-  if (!transformationOld.is_empty() && iteration++ % range == 0)
+  ElemType transformationDiff = 0;
+  if (!transformationOld.is_empty() && iteration++ % updateInterval == 0)
   {
     // Calculate norm of change in transformation.
     transformationDiff = arma::norm(transformation - transformationOld);
@@ -506,7 +493,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
           norm);
     }
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
       // Re-calculate impostors on transformed dataset.
       constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -516,10 +503,10 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
   gradient.zeros(transformation.n_rows, transformation.n_cols);
 
   // Calculate gradient due to target neighbors.
-  arma::mat cij = pCij;
+  MatType cij = pCij;
 
   // Calculate gradient due to impostors.
-  arma::mat cil = zeros(dataset.n_rows, dataset.n_rows);
+  MatType cil = zeros<MatType>(dataset.n_rows, dataset.n_rows);
 
   for (size_t i = 0; i < dataset.n_cols; ++i)
   {
@@ -530,7 +517,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (!transformationOld.is_empty() && evalOld(l, j, i) < -1)
@@ -546,7 +533,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -581,7 +568,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
         }
 
         // Caculate gradient due to impostors.
-        arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+        VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
         cil += diff * trans(diff);
 
         diff = dataset.col(i) - dataset.col(impostors(l, i));
@@ -598,21 +585,22 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
 }
 
 //! Compute gradient over a batch of data points.
-template<typename DistanceType>
+template<typename MatType, typename LabelsType, typename DistanceType>
 template<typename GradType>
-void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
-                                          const size_t begin,
-                                          GradType& gradient,
-                                          const size_t batchSize)
+void LMNNFunction<MatType, LabelsType, DistanceType>::Gradient(
+    const MatType& transformation,
+    const size_t begin,
+    GradType& gradient,
+    const size_t batchSize)
 {
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
   // Calculate norm of change in transformation.
-  std::map<size_t, double> transformationDiffs;
+  std::unordered_map<size_t, ElemType> transformationDiffs;
   TransDiff(transformationDiffs, transformation, begin, batchSize);
 
-  if (impBounds && iteration++ % range == 0)
+  if (impBounds && iteration++ % updateInterval == 0)
   {
     // Track number of data points to use for impostors calculatiom.
     size_t numPoints = 0;
@@ -638,7 +626,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
     constraint.Impostors(impostors, distanceMat,
         transformedDataset, labels, norm, points, numPoints);
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
     // Re-calculate impostors on transformed dataset.
     constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -647,15 +635,15 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
 
   gradient.zeros(transformation.n_rows, transformation.n_cols);
 
-  arma::mat cij = zeros(dataset.n_rows, dataset.n_rows);
-  arma::mat cil = zeros(dataset.n_rows, dataset.n_rows);
+  MatType cij = zeros<MatType>(dataset.n_rows, dataset.n_rows);
+  MatType cil = zeros<MatType>(dataset.n_rows, dataset.n_rows);
 
   for (size_t i = begin; i < begin + batchSize; ++i)
   {
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate gradient due to target neighbors.
-      arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+      VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
       cij += diff * trans(diff);
     }
 
@@ -666,7 +654,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (lastTransformationIndices(i) && evalOld(l, j, i) < -1)
@@ -682,7 +670,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -719,7 +707,7 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
         }
 
         // Caculate gradient due to impostors.
-        arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+        VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
         cil += diff * trans(diff);
 
         diff = dataset.col(i) - dataset.col(impostors(l, i));
@@ -736,25 +724,26 @@ void LMNNFunction<DistanceType>::Gradient(const arma::mat& transformation,
 }
 
 //! Compute cost & gradient over whole dataset.
-template<typename DistanceType>
+template<typename MatType, typename LabelsType, typename DistanceType>
 template<typename GradType>
-double LMNNFunction<DistanceType>::EvaluateWithGradient(
-    const arma::mat& transformation,
+typename MatType::elem_type
+LMNNFunction<MatType, LabelsType, DistanceType>::EvaluateWithGradient(
+    const MatType& transformation,
     GradType& gradient)
 {
-  double cost = 0;
+  ElemType cost = 0;
 
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
-  double transformationDiff = 0;
+  ElemType transformationDiff = 0;
   if (!transformationOld.is_empty())
   {
     // Calculate norm of change in transformation.
     transformationDiff = arma::norm(transformation - transformationOld);
   }
 
-  if (!transformationOld.is_empty() && iteration++ % range == 0)
+  if (!transformationOld.is_empty() && iteration++ % updateInterval == 0)
   {
     if (impBounds)
     {
@@ -781,7 +770,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
           norm);
     }
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
       // Re-calculate impostors on transformed dataset.
       constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -791,17 +780,17 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
   gradient.zeros(transformation.n_rows, transformation.n_cols);
 
   // Calculate gradient due to target neighbors.
-  arma::mat cij = pCij;
+  MatType cij = pCij;
 
   // Calculate gradient due to impostors.
-  arma::mat cil = zeros(dataset.n_rows, dataset.n_rows);
+  MatType cil = zeros<MatType>(dataset.n_rows, dataset.n_rows);
 
   for (size_t i = 0; i < dataset.n_cols; ++i)
   {
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate cost due to distance between target neighbors & data point.
-      double eval = distance.Evaluate(transformedDataset.col(i),
+      ElemType eval = distance.Evaluate(transformedDataset.col(i),
                         transformedDataset.col(targetNeighbors(j, i)));
       cost += (1 - regularization) * eval;
     }
@@ -813,7 +802,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (!transformationOld.is_empty() && evalOld(l, j, i) < -1)
@@ -829,7 +818,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -858,7 +847,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
         cost += regularization * (1 + eval);
 
         // Caculate gradient due to impostors.
-        arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+        VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
         cil += diff * trans(diff);
 
         diff = dataset.col(i) - dataset.col(impostors(l, i));
@@ -877,24 +866,25 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
 }
 
 //! Compute cost & gradient over a batch of data points.
-template<typename DistanceType>
+template<typename MatType, typename LabelsType, typename DistanceType>
 template<typename GradType>
-double LMNNFunction<DistanceType>::EvaluateWithGradient(
-    const arma::mat& transformation,
+typename MatType::elem_type
+LMNNFunction<MatType, LabelsType, DistanceType>::EvaluateWithGradient(
+    const MatType& transformation,
     const size_t begin,
     GradType& gradient,
     const size_t batchSize)
 {
-  double cost = 0;
+  ElemType cost = 0;
 
   // Calculate norm of change in transformation.
-  std::map<size_t, double> transformationDiffs;
+  std::unordered_map<size_t, ElemType> transformationDiffs;
   TransDiff(transformationDiffs, transformation, begin, batchSize);
 
   // Apply distance metric over dataset.
   transformedDataset = transformation * dataset;
 
-  if (impBounds && iteration++ % range == 0)
+  if (impBounds && iteration++ % updateInterval == 0)
   {
     // Track number of data points to use for impostors calculatiom.
     size_t numPoints = 0;
@@ -920,7 +910,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
     constraint.Impostors(impostors, distanceMat,
         transformedDataset, labels, norm, points, numPoints);
   }
-  else if (iteration++ % range == 0)
+  else if (iteration++ % updateInterval == 0)
   {
     // Re-calculate impostors on transformed dataset.
     constraint.Impostors(impostors, distanceMat, transformedDataset, labels,
@@ -929,20 +919,20 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
 
   gradient.zeros(transformation.n_rows, transformation.n_cols);
 
-  arma::mat cij = zeros(dataset.n_rows, dataset.n_rows);
-  arma::mat cil = zeros(dataset.n_rows, dataset.n_rows);
+  MatType cij = zeros<MatType>(dataset.n_rows, dataset.n_rows);
+  MatType cil = zeros<MatType>(dataset.n_rows, dataset.n_rows);
 
   for (size_t i = begin; i < begin + batchSize; ++i)
   {
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate cost due to distance between target neighbors & data point.
-      double eval = distance.Evaluate(transformedDataset.col(i),
+      ElemType eval = distance.Evaluate(transformedDataset.col(i),
                         transformedDataset.col(targetNeighbors(j, i)));
       cost += (1 - regularization) * eval;
 
       // Calculate gradient due to target neighbors.
-      arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+      VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
       cij += diff * trans(diff);
     }
 
@@ -953,7 +943,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
       {
         // Calculate cost due to {data point, target neighbors, impostors}
         // triplets.
-        double eval = 0;
+        ElemType eval = 0;
 
         // Bounds for eval.
         if (lastTransformationIndices(i) && evalOld(l, j, i) < -1)
@@ -969,7 +959,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
         // Calculate exact eval value.
         if (eval > -1)
         {
-          if (iteration - 1 % range == 0)
+          if (iteration - 1 % updateInterval == 0)
           {
             eval = distance.Evaluate(transformedDataset.col(i),
                      transformedDataset.col(targetNeighbors(j, i))) -
@@ -998,7 +988,7 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
         cost += regularization * (1 + eval);
 
         // Caculate gradient due to impostors.
-        arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+        VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
         cil += diff * trans(diff);
 
         diff = dataset.col(i) - dataset.col(impostors(l, i));
@@ -1016,8 +1006,8 @@ double LMNNFunction<DistanceType>::EvaluateWithGradient(
   return cost;
 }
 
-template<typename DistanceType>
-inline void LMNNFunction<DistanceType>::Precalculate()
+template<typename MatType, typename LabelsType, typename DistanceType>
+inline void LMNNFunction<MatType, LabelsType, DistanceType>::Precalculate()
 {
   pCij.zeros(dataset.n_rows, dataset.n_rows);
 
@@ -1026,7 +1016,7 @@ inline void LMNNFunction<DistanceType>::Precalculate()
     for (size_t j = 0; j < k ; ++j)
     {
       // Calculate gradient due to target neighbors.
-      arma::vec diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
+      VecType diff = dataset.col(i) - dataset.col(targetNeighbors(j, i));
       pCij += diff * trans(diff);
     }
   }
diff --git a/src/mlpack/methods/lmnn/lmnn_impl.hpp b/src/mlpack/methods/lmnn/lmnn_impl.hpp
index 740a3a6c0c..5dfc87bb78 100644
--- a/src/mlpack/methods/lmnn/lmnn_impl.hpp
+++ b/src/mlpack/methods/lmnn/lmnn_impl.hpp
@@ -21,27 +21,83 @@ namespace mlpack {
  * Takes in a reference to the dataset. Copies the data, initializes
  * all of the member variables and constraint object and generate constraints.
  */
-template<typename DistanceType, typename OptimizerType>
-LMNN<DistanceType, OptimizerType>::LMNN(const arma::mat& dataset,
-                       const arma::Row<size_t>& labels,
-                       const size_t k,
-                       const DistanceType distance) :
-    dataset(dataset),
-    labels(labels),
+template<typename DistanceType, typename DeprecatedOptimizerType>
+LMNN<DistanceType, DeprecatedOptimizerType>::LMNN(
+    const arma::mat& dataset,
+    const arma::Row<size_t>& labels,
+    const size_t k,
+    const DistanceType distance) :
+    dataset(&dataset),
+    labels(&labels),
     k(k),
     regularization(0.5),
-    range(1),
+    updateInterval(1),
     distance(distance)
 { /* nothing to do */ }
 
-template<typename DistanceType, typename OptimizerType>
-template<typename... CallbackTypes>
-void LMNN<DistanceType, OptimizerType>::LearnDistance(arma::mat& outputMatrix,
+template<typename DistanceType, typename DeprecatedOptimizerType>
+LMNN<DistanceType, DeprecatedOptimizerType>::LMNN(
+    const size_t k,
+    const double regularization,
+    const size_t updateInterval,
+    const DistanceType distance) :
+    k(k),
+    regularization(regularization),
+    updateInterval(updateInterval),
+    distance(distance)
+{ /* nothing to do */ }
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename... CallbackTypes, typename, typename>
+void LMNN<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    arma::mat& outputMatrix,
     CallbackTypes&&... callbacks)
+{
+  if (!dataset || !labels)
+  {
+    throw std::runtime_error("LMNN::LearnDistance(): cannot call without a "
+        "dataset!");
+  }
+
+  LearnDistance(*dataset, *labels, outputMatrix, optimizer,
+      std::forward<CallbackTypes>(callbacks)...);
+}
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename MatType,
+         typename LabelsType,
+         typename... CallbackTypes,
+         typename /* SFINAE check that first callback is not an optimizer */,
+         typename /* callback SFINAE check */>
+void LMNN<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    const MatType& dataset,
+    const LabelsType& labels,
+    MatType& outputMatrix,
+    CallbackTypes&&... callbacks) const
+{
+  // This should be replaced with ens::StandardSGD when the deprecated members
+  // are removed for mlpack 5.0.0.
+  DeprecatedOptimizerType opt;
+  LearnDistance(dataset, labels, outputMatrix, opt,
+      std::forward<CallbackTypes>(callbacks)...);
+}
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename MatType,
+         typename LabelsType,
+         typename OptimizerType,
+         typename... CallbackTypes,
+         typename /* SFINAE check that opt is an ensmallen optimizer */>
+void LMNN<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    const MatType& dataset,
+    const LabelsType& labels,
+    MatType& outputMatrix,
+    OptimizerType& opt,
+    CallbackTypes&&... callbacks) const
 {
   // LMNN objective function.
-  LMNNFunction<DistanceType> objFunction(dataset, labels, k,
-      regularization, range);
+  LMNNFunction<MatType, LabelsType, DistanceType> objFunction(dataset, labels,
+      k, regularization, updateInterval);
 
   // See if we were passed an initialized matrix. outputMatrix (L) must be
   // having r x d dimensionality.
@@ -49,15 +105,23 @@ void LMNN<DistanceType, OptimizerType>::LearnDistance(arma::mat& outputMatrix,
       (outputMatrix.n_rows > dataset.n_rows) ||
       !(arma::is_finite(outputMatrix)))
   {
-    Log::Info << "Initial learning point have invalid dimensionality.  "
-        "Identity matrix will be used as initial learning point for "
-         "optimization." << std::endl;
     outputMatrix.eye(dataset.n_rows, dataset.n_rows);
   }
 
-  optimizer.Optimize(objFunction, outputMatrix, callbacks...);
+  opt.Optimize(objFunction, outputMatrix, callbacks...);
 }
 
+// Serialize the LMNN object.
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename Archive>
+void LMNN<DistanceType, DeprecatedOptimizerType>::serialize(
+    Archive& ar, const unsigned int /* version */)
+{
+  ar(CEREAL_NVP(k));
+  ar(CEREAL_NVP(regularization));
+  ar(CEREAL_NVP(updateInterval));
+  ar(CEREAL_NVP(distance));
+}
 
 } // namespace mlpack
 
diff --git a/src/mlpack/methods/lmnn/lmnn_main.cpp b/src/mlpack/methods/lmnn/lmnn_main.cpp
index 821c4bd62d..0c270624e6 100644
--- a/src/mlpack/methods/lmnn/lmnn_main.cpp
+++ b/src/mlpack/methods/lmnn/lmnn_main.cpp
@@ -57,7 +57,7 @@ BINDING_LONG_DESC(
     PRINT_PARAM_STRING("regularization") + "), In addition, this "
     "implementation of LMNN includes a parameter to decide the interval "
     "after which impostors must be re-calculated (specified with " +
-    PRINT_PARAM_STRING("range") + ")."
+    PRINT_PARAM_STRING("update_interval") + ")."
     "\n\n"
     "Output can either be the learned distance matrix (specified with " +
     PRINT_PARAM_STRING("output") +"), or the transformed dataset "
@@ -124,11 +124,11 @@ BINDING_EXAMPLE(
     PRINT_CALL("lmnn", "input", "iris", "labels", "iris_labels", "k", 3,
     "optimizer", "bbsgd", "output", "output") +
     "\n\n"
-    "An another program call making use of range & regularization parameter "
-    "with dataset having labels as last column can be made as: "
+    "Another program call making use of update interval & regularization "
+    "parameter with dataset having labels as last column can be made as: "
     "\n\n" +
     PRINT_CALL("lmnn", "input", "letter_recognition", "k", 5,
-    "range", 10, "regularization", 0.4, "output", "output"));
+    "update_interval", 10, "regularization", 0.4, "output", "output"));
 
 // See also...
 BINDING_SEE_ALSO("@nca", "#nca");
@@ -174,8 +174,8 @@ PARAM_DOUBLE_IN("step_size", "Step size for AMSGrad, BB_SGD and SGD (alpha).",
 PARAM_FLAG("linear_scan", "Don't shuffle the order in which data points are "
     "visited for SGD or mini-batch SGD.", "L");
 PARAM_INT_IN("batch_size", "Batch size for mini-batch SGD.", "b", 50);
-PARAM_INT_IN("range", "Number of iterations after which impostors needs to be "
-    "recalculated", "R", 1);
+PARAM_INT_IN("update_interval", "Number of iterations after which impostors "
+    "need to be recalculated.", "R", 1);
 PARAM_INT_IN("seed", "Random seed.  If 0, 'std::time(NULL)' is used.", "s", 0);
 
 using namespace mlpack;
@@ -264,8 +264,8 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers)
 
   RequireParamValue<int>(params, "k", [](int x) { return x > 0; }, true,
       "number of targets must be positive");
-  RequireParamValue<int>(params, "range", [](int x) { return x > 0; }, true,
-      "range must be positive");
+  RequireParamValue<int>(params, "update_interval", [](int x) { return x > 0; },
+      true, "update interval must be positive");
   RequireParamValue<int>(params, "batch_size", [](int x) { return x > 0; }, true,
       "batch size must be positive");
   RequireParamValue<double>(params, "regularization", [](double x)
@@ -294,7 +294,7 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers)
   const bool printAccuracy = params.Has("print_accuracy");
   const bool shuffle = !params.Has("linear_scan");
   const size_t batchSize = (size_t) params.Get<int>("batch_size");
-  const size_t range = (size_t) params.Get<int>("range");
+  const size_t updateInterval = (size_t) params.Get<int>("update_interval");
   const size_t rank = (size_t) params.Get<int>("rank");
 
   // Load data.
@@ -359,56 +359,49 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers)
 
   // Now create the LMNN object and run the optimization.
   timers.Start("lmnn_optimization");
+  LMNN lmnn(k, regularization, updateInterval);
   if (optimizerType == "amsgrad")
   {
-    LMNN<LMetric<2>> lmnn(data, labels, k);
-    lmnn.Regularization() = regularization;
-    lmnn.Range() = range;
-    lmnn.Optimizer().StepSize() = stepSize;
-    lmnn.Optimizer().MaxIterations() = passes * data.n_cols;
-    lmnn.Optimizer().Tolerance() = tolerance;
-    lmnn.Optimizer().Shuffle() = shuffle;
-    lmnn.Optimizer().BatchSize() = batchSize;
+    ens::AMSGrad opt;
+    opt.StepSize() = stepSize;
+    opt.MaxIterations() = passes * data.n_cols;
+    opt.Tolerance() = tolerance;
+    opt.Shuffle() = shuffle;
+    opt.BatchSize() = batchSize;
 
-    lmnn.LearnDistance(distance);
+    lmnn.LearnDistance(data, labels, distance, opt);
   }
   else if (optimizerType == "bbsgd")
   {
-    LMNN<LMetric<2>, ens::BBS_BB> lmnn(data, labels, k);
-    lmnn.Regularization() = regularization;
-    lmnn.Range() = range;
-    lmnn.Optimizer().StepSize() = stepSize;
-    lmnn.Optimizer().MaxIterations() = passes * data.n_cols;
-    lmnn.Optimizer().Tolerance() = tolerance;
-    lmnn.Optimizer().Shuffle() = shuffle;
-    lmnn.Optimizer().BatchSize() = batchSize;
+    ens::BBS_BB opt;
+    opt.StepSize() = stepSize;
+    opt.MaxIterations() = passes * data.n_cols;
+    opt.Tolerance() = tolerance;
+    opt.Shuffle() = shuffle;
+    opt.BatchSize() = batchSize;
 
-    lmnn.LearnDistance(distance);
+    lmnn.LearnDistance(data, labels, distance, opt);
   }
   else if (optimizerType == "sgd")
   {
     // Using SGD is not recommended as the learning matrix can
     // diverge to inf causing serious memory problems.
-    LMNN<LMetric<2>, ens::StandardSGD> lmnn(data, labels, k);
-    lmnn.Regularization() = regularization;
-    lmnn.Range() = range;
-    lmnn.Optimizer().StepSize() = stepSize;
-    lmnn.Optimizer().MaxIterations() = passes * data.n_cols;
-    lmnn.Optimizer().Tolerance() = tolerance;
-    lmnn.Optimizer().Shuffle() = shuffle;
-    lmnn.Optimizer().BatchSize() = batchSize;
+    ens::StandardSGD opt;
+    opt.StepSize() = stepSize;
+    opt.MaxIterations() = passes * data.n_cols;
+    opt.Tolerance() = tolerance;
+    opt.Shuffle() = shuffle;
+    opt.BatchSize() = batchSize;
 
-    lmnn.LearnDistance(distance);
+    lmnn.LearnDistance(data, labels, distance, opt);
   }
   else if (optimizerType == "lbfgs")
   {
-    LMNN<LMetric<2>, ens::L_BFGS> lmnn(data, labels, k);
-    lmnn.Regularization() = regularization;
-    lmnn.Range() = range;
-    lmnn.Optimizer().MaxIterations() = maxIterations;
-    lmnn.Optimizer().MinGradientNorm() = tolerance;
+    ens::L_BFGS opt;
+    opt.MaxIterations() = maxIterations;
+    opt.MinGradientNorm() = tolerance;
 
-    lmnn.LearnDistance(distance);
+    lmnn.LearnDistance(data, labels, distance, opt);
   }
   timers.Stop("lmnn_optimization");
 
diff --git a/src/mlpack/methods/nca/nca.hpp b/src/mlpack/methods/nca/nca.hpp
index 4bf817124c..fd66154b28 100644
--- a/src/mlpack/methods/nca/nca.hpp
+++ b/src/mlpack/methods/nca/nca.hpp
@@ -42,7 +42,7 @@ namespace mlpack {
  * @endcode
  */
 template<typename DistanceType = SquaredEuclideanDistance,
-         typename OptimizerType = ens::StandardSGD>
+         typename DeprecatedOptimizerType = ens::StandardSGD>
 class NCA
 {
  public:
@@ -55,10 +55,18 @@ class NCA
    * @param labels Input dataset labels.
    * @param distance Instantiated distance metric to use.
    */
+  [[deprecated("Will be removed in mlpack 5.0.0.  Pass the dataset directly to "
+               "LearnDistance() instead.")]]
   NCA(const arma::mat& dataset,
       const arma::Row<size_t>& labels,
       DistanceType distance = DistanceType());
 
+  /**
+   * Construct the Neighborhood Components Analysis object, optionally with an
+   * instantiated distance metric.
+   */
+  NCA(DistanceType distance = DistanceType());
+
   /**
    * Perform Neighborhood Components Analysis.  The output distance learning
    * matrix is written into the passed reference.  If LearnDistance() is called
@@ -71,32 +79,109 @@ class NCA
    * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
    *      See https://www.ensmallen.org/docs.html#callback-documentation.
    */
-  template<typename... CallbackTypes>
+  template<typename... CallbackTypes,
+           typename = typename std::enable_if<IsEnsCallbackTypes<
+               CallbackTypes...
+           >::value>::type,
+           typename = typename std::enable_if<
+               !FirstElementIsArma<CallbackTypes...>::value
+           >::type>
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the version that takes a "
+               "dataset as a parameter.")]]
   void LearnDistance(arma::mat& outputMatrix, CallbackTypes&&... callbacks);
 
+  /**
+   * Perform Neighborhood Components Analysis.  The output distance learning
+   * matrix is written into the passed reference.  If LearnDistance() is called
+   * with an outputMatrix which has the correct size (dataset.n_rows x
+   * dataset.n_rows), that matrix will be used as the starting point for
+   * optimization.
+   *
+   * @param dataset Dataset to learn distance metric on.
+   * @param labels Labels for dataset.
+   * @param outputMatrix Covariance matrix of Mahalanobis distance.
+   * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
+   *      See https://www.ensmallen.org/docs.html#callback-documentation.
+   */
+  template<typename MatType,
+           typename LabelsType,
+           typename... CallbackTypes,
+           typename = typename std::enable_if<!IsEnsOptimizer<
+               typename First<CallbackTypes...>::type,
+               SoftmaxErrorFunction<MatType, LabelsType, DistanceType>,
+               MatType
+           >::value>::type,
+           typename = typename std::enable_if<IsEnsCallbackTypes<
+               CallbackTypes...
+           >::value>::type>
+  void LearnDistance(const MatType& dataset,
+                     const LabelsType& labels,
+                     MatType& outputMatrix,
+                     CallbackTypes&&... callbacks) const;
+
+  /**
+   * Perform Neighborhood Components Analysis.  The output distance learning
+   * matrix is written into the passed reference.  If LearnDistance() is called
+   * with an outputMatrix which has the correct size (dataset.n_rows x
+   * dataset.n_rows), that matrix will be used as the starting point for
+   * optimization.
+   *
+   * @param dataset Dataset to learn distance metric on.
+   * @param labels Labels for dataset.
+   * @param optimizer Instantiated ensmallen optimizer to use for NCA.
+   * @param outputMatrix Covariance matrix of Mahalanobis distance.
+   * @param callbacks Callback function for ensmallen optimizer `OptimizerType`.
+   *      See https://www.ensmallen.org/docs.html#callback-documentation.
+   */
+  template<typename MatType,
+           typename LabelsType,
+           typename OptimizerType,
+           typename... CallbackTypes,
+           typename = typename std::enable_if<IsEnsOptimizer<
+               OptimizerType,
+               SoftmaxErrorFunction<MatType, LabelsType, DistanceType>,
+               MatType
+           >::value>::type>
+  void LearnDistance(const MatType& dataset,
+                     const LabelsType& labels,
+                     MatType& outputMatrix,
+                     OptimizerType& optimizer,
+                     CallbackTypes&&... callbacks) const;
+
   //! Get the dataset reference.
-  const arma::mat& Dataset() const { return dataset; }
+  [[deprecated("Will be removed in mlpack 5.0.0.")]]
+  const arma::mat& Dataset() const { return *dataset; }
   //! Get the labels reference.
-  const arma::Row<size_t>& Labels() const { return labels; }
+  [[deprecated("Will be removed in mlpack 5.0.0.")]]
+  const arma::Row<size_t>& Labels() const { return *labels; }
 
   //! Get the optimizer.
-  const OptimizerType& Optimizer() const { return optimizer; }
-  OptimizerType& Optimizer() { return optimizer; }
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  const DeprecatedOptimizerType& Optimizer() const { return optimizer; }
+  //! Modify the optimizer.
+  [[deprecated("Will be removed in mlpack 5.0.0.  Use the LearnDistance() "
+               "version that takes the optimizer as a parameter instead.")]]
+  DeprecatedOptimizerType& Optimizer() { return optimizer; }
+
+  //! Get the distance.
+  const DistanceType Distance() const { return distance; }
+  //! Modify the distance.
+  DistanceType& Distance() { return distance; }
+
+  template<typename Archive>
+  void serialize(Archive& ar, const unsigned int /* version */);
 
  private:
-  //! Dataset reference.
-  const arma::mat& dataset;
-  //! Labels reference.
-  const arma::Row<size_t>& labels;
+  //! Dataset pointer (will be removed in mlpack 5.0.0).
+  const arma::mat* dataset;
+  //! Labels reference (will be removed in mlpack 5.0.0).
+  const arma::Row<size_t>* labels;
+  //! The optimizer to use (will be removed in mlpack 5.0.0).
+  DeprecatedOptimizerType optimizer;
 
   //! Distance to be used.
   DistanceType distance;
-
-  //! The function to optimize.
-  SoftmaxErrorFunction<DistanceType> errorFunction;
-
-  //! The optimizer to use.
-  OptimizerType optimizer;
 };
 
 } // namespace mlpack
diff --git a/src/mlpack/methods/nca/nca_impl.hpp b/src/mlpack/methods/nca/nca_impl.hpp
index c0b9dd7098..6b3f9f9103 100644
--- a/src/mlpack/methods/nca/nca_impl.hpp
+++ b/src/mlpack/methods/nca/nca_impl.hpp
@@ -18,27 +18,89 @@
 namespace mlpack {
 
 // Just set the internal matrix reference.
-template<typename DistanceType, typename OptimizerType>
-NCA<DistanceType, OptimizerType>::NCA(const arma::mat& dataset,
-                                    const arma::Row<size_t>& labels,
-                                    DistanceType distance) :
-    dataset(dataset),
-    labels(labels),
-    distance(distance),
-    errorFunction(dataset, labels, distance)
+template<typename DistanceType, typename DeprecatedOptimizerType>
+NCA<DistanceType, DeprecatedOptimizerType>::NCA(
+    const arma::mat& dataset,
+    const arma::Row<size_t>& labels,
+    DistanceType distance) :
+    dataset(&dataset),
+    labels(&labels),
+    distance(std::move(distance))
 { /* Nothing to do. */ }
 
-template<typename DistanceType, typename OptimizerType>
-template<typename... CallbackTypes>
-void NCA<DistanceType, OptimizerType>::LearnDistance(arma::mat& outputMatrix,
+template<typename DistanceType, typename DeprecatedOptimizerType>
+NCA<DistanceType, DeprecatedOptimizerType>::NCA(DistanceType distance) :
+    distance(std::move(distance))
+{ /* Nothing to do. */ }
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename... CallbackTypes,
+         typename /* callback SFINAE check */,
+         typename /* SFINAE check to disambiguate overloads */>
+void NCA<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    arma::mat& outputMatrix,
     CallbackTypes&&... callbacks)
 {
+  if (!dataset || !labels)
+  {
+    throw std::runtime_error("NCA::LearnDistance(): cannot call without a "
+        "dataset!");
+  }
+
+  LearnDistance(*dataset, *labels, outputMatrix, optimizer,
+      std::forward<CallbackTypes>(callbacks)...);
+}
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename MatType,
+         typename LabelsType,
+         typename... CallbackTypes,
+         typename /* SFINAE check that first callback is not an optimizer */,
+         typename /* callback SFINAE check */>
+void NCA<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    const MatType& dataset,
+    const LabelsType& labels,
+    MatType& outputMatrix,
+    CallbackTypes&&... callbacks) const
+{
+  // This should be replaced with ens::StandardSGD when the deprecated members
+  // are removed for mlpack 5.0.0.
+  DeprecatedOptimizerType opt;
+  LearnDistance(dataset, labels, outputMatrix, opt,
+      std::forward<CallbackTypes>(callbacks)...);
+}
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename MatType,
+         typename LabelsType,
+         typename OptimizerType,
+         typename... CallbackTypes,
+         typename /* SFINAE check that opt is an ensmallen optimizer */>
+void NCA<DistanceType, DeprecatedOptimizerType>::LearnDistance(
+    const MatType& dataset,
+    const LabelsType& labels,
+    MatType& outputMatrix,
+    OptimizerType& opt,
+    CallbackTypes&&... callbacks) const
+{
+  SoftmaxErrorFunction<MatType, LabelsType, DistanceType> errorFunction(
+      dataset, labels, distance);
+
   // See if we were passed an initialized matrix.
   if ((outputMatrix.n_rows != dataset.n_rows) ||
       (outputMatrix.n_cols != dataset.n_rows))
     outputMatrix.eye(dataset.n_rows, dataset.n_rows);
 
-  optimizer.Optimize(errorFunction, outputMatrix, callbacks...);
+  opt.Optimize(errorFunction, outputMatrix,
+      std::forward<CallbackTypes>(callbacks)...);
+}
+
+template<typename DistanceType, typename DeprecatedOptimizerType>
+template<typename Archive>
+void NCA<DistanceType, DeprecatedOptimizerType>::serialize(
+    Archive& ar, const unsigned int /* version */)
+{
+  ar(CEREAL_NVP(distance));
 }
 
 } // namespace mlpack
diff --git a/src/mlpack/methods/nca/nca_main.cpp b/src/mlpack/methods/nca/nca_main.cpp
index cae98ce87c..598efe0e12 100644
--- a/src/mlpack/methods/nca/nca_main.cpp
+++ b/src/mlpack/methods/nca/nca_main.cpp
@@ -240,30 +240,31 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers)
 
   // Now create the NCA object and run the optimization.
   timers.Start("nca_optimization");
+  NCA nca;
   if (optimizerType == "sgd")
   {
-    NCA<LMetric<2> > nca(data, labels);
-    nca.Optimizer().StepSize() = stepSize;
-    nca.Optimizer().MaxIterations() = maxIterations;
-    nca.Optimizer().Tolerance() = tolerance;
-    nca.Optimizer().Shuffle() = shuffle;
-    nca.Optimizer().BatchSize() = batchSize;
+    ens::StandardSGD opt;
+    opt.StepSize() = stepSize;
+    opt.MaxIterations() = maxIterations;
+    opt.Tolerance() = tolerance;
+    opt.Shuffle() = shuffle;
+    opt.BatchSize() = batchSize;
 
-    nca.LearnDistance(distance);
+    nca.LearnDistance(data, labels, distance, opt);
   }
   else if (optimizerType == "lbfgs")
   {
-    NCA<LMetric<2>, ens::L_BFGS> nca(data, labels);
-    nca.Optimizer().NumBasis() = numBasis;
-    nca.Optimizer().MaxIterations() = maxIterations;
-    nca.Optimizer().ArmijoConstant() = armijoConstant;
-    nca.Optimizer().Wolfe() = wolfe;
-    nca.Optimizer().MinGradientNorm() = tolerance;
-    nca.Optimizer().MaxLineSearchTrials() = maxLineSearchTrials;
-    nca.Optimizer().MinStep() = minStep;
-    nca.Optimizer().MaxStep() = maxStep;
+    ens::L_BFGS opt;
+    opt.NumBasis() = numBasis;
+    opt.MaxIterations() = maxIterations;
+    opt.ArmijoConstant() = armijoConstant;
+    opt.Wolfe() = wolfe;
+    opt.MinGradientNorm() = tolerance;
+    opt.MaxLineSearchTrials() = maxLineSearchTrials;
+    opt.MinStep() = minStep;
+    opt.MaxStep() = maxStep;
 
-    nca.LearnDistance(distance);
+    nca.LearnDistance(data, labels, distance, opt);
   }
   timers.Stop("nca_optimization");
 
diff --git a/src/mlpack/methods/nca/nca_softmax_error_function.hpp b/src/mlpack/methods/nca/nca_softmax_error_function.hpp
index e4165c9c12..88237d88b3 100644
--- a/src/mlpack/methods/nca/nca_softmax_error_function.hpp
+++ b/src/mlpack/methods/nca/nca_softmax_error_function.hpp
@@ -40,10 +40,17 @@ namespace mlpack {
  * operate on one point in the dataset.  This is useful for optimizers like
  * stochastic gradient descent (see mlpack::optimization::SGD).
  */
-template<typename DistanceType = SquaredEuclideanDistance>
+template<typename MatType = arma::mat,
+         typename LabelsType = arma::Row<size_t>,
+         typename DistanceType = SquaredEuclideanDistance>
 class SoftmaxErrorFunction
 {
  public:
+  // Convenience typedef for element type of data.
+  typedef typename MatType::elem_type ElemType;
+  // Convenience typedef for column vector of data.
+  typedef typename GetColType<MatType>::type VecType;
+
   /**
    * Initialize with the given kernel; useful when the kernel has some state to
    * store, which is set elsewhere.  If no kernel is given, an empty kernel is
@@ -54,8 +61,8 @@ class SoftmaxErrorFunction
    * @param labels Vector of class labels for each point in the dataset.
    * @param metric Instantiated metric (optional).
    */
-  SoftmaxErrorFunction(const arma::mat& dataset,
-                       const arma::Row<size_t>& labels,
+  SoftmaxErrorFunction(const MatType& dataset,
+                       const LabelsType& labels,
                        DistanceType metric = DistanceType());
 
   /**
@@ -70,7 +77,7 @@ class SoftmaxErrorFunction
    *
    * @param covariance Covariance matrix of Mahalanobis distance.
    */
-  double Evaluate(const arma::mat& covariance);
+  ElemType Evaluate(const MatType& covariance);
 
   /**
    * Evaluate the softmax objective function for the given covariance matrix on
@@ -84,9 +91,9 @@ class SoftmaxErrorFunction
    * @param begin Index of the initial point to use for objective function.
    * @param batchSize Number of points to use for objective function.
    */
-  double Evaluate(const arma::mat& covariance,
-                  const size_t begin,
-                  const size_t batchSize = 1);
+  ElemType Evaluate(const MatType& covariance,
+                    const size_t begin,
+                    const size_t batchSize = 1);
 
   /**
    * Evaluate the gradient of the softmax function for the given covariance
@@ -96,7 +103,7 @@ class SoftmaxErrorFunction
    * @param covariance Covariance matrix of Mahalanobis distance.
    * @param gradient Matrix to store the calculated gradient in.
    */
-  void Gradient(const arma::mat& covariance, arma::mat& gradient);
+  void Gradient(const MatType& covariance, MatType& gradient);
 
   /**
    * Evaluate the gradient of the softmax function for the given covariance
@@ -114,7 +121,7 @@ class SoftmaxErrorFunction
    * @param gradient Matrix to store the calculated gradient in.
    */
   template <typename GradType>
-  void Gradient(const arma::mat& covariance,
+  void Gradient(const MatType& covariance,
                 const size_t begin,
                 GradType& gradient,
                 const size_t batchSize = 1);
@@ -122,7 +129,7 @@ class SoftmaxErrorFunction
   /**
    * Get the initial point.
    */
-  const arma::mat GetInitialPoint() const;
+  const MatType GetInitialPoint() const;
 
   /**
    * Get the number of functions the objective function can be decomposed into.
@@ -132,23 +139,23 @@ class SoftmaxErrorFunction
 
  private:
   //! The dataset.  This is an alias until Shuffle() is called.
-  arma::mat dataset;
+  MatType dataset;
   //! Labels for each point in the dataset.  This is an alias until Shuffle() is
   //! called.
-  arma::Row<size_t> labels;
+  LabelsType labels;
 
   //! The instantiated metric.
   DistanceType distance;
 
   //! Last coordinates.  Used for the non-separable Evaluate() and Gradient().
-  arma::mat lastCoordinates;
+  MatType lastCoordinates;
   //! Stretched dataset.  Kept internal to avoid memory reallocations.
-  arma::mat stretchedDataset;
+  MatType stretchedDataset;
   //! Holds calculated p_i, for the non-separable Evaluate() and Gradient().
-  arma::vec p;
+  VecType p;
   //! Holds denominators for calculation of p_ij, for the non-separable
   //! Evaluate() and Gradient().
-  arma::vec denominators;
+  VecType denominators;
 
   //! False if nothing has ever been precalculated (only at construction time).
   bool precalculated;
@@ -166,7 +173,7 @@ class SoftmaxErrorFunction
    *
    * @param coordinates Coordinates matrix to use for precalculation.
    */
-  void Precalculate(const arma::mat& coordinates);
+  void Precalculate(const MatType& coordinates);
 };
 
 } // namespace mlpack
diff --git a/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp b/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
index 2708e6e247..01e133f338 100644
--- a/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
+++ b/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
@@ -20,10 +20,10 @@
 namespace mlpack {
 
 // Initialize with the given kernel.
-template<typename DistanceType>
-SoftmaxErrorFunction<DistanceType>::SoftmaxErrorFunction(
-    const arma::mat& datasetIn,
-    const arma::Row<size_t>& labelsIn,
+template<typename MatType, typename LabelsType, typename DistanceType>
+SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::SoftmaxErrorFunction(
+    const MatType& datasetIn,
+    const LabelsType& labelsIn,
     DistanceType distance) :
     distance(distance),
     precalculated(false)
@@ -33,11 +33,11 @@ SoftmaxErrorFunction<DistanceType>::SoftmaxErrorFunction(
 }
 
 //! Shuffle the dataset.
-template<typename DistanceType>
-void SoftmaxErrorFunction<DistanceType>::Shuffle()
+template<typename MatType, typename LabelsType, typename DistanceType>
+void SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Shuffle()
 {
-  arma::mat newDataset;
-  arma::Row<size_t> newLabels;
+  MatType newDataset;
+  LabelsType newLabels;
 
   ShuffleData(dataset, labels, newDataset, newLabels);
 
@@ -49,8 +49,10 @@ void SoftmaxErrorFunction<DistanceType>::Shuffle()
 }
 
 //! The non-separable implementation, which uses Precalculate() to save time.
-template<typename DistanceType>
-double SoftmaxErrorFunction<DistanceType>::Evaluate(const arma::mat& coordinates)
+template<typename MatType, typename LabelsType, typename DistanceType>
+typename MatType::elem_type
+SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Evaluate(
+    const MatType& coordinates)
 {
   // Calculate the denominators and numerators, if necessary.
   Precalculate(coordinates);
@@ -61,19 +63,23 @@ double SoftmaxErrorFunction<DistanceType>::Evaluate(const arma::mat& coordinates
 
 //! The separated objective function, which does not use Precalculate(),
 //! for a given batch size and from an initial index.
-template<typename DistanceType>
-double SoftmaxErrorFunction<DistanceType>::Evaluate(const arma::mat& coordinates,
-                                                    const size_t begin,
-                                                    const size_t batchSize)
+template<typename MatType, typename LabelsType, typename DistanceType>
+typename MatType::elem_type
+SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Evaluate(
+    const MatType& coordinates,
+    const size_t begin,
+    const size_t batchSize)
 {
   // Unfortunately each evaluation will take O(N) time because it requires a
   // scan over all points in the dataset.  Our objective is to compute p_i.
-  double denominator = 0;
-  double numerator = 0;
-  double result = 0;
+  ElemType denominator = 0;
+  ElemType numerator = 0;
+  ElemType result = 0;
 
   // It's quicker to do this now than one point at a time later.
   stretchedDataset = coordinates * dataset;
+
+  #pragma omp parallel for reduction(+:result)
   for (size_t i = begin; i < begin + batchSize; ++i)
   {
     for (size_t k = 0; k < dataset.n_cols; ++k)
@@ -83,7 +89,7 @@ double SoftmaxErrorFunction<DistanceType>::Evaluate(const arma::mat& coordinates
         continue;
 
       // We want to evaluate exp(-D(A x_i, A x_k)).
-      double eval = std::exp(-distance.Evaluate(
+      ElemType eval = std::exp(-distance.Evaluate(
           stretchedDataset.unsafe_col(i), stretchedDataset.unsafe_col(k)));
 
       // If they are in the same class, update the numerator.
@@ -108,9 +114,9 @@ double SoftmaxErrorFunction<DistanceType>::Evaluate(const arma::mat& coordinates
 }
 
 //! The non-separable implementation, where Precalculate() is used.
-template<typename DistanceType>
-void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
-                                                  arma::mat& gradient)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Gradient(
+    const MatType& coordinates, MatType& gradient)
 {
   // Calculate the denominators and numerators, if necessary.
   Precalculate(coordinates);
@@ -127,22 +133,22 @@ void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
   //     (((p_i - (1 / p_i)) p_ik) + ((p_k - (1 / p_k)) p_ki)) x_ik x_ik^T
   //   otherwise, add
   //     (p_i p_ik + p_k p_ki) x_ik x_ik^T
-  arma::mat sum;
+  MatType sum;
   sum.zeros(stretchedDataset.n_rows, stretchedDataset.n_rows);
   for (size_t i = 0; i < stretchedDataset.n_cols; ++i)
   {
     for (size_t k = (i + 1); k < stretchedDataset.n_cols; ++k)
     {
       // Calculate p_ik and p_ki first.
-      double eval = std::exp(-distance.Evaluate(
+      ElemType eval = std::exp(-distance.Evaluate(
           stretchedDataset.unsafe_col(i), stretchedDataset.unsafe_col(k)));
-      double p_ik = 0, p_ki = 0;
+      ElemType p_ik = 0, p_ki = 0;
       p_ik = eval / denominators(i);
       p_ki = eval / denominators(k);
 
       // Subtract x_i from x_k.  We are not using stretched points here.
-      arma::vec x_ik = dataset.col(i) - dataset.col(k);
-      arma::mat secondTerm = (x_ik * trans(x_ik));
+      VecType x_ik = dataset.col(i) - dataset.col(k);
+      MatType secondTerm = (x_ik * trans(x_ik));
 
       if (labels[i] == labels[k])
         sum += ((p[i] - 1) * p_ik + (p[k] - 1) * p_ki) * secondTerm;
@@ -156,19 +162,20 @@ void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
 }
 
 //! The separable implementation for a given batch size and an initial index.
-template <typename DistanceType>
+template <typename MatType, typename LabelsType, typename DistanceType>
 template <typename GradType>
-void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
-                                                  const size_t begin,
-                                                  GradType& gradient,
-                                                  const size_t batchSize)
+void SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Gradient(
+    const MatType& coordinates,
+    const size_t begin,
+    GradType& gradient,
+    const size_t batchSize)
 {
   // The gradient involves two matrix terms which are eventually combined into
   // one.
   GradType firstTerm, secondTerm;
   // We will need to calculate p_i before this evaluation is done, so
   // these two variables will hold the information necessary for that.
-  double numerator, denominator;
+  ElemType numerator, denominator;
 
   gradient.zeros(coordinates.n_rows, coordinates.n_rows);
 
@@ -189,7 +196,7 @@ void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
         continue;
 
       // Calculate the numerator of p_ik.
-      double eval = std::exp(-distance.Evaluate(
+      ElemType eval = std::exp(-distance.Evaluate(
           stretchedDataset.unsafe_col(i), stretchedDataset.unsafe_col(k)));
 
       // If the points are in the same class, we must add to the second term of
@@ -210,7 +217,7 @@ void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
     }
 
     // Calculate p_i.
-    double p = 0;
+    ElemType p = 0;
     if (denominator == 0)
     {
       Log::Warn << "Denominator of p_" << i << " is 0!" << std::endl;
@@ -231,15 +238,16 @@ void SoftmaxErrorFunction<DistanceType>::Gradient(const arma::mat& coordinates,
   }
 }
 
-template<typename DistanceType>
-const arma::mat SoftmaxErrorFunction<DistanceType>::GetInitialPoint() const
+template<typename MatType, typename LabelsType, typename DistanceType>
+const MatType
+SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::GetInitialPoint() const
 {
-  return arma::eye<arma::mat>(dataset.n_rows, dataset.n_rows);
+  return arma::eye<MatType>(dataset.n_rows, dataset.n_rows);
 }
 
-template<typename DistanceType>
-void SoftmaxErrorFunction<DistanceType>::Precalculate(
-    const arma::mat& coordinates)
+template<typename MatType, typename LabelsType, typename DistanceType>
+void SoftmaxErrorFunction<MatType, LabelsType, DistanceType>::Precalculate(
+    const MatType& coordinates)
 {
   // Ensure it is the right size.
   if (lastCoordinates.n_rows != coordinates.n_rows ||
@@ -265,22 +273,30 @@ void SoftmaxErrorFunction<DistanceType>::Precalculate(
   // order of O((n * (n + 1)) / 2), which really isn't all that great.
   p.zeros(stretchedDataset.n_cols);
   denominators.zeros(stretchedDataset.n_cols);
+
+  // A collapse(2) would be helpful here, but appears to not be supported fully
+  // until OpenMP 5.0.
+  #pragma omp parallel for
   for (size_t i = 0; i < stretchedDataset.n_cols; ++i)
   {
     for (size_t j = (i + 1); j < stretchedDataset.n_cols; ++j)
     {
       // Evaluate exp(-d(x_i, x_j)).
-      double eval = std::exp(-distance.Evaluate(
+      ElemType eval = std::exp(-distance.Evaluate(
           stretchedDataset.unsafe_col(i), stretchedDataset.unsafe_col(j)));
 
       // Add this to the denominators of both p_i and p_j: K(i, j) = K(j, i).
+      #pragma omp atomic
       denominators[i] += eval;
+      #pragma omp atomic
       denominators[j] += eval;
 
       // If i and j are the same class, add to numerator of both.
       if (labels[i] == labels[j])
       {
+        #pragma omp atomic
         p[i] += eval;
+        #pragma omp atomic
         p[j] += eval;
       }
     }
@@ -290,6 +306,7 @@ void SoftmaxErrorFunction<DistanceType>::Precalculate(
   p /= denominators;
 
   // Clean up any bad values.
+  #pragma omp parallel for
   for (size_t i = 0; i < stretchedDataset.n_cols; ++i)
   {
     if (denominators[i] == 0.0)
@@ -297,7 +314,7 @@ void SoftmaxErrorFunction<DistanceType>::Precalculate(
       Log::Debug << "Denominator of p_{" << i << ", j} is 0." << std::endl;
 
       // Set to usable values.
-      denominators[i] = std::numeric_limits<double>::infinity();
+      denominators[i] = std::numeric_limits<ElemType>::infinity();
       p[i] = 0;
     }
   }
diff --git a/src/mlpack/methods/neighbor_search/neighbor_search.hpp b/src/mlpack/methods/neighbor_search/neighbor_search.hpp
index b8f61c0155..d2470b1b19 100644
--- a/src/mlpack/methods/neighbor_search/neighbor_search.hpp
+++ b/src/mlpack/methods/neighbor_search/neighbor_search.hpp
@@ -222,9 +222,11 @@ class NeighborSearch
    * @param distances Matrix storing distances of neighbors for each query
    *     point.
    */
+  // TODO: templatize further to remove Armadillo type requirement
+  template<typename IndexType = size_t>
   void Search(const MatType& querySet,
               const size_t k,
-              arma::Mat<size_t>& neighbors,
+              arma::Mat<IndexType>& neighbors,
               arma::Mat<ElemType>& distances);
 
   /**
@@ -247,9 +249,11 @@ class NeighborSearch
    * @param sameSet Denotes whether or not the reference and query sets are the
    *      same.
    */
+  // TODO: templatize further to remove Armadillo type requirement
+  template<typename IndexType = size_t>
   void Search(Tree& queryTree,
               const size_t k,
-              arma::Mat<size_t>& neighbors,
+              arma::Mat<IndexType>& neighbors,
               arma::Mat<ElemType>& distances,
               bool sameSet = false);
 
@@ -267,8 +271,10 @@ class NeighborSearch
    * @param distances Matrix storing distances of neighbors for each query
    *      point.
    */
+  // TODO: templatize further to remove Armadillo type requirement
+  template<typename IndexType = size_t>
   void Search(const size_t k,
-              arma::Mat<size_t>& neighbors,
+              arma::Mat<IndexType>& neighbors,
               arma::Mat<ElemType>& distances);
 
   /**
@@ -300,8 +306,10 @@ class NeighborSearch
    *     query point.
    * @return Recall.
    */
-  static double Recall(arma::Mat<size_t>& foundNeighbors,
-                       arma::Mat<size_t>& realNeighbors);
+  // TODO: templatize further to remove Armadillo type requirement
+  template<typename IndexType = size_t>
+  static double Recall(arma::Mat<IndexType>& foundNeighbors,
+                       arma::Mat<IndexType>& realNeighbors);
 
   //! Return the total number of base case evaluations performed during the last
   //! search.
diff --git a/src/mlpack/methods/neighbor_search/neighbor_search_impl.hpp b/src/mlpack/methods/neighbor_search/neighbor_search_impl.hpp
index c73b7855f0..0aa0a7de0d 100644
--- a/src/mlpack/methods/neighbor_search/neighbor_search_impl.hpp
+++ b/src/mlpack/methods/neighbor_search/neighbor_search_impl.hpp
@@ -360,11 +360,12 @@ template<typename SortPolicy,
                   typename TreeMatType> class TreeType,
          template<typename> class DualTreeTraversalType,
          template<typename> class SingleTreeTraversalType>
+template<typename IndexType>
 void NeighborSearch<SortPolicy, DistanceType, MatType, TreeType,
 DualTreeTraversalType, SingleTreeTraversalType>::Search(
     const MatType& querySet,
     const size_t k,
-    arma::Mat<size_t>& neighbors,
+    arma::Mat<IndexType>& neighbors,
     arma::Mat<ElemType>& distances)
 {
   if (k > referenceSet->n_cols)
@@ -385,7 +386,7 @@ DualTreeTraversalType, SingleTreeTraversalType>::Search(
   // indices back to their original indices when this computation is finished.
   // To avoid an extra copy, we will store the neighbors and distances in a
   // separate matrix.
-  arma::Mat<size_t>* neighborPtr = &neighbors;
+  arma::Mat<IndexType>* neighborPtr = &neighbors;
   arma::Mat<ElemType>* distancePtr = &distances;
 
   // Mapping is only necessary if the tree rearranges points.
@@ -394,10 +395,10 @@ DualTreeTraversalType, SingleTreeTraversalType>::Search(
     if (searchMode == DUAL_TREE_MODE)
     {
       distancePtr = new arma::Mat<ElemType>; // Query indices need to be mapped.
-      neighborPtr = new arma::Mat<size_t>;
+      neighborPtr = new arma::Mat<IndexType>;
     }
     else if (!oldFromNewReferences.empty())
-      neighborPtr = new arma::Mat<size_t>; // Reference indices need mapping.
+      neighborPtr = new arma::Mat<IndexType>; // Reference indices need mapping.
   }
 
   // Set the size of the neighbor and distance matrices.
@@ -565,11 +566,12 @@ template<typename SortPolicy,
                   typename TreeMatType> class TreeType,
          template<typename> class DualTreeTraversalType,
          template<typename> class SingleTreeTraversalType>
+template<typename IndexType>
 void NeighborSearch<SortPolicy, DistanceType, MatType, TreeType,
 DualTreeTraversalType, SingleTreeTraversalType>::Search(
     Tree& queryTree,
     const size_t k,
-    arma::Mat<size_t>& neighbors,
+    arma::Mat<IndexType>& neighbors,
     arma::Mat<ElemType>& distances,
     bool sameSet)
 {
@@ -593,10 +595,10 @@ DualTreeTraversalType, SingleTreeTraversalType>::Search(
   const MatType& querySet = queryTree.Dataset();
 
   // We won't need to map query indices, but will we need to map distances?
-  arma::Mat<size_t>* neighborPtr = &neighbors;
+  arma::Mat<IndexType>* neighborPtr = &neighbors;
 
   if (!oldFromNewReferences.empty() && TreeTraits<Tree>::RearrangesDataset)
-    neighborPtr = new arma::Mat<size_t>;
+    neighborPtr = new arma::Mat<IndexType>;
 
   neighborPtr->set_size(k, querySet.n_cols);
   distances.set_size(k, querySet.n_cols);
@@ -644,10 +646,11 @@ template<typename SortPolicy,
                   typename TreeMatType> class TreeType,
          template<typename> class DualTreeTraversalType,
          template<typename> class SingleTreeTraversalType>
+template<typename IndexType>
 void NeighborSearch<SortPolicy, DistanceType, MatType, TreeType,
 DualTreeTraversalType, SingleTreeTraversalType>::Search(
     const size_t k,
-    arma::Mat<size_t>& neighbors,
+    arma::Mat<IndexType>& neighbors,
     arma::Mat<ElemType>& distances)
 {
   if (k > referenceSet->n_cols)
@@ -669,14 +672,14 @@ DualTreeTraversalType, SingleTreeTraversalType>::Search(
   baseCases = 0;
   scores = 0;
 
-  arma::Mat<size_t>* neighborPtr = &neighbors;
+  arma::Mat<IndexType>* neighborPtr = &neighbors;
   arma::Mat<ElemType>* distancePtr = &distances;
 
   if (!oldFromNewReferences.empty() && TreeTraits<Tree>::RearrangesDataset)
   {
     // We will always need to rearrange in this case.
     distancePtr = new MatType;
-    neighborPtr = new arma::Mat<size_t>;
+    neighborPtr = new arma::Mat<IndexType>;
   }
 
   // Initialize results.
@@ -861,10 +864,11 @@ template<typename SortPolicy,
                   typename TreeMatType> class TreeType,
          template<typename> class DualTreeTraversalType,
          template<typename> class SingleTreeTraversalType>
+template<typename IndexType>
 double NeighborSearch<SortPolicy, DistanceType, MatType, TreeType,
 DualTreeTraversalType, SingleTreeTraversalType>::Recall(
-    arma::Mat<size_t>& foundNeighbors,
-    arma::Mat<size_t>& realNeighbors)
+    arma::Mat<IndexType>& foundNeighbors,
+    arma::Mat<IndexType>& realNeighbors)
 {
   if (foundNeighbors.n_rows != realNeighbors.n_rows ||
       foundNeighbors.n_cols != realNeighbors.n_cols)
diff --git a/src/mlpack/methods/neighbor_search/neighbor_search_rules.hpp b/src/mlpack/methods/neighbor_search/neighbor_search_rules.hpp
index ef886663a7..057515ca7a 100644
--- a/src/mlpack/methods/neighbor_search/neighbor_search_rules.hpp
+++ b/src/mlpack/methods/neighbor_search/neighbor_search_rules.hpp
@@ -63,7 +63,10 @@ class NeighborSearchRules
    * @param distances Matrix storing distances of neighbors for each query
    *     point.
    */
-  void GetResults(arma::Mat<size_t>& neighbors, arma::Mat<ElemType>& distances);
+  // TODO: templatize fully to remove requirement of Armadillo matrix
+  template<typename IndexType = size_t>
+  void GetResults(arma::Mat<IndexType>& neighbors,
+                  arma::Mat<ElemType>& distances);
 
   /**
    * Get the distance from the query point to the reference point.
diff --git a/src/mlpack/methods/neighbor_search/neighbor_search_rules_impl.hpp b/src/mlpack/methods/neighbor_search/neighbor_search_rules_impl.hpp
index 24fa48c6d5..cfbd350090 100644
--- a/src/mlpack/methods/neighbor_search/neighbor_search_rules_impl.hpp
+++ b/src/mlpack/methods/neighbor_search/neighbor_search_rules_impl.hpp
@@ -59,8 +59,9 @@ NeighborSearchRules<SortPolicy, DistanceType, TreeType>::NeighborSearchRules(
 }
 
 template<typename SortPolicy, typename DistanceType, typename TreeType>
+template<typename IndexType>
 void NeighborSearchRules<SortPolicy, DistanceType, TreeType>::GetResults(
-    arma::Mat<size_t>& neighbors,
+    arma::Mat<IndexType>& neighbors,
     arma::Mat<ElemType>& distances)
 {
   neighbors.set_size(k, querySet.n_cols);
@@ -71,7 +72,7 @@ void NeighborSearchRules<SortPolicy, DistanceType, TreeType>::GetResults(
     CandidateList& pqueue = candidates[i];
     for (size_t j = 1; j <= k; ++j)
     {
-      neighbors(k - j, i) = pqueue.top().second;
+      neighbors(k - j, i) = (IndexType) pqueue.top().second;
       distances(k - j, i) = pqueue.top().first;
       pqueue.pop();
     }
diff --git a/src/mlpack/tests/callback_test.cpp b/src/mlpack/tests/callback_test.cpp
index bd98dcad85..ae9b912c9d 100644
--- a/src/mlpack/tests/callback_test.cpp
+++ b/src/mlpack/tests/callback_test.cpp
@@ -151,12 +151,13 @@ TEST_CASE("LMNNWithOptimizerCallback", "[CallbackTest]")
                       " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  LMNN<> lmnn(dataset, labels, 1);
+  LMNN<> lmnn(1);
 
   arma::mat outputMatrix;
   std::stringstream stream;
 
-  lmnn.LearnDistance(outputMatrix, ens::ProgressBar(70, stream));
+  lmnn.LearnDistance(dataset, labels, outputMatrix,
+      ens::ProgressBar(70, stream));
   REQUIRE(stream.str().length() > 0);
 }
 
@@ -170,12 +171,11 @@ TEST_CASE("NCAWithOptimizerCallback", "[CallbackTest]")
                    " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  NCA<SquaredEuclideanDistance> nca(data, labels);
-
   arma::mat outputMatrix;
   std::stringstream stream;
 
-  nca.LearnDistance(outputMatrix, ens::ProgressBar(70, stream));
+  NCA nca;
+  nca.LearnDistance(data, labels, outputMatrix, ens::ProgressBar(70, stream));
   REQUIRE(stream.str().length() > 0);
 }
 
diff --git a/src/mlpack/tests/lmnn_test.cpp b/src/mlpack/tests/lmnn_test.cpp
index 44eb0b0d66..6b0c0d254a 100644
--- a/src/mlpack/tests/lmnn_test.cpp
+++ b/src/mlpack/tests/lmnn_test.cpp
@@ -30,25 +30,27 @@ using namespace ens;
  * The target neighbors function should be correct.
  * point.
  */
-TEST_CASE("LMNNTargetNeighborsTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNTargetNeighborsTest", "[LMNNTest]", float, double)
 {
-  // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  typedef TestType ElemType;
 
-  Constraints<> constraint(dataset, labels, 1);
+  // Useful but simple dataset with six points and two classes.
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
+
+  Constraints<arma::Mat<ElemType>, arma::Row<size_t>> constraint(dataset,
+      labels, 1);
 
   // Calculate norm of datapoints.
-  arma::vec norm(dataset.n_cols);
+  arma::Col<ElemType> norm(dataset.n_cols);
   for (size_t i = 0; i < dataset.n_cols; ++i)
   {
     norm(i) = arma::norm(dataset.col(i));
   }
 
   //! Store target neighbors of data points.
-  arma::Mat<size_t> targetNeighbors =
-      arma::Mat<size_t>(1, dataset.n_cols, arma::fill::zeros);
+  arma::umat targetNeighbors(1, dataset.n_cols, arma::fill::zeros);
 
   constraint.TargetNeighbors(targetNeighbors, dataset, labels, norm);
 
@@ -63,25 +65,27 @@ TEST_CASE("LMNNTargetNeighborsTest", "[LMNNTest]")
 /**
  * The impostors function should be correct.
  */
-TEST_CASE("LMNNImpostorsTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNImpostorsTest", "[LMNNTest]", float, double)
 {
-  // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  typedef TestType ElemType;
 
-  Constraints<> constraint(dataset, labels, 1);
+  // Useful but simple dataset with six points and two classes.
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
+
+  Constraints<arma::Mat<ElemType>, arma::Row<size_t>> constraint(dataset,
+      labels, 1);
 
   // Calculate norm of datapoints.
-  arma::vec norm(dataset.n_cols);
+  arma::Col<ElemType> norm(dataset.n_cols);
   for (size_t i = 0; i < dataset.n_cols; ++i)
   {
     norm(i) = arma::norm(dataset.col(i));
   }
 
   //! Store impostors of data points.
-  arma::Mat<size_t> impostors =
-      arma::Mat<size_t>(1, dataset.n_cols, arma::fill::zeros);
+  arma::umat impostors(1, dataset.n_cols, arma::fill::zeros);
 
   constraint.Impostors(impostors, dataset, labels, norm);
 
@@ -101,300 +105,339 @@ TEST_CASE("LMNNImpostorsTest", "[LMNNTest]")
  * The LMNN function should return the identity matrix as its initial
  * point.
  */
-TEST_CASE("LMNNInitialPointTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNInitialPointTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Cheap fake dataset.
-  arma::mat dataset = arma::randu(5, 5);
+  arma::Mat<ElemType> dataset = arma::randu<arma::Mat<ElemType>>(5, 5);
   arma::Row<size_t> labels = "0 1 1 0 0";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.5, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.5, 1);
 
   // Verify the initial point is the identity matrix.
-  arma::mat initialPoint = lmnnfn.GetInitialPoint();
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<ElemType, float>::value ? 1e-4 : 1e-5;
+  arma::Mat<ElemType> initialPoint = lmnnfn.GetInitialPoint();
   for (int row = 0; row < 5; row++)
   {
     for (int col = 0; col < 5; col++)
     {
       if (row == col)
-        REQUIRE(initialPoint(row, col) == Approx(1.0).epsilon(1e-7));
+        REQUIRE(initialPoint(row, col) == Approx(1.0).epsilon(eps));
       else
-        REQUIRE(initialPoint(row, col) == Approx(0.0).margin(1e-5));
+        REQUIRE(initialPoint(row, col) == Approx(0.0).margin(margin));
     }
   }
 }
 
 /***
- * Ensure non-seprable objective function is right.
+ * Ensure non-separable objective function is right.
  */
-TEST_CASE("LMNNInitialEvaluationTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNInitialEvaluationTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  double objective = lmnnfn.Evaluate(arma::eye<arma::mat>(2, 2));
+  ElemType objective = lmnnfn.Evaluate(arma::eye<arma::Mat<ElemType>>(2, 2));
 
   // Result calculated by hand.
-  REQUIRE(objective == Approx(9.456).epsilon(1e-7));
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  REQUIRE(objective == Approx(9.456).epsilon(eps));
 }
 
 /**
- * Ensure non-seprable gradient function is right.
+ * Ensure non-separable gradient function is right.
  */
-TEST_CASE("LMNNInitialGradientTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNInitialGradientTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  arma::mat gradient;
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
+  arma::Mat<ElemType> gradient;
+  arma::Mat<ElemType> coordinates = arma::eye<arma::Mat<ElemType>>(2, 2);
   lmnnfn.Gradient(coordinates, gradient);
 
   // Result calculated by hand.
-  REQUIRE(gradient(0, 0) == Approx(-0.288).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).margin(1e-5));
-  REQUIRE(gradient(0, 1) == Approx(0.0).margin(1e-5));
-  REQUIRE(gradient(1, 1) == Approx(12.0).epsilon(1e-7));
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<ElemType, float>::value ? 1e-4 : 1e-5;
+  REQUIRE(gradient(0, 0) == Approx(-0.288).epsilon(eps));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(12.0).epsilon(eps));
 }
 
 /***
- * Ensure non-seprable EvaluateWithGradient function is right.
+ * Ensure non-separable EvaluateWithGradient function is right.
  */
-TEST_CASE("LMNNInitialEvaluateWithGradientTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNInitialEvaluateWithGradientTest", "[LMNNTest]", float,
+    double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  arma::mat gradient;
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
-  double objective = lmnnfn.EvaluateWithGradient(coordinates, gradient);
+  arma::Mat<ElemType> gradient;
+  arma::Mat<ElemType> coordinates = arma::eye<arma::Mat<ElemType>>(2, 2);
+  ElemType objective = lmnnfn.EvaluateWithGradient(coordinates, gradient);
+
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<ElemType, float>::value ? 1e-4 : 1e-5;
 
   // Result calculated by hand.
-  REQUIRE(objective == Approx(9.456).epsilon(1e-7));
+  REQUIRE(objective == Approx(9.456).epsilon(eps));
   // Check Gradient
-  REQUIRE(gradient(0, 0) == Approx(-0.288).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).margin(1e-5));
-  REQUIRE(gradient(0, 1) == Approx(0.0).margin(1e-5));
-  REQUIRE(gradient(1, 1) == Approx(12.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.288).epsilon(eps));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(12.0).epsilon(eps));
 }
 
 /**
  * Ensure the separable objective function is right.
  */
-TEST_CASE("LMNNSeparableObjectiveTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNSeparableObjectiveTest", "[LMNNTest]", float, double)
 {
-  // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  typedef TestType ElemType;
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  // Useful but simple dataset with six points and two classes.
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
+
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
   // Result calculated by hand.
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
-  REQUIRE(lmnnfn.Evaluate(coordinates, 0, 1) == Approx(1.576).epsilon(1e-7));
-  REQUIRE(lmnnfn.Evaluate(coordinates, 1, 1) == Approx(1.576).epsilon(1e-7));
-  REQUIRE(lmnnfn.Evaluate(coordinates, 2, 1) == Approx(1.576).epsilon(1e-7));
-  REQUIRE(lmnnfn.Evaluate(coordinates, 3, 1) == Approx(1.576).epsilon(1e-7));
-  REQUIRE(lmnnfn.Evaluate(coordinates, 4, 1) == Approx(1.576).epsilon(1e-7));
-  REQUIRE(lmnnfn.Evaluate(coordinates, 5, 1) == Approx(1.576).epsilon(1e-7));
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  arma::Mat<ElemType> coordinates = arma::eye<arma::Mat<ElemType>>(2, 2);
+  REQUIRE(lmnnfn.Evaluate(coordinates, 0, 1) == Approx(1.576).epsilon(eps));
+  REQUIRE(lmnnfn.Evaluate(coordinates, 1, 1) == Approx(1.576).epsilon(eps));
+  REQUIRE(lmnnfn.Evaluate(coordinates, 2, 1) == Approx(1.576).epsilon(eps));
+  REQUIRE(lmnnfn.Evaluate(coordinates, 3, 1) == Approx(1.576).epsilon(eps));
+  REQUIRE(lmnnfn.Evaluate(coordinates, 4, 1) == Approx(1.576).epsilon(eps));
+  REQUIRE(lmnnfn.Evaluate(coordinates, 5, 1) == Approx(1.576).epsilon(eps));
 }
 
 /**
  * Ensure the separable gradient is right.
  */
-TEST_CASE("LMNNSeparableGradientTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNSeparableGradientTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
-  arma::mat gradient(2, 2);
+  arma::Mat<ElemType> coordinates = arma::eye<arma::Mat<ElemType>>(2, 2);
+  arma::Mat<ElemType> gradient(2, 2);
 
   lmnnfn.Gradient(coordinates, 0, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<ElemType, float>::value ? 1e-4 : 1e-5;
+
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   lmnnfn.Gradient(coordinates, 1, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   lmnnfn.Gradient(coordinates, 2, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   lmnnfn.Gradient(coordinates, 3, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   lmnnfn.Gradient(coordinates, 4, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   lmnnfn.Gradient(coordinates, 5, gradient, 1);
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 }
 
 /**
  * Ensure the separable EvaluateWithGradient function is right.
  */
-TEST_CASE("LMNNSeparableEvaluateWithGradientTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNSeparableEvaluateWithGradientTest", "[LMNNTest]", float,
+    double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
-  arma::mat gradient(2, 2);
+  arma::Mat<ElemType> coordinates = arma::eye<arma::Mat<ElemType>>(2, 2);
+  arma::Mat<ElemType> gradient(2, 2);
 
-  double objective = lmnnfn.EvaluateWithGradient(coordinates, 0, gradient, 1);
+  ElemType objective = lmnnfn.EvaluateWithGradient(coordinates, 0, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  const double eps = std::is_same<ElemType, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<ElemType, float>::value ? 1e-4 : 1e-5;
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
+
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   objective = lmnnfn.EvaluateWithGradient(coordinates, 1, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   objective = lmnnfn.EvaluateWithGradient(coordinates, 2, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   objective = lmnnfn.EvaluateWithGradient(coordinates, 3, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   objective = lmnnfn.EvaluateWithGradient(coordinates, 4, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 
   objective = lmnnfn.EvaluateWithGradient(coordinates, 5, gradient, 1);
 
-  REQUIRE(objective == Approx(1.576).epsilon(1e-7));
+  REQUIRE(objective == Approx(1.576).epsilon(eps));
 
-  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(1e-7));
-  REQUIRE(gradient(0, 1) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 0) == Approx(0.0).epsilon(1e-7));
-  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(1e-7));
+  REQUIRE(gradient(0, 0) == Approx(-0.048).epsilon(eps));
+  REQUIRE(gradient(0, 1) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 0) == Approx(0.0).margin(margin));
+  REQUIRE(gradient(1, 1) == Approx(2.0).epsilon(eps));
 }
 
 // Check that final objective value using SGD optimizer is optimal.
-TEST_CASE("LMNNSGDSimpleDatasetTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNSGDSimpleDatasetTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNN<> lmnn(dataset, labels, 1);
+  LMNN<> lmnn(1);
 
-  arma::mat outputMatrix;
-  lmnn.LearnDistance(outputMatrix);
+  arma::Mat<ElemType> outputMatrix;
+  lmnn.LearnDistance(dataset, labels, outputMatrix);
 
   // Ensure that the objective function is better now.
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  double initObj = lmnnfn.Evaluate(arma::eye<arma::mat>(2, 2));
-  double finalObj = lmnnfn.Evaluate(outputMatrix);
+  ElemType initObj = lmnnfn.Evaluate(arma::eye<arma::Mat<ElemType>>(2, 2));
+  ElemType finalObj = lmnnfn.Evaluate(outputMatrix);
 
   // finalObj must be less than initObj.
   REQUIRE(finalObj < initObj);
 }
 
 // Check that final objective value using L-BFGS optimizer is optimal.
-TEST_CASE("LMNNLBFGSSimpleDatasetTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNLBFGSSimpleDatasetTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
-  LMNN<SquaredEuclideanDistance, L_BFGS> lmnn(dataset, labels, 1);
+  LMNN lmnn(1);
 
-  arma::mat outputMatrix;
-  lmnn.LearnDistance(outputMatrix);
+  arma::Mat<ElemType> outputMatrix;
+  ens::L_BFGS lbfgs;
+  lmnn.LearnDistance(dataset, labels, outputMatrix, lbfgs);
 
   // Ensure that the objective function is better now.
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
-  double initObj = lmnnfn.Evaluate(arma::eye<arma::mat>(2, 2));
-  double finalObj = lmnnfn.Evaluate(outputMatrix);
+  ElemType initObj = lmnnfn.Evaluate(arma::eye<arma::Mat<ElemType>>(2, 2));
+  ElemType finalObj = lmnnfn.Evaluate(outputMatrix);
 
   // finalObj must be less than initObj.
   REQUIRE(finalObj < initObj);
 }
 
-double KnnAccuracy(const arma::mat& dataset,
-                   const arma::Row<size_t>& labels,
+template<typename MatType, typename LabelsType>
+double KnnAccuracy(const MatType& dataset,
+                   const LabelsType& labels,
                    const size_t k)
 {
-  arma::Row<size_t> uniqueLabels = arma::unique(labels);
+  typedef typename MatType::elem_type ElemType;
+
+  LabelsType uniqueLabels = arma::unique(labels);
 
   arma::Mat<size_t> neighbors;
-  arma::mat distances;
+  arma::Mat<ElemType> distances;
 
-  KNN knn;
+  NeighborSearch<NearestNeighborSort, EuclideanDistance, MatType> knn;
 
   knn.Train(dataset);
   knn.Search(k, neighbors, distances);
@@ -404,43 +447,44 @@ double KnnAccuracy(const arma::mat& dataset,
 
   for (size_t i = 0; i < dataset.n_cols; ++i)
   {
-    arma::vec Map;
-    Map.zeros(uniqueLabels.n_cols);
+    arma::Col<ElemType> m;
+    m.zeros(uniqueLabels.n_cols);
 
     for (size_t j = 0; j < k; ++j)
-      Map(labels(neighbors(j, i))) +=
-          1 / std::pow(distances(j, i) + 1, 2);
+      m(labels(neighbors(j, i))) += 1 / std::pow(distances(j, i) + 1, 2);
 
-    size_t index = ConvTo<size_t>::From(arma::find(Map
-        == arma::max(Map)));
+    size_t index = ConvTo<size_t>::From(arma::find(m == arma::max(m)));
 
     // Increase count if labels match.
     if (index == labels(i))
         count++;
   }
 
-  // return accuracy.
+  // Return accuracy.
   return ((double) count / dataset.n_cols) * 100;
 }
 
 // Check that final accuracy is greater than initial accuracy on
 // simple dataset.
-TEST_CASE("LMNNAccuracyTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNAccuracyTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
 
   // Taking k = 3 as the case of k = 1 can be easily observed.
   double initAccuracy = KnnAccuracy(dataset, labels, 3);
 
-  LMNN<> lmnn(dataset, labels, 2);
+  LMNN<> lmnn(2);
 
-  arma::mat outputMatrix;
-  lmnn.LearnDistance(outputMatrix);
+  arma::Mat<ElemType> outputMatrix;
+  lmnn.LearnDistance(dataset, labels, outputMatrix);
 
-  double finalAccuracy = KnnAccuracy(outputMatrix * dataset, labels, 3);
+  arma::Mat<ElemType> transformedData = outputMatrix * dataset;
+  double finalAccuracy = KnnAccuracy(transformedData, labels, 3);
 
   // finalObj must be less than initObj.
   REQUIRE(initAccuracy < finalAccuracy);
@@ -452,18 +496,20 @@ TEST_CASE("LMNNAccuracyTest", "[LMNNTest]")
 // Check that accuracy while learning square distance matrix is the same as when
 // we are learning low rank matrix.  I'm ok if this passes only once out of
 // three tries.
-TEST_CASE("LMNNLowRankAccuracyLBFGSTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNLowRankAccuracyLBFGSTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   bool success = false;
   for (size_t trial = 0; trial < 3; ++trial)
   {
-    arma::mat dataPart1;
+    arma::Mat<ElemType> dataPart1;
     dataPart1.randn(5, 50);
 
     arma::Row<size_t> labelsPart1(50);
     labelsPart1.fill(0);
 
-    arma::mat dataPart2;
+    arma::Mat<ElemType> dataPart2;
     dataPart2.randn(5, 50);
 
     arma::Row<size_t> labelsPart2(50);
@@ -473,26 +519,29 @@ TEST_CASE("LMNNLowRankAccuracyLBFGSTest", "[LMNNTest]")
     arma::uvec ordering = arma::shuffle(arma::linspace<arma::uvec>(0, 99, 100));
 
     // Generate datasets.
-    arma::mat dataset = join_rows(dataPart1, dataPart2);
+    arma::Mat<ElemType> dataset = join_rows(dataPart1, dataPart2);
     dataset = dataset.cols(ordering);
 
     // Generate labels.
     arma::Row<size_t> labels = join_rows(labelsPart1, labelsPart2);
     labels = labels.cols(ordering);
 
-    LMNN<SquaredEuclideanDistance, L_BFGS> lmnn(dataset, labels, 1);
+    LMNN<SquaredEuclideanDistance> lmnn(1);
 
     // Learn a square matrix.
-    arma::mat outputMatrix;
-    lmnn.LearnDistance(outputMatrix);
+    arma::Mat<ElemType> outputMatrix;
+    L_BFGS lbfgs;
+    lmnn.LearnDistance(dataset, labels, outputMatrix, lbfgs);
 
-    double acc1 = KnnAccuracy(outputMatrix * dataset, labels, 1);
+    arma::Mat<ElemType> transformedData = outputMatrix * dataset;
+    double acc1 = KnnAccuracy(transformedData, labels, 1);
 
     // Learn a low rank matrix.
-    outputMatrix = arma::randu(4, 5);
-    lmnn.LearnDistance(outputMatrix);
+    outputMatrix = arma::randu<arma::Mat<ElemType>>(4, 5);
+    lmnn.LearnDistance(dataset, labels, outputMatrix, lbfgs);
 
-    double acc2 = KnnAccuracy(outputMatrix * dataset, labels, 1);
+    transformedData = outputMatrix * dataset;
+    double acc2 = KnnAccuracy(transformedData, labels, 1);
 
     // We keep the tolerance very high.  We need to ensure the accuracy drop
     // isn't any more than 10%.
@@ -507,18 +556,20 @@ TEST_CASE("LMNNLowRankAccuracyLBFGSTest", "[LMNNTest]")
 // Check that accuracy while learning square distance matrix is the same as when
 // we are learning low rank matrix.  I'm ok if this passes only once out of
 // three tries.
-TEST_CASE("LMNNLowRankAccuracyTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNLowRankAccuracyTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   bool success = false;
   for (size_t trial = 0; trial < 3; ++trial)
   {
-    arma::mat dataPart1;
+    arma::Mat<ElemType> dataPart1;
     dataPart1.randn(5, 50);
 
     arma::Row<size_t> labelsPart1(50);
     labelsPart1.fill(0);
 
-    arma::mat dataPart2;
+    arma::Mat<ElemType> dataPart2;
     dataPart2.randn(5, 50);
 
     arma::Row<size_t> labelsPart2(50);
@@ -528,26 +579,28 @@ TEST_CASE("LMNNLowRankAccuracyTest", "[LMNNTest]")
     arma::uvec ordering = arma::shuffle(arma::linspace<arma::uvec>(0, 99, 100));
 
     // Generate datasets.
-    arma::mat dataset = join_rows(dataPart1, dataPart2);
+    arma::Mat<ElemType> dataset = join_rows(dataPart1, dataPart2);
     dataset = dataset.cols(ordering);
 
     // Generate labels.
     arma::Row<size_t> labels = join_rows(labelsPart1, labelsPart2);
     labels = labels.cols(ordering);
 
-    LMNN<> lmnn(dataset, labels, 1);
+    LMNN<> lmnn(1);
 
     // Learn a square matrix.
-    arma::mat outputMatrix;
-    lmnn.LearnDistance(outputMatrix);
+    arma::Mat<ElemType> outputMatrix;
+    lmnn.LearnDistance(dataset, labels, outputMatrix);
 
-    double acc1 = KnnAccuracy(outputMatrix * dataset, labels, 1);
+    arma::Mat<ElemType> transformedData = outputMatrix * dataset;
+    double acc1 = KnnAccuracy(transformedData, labels, 1);
 
     // Learn a low rank matrix.
-    outputMatrix = arma::randu(4, 5);
-    lmnn.LearnDistance(outputMatrix);
+    outputMatrix = arma::randu<arma::Mat<ElemType>>(4, 5);
+    lmnn.LearnDistance(dataset, labels, outputMatrix);
 
-    double acc2 = KnnAccuracy(outputMatrix * dataset, labels, 1);
+    transformedData = outputMatrix * dataset;
+    double acc2 = KnnAccuracy(transformedData, labels, 1);
 
     // We keep the tolerance very high.  We need to ensure the accuracy drop
     // isn't any more than 10%.
@@ -621,29 +674,31 @@ TEST_CASE("LMNNLowRankAccuracyBBSGDTest", "[LMNNTest]")
 // Comprehensive gradient tests by Marcus Edel & Ryan Curtin.
 
 // Simple numerical gradient checker.
-template<class FunctionType>
+template<typename FunctionType, typename MatType>
 double CheckGradient(FunctionType& function,
-                     arma::mat& coordinates,
-                     const double eps = 1e-7)
+                     MatType& coordinates,
+                     const typename MatType::elem_type eps = 1e-7)
 {
+  typedef typename MatType::elem_type ElemType;
+
   // Get gradients for the current parameters.
-  arma::mat orgGradient, gradient, estGradient;
+  MatType orgGradient, gradient, estGradient;
   function.Gradient(coordinates, orgGradient);
 
-  estGradient = arma::zeros(orgGradient.n_rows, orgGradient.n_cols);
+  estGradient = arma::zeros<MatType>(orgGradient.n_rows, orgGradient.n_cols);
 
   // Compute numeric approximations to gradient.
   for (size_t i = 0; i < orgGradient.n_elem; ++i)
   {
-    double tmp = coordinates(i);
+    ElemType tmp = coordinates(i);
 
     // Perturb parameter with a positive constant and get costs.
     coordinates(i) += eps;
-    double costPlus = function.Evaluate(coordinates);
+    ElemType costPlus = function.Evaluate(coordinates);
 
     // Perturb parameter with a negative constant and get costs.
     coordinates(i) -= (2 * eps);
-    double costMinus = function.Evaluate(coordinates);
+    ElemType costMinus = function.Evaluate(coordinates);
 
     // Restore the parameter value.
     coordinates(i) = tmp;
@@ -657,74 +712,84 @@ double CheckGradient(FunctionType& function,
       arma::norm(orgGradient + estGradient);
 }
 
-TEST_CASE("LMNNFunctionGradientTest", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNFunctionGradientTest", "[LMNNTest]", float, double)
 {
+  typedef TestType ElemType;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
   // 10 trials with random positions.
   for (size_t i = 0; i < 10; ++i)
   {
-    arma::mat coordinates(2, 2, arma::fill::randn);
+    arma::Mat<ElemType> coordinates(2, 2, arma::fill::randn);
     CheckGradient(lmnnfn, coordinates);
   }
 }
 
-TEST_CASE("LMNNFunctionGradientTest2", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNFunctionGradientTest2", "[LMNNTest]", float, double)
 {
-  // Useful but simple dataset with six points and two classes.
-  arma::mat dataset        = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
-                             " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
-  arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
+  typedef TestType ElemType;
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  // Useful but simple dataset with six points and two classes.
+  arma::Mat<ElemType> dataset = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+                                " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
+  arma::Row<size_t> labels    = " 0    0    0    1    1    1   ";
+
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
   // 10 trials with random positions.
   for (size_t i = 0; i < 10; ++i)
   {
-    arma::mat coordinates(2, 2, arma::fill::randu);
+    arma::Mat<ElemType> coordinates(2, 2, arma::fill::randu);
     CheckGradient(lmnnfn, coordinates);
   }
 }
 
-TEST_CASE("LMNNFunctionGradientTest3", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNFunctionGradientTest3", "[LMNNTest]", float, double)
 {
-  arma::mat dataset;
+  typedef TestType ElemType;
+
+  arma::Mat<ElemType> dataset;
   arma::Row<size_t> labels;
   if (!data::Load("iris.csv", dataset))
     FAIL("Cannot load dataset iris.csv");
   if (!data::Load("iris_labels.txt", labels))
     FAIL("Cannot load dataset iris_labels.txt");
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
   // 10 trials with random positions.
   for (size_t i = 0; i < 10; ++i)
   {
-    arma::mat coordinates(dataset.n_rows, dataset.n_rows, arma::fill::randn);
+    arma::Mat<ElemType> coordinates(dataset.n_rows, dataset.n_rows,
+        arma::fill::randn);
     CheckGradient(lmnnfn, coordinates);
   }
 }
 
-TEST_CASE("LMNNFunctionGradientTest4", "[LMNNTest]")
+TEMPLATE_TEST_CASE("LMNNFunctionGradientTest4", "[LMNNTest]", float, double)
 {
-  arma::mat dataset;
+  typedef TestType ElemType;
+
+  arma::Mat<ElemType> dataset;
   arma::Row<size_t> labels;
   if (!data::Load("iris.csv", dataset))
     FAIL("Cannot load dataset iris.csv");
   if (!data::Load("iris_labels.txt", labels))
     FAIL("Cannot load dataset iris_labels.txt");
 
-  LMNNFunction<> lmnnfn(dataset, labels, 1, 0.6, 1);
+  LMNNFunction<arma::Mat<ElemType>> lmnnfn(dataset, labels, 1, 0.6, 1);
 
   // 10 trials with random positions.
   for (size_t i = 0; i < 10; ++i)
   {
-    arma::mat coordinates(dataset.n_rows, dataset.n_rows, arma::fill::randu);
+    arma::Mat<ElemType> coordinates(dataset.n_rows, dataset.n_rows,
+        arma::fill::randu);
     CheckGradient(lmnnfn, coordinates);
   }
 }
diff --git a/src/mlpack/tests/main_tests/lmnn_test.cpp b/src/mlpack/tests/main_tests/lmnn_test.cpp
index 13b68b3d47..901033856a 100644
--- a/src/mlpack/tests/main_tests/lmnn_test.cpp
+++ b/src/mlpack/tests/main_tests/lmnn_test.cpp
@@ -542,7 +542,7 @@ TEST_CASE_METHOD(LMNNTestFixture, "LMNNDiffRegularizationTest",
 }
 
 /**
- * Ensure that different value of range results in a
+ * Ensure that different value of update interval results in a
  * different output matrix.
  */
 TEST_CASE_METHOD(LMNNTestFixture, "LMNNDiffRangeTest",
@@ -573,7 +573,7 @@ TEST_CASE_METHOD(LMNNTestFixture, "LMNNDiffRangeTest",
   SetInputParam("input", std::move(inputData));
   SetInputParam("labels", std::move(labels));
   SetInputParam("linear_scan",  (bool) true);
-  SetInputParam("range", 100);
+  SetInputParam("update_interval", 100);
 
   RUN_BINDING();
 
@@ -674,9 +674,9 @@ TEST_CASE_METHOD(LMNNTestFixture, "LMNNDiffPassesTest",
 }
 
 /**
- * Ensure that number of targets, range, batch size must be always positive
- * and regularization, step size, max iterations, rank, passes & tolerance are
- * always non-negative 
+ * Ensure that number of targets, update interval, batch size must be always
+ * positive and regularization, step size, max iterations, rank, passes &
+ * tolerance are always non-negative.
  */
 TEST_CASE_METHOD(LMNNTestFixture, "LMNNBoundsTest",
                 "[LMNNMainTest][BindingTests]")
@@ -701,12 +701,12 @@ TEST_CASE_METHOD(LMNNTestFixture, "LMNNBoundsTest",
   // Reset settings.
   ResetSettings();
 
-  // Test for range value.
+  // Test for update interval value.
 
   // Input training data.
   SetInputParam("input", inputData);
   SetInputParam("labels", labels);
-  SetInputParam("range", (int) 0);
+  SetInputParam("update_interval", (int) 0);
 
   REQUIRE_THROWS_AS(RUN_BINDING(), std::runtime_error);
 
diff --git a/src/mlpack/tests/nca_test.cpp b/src/mlpack/tests/nca_test.cpp
index 2dcfc0909d..20df35df87 100644
--- a/src/mlpack/tests/nca_test.cpp
+++ b/src/mlpack/tests/nca_test.cpp
@@ -26,26 +26,31 @@ using namespace ens;
  * The Softmax error function should return the identity matrix as its initial
  * point.
  */
-TEST_CASE("SoftmaxInitialPoint", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxInitialPoint", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Cheap fake dataset.
-  arma::mat data;
+  arma::Mat<eT> data;
   data.randu(5, 5);
   arma::Row<size_t> labels;
   labels.zeros(5);
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
   // Verify the initial point is the identity matrix.
-  arma::mat initialPoint = sef.GetInitialPoint();
+  arma::Mat<eT> initialPoint = sef.GetInitialPoint();
+  const double eps = std::is_same<eT, float>::value ? 1e-4 : 1e-7;
+  const double margin = std::is_same<eT, float>::value ? 1e-4 : 1e-5;
   for (int row = 0; row < 5; row++)
   {
     for (int col = 0; col < 5; col++)
     {
       if (row == col)
-        REQUIRE(initialPoint(row, col) == Approx(1.0).epsilon(1e-7));
+        REQUIRE(initialPoint(row, col) == Approx(1.0).epsilon(eps));
       else
-        REQUIRE(initialPoint(row, col) == Approx(0.0).margin(1e-5));
+        REQUIRE(initialPoint(row, col) == Approx(0.0).margin(margin));
     }
   }
 }
@@ -54,16 +59,19 @@ TEST_CASE("SoftmaxInitialPoint", "[NCATesT]")
  * On a simple fake dataset, ensure that the initial function evaluation is
  * correct.
  */
-TEST_CASE("SoftmaxInitialEvaluation", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxInitialEvaluation", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  double objective = sef.Evaluate(arma::eye<arma::mat>(2, 2));
+  eT objective = sef.Evaluate(arma::eye<arma::Mat<eT>>(2, 2));
 
   // Result painstakingly calculated by hand by rcurtin (recorded forever in his
   // notebook).  As a result of lack of precision of the by-hand result, the
@@ -75,22 +83,28 @@ TEST_CASE("SoftmaxInitialEvaluation", "[NCATesT]")
  * On a simple fake dataset, ensure that the initial gradient evaluation is
  * correct.
  */
-TEST_CASE("SoftmaxInitialGradient", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxInitialGradient", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  arma::mat gradient;
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
+  arma::Mat<eT> gradient;
+  arma::Mat<eT> coordinates(2, 2, arma::fill::eye);
   sef.Gradient(coordinates, gradient);
 
   // Results painstakingly calculated by hand by rcurtin (recorded forever in
   // his notebook).  As a result of lack of precision of the by-hand result, the
   // tolerance is fairly high.
+  //
+  // UPDATE 2024: that notebook definitely got thrown away over a decade ago.  I
+  // don't even remember what it looked like.
   REQUIRE(gradient(0, 0) == Approx(-0.089766).epsilon(0.0005));
   REQUIRE(gradient(1, 0) == Approx(0.0).margin(1e-5));
   REQUIRE(gradient(0, 1) == Approx(0.0).margin(1e-5));
@@ -101,36 +115,43 @@ TEST_CASE("SoftmaxInitialGradient", "[NCATesT]")
  * On optimally separated datasets, ensure that the objective function is
  * optimal (equal to the negative number of points).
  */
-TEST_CASE("SoftmaxOptimalEvaluation", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxOptimalEvaluation", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Simple optimal dataset.
-  arma::mat data           = " 500  500 -500 -500;"
+  arma::Mat<eT> data       = " 500  500 -500 -500;"
                              "   1    0    1    0 ";
   arma::Row<size_t> labels = "   0    0    1    1 ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  double objective = sef.Evaluate(arma::eye<arma::mat>(2, 2));
+  eT objective = sef.Evaluate(arma::eye<arma::Mat<eT>>(2, 2));
 
   // Use a very close tolerance for optimality; we need to be sure this function
   // gives optimal results correctly.
-  REQUIRE(objective == Approx(-4.0).epsilon(1e-12));
+  const double eps = std::is_same<eT, float>::value ? 1e-6 : 1e-12;
+  REQUIRE(objective == Approx(-4.0).epsilon(eps));
 }
 
 /**
  * On optimally separated datasets, ensure that the gradient is zero.
  */
-TEST_CASE("SoftmaxOptimalGradient", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxOptimalGradient", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Simple optimal dataset.
-  arma::mat data           = " 500  500 -500 -500;"
+  arma::Mat<eT> data       = " 500  500 -500 -500;"
                              "   1    0    1    0 ";
   arma::Row<size_t> labels = "   0    0    1    1 ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  arma::mat gradient;
-  sef.Gradient(arma::eye<arma::mat>(2, 2), gradient);
+  arma::Mat<eT> gradient;
+  sef.Gradient(arma::eye<arma::Mat<eT>>(2, 2), gradient);
 
   REQUIRE(gradient(0, 0) == Approx(0.0).margin(1e-5));
   REQUIRE(gradient(0, 1) == Approx(0.0).margin(1e-5));
@@ -141,19 +162,22 @@ TEST_CASE("SoftmaxOptimalGradient", "[NCATesT]")
 /**
  * Ensure the separable objective function is right.
  */
-TEST_CASE("SoftmaxSeparableObjective", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxSeparableObjective", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
   // Results painstakingly calculated by hand by rcurtin (recorded forever in
   // his notebook).  As a result of lack of precision of the by-hand result, the
   // tolerance is fairly high.
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
+  arma::Mat<eT> coordinates = arma::eye<arma::Mat<eT>>(2, 2);
   REQUIRE(sef.Evaluate(coordinates, 0, 1) == Approx(-0.22480).epsilon(0.0001));
   REQUIRE(sef.Evaluate(coordinates, 1, 1) == Approx(-0.30613).epsilon(0.0001));
   REQUIRE(sef.Evaluate(coordinates, 2, 1) == Approx(-0.22480).epsilon(0.0001));
@@ -165,16 +189,19 @@ TEST_CASE("SoftmaxSeparableObjective", "[NCATesT]")
 /**
  * Ensure the optimal separable objective function is right.
  */
-TEST_CASE("OptimalSoftmaxSeparableObjective", "[NCATesT]")
+TEMPLATE_TEST_CASE("OptimalSoftmaxSeparableObjective", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Simple optimal dataset.
-  arma::mat data           = " 500  500 -500 -500;"
+  arma::Mat<eT> data       = " 500  500 -500 -500;"
                              "   1    0    1    0 ";
   arma::Row<size_t> labels = "   0    0    1    1 ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
+  arma::Mat<eT> coordinates = arma::eye<arma::Mat<eT>>(2, 2);
 
   // Use a very close tolerance for optimality; we need to be sure this function
   // gives optimal results correctly.
@@ -187,17 +214,20 @@ TEST_CASE("OptimalSoftmaxSeparableObjective", "[NCATesT]")
 /**
  * Ensure the separable gradient is right.
  */
-TEST_CASE("SoftmaxSeparableGradient", "[NCATesT]")
+TEMPLATE_TEST_CASE("SoftmaxSeparableGradient", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  arma::mat coordinates = arma::eye<arma::mat>(2, 2);
-  arma::mat gradient(2, 2);
+  arma::Mat<eT> coordinates = arma::eye<arma::Mat<eT>>(2, 2);
+  arma::Mat<eT> gradient(2, 2);
 
   sef.Gradient(coordinates, 0, gradient, 1);
 
@@ -250,29 +280,33 @@ TEST_CASE("SoftmaxSeparableGradient", "[NCATesT]")
  * On our simple dataset, ensure that the NCA algorithm fully separates the
  * points.
  */
-TEST_CASE("NCASGDSimpleDataset", "[NCATesT]")
+TEMPLATE_TEST_CASE("NCASGDSimpleDataset", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
   // Huge learning rate because this is so simple.
-  NCA<SquaredEuclideanDistance> nca(data, labels);
-  nca.Optimizer().StepSize() = 1.2;
-  nca.Optimizer().MaxIterations() = 300000;
-  nca.Optimizer().Tolerance() = 0;
-  nca.Optimizer().Shuffle() = true;
+  ens::StandardSGD opt;
+  opt.StepSize() = 1.2;
+  opt.MaxIterations() = 300000;
+  opt.Tolerance() = 0;
+  opt.Shuffle() = true;
 
-  arma::mat outputMatrix;
-  nca.LearnDistance(outputMatrix);
+  arma::Mat<eT> outputMatrix;
+  NCA nca;
+  nca.LearnDistance(data, labels, outputMatrix, opt);
 
   // Ensure that the objective function is better now.
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  double initObj = sef.Evaluate(arma::eye<arma::mat>(2, 2));
-  double finalObj = sef.Evaluate(outputMatrix);
-  arma::mat finalGradient;
+  eT initObj = sef.Evaluate(arma::eye<arma::Mat<eT>>(2, 2));
+  eT finalObj = sef.Evaluate(outputMatrix);
+  arma::Mat<eT> finalGradient;
   sef.Gradient(outputMatrix, finalGradient);
 
   // finalObj must be less than initObj.
@@ -284,33 +318,36 @@ TEST_CASE("NCASGDSimpleDataset", "[NCATesT]")
   REQUIRE(arma::norm(finalGradient, 2) < 1e-4);
 }
 
-TEST_CASE("NCALBFGSSimpleDataset", "[NCATesT]")
+TEMPLATE_TEST_CASE("NCALBFGSSimpleDataset", "[NCATest]", float, double)
 {
+  typedef TestType eT;
+
   // Useful but simple dataset with six points and two classes.
-  arma::mat data           = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
+  arma::Mat<eT> data       = "-0.1 -0.1 -0.1  0.1  0.1  0.1;"
                              " 1.0  0.0 -1.0  1.0  0.0 -1.0 ";
   arma::Row<size_t> labels = " 0    0    0    1    1    1   ";
 
-  // Huge learning rate because this is so simple.
-  NCA<SquaredEuclideanDistance, L_BFGS> nca(data, labels);
-  nca.Optimizer().NumBasis() = 5;
+  L_BFGS lbfgs;
+  lbfgs.NumBasis() = 5;
 
-  arma::mat outputMatrix;
-  nca.LearnDistance(outputMatrix);
+  arma::Mat<eT> outputMatrix;
+  NCA nca;
+  nca.LearnDistance(data, labels, outputMatrix, lbfgs);
 
   // Ensure that the objective function is better now.
-  SoftmaxErrorFunction<SquaredEuclideanDistance> sef(data, labels);
+  SoftmaxErrorFunction<arma::Mat<eT>, arma::Row<size_t>,
+      SquaredEuclideanDistance> sef(data, labels);
 
-  double initObj = sef.Evaluate(arma::eye<arma::mat>(2, 2));
-  double finalObj = sef.Evaluate(outputMatrix);
-  arma::mat finalGradient;
+  eT initObj = sef.Evaluate(arma::eye<arma::Mat<eT>>(2, 2));
+  eT finalObj = sef.Evaluate(outputMatrix);
+  arma::Mat<eT> finalGradient;
   sef.Gradient(outputMatrix, finalGradient);
 
   // finalObj must be less than initObj.
   REQUIRE(finalObj < initObj);
   // Verify that final objective is optimal.
-  REQUIRE(finalObj == Approx(-6.0).epsilon(1e-7));
+  REQUIRE(finalObj == Approx(-6.0).epsilon(0.00001));
   // The solution is not unique, so the best we can do is ensure the gradient
   // norm is close to 0.
-  REQUIRE(arma::norm(finalGradient, 2) < 1e-6);
+  REQUIRE(arma::norm(finalGradient, 2) < 1e-5);
 }