Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix handling of Operator.power for non-unitary operators. #13319

Open
wants to merge 2 commits into
base: main
Choose a base branch
from

Conversation

alexanderivrii
Copy link
Contributor

@alexanderivrii alexanderivrii commented Oct 14, 2024

Summary

Fixes #13307 following @jakelishman's suggestion in the description.

Details and comments

After a quick round of profiling on my laptop (Python 10, Windows 11), the old (Schur-based) way to raise an operator to a power is 3x-5x faster than calling scipy.linalg.fractional_matrix_power, so it makes sense to use it when the operator is known to be unitary.

I have renamed the old method to be called Operator.power_if_unitary, this assumes that the given operator is unitary. The Operator.is_power does not make this assumption (and calls fractional_matrix_power). The gate class uses the power_if_unitary method.

P.S. In addition, for the power method and an integer power, I don't see much runtime difference when calling fractional_matrix_power and matrix_power. For the power_if_unitary method calling matrix_power is faster than calling the Schur-based implementation.

@alexanderivrii alexanderivrii added the Changelog: Bugfix Include in the "Fixed" section of the changelog label Oct 14, 2024
@alexanderivrii alexanderivrii added this to the 1.3.0 milestone Oct 14, 2024
@alexanderivrii alexanderivrii requested a review from a team as a code owner October 14, 2024 12:26
@qiskit-bot
Copy link
Collaborator

One or more of the following people are relevant to this code:

  • @Qiskit/terra-core

@@ -22,6 +22,7 @@
from typing import TYPE_CHECKING

import numpy as np
import scipy.linalg
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question: should we import scipy.linalg at the top-level, or locally within the power/power_if_unitary functions (as before)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We've historically forced scipy imports to be runtime imports because the cost of initialising too many scipy submodules is very high, and it made import qiskit feel very sluggish for users (and developers). We still have a lint against it, so chances are there'll be a lint failure if it's a top-level import.

That said, in the last couple of years, scipy have done a lot to make their imports much lazier and on-demand. That means that import scipy now makes scipy.linalg available (but imported only on demand) whereas previously you had to manually import it. I believe they may also have applied some of the lazy loading to subcomponents of the subpackages as well, so the imports might not be as costly as they used to be.

At any rate, we still have the lint right now, so it'll need to be a runtime import.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! Changed in a82b123.

And this also answers my would be followup question of why we call import scipy on the top-level in qsd.py.

@coveralls
Copy link

coveralls commented Oct 14, 2024

Pull Request Test Coverage Report for Build 11327888022

Details

  • 12 of 12 (100.0%) changed or added relevant lines in 2 files are covered.
  • 5 unchanged lines in 3 files lost coverage.
  • Overall coverage increased (+0.03%) to 88.677%

Files with Coverage Reduction New Missed Lines %
crates/qasm2/src/expr.rs 1 94.02%
qiskit/quantum_info/operators/operator.py 1 94.12%
crates/qasm2/src/lex.rs 3 92.73%
Totals Coverage Status
Change from base Build 11305032303: 0.03%
Covered Lines: 73213
Relevant Lines: 82561

💛 - Coveralls

Copy link
Member

@jakelishman jakelishman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This generally seems sensible to me, thanks for the fix.

A minor question (and I'm fine either way): do you prefer adding a separate power_if_unitary method, or adding an assume_unitary=False flag to the existing power method? Either's fine, I'm just making sure we've thought about the interface we want.

Comment on lines +544 to +547
"""Return the matrix power of the operator, provided that operator is
unitary.

This is faster than the ``power`` method.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we might want to put a bigger warning in the docstring of this method that it'll produce nonsense if the operator wasn't unitary.

Also, super minor, but ``power`` can become :meth:`power` to insert a cross-ref.

Comment on lines 565 to 573
import scipy.linalg

# Experimentally, for fractional powers this seems to be 3x faster than
# calling scipy.linalg.fractional_matrix_power(self.data, n)
# calling scipy.linalg.fractional_matrix_power(self.data, exponent)
decomposition, unitary = scipy.linalg.schur(self.data, output="complex")
decomposition_diagonal = decomposition.diagonal()
decomposition_power = [pow(element, n) for element in decomposition_diagonal]
decomposition_power = [pow(element, exponent) for element in decomposition_diagonal]
unitary_power = unitary @ np.diag(decomposition_power) @ unitary.conj().T
ret._data = unitary_power
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just made #13358 to sort out the SciPy 1.14 issues, which as a side-effect modified all this code to use the fraction_matrix_power form with a little rotation to push the branch cut somewhere less likely to cause issues (Lev's idea). That obsoletes the "bugfix" here in power - sorry, I was just trying to get SciPy 1.14 working quicker.

That said, I think the new method / assume_unitary keyword arg to this function is still a very useful idea, so this note is just to make sure that we apply the same branch-cut protection trick here too, and ideally add some tests that it's stable for various roots.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Changelog: Bugfix Include in the "Fixed" section of the changelog
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Operator.power does not handle non-unitary inputs correctly
4 participants