diff --git a/src/pyrecest/tracking/ellipse_geometry.py b/src/pyrecest/tracking/ellipse_geometry.py
index 163e55f63..489ca7230 100644
--- a/src/pyrecest/tracking/ellipse_geometry.py
+++ b/src/pyrecest/tracking/ellipse_geometry.py
@@ -37,7 +37,7 @@ def symmetrize(matrix):
     """Return the symmetric part of ``matrix``."""
 
     matrix = asarray(matrix)
-    return 0.5 * (matrix + matrix.T)
+    return 0.5 * matrix + 0.5 * matrix.T
 
 
 def project_symmetric_covariance(covariance, minimum_eigenvalue=0.0):
diff --git a/tests/tracking/test_ellipse_geometry_extreme_symmetrization.py b/tests/tracking/test_ellipse_geometry_extreme_symmetrization.py
new file mode 100644
index 000000000..c9210b529
--- /dev/null
+++ b/tests/tracking/test_ellipse_geometry_extreme_symmetrization.py
@@ -0,0 +1,28 @@
+from __future__ import annotations
+
+import numpy as np
+import numpy.testing as npt
+from pyrecest.tracking.ellipse_geometry import (
+    project_symmetric_covariance,
+    symmetrize,
+)
+
+
+def test_symmetrize_preserves_extreme_finite_diagonal() -> None:
+    covariance = np.diag([1.0e308, 2.0e307])
+
+    with np.errstate(over="raise", invalid="raise"):
+        symmetric = np.asarray(symmetrize(covariance))
+
+    assert np.all(np.isfinite(symmetric))
+    npt.assert_array_equal(symmetric, covariance)
+
+
+def test_covariance_projection_preserves_extreme_finite_diagonal() -> None:
+    covariance = np.diag([1.0e308, 2.0e307])
+
+    with np.errstate(over="raise", invalid="raise"):
+        projected = np.asarray(project_symmetric_covariance(covariance))
+
+    assert np.all(np.isfinite(projected))
+    npt.assert_array_equal(projected, covariance)