knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE) options(width = 80, tibble.width = Inf)
quad.form()
et seqIn versions prior to 1.2-19, the emulator package included a serious
bug in the quad.form()
family of functions in which the complex
conjugate of the correct answer was returned (which did not matter in
my usual use-case because my matrices were Hermitian). This short
vignette demonstrates that the bug has been fixed. Note that the fix
was considerably more complicated than simply returning the complex
conjugate of the old functions' value, which would have been terribly
inefficient. The actual fix avoids taking more conjugates than
absolutely necessary. The vignette checks all the functions in the
series, including the ones that have not been changed such as
quad.form.inv()
. First load the package:
library("emulator")
We need a helper function to create random complex matrices (NB: we
cannot use the cmvnorm
package because that depends on the
emulator
package):
rcm <- function(row,col){ matrix(rnorm(row*col)+1i*rnorm(row*col),row,col) }
Then use this function to define a square matrix M
with complex
entries (NB: not Hermitian!), and a couple of rectangular matrices,
also complex:
rcm <- function(row,col){matrix(rnorm(row*col)+1i*rnorm(row*col),row,col)} M <- rcm(2,2) x <- rcm(2,3) y <- rcm(3,2) x1 <- rcm(2,3) y1 <- rcm(3,2)
Set up a numerical tester function:
tester <- function(a,b,TOL=1e-13){stopifnot(all(abs(a-b)< TOL))}
(previous versions used a tolerance of 1e-15
, which was
occasionally not met). Now test each function:
ht(x)
= $x^*$ = $\overline{x'}$ (Hermitian transpose):ht(x)=t(Conj(x))
(jj1 <- Conj(t(x))) (jj2 <- t(Conj(x))) (jj3 <- ht(x)) tester(jj1,jj3) tester(jj2,jj3)
cprod()
= $x^*y$:cprod(x,y)=crossprod(Conj(x),y)
(jj1 <- ht(x) %*% x1) (jj2 <- cprod(x,x1)) tester(jj1,jj2)
tcprod()
= $x y^*$:tcprod(x,y)=crossprod(x,Conj(y))
(jj1 <- ht(x1) %*% x) (jj2 <- cprod(x1,x)) tester(jj1,jj2)
quad.form()
= $x^*Mx$:quad.form(M,x)=crossprod(crossprod(M,Conj(x)),x))
(jj1 <- ht(x) %*% M %*% x) (jj2 <- quad.form(M,x)) tester(jj1,jj2)
quad.form.inv()
= $x^*M^{-1}x$:quad.form.inv(M,x)=cprod(x,solve(M,x))
(jj1 <- ht(x) %*% solve(M) %*% x) (jj2 <- quad.form(solve(M),x)) max(abs(jj1-jj2))
quad.3form()
= $x^*My$:quad.3form(M,l,r)=crossprod(crossprod(M,Conj(l)),r)
(jj1 <- ht(x) %*% M %*% x1) (jj2 <- quad.3form(M,x,x1)) tester(jj1,jj2)
quad.3tform()
= $xMy^*$:quad.3tform(M,l,r)=tcrossprod(left,tcrossprod(Conj(right),M))
(jj1 <- y %*% M %*% ht(y1)) (jj2 <- quad.3tform(M,y,y1)) tester(jj1,jj2)
quad.tform()
= $xMx^*$:quad.tform(M,x)=tcrossprod(x,tcrossprod(Conj(x),M))
(jj1 <- y %*% M %*% ht(y)) (jj2 <- quad.tform(M,y)) tester(jj1,jj2)
quad.tform.inv()
= $xM^{-1}x^*$:quad.tform.inv(M,x)=quad.form.inv(M,ht(x))
(jj1 <- y %*% solve(M) %*% ht(y)) (jj2 <- quad.tform.inv(M,y)) tester(jj1,jj2)
quad.diag()
= $\operatorname{diag}(x^*Mx)$ = diag(quad.form())
:quad.diag(M,x)=colSums(crossprod(M,Conj(x)) * x)
(jj1 <- diag(ht(x) %*% M %*% x)) (jj2 <- diag(quad.form(M,x))) (jj3 <- quad.diag(M,x)) tester(jj1,jj3) tester(jj2,jj3)
quad.tdiag()
= $\operatorname{diag}(xMx^*)$ = diag(quad.tform())
:quad.tdiag(M,x)=rowSums(tcrossprod(Conj(x), M) * x)
(jj1 <- diag(y %*% M %*% ht(y))) (jj2 <- diag(quad.tform(M,y))) (jj3 <- quad.tdiag(M,y)) tester(jj1,jj3) tester(jj2,jj3)
quad.3diag()
= $\operatorname{diag}(x^*My)$quad.3diag(M,l,r)=colSums(crossprod(M, Conj(left)) * right)
(jj1 <- diag(ht(x) %*% M %*% x1)) (jj2 <- diag(quad.3form(M,x,x1))) (jj3 <- quad.3diag(M,x,x1)) tester(jj1,jj3) tester(jj2,jj3)
quad.3tdiag()
= $\operatorname{diag}(xMy^*)$quad.3tdiag(M,l,r)=colSums(t(left) * tcprod(M, right))
(jj1 <- diag(y %*% M %*% ht(y1))) (jj2 <- diag(quad.3tform(M,y,y1))) (jj3 <- quad.3tdiag(M,y,y1)) tester(jj1,jj3) tester(jj2,jj3)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.