Reproducing Set Packing Problem

[seanlaw]

I am trying to reproduce the results from a known set packing problem by modeling it as a QUBO problem. Essentially, we are trying to optimize the following objective function along with 2 constraints:

$$ \begin{align} \text{max }x_1+x_2+x_3+x_4 \\ \text{st } x_1+x_3+x_4 \le 1 \\ x_1+x_2\le 1 \end{align} $$

The constraints can be written as QUBO penalties so:

$$\text{max } x_1+x_2+x_3+x_4-Px_1x_3-Px_1x_4-Px_3x_4-Px_1x_2$$

where the scalar penalty $P$ is arbitrarily chosen to be 3. Thus, our Q matrix is:

Q = torch.tensor([[1, -3, -3, -3],
                  [-3, 1,  0,  0],
                  [-3, 0,  1, -3],
                  [-3, 0, -3,  1]
                 ])

The optimal solution is known to be either x1, x2, x3, x4 = (0, 1, 1, 0) or x1, x2, x3, x4 = (0, 1, 0, 1), both with an objective function evaluated to 2 (i.e., this is the maximum value). However, when I use SB to maximize this Q matrix:

sb.maximize(Q, input_type='binary')

I get x1, x2, x3, x4 = (1, 0, 0, 0) and the objective function evaluates to 1, which is smaller than the expected 2 above. Therefore, this SB solution is not optimal. Have I done something wrong or made a mistake? Any help would be greatly appreciated!

[BusyBeaver-42]

The constraints can be written as QUBO penalties so:

$$\max x_1+x_2+x_3+x_4−Px_1x_3−Px_1x_4−Px_3x_4−Px_1x_2$$

where the scalar penalty $P$ is arbitrarily chosen to be 3. Thus, our $Q$ matrix is:

Q = torch.tensor([[1, -3, -3, -3],
                  [-3, 1,  0,  0],
                  [-3, 0,  1, -3],
                  [-3, 0, -3,  1]
                 ])

Actually your matrix uses a penalty $P = 6$. Indeed, for some $i \neq j$ the coefficient in front of $x_i x_j$ is $Q_{i, j} + Q_{j, i}$. However, this is not where the problem comes from.

I get x1, x2, x3, x4 = (1, 0, 0, 0) and the objective function evaluates to 1, which is smaller than the expected 2 above. Therefore, this SB solution is not optimal. Have I done something wrong or made a mistake? Any help would be greatly appreciated!

SB is an approximation algorithm, which means that it will not always find the best values, just some "good" values. Actually there is no theoretical guarantee on how close to the optimum those values are. I ran your exemple and I also got the same result. I also ran it using another variant of SB and got the expected result ([0, 1, 0, 1] and 2).

sb.maximize(Q, input_type='binary', ballistic=True, best_only=True)

[bqth29]

Another possibility is to modify the hyperparameters of the SB algorithm with the sb.set_env function. If you look at the way the SB algorithm is implemented, a discretization time step is required for the underlying symplectic Euler scheme that mimics the spin dynamics evolution. By default, this time step is set to 0.1 but on a small case like the one you addressed, there is a possibility that the time step is too big resulting in an optimization that lacks finesse.

I tried setting the time step to 0.01 instead and got the expected optimal result:

import simulated_bifurcation as sb

sb.set_env(time_step=0.01)
sb.maximize(Q, input_type="binary", ballistic=True)

[seanlaw]

Thank you both for your insightful answers. This is fascinating stuff! I have a few more clarifying questions based on your responses:

Actually your matrix uses a penalty $P = 6$. Indeed, for some $i \ne j$ the coefficient in front of $x_ix_j$ is $Q_{ij}+Q_{ji}$. However, this is not where the problem comes from.

@BusyBeaver-42 Can you please elaborate a little on this. In most of the examples that I've come across so far, I understand that the $Q$ matrix must be symmetrical (and non-null) but, based on your comment, it sounds like the symmetrical/complementary off-diagonal elements should sum up to the total penalty in this particular example? So, in essence, we need to divide the total penalty in half and split them between $Q_{ij}$ and $Q_{ji}$ in order to enforce the required symmetry within the $Q$. Is that correct?

Now that I've reviewed other examples more closely, I am realizing that I may have overlooked this subtle point! Thank you for bringing this to my attention. I am learning a lot.

SB is an approximation algorithm, which means that it will not always find the best values, just some "good" values. Actually there is no theoretical guarantee on how close to the optimum those values are. I ran your exemple and I also got the same result. I also ran it using another variant of SB and got the expected result ([0, 1, 0, 1] and 2).

