QISKIT: ValueError: too many subscripts in einsum DensityMatrix()

Question

I am trying to compute the entanglement entropy of a partition of a quantum system on qiskit. To do this, I call the function DensityMatrix(). If I go above 10 sites (e.g. 12), I get an error like:

ValueError: too many subscripts in einsum

Is there any workaround for this problem?

Here is the code:

#FIXED THETA STEP, BIPARTITION, varying p of theta/p


#Circuit parameters
n=8
steps = np.array(([1]))

#SPECIFY PARAMETERS
A=1;G=0;

#Define arrays
thetas= np.arange(0.001,10,0.1)    
ent_entropy = np.zeros((len(steps),len(thetas)),dtype=complex)
renyi = np.zeros((len(steps),len(thetas)),dtype=complex)
subsystem = [j for j in range(0,n//2)]

for m,step in enumerate(steps):
    
    for k,T in enumerate(thetas):     
        string = "0" * (n)
        qc = QuantumCircuit(n)
        QuantumCircuit.initialize(qc,string)
        
        #init with Kareljan's state
        qc.unitary(prodstaten(n), range(0,n))
        
        matrix = U(A,T/step,G)
        bwassembler(qc,matrix,1,n)
        
        
    
        
        ######################################################################    
        #access density matrix  
        rho = DensityMatrix(qc)
        reduced_rho=partial_trace(rho, subsystem)
        
        #find the eigenvalues of the reduced density matrix
        eigenvalues, _ = np.linalg.eig(reduced_rho)
        
        #compute the entanglement entropy
        entanglement_entropy = -np.sum(eigenvalues * np.log(eigenvalues))
        renyi[m,k]=-np.log(DensityMatrix.purity(reduced_rho))
        ent_entropy[m,k]=entanglement_entropy
        #print(k,ent_entropy[k])
        ######################################################################

Answer 1

Try to transpile your circuit into 1- and 2-qubit gates:

qc = transpile(qc, basis_gates=['cx', 'u'])
rho = DensityMatrix(qc)

This issue happens because Numpy einsum is used internally in DensityMatrix class for generalised matrix-matrix multiplication on higher-order tensors. When the circuit contains 11 qubits, 33 indices would be used to represent the tensor-reduction summation. Numpy has an upper limit on the number of indices it can use in an array, which is numpy.MAXDIM and is 32 by default.

Similar issue is reported here, where the above root cause is given.