Implement random unitary matrices #40381
Conversation
I'm getting a few warnings about tests in this file that take longer than 5s to run (without --long). Here we add "# long time" tags.
|
@JohnCremona you are mentioned by name in the paper :) The CI, lint, and docs all look broken in the latest develop branch but none of it looks related this this PR. I have not been able to check that the docs look OK however. |
I guess that will be referring to a letter to the editor of the Americal Mathematical Monthly which I wrote on 23 Feb 1987! |
| S = A - A.conjugate_transpose() | ||
| U = (S-I).inverse()*(S+I) | ||
| D = identity_matrix(F,n) | ||
| for i in range(n): |
There was a problem hiding this comment.
Why create D, you could reuse I?
There was a problem hiding this comment.
(This code has been sitting in my personal library for a while.) I was probably in a mood where I considered it a sin to have a variable named I that was no longer the identity matrix.
Regardless, I've replaced the diagonal \pm 1 matrix by a loop that uses U.set_row_to_multiple_of_row, so now it's irrelevant. The new method is written in cython so it should be faster anyway. Thanks for the suggestion.
src/sage/matrix/special.py
Outdated
| U = (S-I).inverse()*(S+I) | ||
| D = identity_matrix(F,n) | ||
| for i in range(n): | ||
| if ZZ.random_element(2).is_zero(): |
There was a problem hiding this comment.
Is this really the best way to toss a coin? random()<0.5 would be as good?
There was a problem hiding this comment.
Probably an irrational fear that my 50/50 probability will be wrong by the unit roundoff. I've replaced it. random() will definitely be faster.
src/sage/matrix/special.py
Outdated
| D = identity_matrix(F,n) | ||
| for i in range(n): | ||
| if ZZ.random_element(2).is_zero(): | ||
| D[i,i] *= F(-1) |
There was a problem hiding this comment.
It seems crazy that field elements are not required to have a negate() mthod, but apparently not. But they do have a method called neg() with double underscores for evaluatio of unary minus.
I don't think it is necessary to construct F(-1) every time through the loop. Perhaps D[i,i] = -D[i,i] is better?
There was a problem hiding this comment.
This is handled within set_row_to_multiple_of_row() now (it casts the -1 to a rational for the flint API).
It matched the paper closely but was definitely not the best way. I've improved it, and have mentioned why we are negating some rows. Thanks for the review! |
|
positive review if the failing tests are irrelevant |
New reference (Liebeck and Osborne, 1991) to be used in a new method for generating random unitary matrices.
Add a new function to generate random unitary matrices based on the paper "The Generation of All Rational Orthogonal Matrices" by Liebeck and Osborne (1991).
Use the lower-level methods Matrix.set_row_to_multiple_of_row() and random.random() to implement random_unitary_matrix(). This should bring with it a small speed increase, and actually makes the code simpler. Thanks to John Cremona for the suggestion.
|
Force push to fix the tab character in the reference file |
Add a new
random_unitary_matrix()function that is called frommatrix.random(..., algorithm='unitary').