The QUBO function cannot yield the correct solution to my problem.

I want use QUBO function (https://ww2.mathworks.cn/help/matlab/math/what-is-a-qubo.html) to solve my question in the following:

min Z=2*x12 + 2*x23 + 10*x25 + x34 + 3*x45 + 3*x46 + 2*x54 + x56 + 2*x67 + 2*x12*x23 + 3*x12*x25 + 2*x23*x34 + x25*x54 + 2*x34*x45 + x34*x46 + x25*x56 + 5*x46*x54 + 10*x45*x56 + 20*x46*x67 + x56*x67 ;

s.t. x12 - 1=0;

x23 - x12 + x25 =0;

x23 - x34=0;

x25 + x45 - x54 - x56 =0;

x34 - x45 - x46 + x54 =0;

x46 + x56 - x67=0;

x67 - 1 =0;

all variables: x12,x23, x25,x34,x45,x46,x54,x56,x67 are binary variables

I already know that the optimal solution to this problem is ：

x12=x25=x56=x67=1, other variables =0 the optimal value =20

however when I try to use QUBO function to solve it, I converted all the constraints into penalty factors and added them to the objective function, thereby transforming the problem into an unconstrained one， the code is：

clc
clear
syms x12 x23 x25 x34 x45 x46 x54 x56 x67
P=50; %Penalty coefficient
Z=2*x12+2*x23+10*x25+1*x34+3*x45+ 3*x46+2*x54+1*x56+2*x67 + ...
    2*x12*x23 + 3*x12*x25 + 2*x23*x34 + x25*x54 + x25*x56 +2*x34*x45+ ...
    x34*x46+10*x45*x56+ 20*x46*x67 +5*x54*x46 +x56*x67 ...
    + P*(x12-1)^2 +P*(x12-x23-x25)^2 +P*(x23-x34)^2 ...
    +P*(x45+x25-x54-x56)^2 + P*(x54+x34-x45-x46)^2 + P*(x46+x56-x67)^2 ...
    + P*(x67-1)^2
Q = hessian(Z, [x12 x23 x25 x34 x45 x46 x54 x56 x67])
Q=double(Q)
d = jacobian(Z, [x12 x23 x25 x34 x45 x46 x54 x56 x67]);
d = subs(d, [x12 x23 x25 x34 x45 x46 x54 x56 x67], [0 0 0 0 0 0 0 0 0]);
d = double(d);
qprob = qubo(Q,d)
result = solve(qprob)

in the code, the Hessian matrix and the coefficients of the linear terms of the objective function are determined using the relevant functions. the result is:

It can be seen that the solution obtained by the QUBO function has all 0-1 variable values as 0, and the objective function value is also 0. This result is clearly wrong. I would like to know why caused this result.

[Chuguang Pan]

@Du Wei. Try increase the Penalty coefficient.

[Du Wei]

It doesn't work，Even if I set the value of P to 5000, I still get the same incorrect result as before.

Answer 1

[Stephen23]

The issue is that MATLAB’s hessian returns twice the quadratic coefficient matrix required by qubo.

For the penalized objective, write

The Hessian is

when Q is symmetric. Therefore, this line:

Q = hessian(Z,vars);

passes 2Q, not Q, to the QUBO problem.

Consequently, MATLAB solves

instead of

The quadratic and linear terms are no longer balanced correctly. Increasing P does not repair this, because the same factor-of-two error remains.

For P = 50, the incorrectly constructed QUBO gives:

• all-zero vector: QUBO value 0
• expected vector [1;0;1;0;0;0;0;1;1]: QUBO value 25

So the reported all-zero solution is consistent with the incorrectly supplied matrix.

Corrected code (untested)

Divide the Hessian by two:

syms x12 x23 x25 x34 x45 x46 x54 x56 x67
vars = [x12 x23 x25 x34 x45 x46 x54 x56 x67];
P = 50;
Z0 = 2*x12 + 2*x23 + 10*x25 + x34 + 3*x45 + 3*x46 + 2*x54 + x56 + 2*x67 + 2*x12*x23 + 3*x12*x25 + 2*x23*x34 + x25*x54 + 2*x34*x45 + x34*x46 + x25*x56 + 5*x46*x54 + 10*x45*x56 + 20*x46*x67 + x56*x67;
penalty = ...
    (x12 - 1)^2 + ...
    (x23 - x12 + x25)^2 + ...
    (x23 - x34)^2 + ...
    (x25 + x45 - x54 - x56)^2 + ...
    (x34 - x45 - x46 + x54)^2 + ...
    (x46 + x56 - x67)^2 + ...
    (x67 - 1)^2;
Z = expand(Z0 + P*penalty);
% For Z = x.'*Q*x + d.'*x + constant:
Q = double(hessian(Z,vars))/2;
d = jacobian(Z,vars);
d = subs(d,vars,zeros(1,numel(vars)));
d = double(d(:));
qprob = qubo(Q,d);
result = solve(qprob);
result.BestX
result.BestFunctionValue

[Du Wei]

Thank you very much for your help！I did indeed overlook the fact that calculating Q using the Hessian function requires dividing by 2; the code you provided successfully yields the optimal solution，but why is the optimal value of the objective function obtained -80, when the correct optimal value is 20?

[Stephen23]

Either

x = result.BestX;
val = ...
    2*x(1)+2*x(2)+10*x(3)+x(4)+3*x(5)+3*x(6)+2*x(7)+x(8)+2*x(9) + ...
    2*x(1)*x(2)+3*x(1)*x(3)+2*x(2)*x(4)+x(3)*x(7)+ ...
    2*x(4)*x(5)+x(4)*x(6)+x(3)*x(8)+5*x(6)*x(7)+ ...
    10*x(5)*x(8)+20*x(6)*x(9)+x(8)*x(9);

or equivalently:

val = x.'*Q*x + d.'*x + 100;

[Du Wei]

Your approach essentially involves finding the optimal solution and then substituting it back into the objective function; what confuses me is why the objective function value obtained via QUBO is -80？

[Torsten]

The objective function value of QUBO is

syms x12 x23 x25 x34 x45 x46 x54 x56 x67
vars = [x12 x23 x25 x34 x45 x46 x54 x56 x67];
P = 50;
Z0 = 2*x12 + 2*x23 + 10*x25 + x34 + 3*x45 + 3*x46 + 2*x54 + x56 + 2*x67 + 2*x12*x23 + 3*x12*x25 + 2*x23*x34 + x25*x54 + 2*x34*x45 + x34*x46 + x25*x56 + 5*x46*x54 + 10*x45*x56 + 20*x46*x67 + x56*x67;
penalty = ...
    (x12 - 1)^2 + ...
    (x23 - x12 + x25)^2 + ...
    (x23 - x34)^2 + ...
    (x25 + x45 - x54 - x56)^2 + ...
    (x34 - x45 - x46 + x54)^2 + ...
    (x46 + x56 - x67)^2 + ...
    (x67 - 1)^2;
Z = expand(Z0 + P*penalty);
% For Z = x.'*Q*x + d.'*x + constant:
Q = double(hessian(Z,vars))/2;
d = jacobian(Z,vars);
d = subs(d,vars,zeros(1,numel(vars)));
d = double(d(:));
x = vars.';
Zopt = subs(x.'*Q*x + d.'*x,x,[1 0 1 0 0 0 0 1 1].')

Zopt = 

Your objective function value is

Zopt_with_constant_term_1 = subs(Z,x,[1 0 1 0 0 0 0 1 1].')

Zopt_with_constant_term_1 = 

20

or

Zopt_with_constant_term_2 = subs(Z0,x,[1 0 1 0 0 0 0 1 1].')

Zopt_with_constant_term_2 = 

20

When expanding Z, the constant term "constant" does not influence the optimal x and is not taken into account for QUBO.

[Paul]

The constant term can be included in the defintion of the qubo problem by using the third input argument. The constant term doesn't affect the optimal solution, only the function value. Also, from a notation standpoint, the OP is using the variable d for the linear part though the doc uses "d" for the constant (and "c" for the linear).

[Stephen23]

@Du Wei: Torsten and Paul identified the cause and the solution. To retain the constant, use the third input argument to qubo. It is also clearer to use c for the linear coefficients and d for the constant, consistent with the qubo documentation (untested):

Z = expand(Z0 + P*penalty);
Q = double(hessian(Z,vars))/2;
c = jacobian(Z,vars);
c = double(subs(c,vars,zeros(size(vars))));
c = c(:);
d = double(subs(Z,vars,zeros(size(vars))));
qprob = qubo(Q,c,d);
result = solve(qprob);
result.BestX
result.BestFunctionValue

[Du Wei]

That is indeed the crux of the issue; I am grateful to every expert who helped me resolve it.