From 019700199a23d935dc48f8077738564922909338 Mon Sep 17 00:00:00 2001 From: "Juan M. Bello-Rivas" Date: Wed, 11 Sep 2019 11:20:41 -0700 Subject: [PATCH] Accelerate APPLY-QUBIT-PERMUTATION. --- dqvm/src/permutation.lisp | 292 ++++++++++++++++---------- dqvm/tests/distributed-qvm-tests.lisp | 2 - dqvm/tests/permutation-tests.lisp | 51 +++++ 3 files changed, 233 insertions(+), 112 deletions(-) diff --git a/dqvm/src/permutation.lisp b/dqvm/src/permutation.lisp index 3704c70f..f9f30185 100644 --- a/dqvm/src/permutation.lisp +++ b/dqvm/src/permutation.lisp @@ -4,33 +4,63 @@ (in-package #:dqvm2) -;;; A simple implementation of a permutation data structure. - -;;; Note that (make-permutation) and NIL both represent the identity. +;;; Permutation classes for permuting sets of qubits. +;;; +;;; The value NIL represents the identity permutation. General permutations +;;; are embodied by the PERMUTATION-GENERAL class. Permutations involving a +;;; single transposition swapping 0 with another qubit are represented by the +;;; PERMUTATION-TRANSPOSITION class. +;;; +;;; The generic function APPLY-QUBIT-PERMUTATION does the heavy lifting of +;;; permuting addresses and the class hierarchy laid out here allows us to +;;; accomplish significant speed-ups (applying a PERMUTATION-TRANSPOSITION is +;;; more than three times faster than the equivalent application of a +;;; PERMUTATION-GENERAL object). (deftype transposition () '(or null (cons alexandria:non-negative-fixnum alexandria:non-negative-fixnum))) (defclass permutation () + () + (:documentation "Base class for permutations.")) + +(defclass permutation-general (permutation) ((number-of-transpositions - :initarg :number-of-transpositions :type alexandria:non-negative-integer + :initarg :number-of-transpositions :documentation "Number of transpositions defining the permutation.") (transpositions - :initarg :transpositions :type list + :initarg :transpositions :reader permutation-transpositions :documentation "Bijective map determined by transpositions, stored as an association list sorted by CAR.")) (:default-initargs :transpositions nil) - (:documentation "Permutation acting on sets of qubit indices.")) + (:documentation "Arbitrary permutation acting on sets of qubit indices.")) + +(defclass permutation-transposition () + ((tau + :type (unsigned-byte 6) ; Implies a maximum of 2⁶ = 64 qubits. + :initarg :tau + :initform (error-missing-initform :tau) + :documentation "Positive value of τ in π = (0 τ).")) + (:documentation "Specialized permutation involving a single transposition of the form π = (0 τ) where τ ≠ 0.")) + +(defmethod permutation-transpositions ((permutation permutation-transposition)) + (let ((tau (slot-value permutation 'tau))) + (list (cons 0 tau) (cons tau 0)))) -(defmethod print-object ((permutation permutation) stream) +(defmethod print-object ((permutation permutation-general) stream) (print-unreadable-object (permutation stream :type t :identity t) (let ((transpositions (permutation-transpositions permutation))) (format stream "~:[~:A~;~{~A~^ ~}~]" transpositions transpositions)))) +(defmethod print-object ((permutation permutation-transposition) stream) + (print-unreadable-object (permutation stream :type t :identity t) + (let ((tau (slot-value permutation 'tau))) + (format stream "(0 . ~D) (~D . 0)" tau tau)))) + (defun-inlinable make-permutation (&optional transpositions) "Allocate a permutation defined by TRANSPOSITIONS. @@ -43,11 +73,10 @@ DQVM2> (make-permutation '((2 . 1) (1 . 0))) # Note that in the example above, the transposition (0 2) was automatically added." - (declare (optimize (speed 3) (safety 0)) + (declare #.qvm::*optimize-dangerously-fast* (type list transpositions)) - (let ((permutation (make-instance 'permutation)) - (transpositions* nil) + (let ((transpositions* nil) (domain nil) (codomain nil)) @@ -59,6 +88,8 @@ Note that in the example above, the transposition (0 2) was automatically added. (error "Malformed permutation. A mapping ~D ↦ ~D already existed." (first z) (rest z)))))) + (declare (inline check-transposition)) + (loop :for (a . b) :in transpositions :do (check-transposition a b) (unless (= a b) @@ -69,93 +100,85 @@ Note that in the example above, the transposition (0 2) was automatically added. (loop :for a :of-type alexandria:non-negative-fixnum :in (set-difference codomain domain) :for b :of-type alexandria:non-negative-fixnum - :in (nset-difference domain codomain) + :in (set-difference domain codomain) :unless (= a b) :do (pushnew (cons a b) transpositions* :test #'equal)) - (setf (slot-value permutation 'number-of-transpositions) (length transpositions*) - (slot-value permutation 'transpositions) (sort transpositions* #'< :key #'first)) - - permutation)) - -(defun-inlinable inverse-permutation (permutation) - "Return the inverse of PERMUTATION." - (declare (optimize (speed 3) (safety 0))) - (when permutation - - (let ((inverse-permutation (make-instance 'permutation)) - (transpositions (permutation-transpositions permutation))) - - (setf (slot-value inverse-permutation 'transpositions) (loop :for (a . b) :in transpositions :collect (cons b a)) - (slot-value inverse-permutation 'number-of-transpositions) (slot-value permutation 'number-of-transpositions)) - - inverse-permutation))) - -(defun is-identity-permutation-p (permutation) - "Return T if PERMUTATION is the identity, NIL otherwise." - (if (or (null permutation) (null (permutation-transpositions permutation))) - t - nil)) - -(defun-inlinable apply-permutation (permutation item) - "Apply PERMUTATION to ITEM. - -Examples --------- - -DQVM2> (apply-permutation (make-permutation) 42) -42 - -DQVM2> (apply-permutation (make-permutation '((2 . 0))) 2) -0 - -DQVM2> (apply-permutation (make-permutation '((2 . 1) (1 . 0))) 2) -1" - (declare (optimize (speed 3) (safety 0)) - (type (or null permutation) permutation) - (type alexandria:non-negative-fixnum item)) - (the alexandria:non-negative-fixnum - (if permutation - (alexandria:if-let ((transposition (assoc item (permutation-transpositions permutation)))) - (rest transposition) - item) - item))) - -(defun-inlinable apply-inverse-permutation (permutation item) - "Apply PERMUTATION⁻¹ to ITEM." - (apply-permutation (inverse-permutation permutation) item)) + (cond + ((and (null domain) (null codomain)) nil) + ((and (= 1 (length domain)) + (zerop (min (the qvm:amplitude-address (first domain)) + (the qvm:amplitude-address (first codomain))))) + (make-instance 'permutation-transposition + :tau (max (the qvm:amplitude-address (first domain)) + (the qvm:amplitude-address (first codomain))))) + ((and (= 2 (length domain)) + (null (set-difference domain codomain)) + (zerop (the qvm:amplitude-address (apply #'min domain)))) + (make-instance 'permutation-transposition :tau (apply #'max domain))) + (t + (make-instance 'permutation-general :number-of-transpositions (length transpositions*) + :transpositions (sort transpositions* #'< :key #'first)))))) + +(defgeneric inverse-permutation (permutation) + (:documentation "Return the inverse of PERMUTATION.") + (declare #.qvm::*optimize-dangerously-fast*)) + +(defmethod inverse-permutation ((permutation (eql nil))) + nil) + +(defmethod inverse-permutation ((permutation permutation-transposition)) + permutation) + +(defmethod inverse-permutation ((permutation permutation-general)) + (make-instance 'permutation-general + :transpositions (loop :for (a . b) :in (permutation-transpositions permutation) :collect (cons b a)) + :number-of-transpositions (slot-value permutation 'number-of-transpositions))) + +(defgeneric is-identity-permutation-p (permutation) + (:documentation "Return T if PERMUTATION is the identity, NIL otherwise.")) + +(defmethod is-identity-permutation-p ((permutation (eql nil))) + t) + +(defmethod is-identity-permutation-p ((permutation permutation-transposition)) + nil) ; By construction PERMUTATION-TRANSPOSITION objects cannot be the identity. + +(defmethod is-identity-permutation-p ((permutation permutation-general)) + (null (permutation-transpositions permutation))) (defun compose-permutations (&rest permutations) "Return a new permutation that is the composition of PERMUTATIONS. If PERMUTATIONS is the list π₁, π₂, ..., πₛ, then the result is the composition π₁ ∘ π₂ ∘ ... ∘ πₛ. In other words, the composition starts from right to left as in standard mathematical notation." - (let (transpositions) - - (let (domain) - ;; Aggregate the domain of the composed permutation to get a list of - ;; all possible relevant inputs. - (loop :for permutation :in permutations :when permutation :do - (loop :for transposition :of-type transposition :in (permutation-transpositions permutation) :do - (let ((a (first transposition))) - (declare (type alexandria:non-negative-fixnum a)) - (pushnew a domain)))) - - ;; Now map each domain element to obtain transpositions. - (loop :with codomain := (coerce domain 'vector) - :for permutation :in (nreverse permutations) :when permutation :do - (loop :for i :from 0 :for b :across codomain :do - (setf (aref codomain i) - (apply-permutation permutation (aref codomain i)))) - :finally - (loop :for a :of-type alexandria:non-negative-fixnum :in domain - :for b :of-type alexandria:non-negative-fixnum :across codomain - :unless (= a b) :do - (pushnew (cons a b) transpositions :test #'equal)))) + (let ((transpositions nil) + (domain nil)) + + ;; Aggregate the domain of the composed permutation to get a list of + ;; all possible relevant inputs. + (loop :for permutation :in permutations :when permutation :do + (loop :for transposition :of-type transposition :in (permutation-transpositions permutation) :do + (let ((a (first transposition))) + (declare (type alexandria:non-negative-fixnum a)) + (pushnew a domain)))) + + ;; Now map each domain element to obtain transpositions. + (loop :with codomain := (coerce domain 'vector) + :for permutation :in (nreverse permutations) :when permutation :do + (loop :for i :from 0 :for b :across codomain :do + (setf (aref codomain i) + (apply-permutation permutation (aref codomain i)))) + :finally + (loop :for a :of-type alexandria:non-negative-fixnum :in domain + :for b :of-type alexandria:non-negative-fixnum :across codomain + :unless (= a b) :do + (pushnew (cons a b) transpositions :test #'equal))) (make-permutation transpositions))) -(defun-inlinable apply-qubit-permutation (permutation address) - "Apply PERMUTATION to an index ADDRESS within a wavefunction. +(defgeneric apply-qubit-permutation (permutation address) + (:documentation + "Apply PERMUTATION to an index ADDRESS within a wavefunction. Examples -------- @@ -165,36 +188,56 @@ DQVM2> (apply-qubit-permutation (make-permutation '((2 . 0))) #b100) DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :base 2) 100 -4" +4") + (declare #.qvm::*optimize-dangerously-fast*)) + +(defmethod apply-qubit-permutation ((permutation (eql nil)) address) + address) + +(defmethod apply-qubit-permutation ((permutation permutation-transposition) address) + (declare #.qvm::*optimize-dangerously-fast* + (type (or null permutation) permutation) + ;; (type qvm:amplitude-address address) + (type (unsigned-byte 64) address) ; Imposed maximum number of qubits. + (values qvm:amplitude-address)) + + (let ((tau (slot-value permutation 'tau))) + (declare (type (unsigned-byte 6) tau)) + + (rotatef (ldb (byte 1 0) address) (ldb (byte 1 tau) address)) + address)) + +(defmethod apply-qubit-permutation ((permutation permutation) address) ;; Alternatively, in-place permutations could be implemented following: ;; ;; F. Fich, J. Munro, and P. Poblete, “Permuting in Place,” SIAM ;; J. Comput., vol. 24, no. 2, pp. 266–278, Apr. 1995. - (declare (optimize (speed 3) (safety 0)) + (declare #.qvm::*optimize-dangerously-fast* (type (or null permutation) permutation) - (type qvm:amplitude-address address)) - - (the qvm:amplitude-address - (if permutation - (let* ((transpositions (slot-value permutation 'transpositions)) - (number-of-transpositions (slot-value permutation 'number-of-transpositions)) - (bit-vector (make-array number-of-transpositions :element-type 'bit))) - ;; (declare (dynamic-extent bit-vector)) - - (loop :for index :from 0 - :for transposition :in transpositions :do - (setf (bit bit-vector index) (ldb (byte 1 (first transposition)) - address))) - - (loop :for index :from 0 - :for transposition :of-type transposition :in transpositions :do - (setf address (dpb (bit bit-vector index) - ;; (byte 1 (the (unsigned-byte 6) (rest transposition))) ; Enable this for speed (assumes a maximum of 64 qubits). - (byte 1 (rest transposition)) - address)) - :finally (return address))) - address))) + ;; (type qvm:amplitude-address address) + (type (unsigned-byte 64) address) ; Imposed maximum number of qubits. + (values qvm:amplitude-address)) + + (let* ((transpositions (slot-value permutation 'transpositions)) + (number-of-transpositions (slot-value permutation 'number-of-transpositions)) + (bit-vector (make-array number-of-transpositions :element-type 'bit))) + (declare (type (integer 0 128) number-of-transpositions) + (dynamic-extent bit-vector)) + + (loop :for index :from 0 + :for transposition :in transpositions :do + (setf (bit bit-vector index) (ldb (byte 1 (first transposition)) + address))) + + (loop :for index :from 0 + :for transposition :of-type transposition :in transpositions :do + (setf address (the qvm:amplitude-address + (dpb (bit bit-vector index) + (byte 1 (the (unsigned-byte 6) (rest transposition))) ; Enable this for speed (assumes a maximum of 64 qubits). + ;; (byte 1 (rest transposition)) + address))) + :finally (return address)))) (defun-inlinable apply-inverse-qubit-permutation (permutation address) (apply-qubit-permutation (inverse-permutation permutation) address)) @@ -213,3 +256,32 @@ DQVM2> (write (apply-qubit-permutation (make-permutation '((2 . 0))) #b001) :bas (dotimes (i1 max-value (values)) (let ((i2 (apply-qubit-permutation permutation i1))) (format stream control-string i1 i1 i2 i2))))) + +(defun-inlinable apply-permutation (permutation item) + "Apply PERMUTATION to ITEM. + +Examples +-------- + +DQVM2> (apply-permutation (make-permutation) 42) +42 + +DQVM2> (apply-permutation (make-permutation '((2 . 0))) 2) +0 + +DQVM2> (apply-permutation (make-permutation '((2 . 1) (1 . 0))) 2) +1" + (declare #.qvm::*optimize-dangerously-fast* + (type (or null permutation) permutation) + (type alexandria:non-negative-fixnum item) + (values alexandria:non-negative-fixnum)) + + (if permutation + (alexandria:if-let ((transposition (assoc item (permutation-transpositions permutation)))) + (rest transposition) + item) + item)) + +(defun-inlinable apply-inverse-permutation (permutation item) + "Apply PERMUTATION⁻¹ to ITEM." + (apply-permutation (inverse-permutation permutation) item)) diff --git a/dqvm/tests/distributed-qvm-tests.lisp b/dqvm/tests/distributed-qvm-tests.lisp index e56adf2c..17a7c473 100644 --- a/dqvm/tests/distributed-qvm-tests.lisp +++ b/dqvm/tests/distributed-qvm-tests.lisp @@ -25,8 +25,6 @@ "Find amplitude addresses to exchange when applying NEXT-PERMUTATION and the rank where the amplitudes are located. Returns four sequences: current addresses, new addresses, and the source and target addresses." - (check-type next-permutation permutation) - (let ((permutation (permutation addresses)) (effective-permutation (dqvm2::get-effective-permutation addresses next-permutation)) diff --git a/dqvm/tests/permutation-tests.lisp b/dqvm/tests/permutation-tests.lisp index 50bf6c99..82a1028f 100644 --- a/dqvm/tests/permutation-tests.lisp +++ b/dqvm/tests/permutation-tests.lisp @@ -26,6 +26,13 @@ (is (= (apply-inverse-permutation permutation 1) 1)) (is (= (apply-inverse-permutation permutation 2) 2))) + (let ((permutation (make-permutation '((2 . 1))))) + (is (eq (type-of permutation) 'dqvm2::permutation-general))) + + (let ((permutation (make-permutation '((2 . 0))))) + (is (eq (type-of permutation) 'dqvm2::permutation-transposition)) + (is (eq permutation (inverse-permutation permutation)))) + (let ((permutation (make-permutation '((2 . 1) (1 . 0))))) (is (not (is-identity-permutation-p permutation))) @@ -72,3 +79,47 @@ (is (= (apply-permutation composition 0) 2)) (is (= (apply-permutation composition 1) 1)) (is (= (apply-permutation composition 2) 0)))) + +(deftest benchmark-apply-qubit-permutation () + (labels ((get-elapsed-time-in-seconds (start stop) + "Compute elapsed time in seconds between START and STOP." + (float (/ (- stop start) internal-time-units-per-second))) + + (time-apply-qubit-permutation (permutation number-of-qubits) + "Measure the time taken by calls to APPLY-QUBIT-PERMUTATION on addresses from 0 to 2^NUMBER-OF-QUBITS." + (let ((start (get-internal-real-time))) + (dotimes (x (expt 2 number-of-qubits)) + (let ((y (apply-qubit-permutation permutation x))) + (values x y))) + (get-elapsed-time-in-seconds start (get-internal-real-time)))) + + ;; (time-map-reordered-amplitudes (permutation number-of-qubits) + ;; "Measure the time taken by calls to MAP-REORDERED-AMPLITUDES on addresses from 0 to 2^NUMBER-OF-QUBITS." + ;; (let ((nat-tuple (apply #'qvm::nat-tuple + ;; (mapcar (lambda (x) + ;; (apply-permutation permutation x)) + ;; (loop :for i :below number-of-qubits :collect i)))) + ;; (start (get-internal-real-time))) + ;; (qvm::map-reordered-amplitudes 0 (lambda (x y) (values x y)) nat-tuple) + ;; (get-elapsed-time-in-seconds start (get-internal-real-time)))) + ) + + (let* ((tau 4) + (number-of-qubits 24) + (permutation-0 (make-instance 'dqvm2::permutation-transposition :tau tau)) + (permutation-1 (make-instance 'dqvm2::permutation-general + :number-of-transpositions 2 + :transpositions (list (cons 0 tau) (cons tau 0))))) + + (loop :for x :below (expt 2 (1+ tau)) :do + (is (= (apply-qubit-permutation permutation-0 x) + (apply-qubit-permutation permutation-1 x)))) + + (is (> (/ (time-apply-qubit-permutation permutation-1 number-of-qubits) + (time-apply-qubit-permutation permutation-0 number-of-qubits)) + 3)) + + ;; (is (> (/ (time-map-reordered-amplitudes permutation-0 number-of-qubits) + ;; (time-apply-qubit-permutation permutation-0 number-of-qubits)) + ;; 7)) + )))