diff --git a/RELEASES.md b/RELEASES.md
index f91b42bc0..febbe4b64 100644
--- a/RELEASES.md
+++ b/RELEASES.md
@@ -1,5 +1,12 @@
 # Releases
 
+## 0.9.8dev
+
+#### Closed issues
+
+- Fix mean centering in `ot.dr.fda` and `ot.dr.wda`: `np.mean(X)` returned a scalar instead of the per-feature mean, so `proj` did not center the data as documented. In `ot.dr.fda` the same pattern in the class means made the between-class scatter matrix independent of which features separate the classes, and FDA returned a non-discriminant direction (PR #840)
+- `ot.dr.fda` and `ot.dr.wda` no longer modify the input array `X` in place (PR #840)
+
 ## 0.9.7.post1
 
 This release is identical to 0.9.7 but will allow the upload of a source distribution to PyPI and release on conda-forge (that requires a source distribution).
diff --git a/ot/dr.py b/ot/dr.py
index 4914e6b84..9cfabff5f 100644
--- a/ot/dr.py
+++ b/ot/dr.py
@@ -101,8 +101,8 @@ def fda(X, y, p=2, reg=1e-16):
         projection function including mean centering
     """
 
-    mx = np.mean(X)
-    X -= mx.reshape((1, -1))
+    mx = np.mean(X, axis=0)
+    X = X - mx.reshape((1, -1))
 
     # data split between classes
     d = X.shape[1]
@@ -119,7 +119,7 @@ def fda(X, y, p=2, reg=1e-16):
     mxc = np.zeros((d, nc))
 
     for i in range(nc):
-        mxc[:, i] = np.mean(xc[i])
+        mxc[:, i] = np.mean(xc[i], axis=0)
 
     mx0 = np.mean(mxc, 1)
     Cb = 0
@@ -218,8 +218,8 @@ def wda(
     else:
         raise ValueError("Unknown Sinkhorn method '%s'." % sinkhorn_method)
 
-    mx = np.mean(X)
-    X -= mx.reshape((1, -1))
+    mx = np.mean(X, axis=0)
+    X = X - mx.reshape((1, -1))
 
     # data split between classes
     d = X.shape[1]
diff --git a/test/test_dr.py b/test/test_dr.py
index 2ac026d23..dcb477717 100644
--- a/test/test_dr.py
+++ b/test/test_dr.py
@@ -39,6 +39,51 @@ def test_fda():
     np.testing.assert_allclose(np.sum(Pfda**2, 0), np.ones(p))
 
 
+@pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)")
+def test_fda_recovers_discriminant_direction():
+    # classes are separated along the first feature only, all others are noise,
+    # so FDA must return a direction aligned with e_0
+    rng = np.random.RandomState(1)
+    n_features = 5
+    xs = np.concatenate(
+        [
+            rng.randn(60, n_features) * 0.3 + shift * np.eye(1, n_features, 0)
+            for shift in [-6.0, 0.0, 6.0]
+        ]
+    )
+    ys = np.repeat([0, 1, 2], 60)
+
+    Pfda, _ = ot.dr.fda(xs, ys, p=1)
+
+    direction = Pfda[:, 0] / np.linalg.norm(Pfda[:, 0])
+    np.testing.assert_array_less(0.95, np.abs(direction[0]))
+
+
+@pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)")
+def test_fda_projection_is_centered():
+    rng = np.random.RandomState(0)
+    xs, ys = ot.datasets.make_data_classif("gaussrot", 90, random_state=rng)
+    xs = xs + 10.0  # off-centered data makes an absent centering visible
+
+    _, projfda = ot.dr.fda(xs, ys, p=1)
+
+    np.testing.assert_allclose(projfda(xs).mean(axis=0), 0.0, atol=1e-10)
+
+
+@pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)")
+def test_fda_wda_do_not_modify_input():
+    rng = np.random.RandomState(0)
+    xs, ys = ot.datasets.make_data_classif("gaussrot", 90, random_state=rng)
+    xs = xs + 10.0
+
+    xs_copy = xs.copy()
+    ot.dr.fda(xs, ys, p=1)
+    np.testing.assert_allclose(xs, xs_copy)
+
+    ot.dr.wda(xs, ys, p=1, reg=1.0, k=5, maxiter=5)
+    np.testing.assert_allclose(xs, xs_copy)
+
+
 @pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)")
 def test_wda():
     n_samples = 100  # nb samples in source and target datasets