From 41e86b020fe2c54c1f1fcd98ae2aea40230af6cd Mon Sep 17 00:00:00 2001
From: stecrotti <sky_96@live.it>
Date: Sat, 18 Oct 2025 17:35:46 +0200
Subject: [PATCH 1/3] Another fix for complex derivative

---
 Project.toml                           |  2 +-
 src/MatrixProductStates/derivatives.jl |  8 ++++----
 src/tensor_train.jl                    | 10 +++++-----
 test/mps.jl                            | 14 +++++---------
 4 files changed, 15 insertions(+), 19 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0e8be88..a344513 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TensorTrains"
 uuid = "89893e69-996d-40b1-ba32-8ff5f34c0dd5"
-authors = ["stecrotti <stefano.crotti@polito.it>", "abraunst <alfredo.braunstein@polito.it"]
 version = "0.13.0"
+authors = ["stecrotti <stefano.crotti@polito.it>", "abraunst <alfredo.braunstein@polito.it"]
 
 [deps]
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
diff --git a/src/MatrixProductStates/derivatives.jl b/src/MatrixProductStates/derivatives.jl
index acba284..6590264 100644
--- a/src/MatrixProductStates/derivatives.jl
+++ b/src/MatrixProductStates/derivatives.jl
@@ -12,7 +12,7 @@ function grad_normalization_canonical(p::MPS, k::Integer)
     @tullio zz = Aᵏconj_[m,n,x] * Aᵏ_[m,n,x]
     z2 = abs2(float(p.ψ.z))
     z = zz / z2
-    gradz = 2 * Aᵏ / z2
+    gradz = 2 * conj(Aᵏ) / z2
     return gradz, z
 end
 
@@ -24,7 +24,7 @@ function grad_loglikelihood(p::MPS, k::Integer, X)
     gA = - Zprime ./ Z
     for x in X
         gr, val = grad_evaluate(p.ψ, k, x)
-        gA[:,:,x[k]...] .+= 2/T * gr / val
+        gA[:,:,x[k]...] .+= 2 * conj(gr) / (T * val)
         ll += 1/T * log(abs2(val))
     end
     return gA, ll
@@ -70,7 +70,7 @@ function grad_loglikelihood_two_site(p::MPS, k::Integer, X;
     weights = ones(length(X)))
 
     Zprime, Z = grad_normalization_two_site_canonical(p, k; Aᵏᵏ⁺¹)
-    ll = -log(Z) * mean(weights)
+    ll = - log(Z) * mean(weights)
     T = length(X)
     gA = - Zprime / Z * mean(weights)
 
