# R/Toeplitz.R In SuperGauss: Superfast Likelihood Inference for Stationary Gaussian Time Series

#### Documented in Toeplitz

#' @title Constructor and methods for Toeplitz matrix objects.
#'
#' @description The \code{Toeplitz} class contains efficient methods for linear algebra with symmetric positive definite (i.e., variance) Toeplitz matrices.
#'
#' @aliases setAcf getAcf traceT2 traceT4 show.Toeplitz %*% determinant solve %*%,ANY,Toeplitz-method %*%,Toeplitz,ANY-method determinant,Toeplitz-method dim,Toeplitz-method ncol,Toeplitz-method nrow,Toeplitz-method show,Toeplitz-method solve,Toeplitz-method
#' @section Methods:
#' If \code{Toep} is a \code{Toeplitz} object with first row/column given by \code{acf}, then:
#' \describe{
#' \item{\code{Toep$setAcf(acf)}}{Sets the autocorrelation of the matrix.} #' \item{\code{Toep$getAcf()}}{Gets the autocorrelation of the matrix.}
#' \item{\code{nrow(Toep)}, \code{ncol(Toep)}, \code{dim(Toep)}}{Selected dimension(s) of the matrix.}
#' \item{\code{Toep \%*\% X}, \code{X \%*\% Toep}}{Toeplitz-Matrix and Matrix-Toeplitz multiplication.  Also works if \code{X} is a vector.}
#' \item{\code{solve(Toep, X)}, \code{solve(Toep)}}{Solves Toeplitz systems of equations.  When second argument is missing, returns the inverse of the Toeplitz matrix.}
#' \item{\code{determinant(Toep)}}{Log-determinant of the Toeplitz matrix, i.e., same thing as \code{log(det(toeplitz(acf)))}.}
#' \item{\code{Toep$traceT2(acf2)}}{If \code{T1 == toeplitz(acf)} and \code{T2 == toeplitz(acf2)}, computes the trace of \code{solve(T1, T2)}. This is used in the computation of the gradient of Gaussian likelihoods with Toeplitz variance matrix.} #' \item{\code{Toep$traceT4(acf2, acf3)}}{If \code{T1 == toeplitz(acf)}, \code{T2 == toeplitz(acf2)}, and \code{T3 == toeplitz(acf3)}, computes the trace of \code{solve(T1, T2) \%*\% solve(T1, T3)}.  This is used in the computation of the Hessian of Gaussian likelihoods with Toeplitz variance matrix.}
#' }
#' @details It is assumed that the autocorrelation of the \code{Toeplitz} object defines a valid (i.e., positive definite) variance matrix.  The multiplication algorithms still work when this is not the case but the other algorithms do not (return values typically contain \code{NaN}s).
#' @examples
#' # construction
#' acf <- exp(-(1:5))
#' Toep <- Toeplitz(acf = acf)
#' # alternatively, can allocate space first
#' Toep <- Toeplitz(n = length(acf))
#' Toep$setAcf(acf = acf) #' #' dim(Toep) # == c(nrow(Toep), ncol(Toep)) #' Toep # show method #' Toep$getAcf() # extract the acf
#'
#' # linear algebra
#' X <- matrix(rnorm(10), 5, 2)
#' Toep %*% X
#' t(X) %*% Toep
#' solve(Toep, X)
#' determinant(Toep) # log-determinant
#' @export
.Toeplitz <- setRefClass("Toeplitz",
fields = list(cpp_ptr = "externalptr",
size = "numeric"))
.Toeplitz$lock("cpp_ptr") # locked fields .Toeplitz$lock("size")
# internal constructor
.Toeplitz$methods(initialize = function(n) { cpp_ptr <<- .Toeplitz_constructor(n) size <<- n }) # exported constructor #' @rdname Toeplitz-class #' @param n Size of the Toeplitz matrix. #' @param acf Autocorrelation vector of Toeplitz matrix. #' @return A \code{Toeplitz} object. #' @export Toeplitz <- function(n, acf) { if(missing(n)){ n <- length(acf) } Tz <- .Toeplitz$new(n)
if(!missing(acf)) {
Tz$setAcf(acf) } Tz } #--- custom methods ------------------------------------------------------------ # setter .Toeplitz$methods(setAcf = function(acf) {
if(length(acf) != size) {
stop("acf has wrong length.")
}
.Toeplitz_setAcf(cpp_ptr, acf)
})

