quil-lang / quilc

The optimizing Quil compiler.
Apache License 2.0
452 stars 73 forks source link

Work on removing FIND-DIAGONALIZER-IN-E-BASIS #845

Closed stylewarning closed 1 year ago

stylewarning commented 1 year ago

This is an attempt to replace F-D-I-E-B with @kilimanjaro's idea of using a Schur decomposition.

This aims to fix #841

CC @genos

genos commented 1 year ago

I think the original Python implementation from #841 fell prey to the note from scipy.linalg.qz's documentation, namely

Notes Q is transposed versus the equivalent function in Matlab.

Specifically, in the comment

# Generalized Schur decomposition writes
#    Re[U] = L D_r R^T,  Im[U] = L D_i R^T

I think it should be L^T and R:

In [1]: from scipy.linalg import norm, qz

In [2]: from scipy.stats import unitary_group

In [3]: rng = unitary_group(dim=2, seed=1729)

In [4]: print(u := rng.rvs())
[[-0.29585687+0.47303725j -0.76750556+0.31565756j]
 [ 0.71123856-0.42760284j -0.55613711-0.04479994j]]

In [5]: diag_r, diag_i, left, right = qz(u.real, u.imag)

In [6]: v = (left.T @ diag_r @ right) + 1j * (left.T @ diag_i @ right)

In [7]: norm(u - v)
Out[7]: 3.0436364816617974e-16
genos commented 1 year ago

I take back my earlier comment; for unitary U, I think orthogonal-decomposition-of-UU^T does work correctly. Defining

(defun gen-gaussian (n)
  (magicl:from-list
    (loop :for i :from 0 :below (* n n) :collect (alexandria:gaussian-random))
    (list n n)))

(defun gen-unitary (n)
  (let* ((m (magicl:.complex (gen-gaussian n) (gen-gaussian n)))
         (z (magicl:./ m (sqrt 2))))
    (multiple-value-bind (q r) (magicl:qr z)
      (let* ((d (magicl:diag r))
             (normalized-d
               (magicl:from-list
                 (loop :for x :in d :collect (/ x (abs x)))
                 (list n 1)))
             (matrix-d (magicl:hstack (loop :for _ :from 0 :below n :collect normalized-d))))
        (magicl:.* q matrix-d)))))

and setting *check-math* to T, the following works without complaining:

* (loop :for _ :from 0 :below 100 :do (orthogonal-decomposition-of-uu^t (gen-unitary 100)))
NIL

As such, I'm not sure what's up with the CI failures.

braised-babbage commented 1 year ago

FWIW there's a pretty serious flaw in my pitch to use gges for this. I detailed a bit more in https://github.com/quil-lang/quilc/issues/841

stylewarning commented 1 year ago

going to close this since it needs re-thinking