unitary eigenvalue decomposition #2736
Conversation
balopat
force-pushed
the
unitary_eig
branch
3 times, most recently
from
491e1a1
to
b54e48b
Compare
There was a problem hiding this comment.
IIUC, the problem with np.linalg.eig is not that the second return value does not contain valid eigenvectors. Instead, we complain the eigenvectors are not orthonormalized (in the case of nearly degenerate eigenvalues). If this is so, wouldn't it be more straightforward to take the eigenvectors from np.linalg.eig and use something like Gram-Schmidt process to orthonormalize them? The advantage of that is that we can limit orthonormalization to the degenerate eigenspaces only.
Oh wow, I'm not sure how I never thought about this. We also already have a similar function: |
cirq/linalg/decompositions.py
Outdated
| atol: float = 1e-8, | ||
| ) -> Tuple[np.array, List[np.ndarray]]: | ||
| """An eigendecomposition that ensures eigenvectors are perpendicular. | ||
| def orthogonal_eigendecompose( |
There was a problem hiding this comment.
nit: Since we're also normalizing, this should be orthonormal_eigendecompose.
There was a problem hiding this comment.
LGTM, thanks!
Before we run ahead and merge it - from the numpy thread there was a recommendation to use the Schur decomposition instead, from scipy.
I tried it and it works really well:
def unitary_eig(matrix: np.ndarray,
check_preconditions: bool = True,
atol: float = 1e-8,
rtol: float = 1e-5) -> Tuple[np.array,
np.ndarray]:
if check_preconditions:
if not predicates.is_unitary(matrix, atol=atol, rtol=rtol):
raise ValueError(
'Input must correspond to a unitary matrix '
f'.Received input:\n{matrix}')
R, V = scipy.linalg.schur(matrix)
eigvals = R.diagonal()
return eigvals, VI guess the only thing is that here we really do have to check for the precondition for the matrix to be unitary otherwise R is going to be upper triangular and not diagonal.
But, this seems to pass on all the tests and is much faster on the example performance tests (for some reason the kak perf tests don't change that much or they get slightly worse).
What do you think about swapping out the orthonormal_eig to this implementation instead?
Below is the output of check/performance: 0004_8612933 is master and NOW is this branch with the Schur decomposition (+ unitary check).
----------------------------------------------------------------------------------------------------------------------- benchmark: 52 tests ------------------------------------------------------------------------------------------------------------------------
Name (time in us) Min Max Mean StdDev Median IQR Outliers OPS Rounds Iterations
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
test_kak_decomposition_perf[target2] (NOW) 453.1560 (1.0) 4,172.8890 (4.10) 520.4031 (1.04) 122.9678 (1.78) 488.3245 (1.04) 69.7980 (2.66) 101;98 1,921.5872 (0.96) 1596 1
test_kak_decomposition_perf[target2] (0004_8612933) 454.0000 (1.00) 1,018.0200 (1.0) 498.4696 (1.0) 76.3795 (1.11) 468.8810 (1.0) 26.2720 (1.0) 166;251 2,006.1403 (1.0) 1695 1
test_kak_decomposition_perf[target5] (NOW) 495.4170 (1.09) 1,308.0510 (1.28) 577.5186 (1.16) 69.0747 (1.0) 566.9970 (1.21) 80.2255 (3.05) 216;50 1,731.5460 (0.86) 1504 1
test_kak_decomposition_perf[target4] (0004_8612933) 497.5530 (1.10) 2,130.5120 (2.09) 592.6066 (1.19) 122.6044 (1.77) 548.9960 (1.17) 104.7110 (3.99) 160;100 1,687.4600 (0.84) 1847 1
test_kak_decomposition_perf[target4] (NOW) 498.6510 (1.10) 1,111.1840 (1.09) 584.0454 (1.17) 70.8838 (1.03) 573.4960 (1.22) 73.3287 (2.79) 220;60 1,712.1957 (0.85) 1441 1
test_kak_decomposition_perf[target3] (NOW) 501.4240 (1.11) 1,182.3750 (1.16) 586.6997 (1.18) 78.5033 (1.14) 573.6480 (1.22) 80.8970 (3.08) 134;57 1,704.4495 (0.85) 1336 1
test_kak_decomposition_perf[target3] (0004_8612933) 508.3360 (1.12) 1,126.6820 (1.11) 576.2538 (1.16) 92.7423 (1.34) 542.7360 (1.16) 63.5370 (2.42) 182;168 1,735.3465 (0.87) 1667 1
test_kak_decomposition_perf[target5] (0004_8612933) 520.2390 (1.15) 1,768.1310 (1.74) 604.9632 (1.21) 128.5244 (1.86) 552.1750 (1.18) 80.6710 (3.07) 40;36 1,652.9932 (0.82) 585 1
test_kak_decomposition_perf[target0] (NOW) 525.7840 (1.16) 2,812.9000 (2.76) 656.8992 (1.32) 177.6252 (2.57) 623.8220 (1.33) 106.1257 (4.04) 95;95 1,522.3035 (0.76) 1283 1
test_kak_decomposition_perf[target0] (0004_8612933) 540.2330 (1.19) 1,055.2170 (1.04) 589.2598 (1.18) 70.0335 (1.01) 568.0630 (1.21) 45.8470 (1.75) 108;89 1,697.0442 (0.85) 1028 1
test_kak_decomposition_perf[target1] (0004_8612933) 544.3260 (1.20) 1,460.0730 (1.43) 618.8983 (1.24) 109.5644 (1.59) 582.7875 (1.24) 48.2155 (1.84) 180;244 1,615.7745 (0.81) 1692 1
test_kak_decomposition_perf[target1] (NOW) 551.6020 (1.22) 3,115.9590 (3.06) 648.0131 (1.30) 126.7628 (1.84) 601.3465 (1.28) 88.4350 (3.37) 124;100 1,543.1787 (0.77) 1554 1
test_kak_decomposition_perf[target7] (0004_8612933) 570.4760 (1.26) 2,703.0570 (2.66) 654.9197 (1.31) 127.1578 (1.84) 620.7400 (1.32) 73.2280 (2.79) 140;149 1,526.9048 (0.76) 1526 1
test_kak_decomposition_perf[target10] (NOW) 577.9560 (1.28) 1,683.4710 (1.65) 676.3860 (1.36) 121.2330 (1.76) 638.6635 (1.36) 62.4930 (2.38) 139;164 1,478.4458 (0.74) 1506 1
test_kak_decomposition_perf[target14] (NOW) 582.6810 (1.29) 1,612.6760 (1.58) 669.4592 (1.34) 110.9669 (1.61) 636.2255 (1.36) 43.1510 (1.64) 119;177 1,493.7430 (0.74) 1522 1
test_kak_decomposition_perf[target13] (0004_8612933) 584.0150 (1.29) 2,817.8430 (2.77) 719.0112 (1.44) 190.4825 (2.76) 648.9170 (1.38) 107.2962 (4.08) 128;139 1,390.7990 (0.69) 1549 1
test_kak_decomposition_perf[target12] (NOW) 584.5760 (1.29) 3,288.4840 (3.23) 703.7022 (1.41) 129.4419 (1.87) 692.1585 (1.48) 93.4710 (3.56) 72;61 1,421.0556 (0.71) 1050 1
test_kak_decomposition_perf[target7] (NOW) 585.0490 (1.29) 2,420.7480 (2.38) 686.2563 (1.38) 122.0377 (1.77) 648.3270 (1.38) 78.5292 (2.99) 101;94 1,457.1814 (0.73) 1337 1
test_kak_decomposition_perf[target15] (NOW) 585.5960 (1.29) 1,759.6880 (1.73) 669.3903 (1.34) 104.6757 (1.52) 644.3340 (1.37) 59.3775 (2.26) 97;102 1,493.8967 (0.74) 1252 1
test_kak_decomposition_perf[target6] (NOW) 589.3650 (1.30) 1,408.6630 (1.38) 675.7858 (1.36) 101.3429 (1.47) 653.8495 (1.39) 80.2770 (3.06) 117;94 1,479.7587 (0.74) 1306 1
test_kak_decomposition_perf[target9] (NOW) 590.7170 (1.30) 1,390.9090 (1.37) 683.9849 (1.37) 106.8238 (1.55) 643.7000 (1.37) 74.2195 (2.83) 127;110 1,462.0205 (0.73) 1344 1
test_kak_decomposition_perf[target10] (0004_8612933) 592.1240 (1.31) 4,040.5500 (3.97) 745.7587 (1.50) 273.4335 (3.96) 693.9795 (1.48) 137.6690 (5.24) 97;116 1,340.9163 (0.67) 1566 1
test_kak_decomposition_perf[target13] (NOW) 592.8710 (1.31) 1,473.6770 (1.45) 672.8973 (1.35) 111.1342 (1.61) 637.4530 (1.36) 68.9900 (2.63) 115;115 1,486.1110 (0.74) 1305 1
test_kak_decomposition_perf[target14] (0004_8612933) 597.7220 (1.32) 2,436.2690 (2.39) 720.2268 (1.44) 156.6939 (2.27) 667.0000 (1.42) 122.7630 (4.67) 117;87 1,388.4515 (0.69) 1506 1
test_kak_decomposition_perf[target8] (NOW) 598.4040 (1.32) 1,723.1220 (1.69) 716.2001 (1.44) 93.5505 (1.35) 720.2300 (1.54) 82.6960 (3.15) 338;78 1,396.2579 (0.70) 1503 1
test_kak_decomposition_perf[target11] (NOW) 600.1940 (1.32) 1,608.9120 (1.58) 697.3785 (1.40) 105.6700 (1.53) 667.8805 (1.42) 96.0940 (3.66) 131;79 1,433.9416 (0.71) 1558 1
test_kak_decomposition_perf[target11] (0004_8612933) 604.9370 (1.33) 2,744.7690 (2.70) 736.5345 (1.48) 220.2540 (3.19) 654.5280 (1.40) 106.6078 (4.06) 83;97 1,357.7096 (0.68) 1029 1
test_kak_decomposition_perf[target9] (0004_8612933) 608.1980 (1.34) 3,661.9660 (3.60) 828.4471 (1.66) 379.5916 (5.50) 728.1485 (1.55) 200.6960 (7.64) 113;128 1,207.0776 (0.60) 1554 1
test_kak_decomposition_perf[target15] (0004_8612933) 609.0550 (1.34) 2,686.7020 (2.64) 725.7042 (1.46) 180.5735 (2.61) 658.1670 (1.40) 119.5748 (4.55) 119;111 1,377.9718 (0.69) 1417 1
test_kak_decomposition_perf[target8] (0004_8612933) 611.0050 (1.35) 3,011.5530 (2.96) 775.5193 (1.56) 261.5512 (3.79) 667.5750 (1.42) 192.8190 (7.34) 133;101 1,289.4586 (0.64) 1437 1
test_kak_decomposition_perf[target12] (0004_8612933) 615.6960 (1.36) 2,744.1190 (2.70) 705.3442 (1.42) 147.6240 (2.14) 657.7530 (1.40) 100.9790 (3.84) 92;82 1,417.7476 (0.71) 1312 1
test_kak_decomposition_perf[target6] (0004_8612933) 615.8580 (1.36) 2,538.9690 (2.49) 708.3957 (1.42) 139.9136 (2.03) 659.1830 (1.41) 77.1545 (2.94) 111;116 1,411.6404 (0.70) 1397 1
test_example_runs_hello_qubit_perf (NOW) 1,015.5640 (2.24) 2,300.7970 (2.26) 1,164.0445 (2.34) 170.5853 (2.47) 1,115.1520 (2.38) 141.6220 (5.39) 75;47 859.0737 (0.43) 741 1
test_example_runs_hello_qubit_perf (0004_8612933) 1,248.0460 (2.75) 6,003.6430 (5.90) 1,980.8137 (3.97) 783.2578 (11.34) 1,663.9245 (3.55) 502.0320 (19.11) 77;77 504.8430 (0.25) 598 1
test_example_runs_bernstein_vazirani_perf (NOW) 3,061.5090 (6.76) 8,235.7070 (8.09) 4,302.6822 (8.63) 720.6611 (10.43) 4,271.4155 (9.11) 764.4850 (29.10) 39;5 232.4132 (0.12) 150 1
test_example_runs_bernstein_vazirani_perf (0004_8612933) 4,106.5720 (9.06) 10,394.2620 (10.21) 7,310.6942 (14.67) 1,725.8788 (24.99) 7,132.9370 (15.21) 3,245.9300 (123.55) 25;0 136.7859 (0.07) 56 1
test_example_runs_bell_inequality_perf (NOW) 4,115.4600 (9.08) 7,079.9930 (6.95) 4,636.8957 (9.30) 407.3006 (5.90) 4,508.1790 (9.61) 402.1525 (15.31) 26;12 215.6615 (0.11) 215 1
test_example_runs_bell_inequality_perf (0004_8612933) 4,237.5410 (9.35) 9,211.6460 (9.05) 5,491.2814 (11.02) 1,366.6121 (19.78) 4,904.7460 (10.46) 1,579.3915 (60.12) 18;8 182.1069 (0.09) 117 1
test_example_runs_quantum_teleportation (0004_8612933) 4,413.9710 (9.74) 33,032.9230 (32.45) 8,475.7128 (17.00) 5,006.5829 (72.48) 6,395.2270 (13.64) 4,689.5430 (178.50) 18;14 117.9842 (0.06) 180 1
test_example_runs_quantum_teleportation (NOW) 4,504.1140 (9.94) 13,085.0870 (12.85) 4,925.2735 (9.88) 654.1715 (9.47) 4,775.0970 (10.18) 296.1650 (11.27) 9;13 203.0344 (0.10) 215 1
test_example_runs_grover_perf (NOW) 5,412.6560 (11.94) 9,628.3770 (9.46) 6,730.1221 (13.50) 758.3356 (10.98) 6,590.1550 (14.06) 932.2530 (35.48) 46;5 148.5857 (0.07) 166 1
test_example_runs_grover_perf (0004_8612933) 5,768.8730 (12.73) 15,042.2260 (14.78) 9,878.8287 (19.82) 2,885.1158 (41.77) 10,256.8345 (21.88) 5,475.0000 (208.40) 53;0 101.2266 (0.05) 120 1
test_example_runs_quantum_fourier_transform_perf (0004_8612933) 5,774.6350 (12.74) 11,258.2760 (11.06) 6,929.1787 (13.90) 1,194.9793 (17.30) 6,454.4910 (13.77) 1,117.9678 (42.55) 23;13 144.3172 (0.07) 131 1
test_example_runs_quantum_fourier_transform_perf (NOW) 5,807.0090 (12.81) 10,545.2530 (10.36) 6,572.6624 (13.19) 697.0730 (10.09) 6,445.7480 (13.75) 556.6220 (21.19) 12;8 152.1453 (0.08) 137 1
test_example_runs_superdense_coding (NOW) 34,812.7420 (76.82) 44,272.0770 (43.49) 37,148.5952 (74.53) 2,198.9678 (31.83) 36,333.3065 (77.49) 1,622.1690 (61.75) 6;3 26.9189 (0.01) 26 1
test_example_runs_superdense_coding (0004_8612933) 35,452.6480 (78.23) 96,629.7320 (94.92) 50,196.1411 (100.70) 16,805.7937 (243.30) 42,732.3185 (91.14) 25,740.1135 (979.75) 5;0 19.9219 (0.01) 28 1
test_example_runs_phase_estimator_perf (NOW) 36,722.3290 (81.04) 45,907.6260 (45.10) 39,443.0159 (79.13) 1,986.9435 (28.77) 39,156.3855 (83.51) 2,009.2300 (76.48) 7;1 25.3530 (0.01) 26 1
test_example_runs_phase_estimator_perf (0004_8612933) 41,181.2260 (90.88) 76,534.4080 (75.18) 50,187.7883 (100.68) 8,872.4523 (128.45) 48,483.3870 (103.40) 10,701.7480 (407.34) 3;1 19.9252 (0.01) 19 1
test_example_runs_bcs_mean_field_perf (0004_8612933) 177,387.9150 (391.45) 296,750.0450 (291.50) 235,800.2110 (473.05) 43,874.2855 (635.17) 230,514.5450 (491.63) 54,515.8460 (>1000.0) 2;0 4.2409 (0.00) 5 1
test_example_runs_bcs_mean_field_perf (NOW) 179,008.9870 (395.03) 191,844.1560 (188.45) 185,138.8793 (371.41) 4,649.1307 (67.31) 185,091.5000 (394.75) 7,129.5350 (271.37) 2;0 5.4014 (0.00) 6 1
test_example_runs_hello_line_perf (NOW) 208,783.8220 (460.73) 221,297.9500 (217.38) 213,784.1302 (428.88) 4,842.1199 (70.10) 214,079.2080 (456.57) 6,182.9333 (235.34) 2;0 4.6776 (0.00) 5 1
test_example_runs_hello_line_perf (0004_8612933) 240,935.8590 (531.68) 345,331.2300 (339.22) 274,591.2658 (550.87) 41,590.1253 (602.10) 256,283.7030 (546.59) 43,548.8955 (>1000.0) 1;0 3.6418 (0.00) 5 1
There was a problem hiding this comment.
Note that for real orthogonal matrices the use of the complex factorization needs to be made explicit by passing output='complex'.
There was a problem hiding this comment.
Also, degenerate eigenvalues are not the only problem case with eig. Even if the are no repeated eigenvalues the eigenvectors may not be orthogonal and the matrix must be diagonalized via a similarity transformation, i.e., you need to use the matrix inverse instead of the Hermitian transpose, So normality is still needed even if there are no degenerate eigenvalues. For example
array([[0, 0, 1],
[0, 1, 0],
[0, 0, 2]])
Has eigenvalues [0, 1, 2] and a complete set of eigenvectors, but the eigenvectors are not orthogonal.
There was a problem hiding this comment.
@balopat I like the approach using Schur decomposition. Also, it seems that KAK decomposition performance tests generally do get a speedup as well, but the output is in weird order. Let's make sure we set output='complex' since many interesting unitaries are real.
@charris Makes sense. In our case we will only use this to compute eigendecomposition of unitary matrices, so we have normality.
cirq/linalg/decompositions.py
Outdated
| the same eigenspace and whether eigenvectors are perpendicular. | ||
| atol: Absolute threshold for determining whether eigenvalues are from | ||
| the same eigenspace and whether eigenvectors are perpendicular. | ||
| matrix: a unitary matrix. If not unitary, this method is not |
There was a problem hiding this comment.
s/unitary/normal/ (twice on this line)
|
hi @Strilanc can you please merge? |
| the same eigenspace and whether eigenvectors are perpendicular. | ||
| matrix: a normal matrix. If not normal, this method is not | ||
| guaranteed to return correct eigenvalues. | ||
| check_preconditions: when true and matrix is not unitary, |
There was a problem hiding this comment.
Should be normal instead of unitary, right?
| for i in range(len(g)): | ||
| vecs[g[i]] = q[:, i] | ||
| R, V = scipy.linalg.schur(matrix, output="complex") | ||
| if check_preconditions and not predicates.is_diagonal(R, atol=atol): |
There was a problem hiding this comment.
Preconditions should be checked before performing the computation, right?
While working on #451, I ran into issues with
numpy.linalg.eig(numpy/numpy#15461) that this method helps with. As discussed with @Strilanc, this might be useful to expose for the Cirq community.