diff --git a/pdme/subspace_simulation/mcmc_costs.py b/pdme/subspace_simulation/mcmc_costs.py
index 2891998..4eac4bb 100644
--- a/pdme/subspace_simulation/mcmc_costs.py
+++ b/pdme/subspace_simulation/mcmc_costs.py
@@ -18,3 +18,13 @@ def proportional_costs_vs_actual_measurement(
 		dot_inputs_array, dipoles_to_test
 	)
 	return proportional_cost(actual_measurement_array, vals)
+
+
+def relative_square_diffs(
+	approx: numpy.ndarray, target: numpy.ndarray
+) -> numpy.ndarray:
+	# Assume that both approx and target are arrays of length m
+	# Approx can broadcast if additional indexes to the left
+	# diffs.shape = [ m ]
+	diffs = (approx - target) ** 2 / (target**2)
+	return diffs.sum(axis=-1)
diff --git a/tests/subspace_simulation/test_costs.py b/tests/subspace_simulation/test_costs.py
index 3625c7b..aef46b9 100644
--- a/tests/subspace_simulation/test_costs.py
+++ b/tests/subspace_simulation/test_costs.py
@@ -1,4 +1,5 @@
 import pdme.subspace_simulation
+import pdme.subspace_simulation.mcmc_costs
 import numpy
 
 
@@ -8,3 +9,23 @@ def test_proportional_costs(snapshot):
 
 	actual_result = pdme.subspace_simulation.proportional_cost(a, b).tolist()
 	assert actual_result == snapshot
+
+
+def test_squared_costs_manual():
+	target = numpy.array([100, 400, 900])
+	approx1 = numpy.array([0, 400, 800])
+	approx2 = numpy.array([200, 400, 600])
+
+	expected1 = 1.0123456790123457
+	expected2 = 1.1111111111111111
+
+	actual1 = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(approx1, target)
+	assert actual1 == expected1
+
+	actual2 = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(approx2, target)
+	assert actual2 == expected2
+
+	combined_actual = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(
+		numpy.array([approx1, approx2]), target
+	)
+	numpy.testing.assert_allclose(combined_actual, [expected1, expected2], rtol=1e-14)