Extend Simulation For Convergence

[seanlaw]

I have a simple $Q$ matrix:

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

and when I do:

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

I see:

Iterations:   0%|                                     | 0/10000 [00:00<?, ?it/s]
Iterations: 100%|██████████████████████| 10000/10000 [00:00<00:00, 48377.88it/s]
Bifurcated agents:   0%|                                | 0/128 [00:00<?, ?it/s]
2023-10-04 20:27:29,642 - simulated_bifurcation_optimizer - WARNING - No agent has converged. Returned final positions' signs instead.

What should/could I do to extend the simulations to hopefully reach convergence? I've tried playing around with variables such as max_steps and sampling_period by neither has led to convergence. Any hints or tips would be greatly appreciated!

Update

Doing sb.set_env(time_step=0.001) allowed things to converge but I'm hoping that there may be other options that I could play with.

[BusyBeaver-42]

$Q$ is an extremely special matrix

Short version

There are many solutions which achieve the maximal value. Eventhough the agents have no trouble reaching and staying in the solution space, the agents have trouble staying close to the same solution long enough for the convergence criterion to be met.

Long version

Let $a$, $b$, $c$, $d$ be four binary variables (that is variables in $\set{0, 1}$). Let $x = (a ~ b ~ c ~ d)^T$ and $y = (b ~ c ~ d)^T$. Define the following two matrices.

$$Q = \begin{pmatrix} 1 & 0 & 0 & 0 \\\ 0 & 0 & -0.5 & 0.5 \\\ 0 & -0.5 & 0 & 1 \\\ 0 & 0.5 & 1 & -2 \end{pmatrix}$$ $$P = \begin{pmatrix} 0 & -0.5 & 0.5 \\\ -0.5 & 0 & 1 \\\ 0.5 & 1 & -2 \end{pmatrix}$$

Matrices $Q$ and $P$ are closely related.

$$Q = \begin{pmatrix} 1 & \begin{matrix} 0 & 0 & 0 \end{matrix} \\\ \begin{matrix} 0 \\ 0 \\ 0 \end{matrix} & \mbox{\Huge $P$ } \end{pmatrix}$$

Thus $x^T Q x = a + y^T P y$ (using $a^2 = a$ because $a$ is a binary variable).

Hence the variable $a$ is not coupled with any other variable and the instance is equivalent to maximizing $y^T P y$.

Here are the all the possible values for $y^T P y$.

$b$	$c$	$d$	$y^T Q y$
0	0	0	1
0	0	1	-1
0	1	0	1
0	1	1	1
1	0	0	1
1	0	1	0
1	1	0	0
1	1	1	1

As you can see, 5 binary vectors out of 8 achieve the maximal value: the energy level of the ground states is highly degenerate. Hence the agents have no incentive to prefer one ground state over another which explains the warning about the convergence.

However, by printing the energy levels of the agents you can see that they all converge quickly in terms of energy level. This means than the agents easily reach the solution space, but they do not stay close to a given solution, instead they alternate between solutions.

Currently, the convergence is checked using the state of the agents which is an issue when the energy level of the ground state is highly degenerate as in this case. I am discussing with @bqth29 and @lorpugliese the implementation of a convergence criterion based on the energy levels.

Note that there is only one solution to the minimization of $x^T Q x$ and when running your example using sb.minimize instead of sb.maximize all the agents do converge quickly.

General remarks about the stopping window

The stopping window is something we introduced to allow for early stopping, Goto et al.'s paper uses only max_steps (also we recently implemented a timeout functionality, it will be added in the next release). There are no theoretical guarantees about actually achieving this convergence criterion.

Doing sb.set_env(time_step=0.001) allowed things to converge but I'm hoping that there may be other options that I could play with.

I believe that drastically reducing time_step like that without increasingsampling_period or convergence_threshold can lead to false positives: the states at step, step + sampling_period, step + 2 * sampling_period, ..., are more likely to be equal because they correspond to instants closer in time, not because the agents did converge.

Also as I mentionned above, the stopping window is currently based on the state of the agents whereas it may be more relevant to base it on their energy.

Thank you for showing us this example, this is a case we overlooked and it will help us improve the package! 😃

[seanlaw]

Also as I mentionned above, the stopping window is currently based on the state of the agents whereas it may be more relevant to base it on their energy.

I am thinking if it makes sense to design the API so that it has a default stopping function (i.e., based on energy/agents) but that also can be swapped out by the user for another user-defined function. Of course, this user function would need to have access to the internal state of the simulator (energies, agents, etc). I've seen this done elsewhere where there is something behind the scenes (using Python inspect) that checks whether the user provided function has the right function signature (using functools.partial can also help a lot here to package up a user defined function nicely). So, rather than hardcoding one or two stopping functions, the user can have a little bit more control and define a more complex criterion that, say, considers energy+agents+max_steps+some other constraint (e.g., total wall clock time).

Additionally, if possible, it would be helpful to have more documentation regarding all of the different parameters and environment settings that are available. Right now, all of those things are mostly opaque to the user and, for the ones that are briefly mentioned, it isn't clear what changing those inputs "really" mean. Again, I understand that documentation isn't glamorous work and y'all have done an amazing job thus far so please do not take this as criticism.

Thank you for showing us this example, this is a case we overlooked and it will help us improve the package! 😃

No problem! I manage a high-traffic Python package myself called STUMPY and so I really appreciate the discourse here as I am getting to learn from true experts. In case it matters, the example above (and others that I am working my way through) came from a max 2-sat problem that I was trying to reproduce (See Section 4.3).

[seanlaw]

Also, it might be more useful to use warnings.warn instead of logging.warning as the former can be caught like an exception and the user can specify some other settings to try that might improve convergence and not require human intervention. Currently, with logging.warning, there's no easy way to programmatically know that the simulation did not converge (there's no way to introspect the agents or the state of the simulator that I am aware of).

Please let me know if I should submit a new issue to request this.

[BusyBeaver-42]

Hi Sean, sorry for the late reply, I am a bit overwhelmed at moment.

I am thinking if it makes sense to design the API so that it has a default stopping function (i.e., based on energy/agents) but that also can be swapped out by the user for another user-defined function.

I think this is a nice feature, we will probably implement it at some point. Some big changes are planned for the inner workings of the package (mostly related to the implementation of SB for arbitrary degree polynomials), so this probably will not be added before since it will cause compatibility issues.

Additionally, if possible, it would be helpful to have more documentation regarding all of the different parameters and environment settings that are available. Right now, all of those things are mostly opaque to the user and, for the ones that are briefly mentioned, it isn't clear what changing those inputs "really" mean.

Next release (should be coming soon), will consist mostly of bug fixes and a few features. This documentation will be added in the following release.

Also, it might be more useful to use warnings.warn instead of logging.warning as the former can be caught like an exception and the user can specify some other settings to try that might improve convergence and not require human intervention. Currently, with logging.warning, there's no easy way to programmatically know that the simulation did not converge (there's no way to introspect the agents or the state of the simulator that I am aware of).

Please let me know if I should submit a new issue to request this.

Thanks for the feedback, this will (probably) be fixed in the next release. Yes please open an issue!