How does Qiskit's VQE encode binary digits into a solution from a given Hamiltonian?

So far I have been working with the VQE on different Hamiltonians that happened to have degeneracies, that is there were always at least two different global minima, because the Hamiltonians I was using were describing a problem with "symmetrical" solutions (for example number partitioning problems, see this paper for the corresponding Hamiltonian). So no matter if I got "1010" or "0101" as a result, I could verify that both of these results are equivalent and symmetrical just by looking at the Hamiltonian.

I was trying to implement the exact cover problem from the paper linked above, when I realized that this specific Hamiltonian does not have "symmetrical" results, because of the way it encodes the problem.

$$ H_A = Asum^n_{alpha=1}left( 1 – sum_{i: alphain V_i} x_i right)^2$$

I have written the Hamiltonian for a specific example, and I am confident it is correct, because of the consistent results I get across different examples, even when I try to solve 6-qubit problems. The problem is, I cannot understand why the solutions do not seem to be encoded so that "1" equals to a subset that is chosen to be in the solution and "0" to a subset that is left out. I will show you an example below.

U=[1,2,3,4]
V = np.empty(3,dtype=object)
V[0]=[1,3]
V[1]=[2,3]
V[2]=[2,4]

n=len(U)
N=len(V)

A=n
B=0.5

#encoding spins to binary
Id3 = (I^I^I)
x=np.empty(n, dtype=object)
x[0] = ((Z^I^I) + Id3)/2
x[1] = ((I^Z^I) + Id3)/2
x[2] = ((I^I^Z) + Id3)/2

HA = A*((Id3-x[0])**2 + (Id3-x[1]-x[2])**2 + (Id3-x[0]-x[1])**2 + (Id3-x[2])**2)

HB = B*(x[0]+x[1]+x[2])

H=HA+HB

Now, if you solve this mathematically by substituting 101, you get an HA that is equal to 0, while if you solve it for 010, you get a value of 2 for HA (I do not add HB here because it doesn’t make a difference and it was introduced for more complex problems that require the smallest exact cover available). But the VQE returns a result of 010, no matter how many times I run it, implying that it somehow encodes the solution to 0s instead of 1s. You can try it for yourself. I personally used the COBYLA optimizer, RealAmplitudes as an ansatz and zeros as an initial point.

optimizer = COBYLA(maxiter=1000)
ansatz = RealAmplitudes(3)
initial_state = np.zeros(3*4)
vqe = VQE(ansatz, optimizer, quantum_instance=backend, initial_point=initial_state)
result = vqe.compute_minimum_eigenvalue(H)
print(result)

returns

OrderedDict([   ('aux_operator_eigenvalues', None),
                ('cost_function_evals', 130),
                (   'eigenstate',
                    {   '000': 0.03125,
                        '001': 0.03125,
                        '010': 0.9965761699438734,
                        '011': 0.03125,
                        '101': 0.04419417382415922,
                        '110': 0.03125,
                        '111': 0.03125}),
                ('eigenvalue', 1.06103515625),
                (   'optimal_parameters',
                    {   ParameterVectorElement(θ[1]): -0.2905782902447393,
                        ParameterVectorElement(θ[2]): -0.09217189958775499,
                        ParameterVectorElement(θ[0]): 0.09957547756267641,
                        ParameterVectorElement(θ[3]): 0.039799150529881595,
                        ParameterVectorElement(θ[4]): 1.7026848531510135,
                        ParameterVectorElement(θ[5]): -0.2889546445425341,
                        ParameterVectorElement(θ[6]): -0.3138219479243887,
                        ParameterVectorElement(θ[7]): -0.046670861279977015,
                        ParameterVectorElement(θ[8]): -0.32390403545338103,
                        ParameterVectorElement(θ[9]): 0.3606175535199046,
                        ParameterVectorElement(θ[10]): 1.3037664385836476,
                        ParameterVectorElement(θ[11]): 0.23804467144206667}),
                (   'optimal_point',
                    array([ 0.09957548,  1.30376644,  0.23804467, -0.29057829, -0.0921719 ,
        0.03979915,  1.70268485, -0.28895464, -0.31382195, -0.04667086,
       -0.32390404,  0.36061755])),
                ('optimal_value', 1.06103515625),
                ('optimizer_evals', 130),
                ('optimizer_time', 2.3779969215393066)])

