# Outer function errors, Mutual Information

**URL:** <https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mps, mpo\
**Created:** [June 11, 2025, 7:36pm UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422 "2025-06-11T19:36:40Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![PZuiuc](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/pzuiuc/32/778_2.png) [@PZuiuc](https://itensor.discourse.group/u/PZuiuc)\
**Post date:** [June 11, 2025, 7:36pm UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/1 "2025-06-11T19:36:40Z")

</div>

In the following code, after obtaining the ground state `psi` of the Hamiltonian `H` using DMRG, I verify the normalization by checking `inner(psi, psi)` and `trace(rho)`, where `rho = outer(dag(psi), psi)`. I find that `inner(psi, psi)` equals 1 up to an error on the order of 10⁻¹⁴, while `trace(rho)` deviates from 1 by about 10⁻⁶. I would like to understand the source of this discrepancy. Specifically, does the `outer` function automatically introduce truncation or compression that could explain the larger error in `trace(rho)`? Or might there be a mistake in how I’m computing the trace of `rho`? Thanks!

```
Nuc=60
sites = siteinds("S=1/2",Nuc)

J=1.0
os = OpSum()
for j=1:Nuc-1
      if j==Nuc
          jp=1
      else
          jp=j+1
      end
    os += -4J,"Sz",j,"Sz",jp
    os += -2,"Sx",j
 end
 os += -2,"Sx",Nuc

 H = MPO(os,sites)

 psi0 = productMPS(sites, n -> isodd(n) ? "Up" : "Dn")

 nsweeps = 15
 maxdim = [10,20,100,100,200,200,300,300,400,500]
 noise=[1E-3,1E-5,1E-7,1E-7,1E-9,1E-9,1E-11,1E-11,1E-13,1E-13,1E-14,1E-15]
 cutoff = [1E-12]

 energy0, psi = dmrg(H,psi0; nsweeps, maxdim, noise, cutoff)
 orthogonalize!(psi, 1)
 rho0=outer(psi',psi)

 println(inner(psi,psi))

 aaa=ITensor(1.0)
 for i=1:Nuc
     rho0[i]=rho0[i]*delta(sites[i],sites[i]')
     aaa*=rho0[i]
 end
 println(aaa)

```

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 12, 2025, 1:53am UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/2 "2025-06-12T01:53:02Z")

</div>

I’m not 100% sure, but your guess about the role of `outer` here is likely what’s going on. When I ran your code I found that `outer` gives an MPO with a bond dimension that is less than the square of the MPS bond dimension than is inputted into it. Perhaps this is ok if there are a lot of small or zero singular values of `rho0` but it could be that these are being overly truncated.

In general I would not recommend the use of `outer` except for basic testing on small systems. A much more accurate and efficient way to work with pure density matrices (outer product of an MPS with itself) is to think of forming such a density matrix implicitly, then computing properties like the trace using an algorithm similar to the [algorithm for computing an inner product of two MPS](https://tensornetwork.org/mps/index.html#innerprod).

If you first right-orthogonalize the MPS, as you have done in your code, then actually all the tensors except the first one cancel out if you take a partial trace of \psi \otimes \psi^\dagger on all the sites except site 1. So then the trace of the entire density matrix reduces to just the inner product of the first MPS tensor with itself.

---

<div class="post-metadata">

**Author:** ![PZuiuc](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/pzuiuc/32/778_2.png) [@PZuiuc](https://itensor.discourse.group/u/PZuiuc)\
**Post date:** [June 12, 2025, 2:17am UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/3 "2025-06-12T02:17:41Z")

</div>

Hi Miles,  
Thank you very much for your swift response. Since my goal is to compute the mutual information for a model more complex than the Ising model, I think constructing the full density matrix is indeed necessary. In this case, do you know of any reliable method within ITensor to build the density matrix accurately?

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 12, 2025, 2:34am UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/4 "2025-06-12T02:34:14Z")

</div>

There isn’t a function that you can simply call to obtain the density matrix like this. One reason is that it just doesn’t have an efficient form as an MPS or an MPO. It’s only efficient form is as two MPS in a (“lazy” or formal or implicit) outer product with each other.

What sort of regions do you want to compute the mutual information between? There are efficient ways to compute that if the two regions are sufficiently small or have other structure like if they reach to the edges of the system.

There are also some previous discussions on here asking about mutual information though the answers may not be totally satisfactory.

Finally there is an experimental package ITensorEntropyTools which may help you.  
\< [GitHub - ryanlevy/ITensorEntropyTools.jl: Tools for calculating Entanglement Entropy w/ ITensors.jl](https://github.com/ryanlevy/ITensorEntropyTools.jl)\>

---

<div class="post-metadata">

**Author:** ![PZuiuc](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/pzuiuc/32/778_2.png) [@PZuiuc](https://itensor.discourse.group/u/PZuiuc)\
**Post date:** [June 12, 2025, 3:18am UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/5 "2025-06-12T03:18:10Z")

</div>

Thanks a lot for sharing this resource. The mutual information I want to calculate is between two spins inside a chain, which are far away from each other. I will try to look at the ITensorEntropyTools to see if I can find solution.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 13, 2025, 3:03pm UTC](https://itensor.discourse.group/t/outer-function-errors-mutual-information/2422/6 "2025-06-13T15:03:18Z")

</div>

Yes, the mutual information between two individual sites is something that is efficiently computable with MPS. But the best approach is a more “custom” or “targeted” method for computing the reduced density matrices involved, rather than ever forming the entire (pure) density matrix of the whole system as a single tensor network like an MPO.

Here is what I mean. The diagrams shown below depict how the three reduced density matrices involved in the mutual information can be written as efficient tensor network contractions (efficient because you can contract all of these from the ‘side’ i.e. from left to right, alternating between the top and bottom layers).

 ![Screenshot 2025-06-13 at 10.56.22 AM](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/0e4ba63702266371982e3bbc77258ecac001ae40.jpeg)  
 ![Screenshot 2025-06-13 at 10.56.26 AM](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/c6811b637ac67fc74fb1bd2cacca8aa2de127b4f.png)

As a possible optimization, if the MPS is “gauged” so that its orthogonality center is at i for the first diagram, j for the second diagram, and somewhere in the range [i,j] for the third diagram, then the outer MPS tensors ‘cancel’ and can be omitted from the computations. Then for \rho\_i for example, only the i th MPS tensor is needed for the computation.

So some pseudocode (i.e. to be tested and checked) for the first RDM is:

```auto
psi = orthogonalize(psi,i)
rhoi = prime(psi[i],"Site") * dag(psi[i])

```

Finally, once each RDM is computed, they can all be diagonalized using `eigen` to find their spectra.