Supporting floating-point exponents in Operator.power

[alexanderivrii]

Summary

Fixes #11454.

Details and comments

In operator.py we simply switch from np.linalg.matrix_power to scipy.linalg.fractional_matrix_power. Thanks to @sadbro for this one-liner. This automatically supports floating-point powers.

In gate.py we switch to using Operator.power method.

Benchmarking

Initially I was sure that the new method (scipy.linalg.fractional_matrix_power, see https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.fractional_matrix_power.html) is faster than the previously implemented method in Gate (based on scipy.linalg.schur), but after running a few more tests (on my windows laptop), I am no longer convinced of this, and this seems to depend on the exponents. I may need to remove the sentence in release notes about the increased performance.

@jakelishman, I know that you have done quite a bit of performance tuning, and so you can probably explain the results I am seeing.

Here is the code that I am running:

def method1(M, t):
    # NEW METHOD
    return scipy.linalg.fractional_matrix_power(M, t)


def method2(M, t):
    # OLD METHOD (adapted from gate.py)
    decomposition, unitary = schur(M, output="complex")
    decomposition_diagonal = decomposition.diagonal()
    decomposition_power = [pow(element, t) for element in decomposition_diagonal]
    unitary_power = unitary @ np.diag(decomposition_power) @ unitary.conj().T
    return unitary_power

time1 = 0
time2 = 0

for seed in range(1000):
	for n in [8, 16, 32, 64]:
		for p in [3, 13, 23, 33, 43]:
			U = random_unitary(n, seed=seed)

			A = U.copy()
			B = U.copy()

			time_start = time.time()
			ret1 = method1(A, p)
			time_end = time.time()
			time1 += time_end - time_start

			time_start = time.time()
			ret2 = method2(B, p)
			time_end = time.time()
			time2 += time_end - time_start

print(f"Method1: time elapsed: {time1:.2f}")
print(f"Method2: time elapsed: {time2:.2f}")

When I run this, I get the time for method1 is 8.36 seconds, while the time for method2 is 25.05 seconds, so we have ~3x performance improvement. However, when I replace the exponents to be floating-point, and explicitly by

[3.1415, 13.5767, 23.4545, 33.4544, 43.122]:

then I get that method 1 takes 119.71 seconds and method 2 26.63 seconds, we we get ~4.5 slowdown.

I still think that whichever method is faster should be implemented in operator.py and gate.py should call that method, but I don't understand which method is better.

[qiskit-bot]

One or more of the the following people are requested to review this:

• @Qiskit/terra-core
• @ikkoham

[jakelishman]

I wouldn't be surprised if the slowdown you're seeing is SciPy is doing a bit more polishing of the final result when asked for fractional_matrix_power instead of schur. You might want to compare the norms $\lVert f^{-1}(f(A)) - A\rVert$ for the two fractional-power functions $f$ and see if one is noticeably closer - that could explain performance differences. Integer powers will likely need less polishing.

Alternatively, integer powers might be being done in a totally different way in the first way: you can evaluate an integer power $A^n$ in $\mathcal O\bigl(\lg(n)\bigr)$ matrix multiplications with a divide-and-conquer strategy - the code looks something like:

def integer_power(matrix, n):
    if n < 0:
        n = -n
        matrix = matrix.inv()
    prev_power, power_of_two = out, matrix.copy()
    out = np.eye(matrix.shape, matrix.dtype)
    while n:
        if n & 1:
            out @= power_of_two
        n >>= 1
        power_of_two @= power_of_two
    return out

That might be significantly faster than calculating the Schur decomposition.

Your results potentially suggest that both of the two things I mentioned are happening.

[coveralls]

Pull Request Test Coverage Report for Build 7539466423

Warning: This coverage report may be inaccurate.

We've detected an issue with your CI configuration that might affect the accuracy of this pull request's coverage report.
To ensure accuracy in future PRs, please see these guidelines.
A quick fix for this PR: rebase it; your next report should be accurate.

• 0 of 0 changed or added relevant lines in 0 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage increased (+1.8%) to 89.359%

---

Totals
Change from base Build 7469929160:	1.8%
Covered Lines:	59412
Relevant Lines:	66487

---

💛 - Coveralls

[alexanderivrii]

Thanks, @jakelishman. I have updated the code to use gate's power method for fractional powers (and removed the release note about the speed-up).

[sbrandhsn]

Excellent work, thank you! :-)