Replicate scipy.linalg.qz

[genos]

Hello!

In trying to deal with a violent error arising in QUILC, @stylewarning (quoting @kilimanjaro) proposed simplifying find-diagonalizer-in-e-basis by using generalized Schur decomposition. In order to do this,

we'll need to get qz into MAGICL.

So, my request: could we please have scipy.linalg.qz, viz. DGGES, replicated in MAGICL?

Thank you!

[stylewarning]

FWIW, MAGICL currently ships with

• MAGICL.LAPACK-CFFI::%%DGGES, the raw CFFI binding to the fortran function
• MAGICL.LAPACK-CFFI:%DGGES, an attempt to make the call easy with Lisp-typed objects.

(There's the ZGGES variant as well available.)

Next step would be to determine if we can successfully call %DGGES

[genos]

My Lips chops are somewhat lacking, but it does look like we can at least reach %DGGES:

~ ∃ rlwrap sbcl --dynamic-space-size 4096
This is SBCL 2.2.6, an implementation of ANSI Common Lisp.
More information about SBCL is available at <http://www.sbcl.org/>.

SBCL is free software, provided as is, with absolutely no warranty.
It is mostly in the public domain; some portions are provided under
BSD-style licenses.  See the CREDITS and COPYING files in the
distribution for more information.
* (ql:quickload :magicl)
To load "magicl":
  Load 1 ASDF system:
    magicl
; Loading "magicl"
...............
(:MAGICL)
* (magicl.lapack-cffi:%dgges)

debugger invoked on a SB-INT:SIMPLE-PROGRAM-ERROR @53A4B0A1 in thread
#<THREAD "main thread" RUNNING {100E9C0003}>:
  invalid number of arguments: 0

Type HELP for debugger help, or (SB-EXT:EXIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
  0: [REPLACE-FUNCTION] Call a different function with the same arguments
  1: [CALL-FORM       ] Call a different form
  2: [ABORT           ] Exit debugger, returning to top level.

(MAGICL.LAPACK-CFFI:%DGGES) [external]
   source: (DEFUN %DGGES
                  (JOBVSL JOBVSR SORT SELCTG N A LDA B LDB SDIM ALPHAR
                   ALPHAI ...)
             (DECLARE (INLINE %%DGGES)
                      (TYPE STRING JOBVSL)
                      (TYPE STRING JOBVSR)
                      (TYPE STRING SORT)
                      (TYPE (SIGNED-BYTE 32) SELCTG)
                      (TYPE (SIGNED-BYTE 32) N)
                      (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT) A)
                      (TYPE (SIGNED-BYTE 32) LDA)
                      (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT) B)
                      (TYPE (SIGNED-BYTE 32) LDB)
                      (TYPE (SIGNED-BYTE 32) SDIM)
                      ...)
             (CFFI:WITH-FOREIGN-OBJECTS ((SELCTG-REF6718 ':INT32)
                                         (N-REF6719 ':INT32)
                                         (LDA-REF6721 ':INT32)
                                         (LDB-REF6723 ':INT32)
                                         (SDIM-REF6724 ':INT32)
                                         (LDVSL-REF6729 ':INT32)
                                         (LDVSR-REF6731 ':INT32)
                                         (LWORK-REF6733 ':INT32)
                                         (INFO-REF6735 ':INT32))
               (SETF (CFFI:MEM-REF SELCTG-REF6718 :INT32) SELCTG)
               (SETF (CFFI:MEM-REF N-REF6719 :INT32) N)
               (SETF (CFFI:MEM-REF LDA-REF6721 :INT32) LDA)
               (SETF (CFFI:MEM-REF LDB-REF6723 :INT32) LDB)
               (SETF (CFFI:MEM-REF SDIM-REF6724 :INT32) SDIM)
               (SETF (CFFI:MEM-REF LDVSL-REF6729 :INT32) LDVSL)
               (SETF (CFFI:MEM-REF LDVSR-REF6731 :INT32) LDVSR)
               (SETF (CFFI:MEM-REF LWORK-REF6733 :INT32) LWORK)
               (SETF (CFFI:MEM-REF INFO-REF6735 :INT32) INFO)
               (MAGICL.CFFI-TYPES:WITH-ARRAY-POINTERS ((A-REF6720 A)
                                                       (B-REF6722 B)
                                                       (ALPHAR-REF6725 ALPHAR)
                                                       (ALPHAI-REF6726 ALPHAI)
                                                       (BETA-REF6727 BETA)
                                                       (VSL-REF6728 VSL)
                                                       (VSR-REF6730 VSR)
                                                       (WORK-REF6732 WORK)
                                                       (BWORK-REF6734 BWORK))
                 (%%DGGES JOBVSL JOBVSR SORT SELCTG-REF6718 N-REF6719 A-REF6720
                  LDA-REF6721 B-REF6722 LDB-REF6723 SDIM-REF6724 ALPHAR-REF6725
                  ...))))

[stylewarning]

(side note: You can use

(describe 'foo)

to describe the symbol foo, or

(describe (function foo))

to describe the function named foo.

Similarly

(apropos "dgges")

will search for any and all symbols with dgges in the name.)

[genos]

Thanks, @stylewarning, it's clearly been a long time since I've done anything serious with Common Lisp 😅

I've been trying to replicate the scipy.linalg.qz example but am running into trouble, since (it looks like?) a matrix isn't of the correct type. The following complains

  The value
    #S(MAGICL:MATRIX/DOUBLE-FLOAT
       :NROWS 3
       :NCOLS 3
       :SIZE 9
       :LAYOUT :COLUMN-MAJOR
       :STORAGE #(1.0d0 4.0d0 7.0d0 2.0d0 5.0d0 8.0d0 3.0d0 6.0d0 9.0d0))

  is not of type
    (SIMPLE-ARRAY DOUBLE-FLOAT)
  when binding MAGICL.LAPACK-CFFI::A

What's more, I think we'll probably need to set SORT and SELCTG appropriately, and I could some help passing a pointer for SELCTG.

(ql:quickload '(:alexandria :magicl))

(setf a (magicl:from-list (alexandria:iota 9 :start 1.0d0) '(3 3))
      b (magicl:from-list
          (loop :for _ :from 0 :below 9 :collect (alexandria:gaussian-random))
          '(3 3))
      alphar (magicl:zeros '(9))
      alphai (magicl:zeros '(9))
      beta   (magicl:zeros '(9))
      work   (magicl:zeros '(100))
      info   -1)

(magicl.lapack-cffi:%dgges
  "N"                ; JOBSVL
  "N"                ; JOBSVR
  "N"                ; SORT
  0                  ; SELCTG [not referenced if SORT = "N"]
  (magicl:size a)    ; N
  a                  ; A
  (magicl:nrows a)   ; LDA
  b                  ; B
  (magicl:nrows b)   ; LDB
  0                  ; SDIM [0 if SORT = "N"]
  alphar             ; ALPHAR
  alphai             ; ALPHAI
  beta               ; BETA
  nil                ; VSL [not referenced if JOBVSL = "N"]
  1                  ; LDVSL
  nil                ; VSR [not referenced if JOBVSR = "N"]
  1                  ; LDVSR
  work               ; WORK
  (magicl:size work) ; LWORK
  nil                ; BWORK [not referenced if SORT = "N"]
  info)              ; INFO

[stylewarning]

These % functions won't be able to take MAGICL objects directly. You'd need to pass naked Lisp arrays constructed by

(make-array n :element-type 'double-float :initial-element 0.0d0)

where n is one-dimensional (even for matrices) and stored in column-major format.

Once an API example works with a % function, then we can get the raw Lisp array properly wrapped in MAGICL's clothes

[stylewarning]

(The naked array from a MAGICL object is in its "storage". That can be accessed with magicl::storage, but that's not a public API.)

[stylewarning]

(last point: since Lisp is strongly typed, it may complain about nil being passed in. You may need to use a dummy value of the right type)

[genos]

I thought I had tried make-array, but still got the "it's not a simple array!" complaint; sorry, will try again.

[genos]

Here's an attempt at function-ifying %dgges, absent any sort of argument checking (are a and b matrices of the appropriate type and same shape?) or error handling (did info remain at zero or did DGGES use it to point to a problem?):

(defun dgges (a b)
  (let*
      ((aa     (magicl:deep-copy-tensor a))
       (bb     (magicl:deep-copy-tensor b))
       (n      (magicl:nrows a))
       (nn     (magicl:size a))
       (alphar (make-array nn :element-type 'double-float :initial-element 0d0))
       (alphai (make-array nn :element-type 'double-float :initial-element 0d0))
       (beta   (make-array nn :element-type 'double-float :initial-element 0d0))
       (vsl    (make-array nn :element-type 'double-float :initial-element 0d0))
       (vsr    (make-array nn :element-type 'double-float :initial-element 0d0))
       (lwork  (* 8 (+ n 2)))
       (work   (make-array lwork :element-type 'double-float :initial-element 0d0))
       (bwork  (make-array 0 :element-type '(signed-byte 32)))
       (info   0))
    (progn
      (magicl.lapack-cffi:%dgges
       "V"                  ; JOBSVL
       "V"                  ; JOBSVR
       "N"                  ; SORT
       0                    ; SELCTG [not referenced if SORT = "N"]
       n                    ; N
       (magicl::storage aa) ; A
       n                    ; LDA
       (magicl::storage bb) ; B
       n                    ; LDB
       0                    ; SDIM [0 if SORT = "N"]
       alphar               ; ALPHAR
       alphai               ; ALPHAI
       beta                 ; BETA
       vsl                  ; VSL
       n                    ; LDVSL
       vsr                  ; VSR
       n                    ; LDVSR
       work                 ; WORK
       lwork                ; LWORK
       bwork                ; BWORK [not referenced if SORT = "N"]
       info)                ; INFO
      (values aa bb (magicl:from-array vsl (list n n)) (magicl:from-array vsr (list n n))))))

Attempting the example from scipy.linalg.qz:

(let ((a (magicl:from-list (alexandria:iota 9 :start 1.0d0) '(3 3)))
      (b (magicl:from-list
           (loop :for _ :from 0 :below 9 :collect (alexandria:gaussian-random))
           '(3 3))))
  (multiple-value-bind (aa bb q z) (dgges a b)
    (format t "AA: ~A~%BB: ~A~%Q: ~A~%Z: ~A~%" aa bb q z)))

; AA: #<MATRIX/DOUBLE-FLOAT (3x3):
;      13.548    -3.855    -9.211
;       0.000    -0.375    -1.270
;       0.000     0.000    -0.000>
; BB: #<MATRIX/DOUBLE-FLOAT (3x3):
;       0.355     0.222     0.011
;       0.000     0.638     0.378
;       0.000     0.000     0.958>
; Q: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.876    -0.439    -0.198
;      -0.422     0.898    -0.123
;      -0.232     0.024     0.972>
; Z: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.173    -0.508    -0.843
;      -0.896    -0.274     0.349
;      -0.408     0.816    -0.408>

If this passes muster, could someone with more expertise and experience with MAGICL help with putting together a PR? Ultimately, I'd love it if this could be merged and then this QUILC problem could also be solved.

[stylewarning]

Great progress! Let me investigate and look into making it integrated. I might want your help to write a test or two though. (:

[genos]

I switched the order of vsl and vsr in the return of the above function; edited, with apologies.

[genos]

Also, contrary to the scipy docs (though perhaps similar to their warning about differing from Matlab?), it looks like it should be dgges(a, b) ↦ (aa, bb, q, z) ∋ a = q† aa z & b = q† bb z:

(let ((a (magicl:from-list (alexandria:iota 9 :start 1.0d0) '(3 3)))
      (b (magicl:from-list
          (loop :for _ :from 0 :below 9 :collect (alexandria:gaussian-random))
          '(3 3))))
  (multiple-value-bind (aa bb q z) (dgges a b)
    (format t "A: ~A~%B: ~A~%AA: ~A~%BB: ~A~%Q: ~A~%Z: ~A~%" a b aa bb q z)
    (format t "Q† @ AA @ Z: ~A~%" (magicl:@ (magicl:dagger q) aa z))
    (format t "Q† @ BB @ Z: ~A~%" (magicl:@ (magicl:dagger q) bb z))))

; A: #<MATRIX/DOUBLE-FLOAT (3x3):
;       1.000     2.000     3.000
;       4.000     5.000     6.000
;       7.000     8.000     9.000>)))
; B: #<MATRIX/DOUBLE-FLOAT (3x3):
;       0.481    -0.539    -0.624
;       0.124    -0.163     0.726
;       0.310     0.156    -0.206>
; AA: #<MATRIX/DOUBLE-FLOAT (3x3):
;       13.548    -3.855    -9.211
;        0.000    -0.375    -1.270
;        0.000     0.000    -0.000>
; BB: #<MATRIX/DOUBLE-FLOAT (3x3):
;        0.355     0.222     0.011
;        0.000     0.638     0.378
;        0.000     0.000     0.958>
; Q: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.173    -0.508    -0.843
;      -0.896    -0.274     0.349
;      -0.408     0.816    -0.408>
; Z: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.876    -0.439    -0.198
;      -0.422     0.898    -0.123
;      -0.232     0.024     0.972>
; Q† @ AA @ Z: #<MATRIX/DOUBLE-FLOAT (3x3):
;                 1.000     2.000     3.000
;                 4.000     5.000     6.000
;                 7.000     8.000     9.000>
; Q† @ BB @ Z: #<MATRIX/DOUBLE-FLOAT (3x3):
;                 0.481    -0.539    -0.624
;                 0.124    -0.163     0.726
;                 0.310     0.156    -0.206>

[genos]

Slightly streamlined (leave anonymous that which we don't need) and beefed up (at least check that info is still zero) version. What's more, this follows the same convention as scipy.linalg.qz with regards to q vs. (dagger q) and z vs. (dagger z). I think we can get away with making bwork an array of length zero since the LAPACK docs claim it's unused if SORT = "N", but I didn't want to chance it.

(defun dgges (a b)
  (let*
      ((aa     (magicl:deep-copy-tensor a))
       (bb     (magicl:deep-copy-tensor b))
       (n      (magicl:nrows a))
       (nn     (magicl:size a))
       (q      (magicl:zeros (magicl:shape a)))
       (z      (magicl:zeros (magicl:shape a)))
       (lwork  (* 8 (+ n 2)))
       (info   0)
       (arr (lambda (i &optional (ty 'double-float)) (make-array i :element-type ty))))
    (progn
      (magicl.lapack-cffi:%dgges
       "V"                                ; JOBSVL
       "V"                                ; JOBSVR
       "N"                                ; SORT
       0                                  ; SELCTG [not referenced if SORT = "N"]
       n                                  ; N
       (magicl::storage aa)               ; A
       n                                  ; LDA
       (magicl::storage bb)               ; B
       n                                  ; LDB
       0                                  ; SDIM [0 if SORT = "N"]
       (funcall arr nn)                   ; ALPHAR
       (funcall arr nn)                   ; ALPHAI
       (funcall arr nn)                   ; BETA
       (magicl::storage q)                ; VSL
       n                                  ; LDVSL
       (magicl::storage z)                ; VSR
       n                                  ; LDVSR
       (funcall arr lwork)                ; WORK
       lwork                              ; LWORK
       (funcall arr nn '(signed-byte 32)) ; BWORK [not referenced if SORT = "N"]
       info)                              ; INFO
      (when (not (zerop info)) (error "Nonzero info value; ~A" info))
      (values aa bb q z))))

;;; example from https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.qz.html
(let ((a (magicl:from-list (alexandria:iota 9 :start 1.0d0) '(3 3)))
      (b (magicl:from-list
          (loop :for _ :from 0 :below 9 :collect (alexandria:gaussian-random))
          '(3 3))))
  (multiple-value-bind (aa bb q z) (dgges a b)
    (format t "A: ~A~%B: ~A~%AA: ~A~%BB: ~A~%Q: ~A~%Z: ~A~%" a b aa bb q z)
    (format t "Q @ AA @ Z†∶ ~A~%" (magicl:@ q aa (magicl:dagger z)))
    (format t "Q @ BB @ Z†∶ ~A~%" (magicl:@ q bb (magicl:dagger z)))
    (format t "Q @ AA @ Z† = A? ~A~%" (magicl:= a (magicl:@ q aa (magicl:dagger z))))
    (format t "Q @ BB @ Z† = B? ~A~%" (magicl:= b (magicl:@ q bb (magicl:dagger z))))))

; A: #<MATRIX/DOUBLE-FLOAT (3x3):
;       1.000     2.000     3.000
;       4.000     5.000     6.000
;       7.000     8.000     9.000>)))
; B: #<MATRIX/DOUBLE-FLOAT (3x3):
;       0.481    -0.539    -0.624
;       0.124    -0.163     0.726
;       0.310     0.156    -0.206>
; AA: #<MATRIX/DOUBLE-FLOAT (3x3):
;       13.548    -3.855    -9.211
;        0.000    -0.375    -1.270
;        0.000     0.000    -0.000>
; BB: #<MATRIX/DOUBLE-FLOAT (3x3):
;        0.355     0.222     0.011
;        0.000     0.638     0.378
;        0.000     0.000     0.958>
; Q: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.173    -0.508    -0.843
;      -0.896    -0.274     0.349
;      -0.408     0.816    -0.408>
; Z: #<MATRIX/DOUBLE-FLOAT (3x3):
;      -0.876    -0.439    -0.198
;      -0.422     0.898    -0.123
;      -0.232     0.024     0.972>
; Q @ AA @ Z†: #<MATRIX/DOUBLE-FLOAT (3x3):
;                 1.000     2.000     3.000
;                 4.000     5.000     6.000
;                 7.000     8.000     9.000>
; Q @ BB @ Z†: #<MATRIX/DOUBLE-FLOAT (3x3):
;                 0.481    -0.539    -0.624
;                 0.124    -0.163     0.726
;                 0.310     0.156    -0.206>
; Q @ AA @ Z† = A? T
; Q @ BB @ Z† = B? T

[genos]

As for tests, would a few examples checked via fiasco do the trick? Given my druthers, I'd go for property-based testing of

(and
  (= a (@ q aa (dagger z))
  (= b (@ q bb (dagger z)))

for suitably floating-point-fuzzy versions of =, but I don't know the standard operating procedure for MAGICL.

[stylewarning]

As for tests, would a few examples checked via fiasco do the trick? Given my druthers, I'd go for property-based testing of

(and

  (= a (@ q aa (dagger z))

  (= b (@ q bb (dagger z)))

for suitably floating-point-fuzzy versions of =, but I don't know the standard operating procedure for MAGICL.

Yes, totally works.

[stylewarning]

@genos I should be able to just turn your example function into a test. Don't worry about it. THat's good enough for me.

If you have additional tests you'd like to write, though, by all means!

[genos]

Thanks @stylewarning, I really appreciate your responsiveness on this!