That makes sense. In general, is there a heuristic/estimated ranking/evidence of which variant tends to work better than the others? Is it common to execute the same optimization with different variants and then take the best solution? It seems that ballistic=True is a good place to start if speed is important.

Another possibility is to modify the hyperparameters of the SB algorithm with the sb.set_env function. If you look at the way the SB algorithm is implemented, a discretization time step is required for the underlying symplectic Euler scheme that mimics the spin dynamics evolution. By default, this time step is set to 0.1 but on a small case like the one you addressed, there is a possibility that the time step is too big resulting in an optimization that lacks finesse.

@bqth29 Thank you for explaining this. After an optimization has completed, is there a way to continue from the last final state of a run BUT use a smaller time_step? For example, maybe I use the default time_step=0.1 to generate a pool of top-K candidates, then starting with each of the top-K candidates, I refine the optimization for each candidate using a smaller time_step=0.01?

On a semi-related note, I came across this new research finding recently that is contrary to conventional wisdom for gradient descent, Risky Giant Steps Solve Optimization Problems Faster. Of course, it's a cute discovery but I wonder if we'll see some sort of dynamic time_step (rather than constant) in the future?

[BusyBeaver-42]

@BusyBeaver-42 Can you please elaborate a little on this. In most of the examples that I've come across so far, I understand that the $Q$ matrix must be symmetrical (and non-null) but, based on your comment, it sounds like the symmetrical/complementary off-diagonal elements should sum up to the total penalty in this particular example? So, in essence, we need to divide the total penalty in half and split them between $Q_{ij}$ and $Q_{ji}$ in order to enforce the required symmetry within the $Q$. Is that correct?

Actually $Q$ can be any square real-valued matrix (it does not need to be symmetrical). The quantity which is minimized or maximized is the following.

$$S = \displaystyle \sum_{i=1}^N \displaystyle \sum_{j=1}^N Q_{i,j} ~ x_i x_j$$

Let's say we want to have a coefficient of $-8$ in front of $x_3 x_5$. Note that $x_3 x_5 = x_5 x_3$, hence in the above sum there are two values of $i, j$ we need to take into account: $i=3, j=5$ and $i=5, j=3$. Thus the sum has the following form.

$$S = ... + Q_{3,5} ~ x_3 x_5 + ... + Q_{5,3} ~ x_5 x_3 + ... = (Q_{3,5} + Q_{5,3}) ~ x_3 x_5 + ...$$

Hence we want to have $Q_{3,5} + Q_{5,3} = -8$. A natural choice for that is to set $Q_{3,5} = Q_{5,3} = -4$ (which corresponds to a symmetrical matrix) but if for some reason it is more convenient to set $Q_{3,5} = 100$ and $Q_{5,3} = -108$ when generating the matrix, the energy function would still be the same because $100 + (-108) = (-4) + (-4) = -8$.

Note there is exactly one way to obtain something like $3 x_2$ in the cost function. Indeed first notice that whether $x_2 = 0$ or $x_2 = 1$, $x_2 ^ 2 = x_2$ holds. And to obtain a $3 x_2 x_2$ in $S$, necessarily $i=j=2$.

That makes sense. In general, is there a heuristic/estimated ranking/evidence of which variant tends to work better than the others? Is it common to execute the same optimization with different variants and then take the best solution? It seems that ballistic=True is a good place to start if speed is important.

According to Goto et. al (the physicists who wrote the scientific papers on which this repository is based), for degree 2 cost functions, ballistic SB (ballistic=True) is faster and discrete SB (ballistic=False) is more accurate, for higher-order cost functions (not yet implemented in this repository), ballistic SB is both faster and more accurate. The heating term (heated=True) is supposed to help escape local minima so it should also improve accuracy. Those are just general considerations however; for some instances another combination of hyperparameters will work better (it's mostly empirical, it is not really possible to guess).

On a semi-related note, I came across this new research finding recently that is contrary to conventional wisdom for gradient descent, Risky Giant Steps Solve Optimization Problems Faster.

Actually the situation here is quite different. This type of modifications of gradient-based methods rely on the fact that even if we take a big step and make an error, this cannot happen too often and overall the cost function will keep decreasing. On the contrary, under the hood the SB algorithm actually finds the numerical solution of a differential equation to simulate the evolution of physical system. The larger the time_step, the more the numerical approximation differs from its physical counterpart, which is problematic because it is the physical counterpart which guarantees that the end result has a low energy.

[seanlaw]

This is super insightful! Thank you @BusyBeaver-42 for this wonderful discussion. As a novice, it's not easy for me to acquire these perspectives just from reading the academic paper so I really appreciate your willingness to expand on my trivial questions!