version 0.0-1

DESCRIPTION

@@ -0,0 +1,22 @@
+Package: quadform
+Type: Package
+Title: Efficient Evaluation of Quadratic Forms
+Version: 0.0-1
+Authors@R: person(given=c("Robin", "K. S."), family="Hankin", role = c("aut","cre"), email="hankin.robin@gmail.com", comment = c(ORCID = "0000-0001-5982-0415"))
+Depends: R (>= 3.0.1)
+Imports: mathjaxr
+Suggests: testthat
+Maintainer: Robin K. S. Hankin <hankin.robin@gmail.com>
+Description: 
+ A range of quadratic forms are evaluated, using efficient methods.
+ Unnecessary transposes are not performed.  Complex values are handled
+ consistently.
+License: GPL
+URL: https://github.com/RobinHankin/quadform
+BugReports: https://github.com/RobinHankin/quadform/issues
+RdMacros: mathjaxr
+NeedsCompilation: no
+Packaged: 2024-04-09 18:16:01 UTC; rhankin
+Author: Robin K. S. Hankin [aut, cre] (<https://orcid.org/0000-0001-5982-0415>)
+Repository: CRAN
+Date/Publication: 2024-04-10 16:00:02 UTC

MD5

@@ -0,0 +1,12 @@
+c108da68ce801fa9d1725ac9507d5541 *DESCRIPTION
+78b35acdd842def4b9d53e591ce8f14d *NAMESPACE
+074e49250a3df0fb9c1c08707bc48feb *R/quadform.R
+c31e38ebe173747c3c0cd8afb5ef1dc6 *README.md
+058d288c653d6b965d1d64085d40b24b *build/quadform.pdf
+8efd13461874a09c24a2de9f1d079f54 *build/stage23.rdb
+5ce56307342698572faee22b10046a4f *inst/quad_form_test.Rmd
+889abde9e951a0ccad03cbeed593e9fd *inst/quad_form_test.html
+4fc4335d77e243224211bf4bf7c2f708 *man/figures/quadform.png
+f1a1034c941c027d4a09b32c7888d640 *man/quad.form.Rd
+5ab6fe48d5c16b3b7476ab646d2d2241 *tests/testthat.R
+40fde0ef8ed23e13cbef7167c783646b *tests/testthat/test_aaa.R

NAMESPACE

@@ -0,0 +1,2 @@
+exportPattern("^[^\\.]")
+import("mathjaxr")

R/quadform.R

@@ -0,0 +1,73 @@
+`ht` <- function(x){
+    if(is.complex(x)){
+        return(t(Conj(x)))
+    } else {
+        return(t(x))
+    }
+}
+
+`cprod` <- function(x,y=NULL){
+    if(is.complex(x) | is.complex(y)){
+        if(is.null(y)){
+            return(crossprod(Conj(x),x))
+        } else {
+            return(crossprod(Conj(x),y))
+        }
+    } else {
+        return(crossprod(x,y))
+    }
+}
+
+`tcprod` <- function(x,y=NULL){
+    if(is.complex(x) | is.complex(y)){
+        if(is.null(y)){
+            return(tcrossprod(x,Conj(x)))
+        } else {
+            return(tcrossprod(x,Conj(y)))
+        }
+    } else {
+        return(tcrossprod(x,y))
+    }
+}
+
+`quad.form.chol` <- function (chol, x){
+        jj <- cprod(chol, x)
+        drop(cprod(jj, jj))
+}
+
+`quad.form` <- function (M, x){ drop(crossprod(crossprod(M,Conj(x)),x)) }
+
+`quad.form.inv` <- function (M, x){ drop(cprod(x, solve(M, x))) }
+
+`quad.3form` <- function(M,left,right){ crossprod(crossprod(M, Conj(left)), right) }
+
+`quad.3form.inv` <- function(M,left,right){ drop(cprod(left, solve(M, right))) }
+
+`quad.3tform` <- function(M,left,right){ tcrossprod(left, tcrossprod(Conj(right), M)) }
+
+`quad.tform` <- function(M,x){ tcrossprod(x, tcrossprod(Conj(x), M)) }
+
+`quad.tform.inv` <- function(M,x){ drop(quad.form.inv(M, ht(x))) }
+
+`quad.diag` <- function(M,x){ colSums(crossprod(M, Conj(x)) * x) }
+
+`quad.tdiag` <- function(M,x){ rowSums(tcrossprod(Conj(x), M) * x) }
+
+`quad.3diag` <- function(M,left,right){ colSums(crossprod(M, Conj(left)) * right) }
+
+`quad.3tdiag` <- function(M,left,right){ colSums(t(left) * tcprod(M, right)) }
+
+cp <- cprod
+tcp <- tcprod
+qf <- quad.form
+qfi <- quad.form.inv
+q3 <- quad.3form
+q3i <- quad.3form.inv
+
+q3t <- quad.3tform
+qt <- quad.tform
+q3i <- quad.tform.inv
+qd <- quad.diag
+qtd <- quad.tdiag
+q3d <- quad.3diag
+q3td <- quad.3tdiag

README.md

@@ -0,0 +1,54 @@
+Quadratic forms in R: the `quadform` package
+================
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+# <img src="man/figures/quadform.png" width = "150" align="right" />
+
+<!-- badges: start -->
+
+[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/quadform)](https://cran.r-project.org/package=quadform)
+<!-- badges: end -->
+
+Quadratic forms are polynomials with all terms of degree 2. Given a
+column vector ${\mathbf x}=(x_1,\ldots,x_n)^\top$ and an $n\times n$
+matrix $M$ then the function
+$f\colon\mathbb{R}^n\longrightarrow\mathbb{R}$ given by
+$f({\mathbf x})=x^TMx$ is a quadratic form; we extend to complex vectors
+by mapping ${\mathbf z}=(z_1,\ldots, z_n)^\top$ to
+${\mathbf z}^*M{\mathbf z}$, where $z^*$ means the complex conjugate of
+$z^T$. These are implemented in the package with `quad.form(M,x)` which
+is essentially
+
+`quad.form <- function(M,x){crossprod(crossprod(M, Conj(x)), x)}.`
+
+This is preferable to `t(x) %*% M %*% x` on several grounds. Firstly, it
+streamlines and simplifies code; secondly, it is more efficient; and
+thirdly it handles the complex case consistently. The package includes
+similar functionality for other related expressions.
+
+The main motivation for the package is nicer code. For example, the
+`emulator` package has to manipulate the following expression:
+
+$$
+\left[H_x-H^\top A^{-1}U\right]^\top
+\left[H^\top\left(H^\top A^{-1}H\right)^{-1}H\right]
+\left[H_x-H^\top A^{-1}U\right].
+$$
+
+Direct R idiom would be:
+
+    t(Hx - t(H) %*% solve(A) %*% U) %*%  t(H) %*% solve(t(H) %*% solve(A) %*% H) %*% H %*%  (Hx - t(H) %*% solve(A) %*% U)
+
+But `quadform` idiom is:
+
+    quad.form(quad.form.inv(quad.form.inv(A,H),H), Hx - quad.3form.inv(A,H,U))
+
+and in terse form becomes:
+
+    qf(qfi(qfi(A,H),H), Hx - q3fi(A,H,U))
+
+which is certainly shorter, arguably more elegant, and possibly faster.
+
+The package is maintained on
+[github](https://github.com/RobinHankin/quadform).

inst/quad_form_test.Rmd

@@ -0,0 +1,214 @@
+---
+title: "Quadratic form testing vignette"
+author: "Robin K. S. Hankin"
+vignette: >
+  %\VignetteIndexEntry{quadratic form tests}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r set-options, echo = FALSE}
+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)
+```
+
+```{r out.width='15%', out.extra='style="float:right; padding:10px"',echo=FALSE}
+knitr::include_graphics(system.file("help/figures/emulator.png", package = "emulator"))
+```
+
+
+# Testing for `quad.form()` et seq
+
+In 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:
+
+```{r}
+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:
+
+```{r}
+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:
+
+```{r}
+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:
+
+## Test of `ht(x)` = $x^*$ = $\overline{x'}$ (Hermitian transpose):
+
+### ```ht(x)=t(Conj(x))```
+
+```{r}
+(jj1 <- Conj(t(x)))
+(jj2 <- t(Conj(x)))
+(jj3 <- ht(x))
+tester(jj1,jj3)
+tester(jj2,jj3)
+```
+
+## Test of `cprod()` = $x^*y$:
+
+### `cprod(x,y)=crossprod(Conj(x),y)`
+
+```{r}
+(jj1 <- ht(x) %*% x1)
+(jj2 <- cprod(x,x1))
+tester(jj1,jj2)
+```
+
+## Test of `tcprod()` = $x y^*$:
+
+### `tcprod(x,y)=crossprod(x,Conj(y))`
+
+```{r}
+(jj1 <- ht(x1) %*% x)
+(jj2 <- cprod(x1,x))
+tester(jj1,jj2)
+```
+
+## Test of `quad.form()` = $x^*Mx$:
+
+### `quad.form(M,x)=crossprod(crossprod(M,Conj(x)),x))`
+
+```{r}
+(jj1 <- ht(x) %*% M %*% x)
+(jj2 <- quad.form(M,x))
+tester(jj1,jj2)
+```
+
+##  Test of `quad.form.inv()` = $x^*M^{-1}x$:
+
+### `quad.form.inv(M,x)=cprod(x,solve(M,x))`
+
+
+```{r}
+(jj1 <- ht(x) %*% solve(M) %*% x)
+(jj2 <- quad.form(solve(M),x))
+max(abs(jj1-jj2))
+```
+
+##  Test of `quad.3form()` = $x^*My$:
+
+### `quad.3form(M,l,r)=crossprod(crossprod(M,Conj(l)),r)`
+
+
+```{r}
+(jj1 <- ht(x) %*% M %*% x1)
+(jj2 <- quad.3form(M,x,x1))
+tester(jj1,jj2)
+```
+
+## Test of `quad.3tform()` = $xMy^*$:
+
+### `quad.3tform(M,l,r)=tcrossprod(left,tcrossprod(Conj(right),M))`
+
+
+```{r}
+(jj1 <- y %*% M %*% ht(y1))
+(jj2 <- quad.3tform(M,y,y1))
+tester(jj1,jj2)
+```
+
+## Test of `quad.tform()` = $xMx^*$:
+
+### `quad.tform(M,x)=tcrossprod(x,tcrossprod(Conj(x),M))`
+
+```{r}
+(jj1 <- y %*% M %*% ht(y))
+(jj2 <- quad.tform(M,y))
+tester(jj1,jj2)
+```
+
+
+## Test of `quad.tform.inv()` = $xM^{-1}x^*$:
+
+### `quad.tform.inv(M,x)=quad.form.inv(M,ht(x))`
+
+```{r}
+(jj1 <- y %*% solve(M) %*% ht(y))
+(jj2 <- quad.tform.inv(M,y))
+tester(jj1,jj2)
+```
+
+## Test of `quad.diag()` = $\operatorname{diag}(x^*Mx)$ = `diag(quad.form())`:
+
+### `quad.diag(M,x)=colSums(crossprod(M,Conj(x)) * x)`
+
+```{r}
+(jj1 <- diag(ht(x) %*% M %*% x))
+(jj2 <- diag(quad.form(M,x)))
+(jj3 <- quad.diag(M,x))
+tester(jj1,jj3)
+tester(jj2,jj3)
+```
+
+##  Test of `quad.tdiag()` =  $\operatorname{diag}(xMx^*)$ = `diag(quad.tform())`:
+
+### `quad.tdiag(M,x)=rowSums(tcrossprod(Conj(x), M) * x)`
+
+
+```{r}
+(jj1 <- diag(y %*% M %*% ht(y)))
+(jj2 <- diag(quad.tform(M,y)))
+(jj3 <- quad.tdiag(M,y))
+tester(jj1,jj3)
+tester(jj2,jj3)
+```
+
+## Test of `quad.3diag()` = $\operatorname{diag}(x^*My)$
+
+### `quad.3diag(M,l,r)=colSums(crossprod(M, Conj(left)) * right)`
+
+
+```{r}
+(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)
+```
+
+## Test of `quad.3tdiag()` = $\operatorname{diag}(xMy^*)$
+
+### `quad.3tdiag(M,l,r)=colSums(t(left) * tcprod(M, right))`
+
+```{r}
+(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)
+```