diff --git a/src/algorithms/grassmann.jl b/src/algorithms/grassmann.jl
index 7da1a0133..bba925521 100644
--- a/src/algorithms/grassmann.jl
+++ b/src/algorithms/grassmann.jl
@@ -81,8 +81,7 @@ function precondition(state, g)
     rtolmin = eps(real(scalartype(state)))^(3 / 4)
     tforeach(eachindex(state); scheduler = MPSKit.Defaults.scheduler[]) do i
         rtol = max(rtolmin, norm(g[i]))
-        ρ = rho_inv_regularized(state.C[i]; rtol)
-        g′[i] = rmul(g[i], ρ)
+        g′[i] = rmul_rho_inv_regularized(g[i], state.C[i]; rtol)
         return nothing
     end
     return g′
@@ -234,15 +233,22 @@ function fg(
 end
 
 """
-    rho_inv_regularized(C; rtol = eps(real(scalartype(C)))^(3 / 4))
+    rmul_rho_inv_regularized(Δ, C; rtol = eps(real(scalartype(C)))^(3 / 4))
 
-Compute the (regularized) inverse of the MPS fixed point `ρ = C * C'`.
-Here we use the Tikhonov regularization, i.e. `inv(ρ) = inv(C * C' + δ²1)`,
+Right-multiply the tangent vector `Δ` with the (regularized) inverse of the MPS fixed point
+`ρ = C * C'`. Here we use the Tikhonov regularization, i.e. `inv(ρ) = inv(C * C' + δ²1)`,
 where the regularization parameter is `δ = rtol * norm(C)`.
+
+The inverse is applied in factored form, `((Δ.Z * U) * inv(S² + δ²)) * U'`, rather than formed
+explicitly: its condition number can reach `1 / rtol²`, and once that exceeds `1 / eps` the
+explicit matrix `U * inv(S² + δ²) * U'` is no longer numerically positive definite, such that
+the preconditioned gradient need not be a descent direction.
 """
-function rho_inv_regularized(C; rtol = eps(real(scalartype(C)))^(3 / 4))
+function rmul_rho_inv_regularized(Δ::GrassmannTangent, C; rtol = eps(real(scalartype(C)))^(3 / 4))
     U, S, _ = svd_compact(C)
-    return U * pinv_tikhonov!!(S; rtol) * U'
+    # the diagonal scaling has to happen before mixing back with U' to preserve definiteness
+    Z′ = (Δ.Z * U) * pinv_tikhonov!!(S; rtol)
+    return GrassmannTangent(Δ.W, Z′ * U')
 end
 
 function pinv_tikhonov!!(S::DiagonalTensorMap{<:Real}; rtol = zero(scalartype(S)))
diff --git a/test/groundstate/groundstate.jl b/test/groundstate/groundstate.jl
index b31b1512a..527947e05 100644
--- a/test/groundstate/groundstate.jl
+++ b/test/groundstate/groundstate.jl
@@ -154,9 +154,6 @@ verbosity_conv = 1
                 ψ₀, H, GradientGrassmann(; verbosity = verbosity_full, maxiter = 2)
             )
 
-            # an explicit `tol` keeps the optimizer from overshooting past the point where the
-            # gradient is floating-point noise: pushed further, the CG line search can hit a
-            # 0/0 in its step-size formula and feed a NaN tangent into the Grassmann retraction
             ψ, envs, δ = find_groundstate(
                 ψ, H, GradientGrassmann(; tol, verbosity = verbosity_conv, maxiter = 50), envs
             )
@@ -195,6 +192,26 @@ verbosity_conv = 1
     end
 end
 
+@testset "GradientGrassmann preconditioner is positive definite" begin
+    # nearly rank-deficient bonds, as for the long-range model above, make the regularized
+    # inverse of `ρ = C * C'` ill-conditioned beyond `1 / eps` once the gradient is small
+    Grassmann = MPSKit.GrassmannMPS.Grassmann
+    Random.seed!(123)
+    V, P = ℙ^6, ℙ^2
+    S = DiagonalTensorMap([1.0, 0.06, 2.0e-10, 1.5e-10, 5.0e-11, 1.0e-11], V)
+    @testset "$T" for T in (Float64, ComplexF64)
+        isdescent = map(1:50) do _
+            AL = randisometry(T, V ⊗ P, V)
+            C = randisometry(T, V, V) * S * randisometry(T, V, V)
+            Δ = MPSKit.GrassmannMPS.rmul(Grassmann.project(randn(T, V ⊗ P, V), AL), C')
+            Δ = MPSKit.GrassmannMPS.scale!(Δ, 1.0e-9 / norm(Δ.Z))
+            PΔ = MPSKit.GrassmannMPS.rmul_rho_inv_regularized(Δ, C; rtol = norm(Δ.Z))
+            return real(Grassmann.inner(AL, Δ, PΔ)) > 0
+        end
+        @test all(isdescent)
+    end
+end
+
 @testset "InfiniteMPS ground state" verbose = true begin
     tol = 1.0e-8
     g = 4.0