diff --git a/calculations/PMM.jl b/calculations/PMM.jl
index eee0769..49b1ba2 100644
--- a/calculations/PMM.jl
+++ b/calculations/PMM.jl
@@ -15,26 +15,30 @@ data_E = [quick_pole_E(V_system(c)) for c in data_c]
 N = 9
 
 # initialize random Hamiltonians
-H0 = randn(ComplexF64, N, N)
-H0 = H0 + transpose(H0) # symmetric
-H1 = randn(ComplexF64, N, N)
-H1 = H1 + transpose(H1) # symmetric
+H0 = randn(N, N)
+H1 = randn(N, N)
 
 # training
 Es = ComplexF64[]
-ψs = Vector{ComplexF64}[]
+ψrs = Vector{ComplexF64}[]
+ψls = Vector{ComplexF64}[]
 
 lr = 0.05
 epochs = 100000
+
 for epoch in 1:epochs
     empty!(Es)
-    empty!(ψs)
+    empty!(ψrs)
+    empty!(ψls)
     for (c, E) in zip(data_c, data_E)
         H = H0 + c * H1
-        evals, evecs = eigen(H)
-        i = nearestIndex(evals, E) # TODO: more robust way to identify the eigenvector
-        push!(Es, evals[i])
-        push!(ψs, evecs[:, i])
+        r_evals, r_evecs = eigen(H)
+        l_evals, l_evecs = eigen(transpose(H))
+        @assert all(r_evals .≈ l_evals) "Right/left eigenvalues do not match"
+        i = nearestIndex(r_evals, E) # TODO: more robust way to identify the eigenvector
+        push!(Es, r_evals[i])
+        push!(ψrs, r_evecs[:, i])
+        push!(ψls, l_evecs[:, i])
     end
 
     if epoch % 1000 == 0
@@ -45,8 +49,8 @@ for epoch in 1:epochs
     # gradient of the loss function
     function grad(c_order=0)
         out = zeros(ComplexF64, N, N)
-        for (c, E_target, ψ, E) in zip(data_c, data_E, ψs, Es)
-            out .+= (c^c_order * conj(E - E_target)) .* (ψ * transpose(ψ))
+        for (c, E_target, ψr, ψl, E) in zip(data_c, data_E, ψrs, ψls, Es)
+            out .+= (c^c_order * conj(E - E_target)) .* (ψl * transpose(ψr))
         end
         return 2 .* real.(out)
     end