Calculate complementary polynomial for QSP using flip convolution.

[Epsilon1024]

Calculate complementary polynomial using Eq. 60 of the paper

[Epsilon1024]

very cool! I left some comments (the docstring ones are useful but not urgent).

Could you add some more tests - larger degrees and verifying the actual GQSP unitaries? I think tweaking the precision is okay for now, say if it works with 1e-5 or something. We can use a more powerful optimizer later on to get more precise results, as long as the optimization (loss function) is correct

I'll add in the rest of the tests (I initially forgot to make this PR as a draft - there are a few things I'm planning to add/clean up).

I just noticed that while the 4th order polynomials fit with 1e-5 tolerance, the larger order polynomials have a much looser tolerance which is not great. (A 10th order gave me an error in the 1e-3 range).

We might need to address this sooner rather than later.

[Epsilon1024]

Update: I was able to get the error below 1e-7 by tweaking some parameters. (I used a different miniization method than was used by the authors of the paper, and I set the gradient to be calculated by a 3 point estimation when the default was 2).

The issue that I'm facing now is that the tests that check for all real coefficents are giving a higher error on the unit circle check.

The calculation seems to be a bit different depending on if we want only real coefficients. According to the [authors' github] (https://github.com/Danimhn/GQSP-Code/blob/main/FindingQ.ipynb), the flip convolution is performed differently for real and complex polynomials. It appears that the minimizer they used (and I assume also the minimizer I am using) won't work for complex numbers. As a result, they deconstruct the real and imaginary components (which is why two different methods are being used).

My understanding is that there can be multiple Q polynomials for any given P (which would mean that for an all-real P, we can have both a complex Q and an all-real Q). Please correct me if this is wrong.

Tomorrow, I plan to look into this further and see if I can bring the error down for the real coefficients. If so, then after some minor cleanup, I should be ready to send this out.

[Epsilon1024]

I also found that the test_generalized_real_qsp_with_symbolic_signal_matrix fails with the new qsp method. The test will run 10 random polynomials of each degree. Oddly enough, most of them pass the test (within a tolerance of 1e-5). However, 1 of these polynomials will produce an error in the 1e-1 range. This seems to only happen with 10th degree polynomials. This one does not seem to be a precision issue like the others, but seems to be a bug. However, I can't find anything obvious in the symbolic_qsp_matrix method that could cause the error. The P and Q polynomial pairs pass the unit circle check just before the matrix is made.

[anurudhp]

I also found that the test_generalized_real_qsp_with_symbolic_signal_matrix fails with the new qsp method. The test will run 10 random polynomials of each degree. Oddly enough, most of them pass the test (within a tolerance of 1e-5). However, 1 of these polynomials will produce an error in the 1e-1 range. This seems to only happen with 10th degree polynomials. This one does not seem to be a precision issue like the others, but seems to be a bug. However, I can't find anything obvious in the symbolic_qsp_matrix method that could cause the error. The P and Q polynomial pairs pass the unit circle check just before the matrix is made.

Could you add an issue? You can extract the failing P, Q from your run and report the bug. Also, does this failing polynomial pair pass test_generalized_qsp_with_complex_poly_on_random_unitaries?

The "symbolic gqsp" checker only computes an upper-bound on the error, which may be larger than the actual error (not sure by how much).

Qualtran/qualtran/bloqs/qsp/generalized_qsp_test.py

Lines 296 to 297 in abfbfab

	def upperbound_matrix_norm(M, U) -> float:
	return float(sum(abs(c) for c in sympy.Poly(M, U).coeffs()))

[Epsilon1024]

There are two small followups I will be adding after this PR is submitted:

1. This version does not normalize P before calculating Q. (From my understanding the first version of the default QSP did not normalize P, but this function was added later). I will be adding the normalization in a follow-up PR (that is nearly finished).

2. The main challenge that I faced in this PR was setting a high tolerance for the case with all real polynomials. I found the method in this PR performs comparatively to that in the code provided by the authors of the paper. I will write out my findings in a notebook to go along with this.

[anurudhp]

Looks good! Mostly docstring+type-hint nits. Also please fix the rng comment, best practice is to use np.random.Generator.

[Epsilon1024]

As a note: I ran some checks to see if normalizing P before running the tests on fast_complementary_polynomial would give us more accurate results. From my checks, it performs slightly worse than not being normalized. This tells me that the precision problem is not related to normalization.

[tanujkhattar]

@anurudhp @Epsilon1024 Can we merge this and continue improvements in follow-up PRs? Maybe open an issue to track the potential improvements?