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.

Sign up for GitHub

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

Add BLOCK-DIAG, MAP-BLOCKS, and PARTIAL-TRACE #104

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Add BLOCK-DIAG, MAP-BLOCKS, and PARTIAL-TRACE #104

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions src/high-level/matrix.lisp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,6 +219,11 @@ ELEMENT-TYPE, CAST, COPY-TENSOR, DEEP-COPY-TENSOR, TREF, SETF TREF)"
"Whether MATRIX is a unitary matrix"
(= matrix (conjugate-transpose matrix) epsilon))

(defun block-matrix-p (matrix)
Copy link
Contributor

@braised-babbage braised-babbage Apr 17, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note: the general usage is that a "block matrix" is a matrix with an associated decomposition, e.g.

[[A, B, C],
 [D, E, F]]

would be a valid block decomposition if nrows(A) = nrows(B) = nrows(C), ncols(A) = ncols(D), and so on. In this usage, the matrix doesn't have to be square. I'd suggest changing the name of your defun and adding a docstring if you want to keep the current behavior.

(and (square-matrix-p matrix)
(let ((dim (nrows matrix)))
(cl:= dim (expt (isqrt dim) 2)))))

(defmacro assert-square-matrix (&rest matrices)
`(progn
,@(loop :for matrix in matrices
Expand Down Expand Up @@ -379,6 +384,17 @@ If fast is t then just change layout. Fast can cause problems when you want to m
(loop :for i :below rows
:collect (tref matrix i i))))))

(defgeneric block-diag (matrix)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same thing here -- seems like you have in mind specifically 2x2 block decompositions.

(:documentation "Get a list of the diagonal block elements of MATRIX")
(:method ((matrix matrix))
(policy-cond:with-expectations (> speed safety)
((assertion (block-matrix-p matrix)))
(let* ((dim (floor (sqrt (nrows matrix)))))
(loop :for i :from 0 :below dim
:for from := (list (* i dim) (* i dim))
:for to := (list (* (1+ i) dim) (* (1+ i) dim))
:collect (slice matrix from to))))))

(defgeneric trace (matrix)
(:documentation "Get the trace of MATRIX (sum of diagonals)")
(:method ((matrix matrix))
Expand All @@ -387,6 +403,24 @@ If fast is t then just change layout. Fast can cause problems when you want to m
(loop :for i :below (nrows matrix)
:sum (tref matrix i i)))))

(defgeneric map-blocks (function matrix &key type)
(:documentation "Map FUNCTION over the NxN blocks of the N²xN² MATRIX, returning a fresh matrix of dimension NxN.")
(:method (function (matrix matrix) &key (type (element-type matrix)))
(policy-cond:with-expectations (> speed safety)
((assertion (block-matrix-p matrix)))
(let* ((dim (floor (sqrt (nrows matrix))))
(freiche (zeros (list dim dim) :type type)))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yum

(loop :for i :from 0 :below dim :do
(loop :for j :from 0 :below dim :do
(let ((from (list (* i dim) (* j dim)))
(to (list (* (1+ i) dim) (* (1+ j) dim))))
(setf (tref freiche i j)
(funcall function (slice matrix from to))))))
freiche))))

(defun partial-trace (matrix)
(map-blocks #'trace matrix))

(defgeneric det (matrix)
(:documentation "Compute the determinant of a square matrix MATRIX")
(:method ((matrix matrix))
Expand Down
3 changes: 3 additions & 0 deletions src/packages.lisp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@

#:map
#:map!
#:map-blocks

#:reshape
#:slice
Expand Down Expand Up @@ -136,6 +137,7 @@
#:scale
#:scale!
#:diag
#:block-diag
#:det
#:upper-triangular
#:lower-triangular
Expand All @@ -146,6 +148,7 @@
#:orthonormalize
#:orthonormalize!
#:trace
#:partial-trace
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gotem

#:direct-sum
#:conjugate-transpose
#:conjugate-transpose!
Expand Down