# getter
.Toeplitz$methods(getAcf = function() { .Toeplitz_getAcf(cpp_ptr) }) # traceT2 .Toeplitz$methods(traceT2 = function(acf2) {
if(!.Toeplitz_hasAcf(cpp_ptr)) {
stop("setAcf has not been called yet")
}
if(length(acf2) != size) {
stop("acf2 has wrong length.")
}
.Toeplitz_traceT2(cpp_ptr, acf2)
})

# traceT4
# TODO: fix this so that small1 is inside C++ code
.Toeplitz$methods(traceT4 = function(acf2, acf3) { if(!.Toeplitz_hasAcf(cpp_ptr)) { stop("setAcf has not been called yet") } if(length(acf2) != size) { stop("acf2 has wrong length.") } if(length(acf3) != size) { stop("acf3 has wrong length.") } .Toeplitz_traceT4(cpp_ptr, acf2, acf3) }) ## .traceT4 <- function(acf, acf2, acf3) { ## N <- acf$size
##   ee <- c(1, rep(0, N-1))
##   small1 <- abs(acf3[1]) < .0001
##   if(small1) acf3[1] <- 1 + acf3[1]
##   T4 <- acf$traceT4(acf2, acf3) ## if(small1) T4 <- T4 - acf$traceT4(acf2, ee)
##   T4
## }

#--- generic methods -----------------------------------------------------------

# show
#' @export
setMethod("show", "Toeplitz", function(object) {
if(.Toeplitz_hasAcf(object$cpp_ptr)) { obj.acf <- object$getAcf()[1:min(6, object$size)] obj.acf <- signif(obj.acf, digits = 3) if(object$size > 6) obj.acf <- c(obj.acf, "...")
} else {
obj.acf <- "NULL"
}
cat("Toeplitz matrix of size", object$size, "\n", "acf: ", obj.acf, "\n") }) # ncol #' @export setMethod("ncol", "Toeplitz", function(x){ x$size
})

# nrow
#' @export
setMethod("nrow", "Toeplitz", function(x){
x$size }) # dim #' @export setMethod("dim", "Toeplitz", function(x){ rep(x$size, 2)
})

# Toeplitz-Matrix multiplication
#' @export
setMethod("%*%", signature(x = "Toeplitz", y = "ANY"), function(x, y) {
if(!.Toeplitz_hasAcf(x$cpp_ptr)) { stop("setAcf has not been called yet") } if(is.vector(y)) y <- as.matrix(y) if(!is.matrix(y)) { stop("Second argument should be a matrix or vector.") } if(nrow(y) != x$size) {
stop("Toeplitz and second argument are non-conformable.")
}
.Toeplitz_Multiply(x$cpp_ptr, y) }) # Matrix-Toeplitz multiplication setMethod("%*%", signature(x = "ANY", y = "Toeplitz"), function(x, y) { if(!.Toeplitz_hasAcf(y$cpp_ptr)) {
stop("setAcf has not been called yet")
}
if(is.vector(x)) x <- as.matrix(x)
if(!is.matrix(x)) {
stop("First argument should be a matrix or vector.")
}
if(ncol(x) != y$size) { stop("First argument and Toeplitz are non-conformable.") } t(.Toeplitz_Multiply(y$cpp_ptr, t(x)))
})

# determinant
#' @export
setMethod("determinant", "Toeplitz",
function(x, logarithm = TRUE, ...) {
if(!.Toeplitz_hasAcf(x$cpp_ptr)) { stop("setAcf has not been called yet") } ldT <- .Toeplitz_Determinant(x$cpp_ptr)
if(!logarithm) {
ldT <- exp(ldT)
}
ldT
})

# solve
#' @export
setMethod("solve", "Toeplitz", function(a, b, ...) {
if(!.Toeplitz_hasAcf(a$cpp_ptr)) { stop("setAcf has not been called yet") } if(missing(b)) b <- diag(a$size)
if(is.vector(b)) b <- as.matrix(b)
if(!is.matrix(b)) {
stop("b must be a matrix or vector.")
}
if(nrow(b) != a$size) { stop("a and b are non-conformable.") } .Toeplitz_Solve(a$cpp_ptr, b)
})

# now make methods display arguments
.DollarNames.Toeplitz <- function(x, pattern)
grep(pattern, getRefClass(class(x))\$methods(), value=TRUE)


## Try the SuperGauss package in your browser

Any scripts or data that you put into this service are public.

SuperGauss documentation built on May 1, 2019, 7:58 p.m.