@@ -79,7 +79,7 @@ function grad_loglikelihood_two_site(p::MPS, k::Integer, X;
         gr, val = grad_evaluate_two_site(p.ψ, k, x;
             Ax_left = prodA_left[n][k-1], Ax_right = prodA_right[n][k+2], Aᵏᵏ⁺¹
             )
-        @inbounds gA[:,:,x[k]...,x[k+1]...] .+= 2/T * gr / val * weights[n]
+        @inbounds gA[:,:,x[k]...,x[k+1]...] .+= 2 * conj(gr) / (T * val) * weights[n]
         ll += 1/T * log(abs2(val)) * weights[n]
     end
     return gA, ll
diff --git a/src/tensor_train.jl b/src/tensor_train.jl
index 93d7b03..55c92f9 100644
--- a/src/tensor_train.jl
+++ b/src/tensor_train.jl
@@ -154,19 +154,19 @@ function orthogonalize_left!(C::TensorTrain{F}; svd_trunc=TruncThresh(1e-6),
     return C
 end
 
-function orthogonalize_center!(C::TensorTrain, l::Integer; svd_trunc=TruncThresh(1e-6))
+function orthogonalize_center!(C::TensorTrain, l::Integer; svd_trunc=TruncThresh(0.0))
     orthogonalize_left!(C; svd_trunc, indices = 1:l-1)
     orthogonalize_right!(C; svd_trunc, indices = l+1:length(C))
 end
 
 """
-    orthogonalize_two_site_center!(C::TensorTrain, k::Integer; svd_trunc=TruncThresh(1e-6))
+    orthogonalize_two_site_center!(C::TensorTrain, k::Integer; svd_trunc=TruncThresh(0.0))
 
 Orthogonalize the tensor train for a two-site DMRG update at positions k and k+1.
 This puts sites 1:k-1 in left-canonical form, sites k+2:N in right-canonical form,
 and leaves sites k and k+1 as the non-orthogonal center for merging.
 """
-function orthogonalize_two_site_center!(C::TensorTrain, k::Integer; svd_trunc=TruncThresh(1e-6))
+function orthogonalize_two_site_center!(C::TensorTrain, k::Integer; svd_trunc=TruncThresh(0.0))
     @assert 1 <= k < length(C) "k must be between 1 and length(C)-1 for two-site update"
     orthogonalize_left!(C; svd_trunc, indices = 1:k-1)
     orthogonalize_right!(C; svd_trunc, indices = k+2:length(C))
@@ -246,7 +246,7 @@ function grad_evaluate(A::TensorTrain, l::Integer, X)
     Ax_center = A[l][:,:,X[l]...]
     z = float(A.z)
     val = only(prodA_left * Ax_center * prodA_right) / z
-    gr = (prodA_right * prodA_left)' / z
+    gr = transpose(prodA_right * prodA_left) / z
     return gr, val
 end
 
@@ -281,7 +281,7 @@ function grad_evaluate_two_site(A::TensorTrain, k::Integer, X;
     Ax_center = Aᵏᵏ⁺¹[:,:,X[k]...,X[k+1]...]
     z = float(A.z)
     val = only(Ax_left * Ax_center * Ax_right) / z
-    gr = (Ax_right * Ax_left)' / z
+    gr = transpose(Ax_right * Ax_left) / z
     return gr, val
 end
 
diff --git a/test/mps.jl b/test/mps.jl
index d4e0638..859f5ca 100644
--- a/test/mps.jl
+++ b/test/mps.jl
@@ -247,23 +247,19 @@ end
     @testset "Complex Derivatives" begin
         F = ComplexF64
         tensors = [rand(F, 1,5,2,2), rand(F, 5,4,2,2),
-            rand(F, 4,10,2,2), rand(F, 10,1,2,2)]
+            rand(F, 4,10,2,2), rand(F, 10,3,2,2), rand(F, 3, 1, 2, 2)]
         ψ = TensorTrain(tensors)
         p = MPS(ψ)
 
         @testset "Gradient of loglikelihood - 2-site" begin
-            X = [sample(p)[1] for _ in 1:10]
+            X = [sample(rng, p)[1] for _ in 1:10]
             for l in 1:length(p)-1
-                orthogonalize_two_site_center!(p, l)
                 p_cp = deepcopy(p)
-                A = _merge_tensors(p_cp[l], p_cp[l+1])
-                dlldA, ll = grad_loglikelihood_two_site(p_cp, l, X)
-                @test ll ≈ loglikelihood(p_cp, X)
-                η = 1e-3
-                lls = map(1:100) do _
+                orthogonalize_two_site_center!(p_cp, l)
+                η = 1e-4
+                lls = map(1:100) do it
                     A = _merge_tensors(p_cp[l], p_cp[l+1])
                     dlldA, ll = grad_loglikelihood_two_site(p_cp, l, X)
-                    @test ll ≈ loglikelihood(p_cp, X)
                     p_cp[l], p_cp[l+1] = TensorTrains._split_tensor(A + η*dlldA)
                     ll
                 end

From 631f9996d849084d80aa016338d50146d8e72db3 Mon Sep 17 00:00:00 2001
From: stecrotti <sky_96@live.it>
Date: Sat, 18 Oct 2025 19:57:32 +0200
Subject: [PATCH 2/3] maybe found the error

---
 experiments/dmrg.jl                    | 4 ++--
 src/MatrixProductStates/derivatives.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/experiments/dmrg.jl b/experiments/dmrg.jl
index 3d27cdb..b5859d2 100644
--- a/experiments/dmrg.jl
+++ b/experiments/dmrg.jl
@@ -35,9 +35,9 @@ function CB()
 end
 
 callback = CB()
-nsweeps = 40
+nsweeps = 20
 two_site_dmrg!(p, X, nsweeps;
-    η=1e-4, ndesc=10, svd_trunc=TruncBond(d_original+2), callback)
+    η=1e-3, ndesc=10, svd_trunc=TruncBond(d_original), callback)
 
 pl1 = plot(callback.nlls, xlabel="it", ylabel="NLL", label="")
 hline!(pl1, [nll], ls=:dash, c=:gray, label="NLL according to generative model")
diff --git a/src/MatrixProductStates/derivatives.jl b/src/MatrixProductStates/derivatives.jl
index 6590264..149a5c3 100644
--- a/src/MatrixProductStates/derivatives.jl
+++ b/src/MatrixProductStates/derivatives.jl
@@ -79,7 +79,7 @@ function grad_loglikelihood_two_site(p::MPS, k::Integer, X;
         gr, val = grad_evaluate_two_site(p.ψ, k, x;
             Ax_left = prodA_left[n][k-1], Ax_right = prodA_right[n][k+2], Aᵏᵏ⁺¹
             )
-        @inbounds gA[:,:,x[k]...,x[k+1]...] .+= 2 * conj(gr) / (T * val) * weights[n]
+        @inbounds gA[:,:,x[k]...,x[k+1]...] .+= 2 * conj(gr) / (T * conj(val)) * weights[n]
         ll += 1/T * log(abs2(val)) * weights[n]
     end
     return gA, ll

From 192dbb02f5690bb2f839374c4301bf48d7a7b66c Mon Sep 17 00:00:00 2001
From: stecrotti <sky_96@live.it>
Date: Sat, 18 Oct 2025 20:15:51 +0200
Subject: [PATCH 3/3] test with experiment

---
 experiments/dmrg.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/experiments/dmrg.jl b/experiments/dmrg.jl
index b5859d2..13cc9e4 100644
--- a/experiments/dmrg.jl
+++ b/experiments/dmrg.jl
@@ -3,15 +3,16 @@ using Plots
 
 d_original = 3
 L = 5
-q = rand_mps(ComplexF64, d_original, L, 4)
+F = ComplexF64
+q = rand_mps(F, d_original, L, 4)
 normalize!(q)
-nsamples = 5*10^3
+nsamples = 10*10^3
 X = [sample(q)[1] for _ in 1:nsamples]
 nll = -loglikelihood(q, X)
 println("Negative Log-Likelihood according to generating distribution q = $nll\n")
 mq = marginals(q)
 
-p = rand_mps(ComplexF64, 2, length(q), 4)
+p = rand_mps(F, 2, length(q), 4)
 
 function CB()
     nlls = zeros(0)
@@ -37,7 +38,7 @@ end
 callback = CB()
 nsweeps = 20
 two_site_dmrg!(p, X, nsweeps;
-    η=1e-3, ndesc=10, svd_trunc=TruncBond(d_original), callback)
+    η=1e-4, ndesc=10, svd_trunc=TruncBond(d_original), callback)
 
 pl1 = plot(callback.nlls, xlabel="it", ylabel="NLL", label="")
 hline!(pl1, [nll], ls=:dash, c=:gray, label="NLL according to generative model")