Dbi scheduling

[Sam-XiaoyueLi]

Checklist:

• Reviewers confirm new code works as expected.
• Tests are passing.
• Coverage does not decrease.
• Documentation is updated.

Context of the PR: given a diagonal operator $D$, we want to know the optimal DBI duration $s$ which gives maximal diagonalization gain in one step. The monotonicity relation implies that a linear change of the cost function $||\sigma(e^{-sW}He^{sW})||$ at the beginning, but at some point it needs to turn around, because $W$ cannot diagonalize $H$ in one step. Generally, we can see a curve such as:

We will use analytical means to fit the first minimum (taylor expansion which is more sensitive to the first minimum). Other methods like hyperopt and grid_search are sensitive to the scheduling search interval and may end up in minima other than the first dip (see Fig below). Taking the first minimum may not lead to the global minimum at that iteration step, however, over many iteration steps, it seems to provide a better stability of diagonalization.

---

This PR is an addition to the DBI model, adding utilities in double_bracket.py to select scheduling (time step) strategies via class DoubleBracketScheduling.

These strategies include: use_hyperopt, use_grid_search, and use_polynomial_approximation Within each double bracket iteration, the step is chosen with function choose_step, which returns the optimal step calculated with a given scheduling strategy.

While hyperopt and grid search have been implemented in previous codes, the polynomial approximation is an analytical strategy newly suggested by @wrightjandrew. However, some modifications have been made:

1. The polynomial order is arbitrary (default to 3)
2. Given option for backup scheduling strategy in case no positive real root is found.

Note that additional parameters for each search (max_evals, etc.) can be given as key arguments in choose_step

[codecov]

Codecov Report

Attention: Patch coverage is 80.58252% with 20 lines in your changes are missing coverage. Please review.

Project coverage is 99.80%. Comparing base (056830f) to head (00f6d8c).
Report is 297 commits behind head on dbf_migrate.

❗ Current head 00f6d8c differs from pull request most recent head 93eb20a. Consider uploading reports for the commit 93eb20a to get more accurate results

Files	Patch %	Lines
src/qibo/models/dbi/utils.py	30.43%	16 Missing ⚠️
src/qibo/models/dbi/double_bracket.py	88.88%	3 Missing ⚠️
src/qibo/models/dbi/utils_scheduling.py	98.11%	1 Missing ⚠️

Additional details and impacted files

@@               Coverage Diff               @@
##           dbf_migrate    #1219      +/-   ##
===============================================
- Coverage       100.00%   99.80%   -0.20%     
===============================================
  Files               70       71       +1     
  Lines            10168    10259      +91     
===============================================
+ Hits             10168    10239      +71     
- Misses               0       20      +20     

Flag	Coverage Δ
unittests	99.80% <80.58%> (-0.20%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Edoardo-Pedicillo]

1. The polynomial order is arbitrary (default to 3)

Why is the default value an odd number ? I will expect to use an even order to approximate better a minimum

[Edoardo-Pedicillo]

We already have this function

qibo/src/qibo/models/dbi/double_bracket.py

Lines 108 to 115 in e38a340

	@property
	def diagonal_h_matrix(self):
	"""Diagonal H matrix."""
	return self.backend.cast(np.diag(np.diag(self.backend.to_numpy(self.h.matrix))))
	@property
	def off_diag_h(self):
	return self.h.matrix - self.diagonal_h_matrix

[Sam-XiaoyueLi]

Yes, but here the h is a general Hamiltonian instead of the property h in class DoubleBracketIteration. Is there a good way to reduce the redundancy such that we can have both the property and a function that works for other hamiltonians?

[marekgluza]

maybe the property in line 115 should actually call 228?

[Edoardo-Pedicillo]

I suggest you define a sigma static method inside DoubleBracketIteration class and call it inside the property off_diag_h and polynomial_step.

[marekgluza]

1. The polynomial order is arbitrary (default to 3)

Why is the default value an odd number ? I will expect to use an even order to approximate better a minimum

@Sam-XiaoyueLi I think this is because if we don't go +1 degree the highest degree will be the error term?

so
$$f(x) = \sum_{n=}^N x^n \partial_x^n f(0)/n! + x^{N+1} R_N$$

We want to assume that $R_N\approx 0$ and so that way somehow have a more stable fit. (I think)

@Edoardo-Pedicillo please close comment if you agree? :)

[Sam-XiaoyueLi]

@Edoardo-Pedicillo @marekgluza
The polynomial here equals the derivative of the loss function, to find the optimal step is to find a real positive root of the polynomial. Hence, the highest order being either odd or even should not be a big issue? Please correct me if I understood wrongly.

[marekgluza]

@Edoardo-Pedicillo @marekgluza The polynomial here equals the derivative of the loss function, to find the optimal step is to find a real positive root of the polynomial. Hence, the highest order being either odd or even should not be a big issue? Please correct me if I understood wrongly.

Can you confirm if starting from 3 or 2 or 4 makes a difference? I thought 3 is minimal?

[wrightjandrew]

@Edoardo-Pedicillo @marekgluza The polynomial here equals the derivative of the loss function, to find the optimal step is to find a real positive root of the polynomial. Hence, the highest order being either odd or even should not be a big issue? Please correct me if I understood wrongly.

Can you confirm if starting from 3 or 2 or 4 makes a difference? I thought 3 is minimal?

From the simulations I made, the third order was always the minimum order necessary (although sometimes higher terms were needed). In any case, second-order never seems to work and I think there must be an underlying reason that I couldn't obtain analytically

[Sam-XiaoyueLi]

@Edoardo-Pedicillo @marekgluza The polynomial here equals the derivative of the loss function, to find the optimal step is to find a real positive root of the polynomial. Hence, the highest order being either odd or even should not be a big issue? Please correct me if I understood wrongly.

Can you confirm if starting from 3 or 2 or 4 makes a difference? I thought 3 is minimal?

Yes they make a difference, in general higher orders are better, but not necessarily the even ones. (Here n refers to the highest order in the loss function derivative, so n=2 has 3 coefficients)

[Edoardo-Pedicillo]

Suggested change

backup_scheduling: DoubleBracketScheduling = None,

For transparency, I prefer that the case of no solutions is handled in choose_step

[Sam-XiaoyueLi]

Thanks for the suggestion, the no solution case has been moved to choose_step. However, would you have a good solution for the problem where no additional key arguments can be supplied to the backup scheduling? Meaning that when the backup scheduling runs, only the default values can be used.

[Edoardo-Pedicillo]

Maybe, the best approach is to remove the backup scheduling option and handle the failure of the optimization in another way, understanding better why it does not converge.
For what I have understood, the optimization algorithm fails because the polynomial has not real roots, so I suppose this could mean two things (correct me if I am wrong): either the function has no minima or the approximation is incorrect.
In the first case you should take one of the range boundaries, in the second one maybe increase the polynomial degree is enough.