After this I tried swapping the sets’ order, just to see how the digits would behave.

U=[1,2,3,4]
V = np.empty(3,dtype=object)
V[0]=[1,3]
V[1]=[2,4]
V[2]=[2,3]

n=len(U)
N=len(V)

A=n
B=0.5

Id3 = (I^I^I)
x=np.empty(n, dtype=object)
x[0] = ((Z^I^I) + Id3)/2
x[1] = ((I^Z^I) + Id3)/2
x[2] = ((I^I^Z) + Id3)/2

HA = A*((Id3-x[0])**2 + (Id3-x[1]-x[2])**2 + (Id3-x[0]-x[2])**2 + (Id3-x[1])**2)
HB = B*(x[0]+x[1]+x[2])

H=HA+HB

backend = QasmSimulator()

optimizer = COBYLA(maxiter=1000)
ansatz = RealAmplitudes(3)
initial_state = np.zeros(3*4)
vqe = VQE(ansatz, optimizer, quantum_instance=backend, initial_point=initial_state)
result = vqe.compute_minimum_eigenvalue(H)
print(result)

groundstates = probability_calculator(result)
bootstraping_point = result.optimal_point.tolist()
plot_histogram(groundstates, figsize=(15, 5))

Confusingly enough, this time I get the results that I would expect. The VQE returns 110 as the solution, which is the result someone would expect after doing the math, with 1s meaning that the subset is chosen to be in the solution.

OrderedDict([   ('aux_operator_eigenvalues', None),
            ('cost_function_evals', 120),
            (   'eigenstate',
                {   '000': 0.03125,
                    '010': 0.06987712429686843,
                    '011': 0.03125,
                    '101': 0.09375,
                    '110': 0.9921567416492215}),
            ('eigenvalue', 8.505859374999998),
            (   'optimal_parameters',
                {   ParameterVectorElement(θ[7]): 0.9180159629622513,
                    ParameterVectorElement(θ[3]): -0.2952549861598156,
                    ParameterVectorElement(θ[4]): 0.4780970008625806,
                    ParameterVectorElement(θ[2]): 0.8850714897234401,
                    ParameterVectorElement(θ[0]): 0.06576169924533613,
                    ParameterVectorElement(θ[6]): 0.17549300683520178,
                    ParameterVectorElement(θ[1]): 1.0728636784578363,
                    ParameterVectorElement(θ[5]): 0.7532584743057651,
                    ParameterVectorElement(θ[11]): 1.0224055378029202,
                    ParameterVectorElement(θ[9]): -0.03195532295070413,
                    ParameterVectorElement(θ[10]): 0.973327271430723,
                    ParameterVectorElement(θ[8]): 0.00735674701360727}),
            (   'optimal_point',
                array([ 0.0657617 ,  0.97332727,  1.02240554,  1.07286368,  0.88507149,
   -0.29525499,  0.478097  ,  0.75325847,  0.17549301,  0.91801596,
    0.00735675, -0.03195532])),
            ('optimal_value', 8.505859374999998),
            ('optimizer_evals', 120),
            ('optimizer_time', 2.728998899459839)])

I just want to make a note here. I use very simple 3-qubit examples to make my point across, but I’ve used the same Hamiltonian for several different 4-qubit and 6-qubit problems and I mostly get results that behave the same way as my first example here, that is returning 0s where I would expect 1s according to how the Hamiltonian is formulated.

My question is why does this happen and if it is related to how the vqe works or there is a logical error in my approach?