Skielex / thinqpbo

A thin QPBO wrapper for Python.
GNU General Public License v3.0
4 stars 0 forks source link

qpbo as qubo-solver? #6

Open RobG-LHind opened 2 weeks ago

RobG-LHind commented 2 weeks ago

Thank you very much for your library.

I try to use it in order to solve qubo formulations ( x^T Q x = min) Sum (i,j, x[i]x[j]qubo[i][j]) qubo being a symmetric matrix, I thought this should do the trick:

import thinqpbo as tq
# import qubogen

qubo = [[-375,  600, -125, -100],
       [ 600, -384, -120,  -96],
       [-125, -120,  225,   20],
       [-100,  -96,   20,  176]]
# qubo = qubogen.qubo_number_partition([25, 24, -5, -4])

graph = tq.QPBOInt() 

nodes_to_add = len(qubo)
first_node_id = graph.add_node(nodes_to_add)

for i in range(0,len(qubo)):
    graph.add_unary_term(i, 0, qubo[i][i])
    for j in range(i + 1, len(qubo)):
        graph.add_pairwise_term(i, j, 0, 0, 0, 2 * qubo[i][j]) # symmetry

graph.solve()
graph.compute_weak_persistencies()
twice_energy = graph.compute_twice_energy()

solution = [graph.get_label(i) for i in range(len(qubo))]
print(f"solution={solution}, energy={twice_energy/2}")

the result is solution=[-1, -1, -1, -1], energy=0.0 should be solution=[1, 0, 1, 0], ... or solution=[0, 1, 0, 1], ...

what am i overlooking?

Skielex commented 2 weeks ago

Honestly, I don't know what a QUBO formulation is, so I can't help you. I think there might be a misunderstanding here. This is a wrapper around an algorithm for solving quadratic pseudo-Boolean optimization problems.

RobG-LHind commented 2 weeks ago

Hello Skielex, thank you for your fast answer. To me it seems that QPBO and QUBO are very similar (esp. it should be possible to map every QUBO to QPBO).

A QUBO is a "quadratic unconstrained binary optimization" H(x) = x^T Q x with no additional terms / constraints with x being a (column-) vector of binary elements b[i] ∈ {0,1}, x^T meaning transposed (row-vector) and Q being a symmetric matrix Q[i][j] ∈ ℝ or Q[i][j] ∈ ℤ} H(x) = Sum_ij( x[i] * Q[i][j] * x[j] ) == Sum_ij (x_i * Q_ij * x_j) == Sum_ij ( Q_ij * x_i * x_j )

so for me it looked like the same Problem as QPBO F(x) = w_0 + Sum_p( w_p(x_p) ) + Sum_pq( w_pq(x_p, x_q) )

if all w_p(0) == 0 and all w_pq(0, x_q) == w_pq(x_p, 0) == 0 this should be the same as F*(x) = w_0 + Sum_p( w_p(1) * x_p ) + Sum_pq( w_pq(1,1) * x_p * x_q )

and with x_p == x_p x_p because x_p ∈ {0,1} ```F(x) = w_0 + Sum_pq( w_pq(1,1) x_p x_q )which now seems similar to H(x) = Sum_ij( Q_ij x_i x_j )```

(After re-reading wikipedia, I notice, that my problem formulations do not guarantee submodularity)

I tried to map H(x) to F*(x) via:

graph = tq.QPBOInt() # can not provide info abour edges
graph.add_note(len(Q))  # Q is quadratic symmetric

for i in range(0, len(Q)):
    for j in range(0, len(Q)):
        if (Q[i][j]):
            graph.add_pairwise_term( i, j, 0, 0, 0, Q[i][j]) # Q_ij(1,1) * x_i * x_j

graph.solve()
graph.compute_weak_persistencies()
twice_energy = graph.compute_twice_energy()

solution = [graph.get_label(i) for i in range(len(qubo))] # solution vector
print(f"solution={solution}, energy={twice_energy/2}") # solution enegergy

so for each non-zero element of Q, I add a pairwise_term.

I get solutions for Q=[[-1, 2], [2, -2]] (solution [0,1], energy=-2) I do not get solutions for Q=[[-8, 6, 2], [6, -9, 3], [2, 3, -5]] (should be [1,0,1] or [0,1,0], energy=-9)

Sorry for the lengthy text. Hope you made it till here and that my point got through.

What am i overlooking?

Thanks in advance Rob

RobG-LHind commented 2 weeks ago

I tried to reduce terms as much as possible for a qubo resulting in

graph.add_pairwise_term(0,1, 24, -48, -40, 80)
graph.add_pairwise_term(0,2, 40,   0, -24,  0)
graph.add_pairwise_term(1,2, 30, -10, -42, 14)

this should result in a solution[1,0,1] or [0,1,0] with energy=-50

                                       E00  E01  E10 E11
[1,0,1] -> graph.add_pairwise_term(0,1, 24, -48, -40, 80) = -40 (E10)
[1,0,1] -> graph.add_pairwise_term(0,2, 40,   0, -24,  0) =   0 (E11)
[1,0,1] -> graph.add_pairwise_term(1,2, 30, -10, -42, 14) = -10 (E01)

do I have to configure sth or use other functions (improve, ...)?

sry for the amount of text - i write it in order to clarify what I try to do.


appendix (deduction of above terms)

Q = qubogen.qubo_number_partition([2,3,1])
array([[-8,  6,  2],
       [ 6, -9,  3],
       [ 2,  3, -5]]) 

-> H(x,y,z) = (x,y,z)^T * Q * (x,y,z)
-> H(x,y,z) = -8xx -9yy -5zz +12xy +4xz +6yz
-> H(x,y,z) =  -8x  -9y  -5z +12xy +4xz +6yz

adding a constant and multiplication with a positive scalar does not change solution. so this can be written as

H(x,y,z) = (-160x + 60)(-120y + 40)  +  (-160x + 100)(-40z + 40)  +   (-96y + 40)(-100z + 75)
pairwise(0,1, 60*40, 60 * -80, -100 * 40, -100 * -80) = pairwise(0,1, 2400, -4800, -4000 -8000)
pairwise(0,2, 100*40, 100 * 0, -60 * 40, -60 * 0)     = pairwise(0,2, 4000,     0, -2400,    0)
pairwise(1,2, 40*75, 40 * -25, -56 * 75, -56 * -25)   = pairwise(1,2, 3000, -1000, -4200, 1400)

pairwise(0,1, 2400, -4800, -4000 -8000)
pairwise(0,2, 4000,     0, -2400,    0)
pairwise(1,2, 3000, -1000, -4200, 1400)
 / 100 
pairwise(0,1, 24, -48, -40 -80)
pairwise(0,2, 40,   0, -24,  0)
pairwise(1,2, 30, -10, -42, 14)
Skielex commented 2 weeks ago

Ok, it's been a while since I worked on this and my knowledge is very limited outside of the graph-based approach to solve these problems. But I'd like to help if I can. However, I won't have time to dive back into this for another week and a half.

Skielex commented 2 weeks ago

Also, I think you need unary terms to get a solution with an energy different from zero.

Skielex commented 2 days ago

Hi @RobG-LHind,

As you mentioned yourself, since you have pairwise terms that are non-submodular, QPBO is not guaranteed to provide a complete solution. I assume this is the reason QPBO can't solve your problem.

Skielex commented 2 days ago

QPBO may reformulate your problem to allow a proper graph representation, see ComputeWeights. You can use the get_twice_pairwise_term and get_twice_unary_term methods to check the terms.