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

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

Closed
wants to merge 1 commit into from

Conversation

stylewarning
Copy link
Member

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
Copy link

genos commented Sep 7, 2022

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
Copy link

genos commented Sep 8, 2022

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
Copy link
Contributor

FWIW there's a pretty serious flaw in my pitch to use gges for this. I detailed a bit more in #841

@stylewarning
Copy link
Member Author

going to close this since it needs re-thinking

@stylewarning stylewarning deleted the diag-to-schur branch September 8, 2022 20:49
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

seeking to simplify find-diagonalizer-in-e-basis
3 participants