Mutual Information matrix

DMRG and Numerical Methods
1

Hi, everyone!
I need to compute a mutual information matrix of a MPS:

I(A,B) = S(A) + S(B) - S(A,B)

where

S(A,B) = -tr( \rho_{AB} log(\rho_{AB}))\\ S(A) = -tr( \rho_{A} log(\rho_{A}))\\ S(B) = -tr( \rho_{B} log(\rho_{B}))

Following this topic, I’ve computed the reduced density matrix for each pair of sites.

Then, I tried to implement this function to compute each entropy term

function MI_matrix(rdm)
    A,B = size(rdm);
    I_AB = zeros(A,B);
    for a in 1:A, b in a+1:B
        rho_AB = copy(rdm[a,b]);
        delta_AA = delta(dag(inds(rho_AB)[1]), dag(inds(rho_AB)[2]));
        delta_BB = delta(dag(inds(rho_AB)[3]), dag(inds(rho_AB)[4]));

        rho_B = delta_AA*rho_AB;
        rho_A = delta_BB*rho_AB;


        S_AB = -tr(reshape(Array(rho_AB,inds(rho_AB)[1], inds(rho_AB)[3], inds(rho_AB)[2], inds(rho_AB)[4]),16,16) 
                * log(reshape(Array(rho_AB,inds(rho_AB)[1], inds(rho_AB)[3], inds(rho_AB)[2], inds(rho_AB)[4]),16,16)));
        S_A = -tr(Array(rho_B, inds(rho_B)[1], inds(rho_B)[2]) * log(Array(rho_B, inds(rho_B)[1], inds(rho_B)[2])));
        S_B = -tr(Array(rho_A, inds(rho_A)[1], inds(rho_A)[2]) * log(Array(rho_A, inds(rho_A)[1], inds(rho_A)[2])));
        
        I_AB[a,b] = I_AB[b,a] = S_A + S_B - S_AB;
    end
    return I_AB
end

rdm is the matrix cotaining all the 2-sites reduced density matrices.

For the part of S(A,B) , it is correct to compute the value in this way?
First, i need to convert each rank-4 tensor \rho_{AB} into a matrix (16x16 because each site has dimension 4) and then I apply the function tr.
For \rho_{A} and \rho_{B} , I’ve used a delta tensor.

For the single site entropy, I obtain something then seems right, but for S(A,B) when I compute log\rho_{AB } the result is NaN.

Do you think that the procedure is correct?

1 Like
2

The main change I would recommend is to diagonalize the density matrices, obtaining their eigenvalues, rather than trying to form \rho \log(\rho). (In fact, the log of a matrix is defined as the log of its eigenvalues. I’m not sure what calling Julia’s log function on an Array might do but I think it’s unlikely to be the correct mathematical operation.)

To diagonalize a positive-definite, Hermitian ITensor, such as a density matrix, you can use the function eigen provided by ITensor. Here is the documentation for it:
https://itensor.github.io/ITensors.jl/dev/ITensorType.html#LinearAlgebra.eigen-Tuple{ITensor,%20Any,%20Any}
which also shows a small example.

A good thing about eigen is that it doesn’t require you to use delta tensors or anything to prepare the density matrix ITensor or turn it into a matrix. You just need to pass in arrays of the indices that you want to be considered collectively as the “row” indices and the “column” indices.

From the diagonal elements of the eigenvalue tensor that is returned, you can do something like:

S = 0.0
for n=1:dim(D, 1)
  p = D[n,n]
  S -= p * log(p)
end

to compute the entropy S.

2 Likes
3

I have also been computing the mutual information for all pairs of sites. I’m curious what is the best way to calculate the RDMs. I have been using the following formulas

\rho^A_{i j} = \braket{\psi | P^A_{i j} | \psi}

and

\rho^{A, B}_{i j k \ell} = \braket{\psi | P^A_{i j} P^B_{k \ell} | \psi}

where P^A_{i j} is an operator that acts on site A and is all zeros except for the entry at (i, j). This can be done fairly cleanly using expect and correlation_matrix, but more importantly if \psi has some conserved quantities then many of the entries in the RDM are known to be zero. For example if \psi is composed of electron sites that conserve particle number and spin only the diagonal entries of \rho^A_{i j} need to be calculated.

I did kind of just hack this together so, so I’d love to hear about a better approach! (I also think it would be a nice addition to ITensor to have a function to calculate RDMs).

1 Like
4

Thank you very much for helping!
Sorry if I didn’t reply sooner, but I’ve had a busy few weeks and only recently have I finally tested the change and it works very well!

I’ve also added a check on the eigenvalues:

#=
Diagonalization of the two-sites reduced density matrix
=#
D_AB, _ = eigen(rho_AB, L_inds, R_inds); 

S_AB = 0.0;
for n=1:dim(D_AB, 1)
    p = D_AB[n,n]
#Added a check if the value is 0.0, I'll skip the calculation
    if p == 0.0
        nothing
    else
        S_AB -= p * log(p)
    end
end

Best,
Fabio.

5

I am new to Julia and ITensors. Could you please tell how do you determine L_inds and R_inds?