Returns the matrix diagonal of the product tensor

def productdiag(tensor: BKTensor) -> BKTensor:
    """Returns the matrix diagonal of the product tensor"""  # DOCME: Explain
    N = rank(tensor)
    tensor = reshape(tensor, [2**(N//2), 2**(N//2)])
    tensor = np.diag(tensor)
    tensor = reshape(tensor, [2]*(N//2))
    return tensor