diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py
index 73c93169..9a549c16 100644
--- a/pyerrors/linalg.py
+++ b/pyerrors/linalg.py
@@ -1,3 +1,4 @@
+from packaging import version
 import numpy as np
 import autograd.numpy as anp  # Thinly-wrapped numpy
 from .obs import derived_observable, CObs, Obs, import_jackknife
@@ -260,6 +261,8 @@ def _mat_mat_op(op, obs, **kwargs):
 
 def eigh(obs, **kwargs):
     """Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh."""
+    if version.parse(np.__version__) >= version.parse("1.25.0"):
+        raise  NotImplementedError("eigh error propagation is not working with numpy>=1.25 and autograd==1.5.")
     w = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[0], obs)
     v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs)
     return w, v
@@ -278,6 +281,8 @@ def pinv(obs, **kwargs):
 
 def svd(obs, **kwargs):
     """Computes the singular value decomposition of a matrix of Obs."""
+    if version.parse(np.__version__) >= version.parse("1.25.0"):
+        raise  NotImplementedError("svd error propagation is not working with numpy>=1.25 and autograd==1.5.")
     u = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[0], obs)
     s = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[1], obs)
     vh = derived_observable(lambda x, **kwargs: anp.linalg.svd(x, full_matrices=False)[2], obs)
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
index babf3797..149bf92a 100644
--- a/tests/linalg_test.py
+++ b/tests/linalg_test.py
@@ -1,3 +1,4 @@
+from packaging import version
 import numpy as np
 import autograd.numpy as anp
 import math
@@ -291,23 +292,26 @@ def test_matrix_functions():
         diff = entry - sym[i, j]
         assert diff.is_zero()
 
-    # Check eigh
-    e, v = pe.linalg.eigh(sym)
-    for i in range(dim):
-        tmp = sym @ v[:, i] - v[:, i] * e[i]
-        for j in range(dim):
-            assert tmp[j].is_zero()
+    # These linalg functions don't work with numpy>=1.25 and autograd==1.5.
+    # Remove this guard once this is fixed in autograd.
+    if version.parse(np.__version__) < version.parse("1.25.0"):
+        # Check eigh
+        e, v = pe.linalg.eigh(sym)
+        for i in range(dim):
+            tmp = sym @ v[:, i] - v[:, i] * e[i]
+            for j in range(dim):
+                assert tmp[j].is_zero()
 
-    # Check eig function
-    e2 = pe.linalg.eig(sym)
-    assert np.all(np.sort(e) == np.sort(e2))
+        # Check eig function
+        e2 = pe.linalg.eig(sym)
+        assert np.all(np.sort(e) == np.sort(e2))
 
-    # Check svd
-    u, v, vh = pe.linalg.svd(sym)
-    diff = sym - u @ np.diag(v) @ vh
+        # Check svd
+        u, v, vh = pe.linalg.svd(sym)
+        diff = sym - u @ np.diag(v) @ vh
 
-    for (i, j), entry in np.ndenumerate(diff):
-        assert entry.is_zero()
+        for (i, j), entry in np.ndenumerate(diff):
+            assert entry.is_zero()
 
     # Check determinant
     assert pe.linalg.det(np.diag(np.diag(matrix))) == np.prod(np.diag(matrix))
diff --git a/tests/mpm_test.py b/tests/mpm_test.py
index 165f09bc..9ab85968 100644
--- a/tests/mpm_test.py
+++ b/tests/mpm_test.py
@@ -1,3 +1,4 @@
+from packaging import version
 import numpy as np
 import pyerrors as pe
 import pytest
@@ -5,10 +6,11 @@ import pytest
 np.random.seed(0)
 
 
-def test_mpm():
-    corr_content = []
-    for t in range(8):
-        f = 0.8 * np.exp(-0.4 * t)
-        corr_content.append(pe.pseudo_Obs(np.random.normal(f, 1e-2 * f), 1e-2 * f, 't'))
+if version.parse(np.__version__) < version.parse("1.25.0"):
+    def test_mpm():
+        corr_content = []
+        for t in range(8):
+            f = 0.8 * np.exp(-0.4 * t)
+            corr_content.append(pe.pseudo_Obs(np.random.normal(f, 1e-2 * f), 1e-2 * f, 't'))
 
-    res = pe.mpm.matrix_pencil_method(corr_content)
+        res = pe.mpm.matrix_pencil_method(corr_content)