diff --git a/pylinalg/matrix.py b/pylinalg/matrix.py
index 65af6f6..f0d620b 100644
--- a/pylinalg/matrix.py
+++ b/pylinalg/matrix.py
@@ -379,7 +379,11 @@ def mat_decompose(matrix, /, *, scaling_signs=None, dtype=None, out=None):
     scaling *= scaling_signs
 
     rotation = out[1] if out is not None else None
-    rotation_matrix = matrix[:-1, :-1] * (1 / scaling)[None, :]
+
+    rotation_matrix = matrix[:-1, :-1].copy().astype(float)
+    mask = scaling != 0
+    rotation_matrix[:, mask] /= scaling[mask][None, :]
+    rotation_matrix[:, ~mask] = 0.0
     rotation = quat_from_mat(rotation_matrix, out=rotation, dtype=dtype)
 
     return translation, rotation, scaling
diff --git a/tests/conftest.py b/tests/conftest.py
index 6436b93..84cc308 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -10,7 +10,7 @@
 import pylinalg as la
 
 
-def pytest_report_header(config, start_path, startdir):
+def pytest_report_header(config):
     # report the CPU model to allow detecting platform-specific problems
     if platform.system() == "Windows":
         name = (
diff --git a/tests/test_matrix.py b/tests/test_matrix.py
index 1fe6647..bfcaa73 100644
--- a/tests/test_matrix.py
+++ b/tests/test_matrix.py
@@ -192,6 +192,22 @@ def test_mat_decompose():
     npt.assert_array_almost_equal(rotation, [0, 0, np.sqrt(2) / 2, np.sqrt(2) / 2])
 
 
+def test_mat_decompose_scaling_0():
+    """Test that the matrices are decomposed correctly when scaling is 0."""
+
+    scaling = [0, 0, 2]
+    rotation = [0, 0, np.sqrt(2) / 2, np.sqrt(2) / 2]
+    translation = [2, 2, 2]
+
+    matrix = la.mat_compose(translation, rotation, scaling)
+    translation_, rotation_, scaling_ = la.mat_decompose(matrix)
+
+    npt.assert_array_almost_equal(translation_, translation)
+    npt.assert_array_almost_equal(scaling_, scaling)
+    # rotation is not uniquely defined when scaling is 0, but it should not be NaN
+    assert not np.isnan(rotation_).any()
+
+
 @pytest.mark.parametrize(
     "signs",
     [