diff --git a/diffpy/snmf/subroutines.py b/diffpy/snmf/subroutines.py
new file mode 100644
index 0000000..093e3b1
--- /dev/null
+++ b/diffpy/snmf/subroutines.py
@@ -0,0 +1,51 @@
+import numpy as np
+
+
+# import scipy.interpolate
+
+
+def objective_function(residual_matrix, stretching_factor_matrix, smoothness, smoothness_term, component_matrix,
+                       sparsity):
+    """Defines the objective function of the algorithm and returns its value.
+
+    Calculates the value of '(||residual_matrix||_F) ** 2 + smoothness * (||smoothness_term *
+    stretching_factor_matrix.T||)**2 + sparsity * sum(component_matrix ** .5)' and returns its value.
+
+    Parameters
+    ----------
+    residual_matrix: 2d array like
+      The matrix where each column is the difference between an experimental PDF/XRD pattern and a calculated PDF/XRD
+      pattern at each grid point. Has dimensions R x M where R is the length of each pattern and M is the amount of
+      patterns.
+
+    stretching_factor_matrix: 2d array like
+      The matrix containing the stretching factors of the calculated component signal. Has dimensions K x M where K is
+      the amount of components and M is the number of experimental PDF/XRD patterns.
+
+    smoothness: float
+      The coefficient of the smoothness term which determines the intensity of the smoothness term and its behavior.
+      It is not very sensitive and is usually adjusted by multiplying it by ten.
+
+    smoothness_term: 2d array like
+      The regularization term that ensures that smooth changes in the component stretching signals are favored.
+      Has dimensions (M-2) x M where M is the amount of experimentally obtained PDF/XRD patterns, the moment amount.
+
+    component_matrix: 2d array like
+      The matrix containing the calculated component signals of the experimental PDF/XRD patterns. Has dimensions R x K
+      where R is the signal length and K is the number of component signals.
+
+    sparsity: float
+      The parameter determining the intensity of the sparsity regularization term which enables the algorithm to
+      exploit the sparse nature of XRD data. It is usually adjusted by doubling.
+
+    Returns
+    -------
+    float
+      The value of the objective function.
+
+    """
+    residual_matrix = np.asarray(residual_matrix)
+    stretching_factor_matrix = np.asarray(stretching_factor_matrix)
+    component_matrix = np.asarray(component_matrix)
+    return .5 * np.linalg.norm(residual_matrix, 'fro') ** 2 + .5 * smoothness * np.linalg.norm(
+        smoothness_term @ stretching_factor_matrix.T, 'fro') ** 2 + sparsity * np.sum(np.sqrt(component_matrix))
diff --git a/diffpy/snmf/tests/test_subroutines.py b/diffpy/snmf/tests/test_subroutines.py
new file mode 100644
index 0000000..fffaada
--- /dev/null
+++ b/diffpy/snmf/tests/test_subroutines.py
@@ -0,0 +1,16 @@
+import pytest
+from diffpy.snmf.subroutines import objective_function
+
+to = [
+    ([[[1, 2], [3, 4]], [[5, 6], [7, 8]], 1e11, [[1, 2], [3, 4]], [[1, 2], [3, 4]], 1], 2.574e14),
+    ([[[11, 2], [31, 4]], [[5, 63], [7, 18]], .001, [[21, 2], [3, 4]], [[11, 22], [3, 40]], 1], 650.4576),
+    ([[[1, 2], [3, 4]], [[5, 6], [7, 8]], 1e11, [[1, 2], [3, 4]], [[1, 2], [3, 4]], 0], 2.574e14),
+
+]
+
+
+@pytest.mark.parametrize("to", to)
+def test_objective_function(to):
+    actual = objective_function(to[0][0], to[0][1], to[0][2], to[0][3], to[0][4], to[0][5])
+    expected = to[1]
+    assert actual == pytest.approx(expected)