Movatterモバイル変換

Type:	Package
Title:	SPArse Matrix
Version:	2.11-1
Date:	2025-01-20
Depends:	R (≥ 4.0)
Imports:	dotCall64, grid, methods, Rcpp (≥ 1.0.8.3)
LinkingTo:	Rcpp
Suggests:	spam64, fields, Matrix, testthat, R.rsp, truncdist, knitr,rmarkdown
VignetteBuilder:	R.rsp, knitr
Description:	Set of functions for sparse matrix algebra. Differences with other sparse matrix packages are: (1) we only support (essentially) one sparse matrix format, (2) based on transparent and simple structure(s), (3) tailored for MCMC calculations within G(M)RF. (4) and it is fast and scalable (with the extension package spam64). Documentation about 'spam' is provided by vignettes included in this package, see also Furrer and Sain (2010) <doi:10.18637/jss.v036.i10>; see 'citation("spam")' for details.
LazyData:	true
License:	LGPL-2 |BSD_3_clause + file LICENSE
URL:	https://www.math.uzh.ch/pages/spam/
BugReports:	https://git.math.uzh.ch/reinhard.furrer/spam/-/issues
NeedsCompilation:	yes
Packaged:	2025-01-20 19:24:40 UTC; furrer
Author:	Reinhard Furrer [aut, cre], Florian Gerber [aut], Roman Flury [aut], Daniel Gerber [ctb], Kaspar Moesinger [ctb], Annina Cincera [ctb], Youcef Saad [ctb] (SPARSEKIT http://www-users.cs.umn.edu/~saad/software/SPARSKIT/), Esmond G. Ng [ctb] (Fortran Cholesky routines), Barry W. Peyton [ctb] (Fortran Cholesky routines), Joseph W.H. Liu [ctb] (Fortran Cholesky routines), Alan D. George [ctb] (Fortran Cholesky routines), Lehoucq B. Rich [ctb] (ARPACK), Maschhoff Kristi [ctb] (ARPACK), Sorensen C. Danny [ctb] (ARPACK), Yang Chao [ctb] (ARPACK)
Maintainer:	Reinhard Furrer <reinhard.furrer@math.uzh.ch>
Repository:	CRAN
Date/Publication:	2025-01-20 23:50:02 UTC

SPArse Matrix Package

Description

spam is a collection of functions for sparse matrixalgebra.

Gereral overview

What is spam and what is it not:

WhileMatrix seems an overshoot of classes andSparseMfocuses mainly on regression type problem, we provide a minimal set ofsparse matrix functions fully functional for everyday spatial statisticslife. There is however some emphasize on Markov chain Monte Carlo typecalculations within the framework of (Gaussian) Markov random fields.

Emphasis is given on a comprehensive, simple, tutorial structure of thecode. The code is S4 based but (in a tutorial spirit) the functions arein a S3 structure visible to the user (exported viaNAMESPACE).

There exist many methods for sparse matrices that work identically as inthe case of ordinary matrices. All the methods are discussed in the helpand can be accessed directly via a*.spam concatenation to thefunction. For example,help(chol.spam) calls the helpdirectly. We deliberately avoided aliases according to analogue helps from thebase package.

Sparseness is used when handling large matrices. Hence, care has beenused to provide efficient and fast routines. Essentially, the functionsdo not transform the sparse structure into full matrices to use standard(available) functionality, followed by a back transform. We agree, moreoperators, functions, etc. should eventually be implemented.

The packagesfields andspam are closely linked.

Author(s)

Reinhard Furrer, with the help of Florian Gerber, KasparMoesinger and many others.
Some Fortran routines were written by Youcef Saad, Esmond G. Ng, Barry W. Peyton, Joseph W.H. Liu, Alan D. George.

References

Reinhard Furrer, Stephan R. Sain (2010)."spam: A Sparse Matrix R Package with Emphasis on MCMCMethods for Gaussian Markov Random Fields.",Journal of Statistical Software, 36(10), 1-25,doi:10.18637/jss.v036.i10.
Florian Gerber, Reinhard Furrer (2015)."Pitfalls in the Implementation of Bayesian Hierarchical Modeling of Areal Count Data: An Illustration Using BYM and Leroux Models.",Journal of Statistical Software, Code Snippets, 63(1), 1-32,doi:10.18637/jss.v063.c01.
F. Gerber, K. Moesinger, R. Furrer (2017),"Extending R packages to support 64-bit compiled code: An illustration with spam64 and GIMMS NDVI3g data.",Computer & Geoscience 104, 109-119,doi:10.1016/j.cageo.2016.11.015."

See Also

Seespam.class for a detailed class description,spam andspam.ops for creation,coercion and algebraic operations.options.

Examples

## Citations:citation('spam')citation('spam', auto=TRUE)## History of changes## Not run: file.show(system.file("NEWS.md", package = "spam"))## End(Not run)

Coercion to a Vector

Description

Coercion ofspam matrices to proper vector objects

Usage

## S4 method for signature 'spam'as.vector(x, mode = "any")

Arguments

Details

This coercion allows smooth transitions between differentmatrix formats, see example below.
The Cholesky factors are first transformed to aspam object.

Value

Ifstructurebased=TRUE, the vectorx@entries.
Conversely, ifstructurebased=FALSE, the result is identical toone withas.vector(as.matrix(x)).

Author(s)

Reinhard Furrer

See Also

spam.options

Examples

x <- diag(2)ifelse( x, x, 1-x)ifelse( x, as.vector(x), 1-as.vector(x))x <- diag.spam(2)options(spam.structurebased=FALSE)ifelse( x, as.vector(x), 1-as.vector(x))options(spam.structurebased=TRUE)ifelse( x, as.vector(x), 1-as.vector(x))

Mathematical Functions

Description

Applies theMath group functions tospam objects

Usage

# ceiling(x)# floor(x)# exp(x, base = exp(1))# log(x, base = exp(1))# sqrt(x)# abs(x)# cumprod(x)# cumsum(x)# cos(x)# sin(x)# tan(x)# acosh(x)

Arguments

Details

It is important to note that the zero entries do not enter theevaluation whenstructurebased=FALSE. The operations are performed on the stored non-zeroelements. This may lead to differences if compared with the sameoperation on a full matrix.

Value

If the optionspam.structurebased=TRUE, all functions operate on the vectorx@entries and return theresult thereof.
Conversely, ifstructurebased=FALSE, the result is identical toone withas.matrix(x) input and anas.spam purger.

Author(s)

Reinhard Furrer

See Also

Summary.spam,Ops.spam andMath2.spam

Examples

getGroupMembers("Math")mat <- matrix(c( 1,2,0,3,0,0,0,4,5),3)smat <- as.spam( mat)cos( mat)cos( smat)options(spam.structurebased=FALSE)cos( smat)sqrt( smat)

Rounding of Numbers

Description

Applies theMath2 group functions to 'spam' objects

Usage

## S4 method for signature 'spam'round(x, digits = 0)## S4 method for signature 'spam'signif(x, digits = 6)

Arguments

Value

All functions operate on the vectorx@entries and return theresult thereof.

Author(s)

Reinhard Furrer

See Also

Ops.spam andMath.spam

Examples

getGroupMembers("Math2")set.seed(12)smat <- diag.spam( rnorm(15))round(smat, 3)

Oral Cavity Cancer

Description

Oral cavity cancer counts in 544 districts in Germany over 1986-1990.

Format

Oral is a dataframe with 3 columns.

Y

observed counts

E

expected counts

SMR

standardized mortality ratios

germany is alist of 544 elements, each describing an individual polygon of the district.

Details

The expected counts depend on the number of people in theregion and their age distribution.
The regions are ordered according the supplied polygon description andadjacency graph.

There is a similar datasetdata(Germany) with larynx cancer cases from the packageINLA as well, with an additional smoking covariate.

Source

The data is available from the packageINLAdistributed fromhttps://www.r-inla.org.

References

Knorr-Held, L. and Rasser, G. (2000) Bayesian Detection of Clusters and Discontinuities in Disease Maps,Biometrics,56,13–21.

See Also

germany.plot.

Rounding of Numbers

Description

Applies theMath2 group functions tospam objects

Usage

# max(x,..., na.rm = FALSE)

Arguments

Details

Thena.rm argument is only meaninful ifNAOK=TRUE.

Value

Ifstructurebased=TRUE, all functions operate on the vectorx@entries and return theresult thereof.
Conversely, ifstructurebased=FALSE, the result is identical toone withas.matrix(x) input.

Author(s)

Reinhard Furrer

See Also

Math.spam andMath2.

Examples

getGroupMembers("Summary")smat <- diag.spam( runif(15))range(smat)options(spam.structurebased=FALSE)range(smat)## Not run: max( log(spam(c(1,-1))), na.rm=TRUE)## End(Not run)# allow 'NA's first:# TODO# options(spam.NAOK=TRUE)# max( log(spam(c(1,-1))), na.rm=TRUE)

Adjacency Structure of the Counties in the Contiguous United States

Description

First and second order adjacency structure of the countiesin the contiguous United States. We consider that two countiesare neighbors if they share atleast one edge of their polygon description inmaps.

Format

Two matrices of classspam

UScounties.storder

Contains a one in thei andj element if countyi is a neighbor of countyj.

UScounties.ndorder

Contains a one in thei andj element if countiesi andj are a neighbors of countyk and countiesi andj are not neighbors.

See Also

map, frommaps.

Examples

# number of counties:n  <- nrow( UScounties.storder)## Not run: # make a precision matrixQ <- diag.spam( n) + .2 * UScounties.storder + .1 * UScounties.ndorderdisplay( as.spam( chol( Q)))## End(Not run)

Monthly Total Precipitation (mm) for April 1948 in the Contiguous United States

Description

This is a useful spatial data set of moderate to large size consisting of 11918locations. Seehttps://www.image.ucar.edu/GSP/Data/US.monthly.met/ for the source of these data.

Format

This data set is an array containing the following columns:

lon,lat

Longitude-latitude position of monitoring stations.

raw

Monthly total precipitation in millimeters for April 1948.

anomaly

Preipitation anomaly for April 1948.

infill

Indicator, which station values were observed (5906 out of the 11918)compared to which were estimated.

Source

https://www.image.ucar.edu/GSP/Data/US.monthly.met/

References

Johns, C., Nychka, D., Kittel, T., and Daly, C. (2003)Infilling sparse records of spatial fields.Journal of the American Statistical Association,98, 796–806.

See Also

RMprecip

Examples

# plot## Not run: library(fields)data(USprecip)par(mfcol=c(2,1))quilt.plot(USprecip[,1:2],USprecip[,3])US( add=TRUE, col=2, lty=2)quilt.plot(USprecip[,1:2],USprecip[,4])US( add=TRUE, col=2, lty=2)## End(Not run)

Administrative Districts of Germany

Description

Constructing the adjacency graphof the administrative districts of Germany

Usage

adjacency.landkreis(loc)

Arguments

Details

The function is included as an example on how toconstruct adjacency matrices form a (common) adjacency structure.For the particular example, note that the nodes are not numberedconsecutively and that they start from zero.

Value

a sparse matrix inspam format.

Author(s)

Reinhard Furrer

References

The adjacency data has been provided by Havard Rue and isalso available inINLA.

See Also

germany.plot super-seedingmap.landkreisfor plotting.
Oral.

Examples

## Not run: loc <- system.file("demodata/germany.adjacency", package="spam")display( adjacency.landkreis( loc))## End(Not run)

Test if Two Sparse Matrices are (Nearly) Equal

Description

Utility to compare twospam objectstesting 'near equality'. Depending on the type of difference, comparison isstill made to some extent, and a report of the differences isreturned.

Usage

## S3 method for class 'spam'all.equal(target, current, tolerance = .Machine$double.eps^0.5,    scale = NULL, check.attributes = FALSE,...)

Arguments

Details

Numerical comparisons forscale = NULL (the default) aretypically on arelative difference scale unless thetarget values are close to zero or infinite. Specifically,the scale is computed as the average absolute value oftarget.If this scale is finite and exceedstolerance, differencesare expressed relative to it; otherwise, absolute differences are used.

Ifscale is numeric (and positive), absolute comparisons aremade after scaling (dividing) byscale. Note that if all ofscale is sufficiently close to 1 (specifically, withintolerance),the difference is still reported as being on an absolute scale.

Do not useall.equal.spam directly inifexpressions: either useisTRUE( all.equal.spam(...)) oridentical if appropriate.

Cholesky decomposition routines use this function to test forsymmetry.

A method formatrix-spam objects is defined as well.

There is the additional catch of a zero matrix being represented byone zero element, see ‘Examples’ below.

Value

EitherTRUE or a vector of 'mode'"character" describing thedifferences betweentarget andcurrent.

Author(s)

Reinhard Furrer

See Also

isSymmetric.spam andcleanup.

Examples

obj <- diag.spam(2)obj[1,2] <- .Machine$double.epsall.equal( diag.spam(2), obj)all.equal( t(obj), obj)all.equal( t(obj), obj*1.1)# We can compare a spam to a matrixall.equal(diag(2),diag.spam(2))# the opposite does often not make sense,# hence, it is not implemented.all.equal(diag.spam(2),diag(2))# A zero matrix contains one element:str(spam(0))# henceall.equal.spam(spam(0,3,3), diag.spam(0,3) )norm(spam(0,3,3) - diag.spam(0,3) )

Apply Functions Over Sparse Matrix Margins

Description

Returns a vector or array or list of values obtained by applying afunction to margins of a sparse matrix.

Usage

apply.spam(X, MARGIN=NULL, FUN, ...)

Arguments

Details

This is a handy wrapper to apply a function to the (nonzero)elements of a sparse matrix. For example, it is possible to apply a covariance matrix to a distancematrix obtained bynearest.dist, see Examples.

A call toapply only coerces the sparse matrix to a regular one.

The basic principle is applying the function to@entries, or tothe extracted columns or rows ([,i,drop=F] or[i,,drop=F]). It is important to note that an empty columncontains at least one zero value and may lead to non intuitiveresults.

This function may evolve over the next few releases.

Value

Similar as a call toapply with a regular matrix. The mostimportant cases are as follows. Theresult is a vector (MARGIN is length 1 andFUN isscalar) or a matrix (MARGIN is length 1 andFUN returnsfixed length vectors, orMARGIN is length 2 andFUN isscalar) or a list (ifFUN returns vectors of different lengths).

Author(s)

Reinhard Furrer

See Also

base:apply for more details on Value.

Examples

S <- as.spam(dist(1:5))S <- apply.spam(S/2, NULL, exp)# instead of # S@entries <- exp( S@entries/2) # Technical detail, a null matrix consists# of one zero element.apply.spam(S,c(1,2),pmax)apply.spam(S,1,range)# A similar example as for the base apply.# However, no dimnames else we would get warnings. x <- as.spam(cbind(x1 = 3, x2 = c(0,0,0, 5:2)))apply.spam(x, 2, mean, trim = .2)col.sums <- apply.spam(x, 2, sum)row.sums <- apply.spam(x, 1, sum)rbind(cbind(x, row.sums), c(col.sums, sum(col.sums)))apply.spam(x, 2, is.vector)# Sort the columns of a matrix# Notice that the result is a list due to the different# lengths induced by the nonzero elementsapply.spam(x, 2, sort)# Function with extra args:cave <- function(x, c1, c2) c(mean(x[c1]), mean(x[c2]))apply(x,1, cave,  c1=1, c2=c(1,2))ma <- spam(c(1:4, 0, 0,0, 6), nrow = 2)maapply.spam(ma, 1, table)  #--> a list of length 2apply.spam(ma, 1, stats::quantile)# 5 x n matrix with rownames

Bandwidth of a Sparse Matrix

Description

Returns the lower and upper bandwidth of a sparse matrix

Usage

bandwidth(A)

Arguments

Details

The matrix does not need to be diagonal. Values can be negativeindicating the the matrix contains a band cinfined in theupper or lower triangular part.

Value

Integer vector containing the lower and upper bandwidth

Author(s)

Reinhard Furrer

See Also

diag.spam.

Examples

bandwidth(spam(c(0, 1), 3, 2))bandwidth(spam(c(0, 0, 1, rep(0, 9)), 4, 3))

Binds Arrays Corner-to-Corner

Description

Creates a sparse block-diagonal matrix.

Usage

bdiag.spam(...)

Arguments

Details

This is a small helper function to create block diagonal sparse matrices. In the two matrix case,bdiag.spam(A,B), this is equivalent to a complicatedrbind(cbind(A, null), cbind(B, t(null))),wherenull is a null matrix of appropriate dimension.

It is recursively defined.

The arrays are coerced to sparse matrices first.

This function is similar to the functionbdiag from the packageMatrix. It is also similar to the functionadiag from the packagemagic. However, here no padding is done and all the dimnames arestripped.

Value

Returns aspam matrix as described above.

Author(s)

Reinhard Furrer

See Also

diag.spam.

Examples

A <- diag.spam(2, 4)           # 2*I4B <- matrix(1,3,3)AB <- bdiag.spam(A,B)# equivalent to:ABalt <- rbind(cbind( A, matrix(0,nrow(A),ncol(B))),               cbind( matrix(0,nrow(B),ncol(A)), B))         norm(AB-ABalt)# Matrices do not need to be square:bdiag.spam(1,2:5,6)

Combine Sparse Matrices by Rows or Columns

Description

Take a sequence of vector, matrix orspam object arguments andcombine bycolumns orrows, respectively.

Usage

# cbind(\dots, force64 = getOption("spam.force64"), deparse.level = 0)# rbind(\dots, deparse.level = 0)

Arguments

Details

rbind andcbind are not exactly symmetric in howthe objects are processed.cbind calls aFortran routine after the input has been coerced tospamobjects. Whereasrbind calls a Fortran routine only in the caseof tospam matrices. Note that row binding is essentially an concatenationof the slots due to the sparse storage format.

Only two objects at a time are processed. If more than two arepresent, a loop concatenates them successively.

A method is defined for aspam object as first argument.

Value

aspam object combining the... argumentscolumn-wise or row-wise. (Exception: if there are no inputs or allthe inputs areNULL, the value isNULL.)

Author(s)

Reinhard Furrer

Examples

x <- cbind.spam(1:5,6)y <- cbind(x, 7)rbind( x, x)# for some large matrices   t( cbind( t(x), t(x)))# might be slightly faster:

Cholesky Factorization for Sparse Matrices

Description

chol performs a Choleskydecomposition of a symmetric positive definite sparse matrixxof classspam.

Usage

# chol(x, \dots)## S4 method for signature 'spam'chol(x, pivot = "MMD", method = "NgPeyton",       memory = list(), eps = getOption("spam.eps"), Rstruct=NULL,       ..., verbose=FALSE)# update.spam.chol.NgPeyton(object, x,...)## S4 method for signature 'spam.chol.NgPeyton'update(object, x,...)

Arguments

Details

chol performs a Cholesky decomposition of a symmetricpositive definite sparse matrixx of classspam. Currently, there is only the block sparse Choleskyalgorithm of Ng and Peyton (1993) implemented (method="NgPeyton").

To pivot/permute the matrix, you can choose between the multiple minimumdegree (pivot="MMD") or reverse Cuthill-Mckee (pivot="RCM")from George and Lui (1981). It is also possible to furnish a specificpermutation in which casepivot is a vector. For compatibilityreasons,pivot can also take a logical in which forFALSEno permutation is done and forTRUE is equivalent toMMD.

Often the sparsity structure is fixed and does not change, but theentries do. In those cases, we can update the Cholesky factor withupdate.spam.chol.NgPeyton by suppling a Cholesky factor and theupdated matrix. ForU <- chol(A),update(U, Anew) andchol(Anew, Rstruct=U) are equivalent.

The optioncholupdatesingular determines how singular matricesare handled byupdate. The function hands back an error("error"), a warning ("warning") or the valueNULL("null").

The Cholesky decompositions requires parameters, linked to memoryallocation. If the default values are too small the Fortran routinereturns an error toR, which allocates more space and calls the Fortranroutine again. The user can also pass better estimates of the allocationsizes tochol with the argumentmemory=list(nnzR=..., nnzcolindices=...). The minimal sizes for a fixed sparsitystructure can be obtained from asummary call, see ‘Examples’.

The output ofchol can be used withforwardsolve andbacksolve to solve a system of linear equations.

Notice that the Cholesky factorization of the packageSparseM is alsobased on the algorithm of Ng and Peyton (1993). Whereas the Choleskyroutine of the packageMatrix are based onCHOLMOD by Timothy A. Davis (C code).

Value

The function returns the Cholesky factor in an object of classspam.chol.method. Recall that the latter is the Choleskyfactor of a reordered matrixx, see alsoordering.

Note

Although the symmetric structure ofx is needed, only the upperdiagonal entries are used. By default, the code does check forsymmetry (contrarily tobase:::chol). However,depending on the matrix size, this is a time consuming test.A test is ignored ifoptions("spam.cholsymmetrycheck") is set toFALSE.

If a permutation is supplied withpivot,options("spam.cholpivotcheck") determines if the permutation istested for validity (defaults toTRUE).

Author(s)

Reinhard Furrer, based on Ng and Peyton (1993) Fortran routines

References

Ng, E. G. and Peyton, B. W. (1993) Block sparse Cholesky algorithms onadvanced uniprocessor computers,SIAM J. Sci. Comput.,14,1034–1056.

Gilbert, J. R., Ng, E. G. and Peyton, B. W. (1994) An efficientalgorithm to compute row and column counts for sparse Choleskyfactorization,SIAM J. Matrix Anal. Appl.,15,1075–1091.

George, A. and Liu, J. (1981)Computer Solution of Large Sparse Positive Definite Systems,Prentice Hall.

See Also

det.spam,solve.spam,forwardsolve.spam,backsolve.spam andordering.

Examples

# generate multivariate normals:set.seed(13)n <- 25    # dimensionN <- 1000  # sample sizeSigma <- .25^abs(outer(1:n,1:n,"-"))Sigma <- as.spam( Sigma, eps=1e-4)cholS <- chol( Sigma)# cholS is the upper triangular part of the permutated matrix Sigmaiord <- ordering(cholS, inv=TRUE)R <- as.spam(cholS)mvsample <- ( array(rnorm(N*n),c(N,n)) %*% R)[,iord]# It is often better to order the sample than the matrix# R itself.# 'mvsample' is of class 'spam'. We need to transform it to a# regular matrix, as there is no method 'var' for 'spam' (should there?).norm( var( as.matrix( mvsample)) - Sigma, type='m')norm( t(R) %*% R - Sigma)# To speed up factorizations, memory allocations can be optimized:opt <- summary(cholS)# here, some elements of Sigma may be changed...cholS <- chol( Sigma, memory=list(nnzR=opt$nnzR,nnzcolindices=opt$nnzc))

Create Circulant Matrices

Description

Creates a circulant matrix inspam format.

Usage

circulant.spam(x, n = NULL, eps = getOption("spam.eps"))

Arguments

Value

The circulant matrix inspam format.

Author(s)

Reinhard Furrer

See Also

circulant from packagemagic,toeplitz.spam

Examples

circulant.spam(c(1,.25,0,0,0))

Cleaning up sparse matrices

Description

Eliminates an zeros in a sparse matrix.

Usage

cleanup(x, eps = getOption("spam.eps"))

Arguments

Details

A sparse matrix may still contain zeros. This function (aliasedtoas.spam) filters these values.
This often causes confusion when testing such matrices for symmetryor comparing apparently equal matrices withall.equal(see ‘Examples’ below.

Author(s)

Reinhard Furrer

See Also

isSymmetric.spam andall.equal.spam.

Examples

A <- diag.spam(2)A[1,2] <- 0all.equal(A, t(A))isSymmetric.spam(A)all.equal(cleanup(A), diag.spam(2))

Force aspam Object to Belong to a Class

Description

These functions manage the relations that allow coercing aspam objectto a given class.

Methods

signature(from = "spam", to = "matrix")

this is essentially equivalent toas.matrix(object).

signature(from = "spam", to = "list")

this is essentially equivalent totriplet(object).

signature(from = "spam", to = "vector")

this is essentially equivalent toobject@entries(structurebased=TRUE) orc(object).

signature(from = "spam", to = "logical")

the entries are forced to logicals (nonzeros only in case ofstructurebased=TRUE).

signature(from = "spam", to = "integer")

the entries are forced to integers (nonzeros only in case ofstructurebased=TRUE).

Examples

ifelse( diag.spam(2)*c(0,1), TRUE, FALSE)

Complexity for Sparse Matrices

Description

A few results of computational complexities for selectedsparse algoritms inspam

Details

A Cholesky factorization of an n-matrix requires n^3/3 flops. In case ofbanded matrices (bandwidth p, p<<n) a factorization requires about 2np^2flops. Forward- and backsolves for banded matrices require essentially2np flops.

George and Liu (1981) proves that any reordering would require atleast O(n^3/2) flops for the factorization and produce at least O(nlog(n)) fill-ins for square lattices with a local neighbor hood.
They also show that algorithms based on nested dissection are optimalin the order of magnitude sense.

More to follow.

References

George, A. and Liu, J. (1981)Computer Solution of Large Sparse Positive Definite Systems,Prentice Hall.

See Also

det,solve,forwardsolve,backsolve andordering.

Slot Modification

Description

Modify slots ofspam objects

Usage

rowpointers( x) <- valuecolindices( x) <- valueentries( x) <- value

Arguments

Details

Various tests are performed. Thus much slower than directassignment.
Slotdimension should be changed throughpad ordim

Value

Modifiedspam object.

Author(s)

Reinhard Furrer

Examples

x <- diag.spam( 2)  rowpointers( x) <- c(1,1,3)# The last line is equivalent to x@rowpointers <- as.integer( c(1,1,3))

Covariance Functions

Description

Evaluate a covariance function.

Usage

covmat(h, theta, ... , type="sph")cov.exp(h, theta, ... , eps= getOption("spam.eps"))cov.sph(h, theta, ... , eps= getOption("spam.eps"))cov.nug(h, theta, ... , eps= getOption("spam.eps"))cov.wend1(h, theta, ... , eps= getOption("spam.eps"))cov.wend2(h, theta, ... , eps= getOption("spam.eps"))cov.wu1(h, theta, ... , eps= getOption("spam.eps"))cov.wu2(h, theta, ... , eps= getOption("spam.eps"))cov.wu3(h, theta, ... , eps= getOption("spam.eps"))cov.mat(h, theta, ... , eps= getOption("spam.eps"))cov.finnmat(h, theta, ... , eps= getOption("spam.eps"))cov.mat12(h, theta, ... , eps= getOption("spam.eps"))cov.mat32(h, theta, ... , eps= getOption("spam.eps"))cov.mat52(h, theta, ... , eps= getOption("spam.eps"))cor.sph(h, range, ... , eps= getOption("spam.eps"))

Arguments

Details

covmat is a wrapper that calls the other functionsaccording to the argumenttype. The nomenclature is similar toprecmat.
The parametrization is (range, [partial-sill = 1], [smoothness = 1], [nugget = 0]), whereonly the range needs to be specified.In case of negative parameter values, a warning is issued and theabsolute value is retained.Although more cryptic, having all arguments as a single vectorsimplifies optimization withoptim.
The parameters are and locations are up to precisionepsilon.That means that all distances smaller thaneps are considered zero; a nugget smaller thaneps is ignored; a range smaller thaneps represents a nugget model; etc.
cov.finnmat() is similar tocov.mat() but with thesqrt(8*smoothness)/range argument in the Bessel function(instead of1/range).cov.mat12() is a wrapper tocov.exp()cov.mat32(), andcov.mat52() are fast version ofcov.mat() with smoothness 3/2 and 5/2, respectively (factor 10).
cor.sph(,range) is a fast version ofcov.sph(,c(range,1,0)).
Currently, the functions distinguish between a sparsespamobjecth and any other numeric type. In the future, this mightchange and appropriate methods will be implemented.

Value

Covariance function evaluated onh.

Author(s)

Reinhard Furrer

References

Any classical book about geostatistics.

See Also

precmat.

Examples

set.seed(123)n <- 200locs <- cbind(runif(n),runif(n))h <- nearest.dist(locs, delta=sqrt(2), upper = NULL)Sigma <- cov.sph(h, c(.3, 1, .1))iidsample <- rnorm(n)cholS <- chol.spam(as.spam(Sigma))iorder <- iord <- ordering(cholS, inv = TRUE)sample <- (iidsample %*% as.spam(cholS))[iorder]plot(locs, col = fields::tim.colors(n = 256)[cut(sample, n)], pch = 20)## Not run: h <- seq(0, to=1, length.out=100)plot( h, cov.exp(h, c(1/3,1)), type='l', ylim=c(0,1))type <- c("sph","wendland1","wendland2","wu1","wu2","wu3")for (i in 1:6)  lines( h, covmat(h, 1, type=type[i]), col=i+1)legend('topright',legend=type, col=2:7, lty=1)## End(Not run)

Spam Matrix Crossproduct

Description

Given matricesx andy as arguments, return a matrixcross-product. This is formally equivalent to (but usuallyslightly faster than) the callt(x) %*% y (crossprod.spam) orx %*% t(y) (tcrossprod.spam).

Usage

     crossprod.spam(x, y = NULL, ...)          tcrossprod.spam(x, y = NULL, ...)

Arguments

Value

A double matrix

Note

Whenx ory are not matrices, they are treated as column orrow matrices.

Author(s)

Reinhard Furrer

Examples

crossprod.spam(diag.spam(2),1:2)

Determinant of a Symmetric Positive Definite Sparse Matrix

Description

det anddeterminant calculate the determinant of asymmetric, positive definite sparse matrix.determinant returnsseparately the modulus of the determinant, optionally on the logarithm scale,and the sign of the determinant.

Usage

det(x, ...)determinant(x, logarithm = TRUE, ...)

Arguments

Details

If the matrix is not positive definite, the function issues awarning and returnsNA.

The determinant is based on the product of the diagonal entries of aCholesky factor, i.e. internally, a Cholesky decomposition isperformed. By default, the NgPeyton algorithm with minimal degreeordering us used. To change the methods or supply additonal parametersto the Cholesky factorization function, it is possible to pass viachol.

The determinant of a Cholesky factor is also defined.

Value

Fordet, the determinant ofx. Fordeterminant, alist with components

Author(s)

Reinhard Furrer

References

Ng, E. G. and B. W. Peyton (1993) Block sparse Cholesky algorithmson advanced uniprocessor computers,SIAM J. Sci. Comput.,14,1034–1056.

See Also

chol.spam

Examples

x <- spam( c(4,3,0,3,5,1,0,1,4), 3)det( x)determinant( x)det( chol( x))

Sparse Matrix diagonals

Description

Extract or replace the diagonal of a matrix, or construct adiagonal matrix.

Usage

## diag(x)## diag(x=1, nrow, ncol, names = TRUE)diag(x) <- valuediag.spam(x=1, nrow, ncol)spam_diag(x=1, nrow, ncol)diag.spam(x) <- value

Arguments

Details

Usingdiag(x) can have unexpected effects ifx is a vectorthat could be of length one. Usediag(x, nrow = length(x)) forconsistent behaviour.

Value

Ifx is a spam matrix thendiag(x) returns the diagonal ofx.

The assignment form sets the diagonal of the sparse matrixx to thegiven value(s).

diag.spam works asdiag for spam matrices:Ifx is a vector (or 1D array) of length two or more, thendiag.spam(x) returns a diagonal matrix whose diagonal isx.spam_diag is an alias fordiag.spam and in the spiritof the result ofdiag is aspam object.

Ifx is a vector of length one thendiag.spam(x) returns anidentity matrix of order the nearest integer tox. Thedimension of the returned matrix can be specified bynrow andncol (the default is square).

The assignment form sets the diagonal of the matrixx to thegiven value(s).

Author(s)

Reinhard Furrer

See Also

upper.tri,lower.tri.

Examples

diag.spam(2, 4)           # 2*I4smat <- diag.spam(1:5)diag( smat)diag( smat) <- 5:1# The last line is equivalent to diag.spam( smat) <- 5:1# Note that diag.spam( 1:5) <- 5:1 not work of course.

Lagged Differences

Description

Returns suitably lagged and iterated differences.

Usage

# diff.spam(x, lag = 1, differences = 1, ...)## S4 method for signature 'spam'diff(x, lag = 1, differences = 1, ...)

Arguments

Value

Aspam matrix with elements similar toas.spam(diff(as.matrix(x), ...)).

Author(s)

Reinhard Furrer

See Also

diff inbase,precmat.

Examples

# incidence matrix for a RW(3) modelD <- diff.spam(diag.spam(10), lag=1, differences=3)t(D)%*%D

Dimensions of an Object

Description

Retrieve or set the dimension of anspam object.

Usage

# dim(x)# dim(x) <- value

Arguments

Details

In older version ofspam, the behavior of the replacementmethod was different and is now implemented inpad.spam.

Value

dim retrievesthedimension slot of the object. It is a vectorof modeinteger.

The replacemnt method changes the dimension of the object by rearranging.

Author(s)

Reinhard Furrer

See Also

pad.spam.

Examples

x <- diag(4)dim(x)<-c(2,8)xs <- diag.spam(4)dim(s) <- c(2,8)  # result is different than xs <- diag.spam(4)pad(s) <- c(7,3)  # any positive value can be used

Graphially Represent the Nonzero Entries

Description

The function represents the nonzero entries in a simplebicolor plot.

Usage

display(x, ...)

Arguments

Details

spam.getOption("imagesize") determines if the sparse matrix iscoerced into a matrix and the plotted withimage.default or ifthe matrix is simply represented as a scatterplot withpch=".". The points are scaled according tocex*getOption("spam.cex")/(nrow + ncol).For some devices or for non-square matrices,cex needs probably some adjustment.

Author(s)

Reinhard Furrer

See Also

image,spam.options

Examples

set.seed(13)smat <- spam_random(8)par(mfcol=c(1,2), pty='s')options(spam.imagesize = 1000)display(smat)options(spam.imagesize = 10)display(smat, cex=.25)# very large but very sparse matrixsmat <- spam_random(2^14, distribution=rnorm, density=1e-5, verbose=TRUE)par(mfcol=c(1, 1), mai=c(.4,.4,.1,.1), pty='s')display(smat)

Eigenvalues for Sparse Matrices

Description

Functions to calculate eigenvalues and eigenvectors ofsparse matrices.It uses the value ofspam.options("inefficiencywarning") to dispatch betweenbase::eigen() or the Implicitly Restarted Arnoldi Process, using 'ARPACK'.

eigen.spam is a wrapper function ofeigen_approx and transforms its output tobase::eigen like.

Usage

eigen.spam(x, nev = 10, symmetric, only.values = FALSE, control = list())eigen_approx(x, nev, ncv, nitr, mode, only.values = FALSE, verbose = FALSE, f_routine)

Arguments

Details

mode = " ":

there are different modes available for this function, each mode returns a different range of eigenvalues.Also the available modes are dependent, whether the input matrix is symmetric or not:

"LM":

Eigenvalues with largest magnitude (sym, non sym), that is, largest eigenvalues in the Euclidean norm of complex numbers.

"SM":

Eigenvalues with smallest magnitude (sym, non sym), that is, smallest eigenvalues in the Euclidean norm of complex numbers.

"LR":

Eigenvalues with largest real part (non sym).

"SR":

Eigenvalues with smallest real part (non sym).

"LI":

Eigenvalues with largest imaginary part (non sym).

"SI":

Eigenvalues with smallest imaginary part (non sym).

"LA":

Eigenvalues with largest algebraic value (sym), that is, largest eigenvalues inclusive of any negative sign.

"SA":

Eigenvalues with smallest algebraic value (syn), that is, smallest eigenvalues inclusive of any negative sign.

ncv:

the largest number of basis vectors that will be used in the Implicitly Restarted Arnoldi Process.Work per major iteration is proportional to x@dimension[1]*ncv*ncv.The default is set ifsymmetric to min(x@dimension[1] + 1, max(2 * nev + 1, 200)) or else to min(x@dimension[1] - 1, max(2 * nev + 1, 100)).Note, this value should not be chosen arbitrary large, but slightly larger thannev.Otherwise it could lead to memory allocation problems.

nitr:

the maximum number of iterations.The default is set toncv + 1000

spamflag = FALSE:

if TRUE, the Implicitly Restarted Arnoldi Process is used,independent of the dimension of the respective matrix (provided matrixis larger than 10x10).

verbose = FALSE:

print additional information.

cmplxeps:

threshold to determine whether a double value is zero, while transforming the ARPACK output to R class complex.The default is set to.Machine$double.eps.

Value

A vector of the length corresponding to the dimension of the input matrix.Containing the requirednev eigenvalues.If requested also the corresponding eigenvectors.In the non symmetric case, the eigenvalues are returned in a matrix with a column containing the real parts and a column containing the imaginary parts of the eigenvalues.The eigenvectors are then returned in two matrices.

Note

The user is advised to choose thecontrol options carefully, see ‘Details’ for more information.

Author(s)

Roman Flury, Reinhard Furrer

References

Lehoucq, R. B. and Sorensen, D. C. and Yang, C. (1997)ARPACK Users Guide: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods.

See Also

Option"inefficiencywarning" inspam.options andspam_random.

Examples

set.seed(81)rspam <- spam_random(42^2, density = .0001, spd = TRUE)SPD <- eigen.spam(rspam, nev = 18, control = list(mode = "SM"),                  only.values = TRUE)any(SPD$values <= 0, na.rm = TRUE)isSymmetric(rspam)# hence the matrix is symmetric positiv definitrspam2 <- spam_random(50^2, density = .0001, spd = FALSE, sym = TRUE,                      distribution = rpois, lambda = 2)SNPD <- eigen.spam(rspam2, nev = 18, control = list(mode = "SM"),                    only.values = TRUE)any(SNPD$values <= 0, na.rm = TRUE)isSymmetric(rspam2)# hence the matrix is symmetric but not positiv definit

Wrapper for Distance Matrix Computation

Description

These functions are simple wrappers tonearest.distto be used infields.

Usage

spam_rdist( x1, x2, delta = 1)spam_rdist.earth( x1, x2, delta = 1, miles=TRUE, R=NULL)

Arguments

Details

These functions are wrappers tordist andrdist.earth infields. They are usedto simplify the use of sparse matrices in functions likemKrig.

For great circle distance, the matricesx1 andx2contain the degrees longitudes in the first and the degrees latitudesin the second column.delta is in degrees. Hence to restrictto distances smaller thandelta.km, one has to specifydelta=delta.km*360/(6378.388*2*pi).

Value

Aspam object containing the distances spanned betweenzero anddelta. The sparse matrix may contain many zeros(e.g., collocated data). However, to calculate covariances, these zerosare essential.

Author(s)

Reinhard Furrer

See Also

nearest.dist

Examples

## Not run: require(fields)look <- mKrig(x,Y, Covariance="Wendland", dimension=2, k=1,    cov.args=list( Distance='spam_rdist'))## End(Not run)

Transformation to Other Sparse Formats

Description

Transform between thespam sparse format to thematrix.csr format ofSparseM anddgRMatrix format ofMatrix

Usage

as.spam.matrix.csr(x)as.dgRMatrix.spam(x)as.dgCMatrix.spam(x)as.spam.dgRMatrix(x)as.spam.dgCMatrix(x)

Arguments

Details

We do not provide anyS4 methods and because of the existingmechanism a standardS3 does not work.

The functions are based onrequire.

Notice thatas.matrix.csr.spam should read asas."matrix.csr".spam.

Value

According to the call, a sparse matrix of classspam,matrix.csr,dgRMatrix ordgCMatrix.

Author(s)

Reinhard Furrer

See Also

triplet,Matrix ormatrix.csr from packageSparseM.

Examples

## Not run: S <- diag.spam(4)R <- as.dgRMatrix.spam( S)C <- as.dgCMatrix.spam( S)as.spam.dgCMatrix(C)slotNames(C)slotNames(R)# For column oriented sparse formats a transpose does not the job,# as the slot names change.# as.spam(R) does not work.## End(Not run)## Not run: # for transformations between SparseM and spam:as.matrix.csr.spam <- function(x,...) {  if (new("matrix.csr")) {    newx <- new("matrix.csr")    slot(newx,"ra",check=FALSE) <- x@entries    slot(newx,"ja",check=FALSE) <- x@colindices    slot(newx,"ia",check=FALSE) <- x@rowpointers    slot(newx,"dimension",check=FALSE) <- x@dimension    return(newx)    }  stop("function requires 'SparseM' package")}# then with `SparseM`:  as.matrix.csr.spam( spamobject )## End(Not run)## Not run: # a dataset contained in Matrixdata(KNex, package='Matrix')summary( KN <- as.spam.dgCMatrix(KNex$mm) )## End(Not run)

Meta-data About Administrative Districts of Germany

Description

Supplementary data used for the display of datafrom the administrative districts of Germany

Format

germany.info is a list with elements

n

544 (number of districts around 1990).

xrep,yrep

representative coordinates of the districts (vectorsof length 544)

xlim,ylim

2-vectors defining the limits of the districts.

polyid

linking the polygons to the districts (599 vector).

id

linking the districts to Community Identification Number.

germany.poly defines the polygons. It is a 17965 by two matrix,each polygon separated by a row ofNAs, each district by two rows.
germany defines the polygons in form of a list (backwards compatibility).

Details

The representative coordinates are calculated based on the meanvalue of the polygon coordinates. This creates sometimes strangevalues, e.g., district Leer.

Author(s)

Reinhard Furrer

References

The meta-data has been constructed based on (essentially)files from the packageINLA, seedemo(Bym).

See alsohttps://de.wikipedia.org/wiki/Amtlicher_Gemeindeschl%C3%BCsselandhttps://en.wikipedia.org/wiki/Districts_of_Germany

See Also

germany.plotOral.

Examples

# Plot the Bundeslaender:germany.plot(germany.info$id%/%1000,col=rep(2:8,3), legend=FALSE)

Plot Administrative Districts of Germany

Description

Displaying dataover the administrative districts of Germany

Usage

germany.plot(vect,  col=NULL, zlim=range(vect), legend=TRUE,              main=NULL, cex.axis=1, cex.main=1.5, add=FALSE, ... )

Arguments

Details

The legend is only added, provided (a) functionimage.plot is available.

The perfect position of the legend is an art per se and depends onvariouspar parameters. One possiblity for finer control is tonot plot it and to manually call the functionimage.plot offields.

Author(s)

Reinhard Furrer

References

See alsohttps://de.wikipedia.org/wiki/Amtlicher_Gemeindeschl%C3%BCsselandhttps://de.wikipedia.org/wiki/Liste_der_Landkreise_in_Deutschland

See Also

Oral.

Examples

data( Oral)germany.plot( Oral$Y/Oral$E)# Plot the Bundeslaender:germany.plot(germany.info$id%/%1000,col=rep(2:8,3), legend=FALSE)

Generalized Multiplication

Description

Multiplying specific submatrices of a spam matrix with different factors.

Usage

gmult(x, splits, fact)

Arguments

Value

Spam matrix, where each specified submatrix is multiplied with a factor.

Author(s)

Florian Gerber, Roman Flury

Examples

x <- spam(1, 15, 15)print(x, minimal = FALSE)splits <- c(1,2,8,ncol(x)+1) # divide matrix into 9 submatricesfact <- array(1:9, c(3,3))   # multiply each submatrix with a different factoroF <- gmult(x, splits, fact)print(oF, minimal = FALSE)

Two trace plots and a scatter plot.

Description

For two (MCMC) chains of the same length trace plotsand a scatter plot are drawn.

Usage

grid_trace2(chain1, chain2 = NULL,            xlim = NULL, ylim1 = NULL, ylim2=NULL,            chain1_lab = "", chain2_lab = "", main = "",            chain1_yaxis_at = NULL, chain2_yaxis_at = NULL,            log = FALSE,            cex_points = unit(0.2, "mm"),            cex_main = unit(1.2, "mm"),            lwd_lines = unit(0.1, "mm"),            lwd_frame = unit(0.8, "mm"),            draw = TRUE)

Arguments

Details

The figure is optimized for a plot of the formatx11(width = 10, height = 3).

Value

A grob object.

Author(s)

Florian Gerber <florian.gerber@math.uzh.ch>

See Also

grid_zoom

Examples

grid_trace2(runif(500), runif(500),            chain1_yaxis_at = seq(.2, 1, by = .2),            chain1_lab = "chain1", chain2_lab = "chain2",            main = "Uniform",            lwd_lines = grid::unit(.5, "mm"))

grid_zoom

Description

This function takes a grob object (e.g. created withpackage grid) and adds a zoom window.

Usage

grid_zoom(inputGrob = pointsGrob(runif(200),runif(200)),          inputViewport = viewport(name='main'),          x = 'topleft', y, just,          ratio = c(.3,.4), zoom_xlim, zoom_ylim,          rect = TRUE, rect_lwd = 1, rect_fill = 'gray92',          draw =TRUE, zoom_fill = 'white',          zoom_frame_gp = gpar(lwd = 1),          zoom_gp = NULL, zoom_xaxis = xaxisGrob(main = FALSE),          zoom_yaxis = NULL)

Arguments

Details

A zoom plot does only make sense if all objects of theinputGrobare specified innative units. Additional caution me be requirefor certain grobs: e.g. a zoom of a circleGrob() is problematic if the xand y axis are stretched by a different amount.

Value

A grob object.

Author(s)

Florian Gerber <florian.gerber@math.uzh.ch>

See Also

grid_trace2

Examples

require(grid)## -- Example 1 --set.seed(133)grid_zoom(inputGrob = pointsGrob(runif(200), runif(200)),          inputViewport = viewport(name = 'main'),          zoom_xlim = c(.2, .3), zoom_ylim = c(.2, .3))## -- Example 2 --## initial plotgrid.newpage()vp <- viewport(width=.8, height=.8, clip='on')gt <- gTree(children=gList(polylineGrob(x=c((0:4)/10, rep(.5, 5), (10:6)/10, rep(.5, 5)),              y=c(rep(.5, 5), (10:6/10), rep(.5, 5), (0:4)/10),              id=rep(1:5, 4), default.units='native',              gp=gpar(col=1:5, lwd=3)),              pointsGrob(runif(1000), runif(1000),pch='.', gp=gpar(cex=3)),              rectGrob(gp=gpar(lwd=3))))pushViewport(vp)grid.draw(gt)## plot with zoom windowgrid.newpage()grid_zoom(inputGrob = gt,          inputViewport = vp,          x='topright', zoom_xlim=c(.6,.73), zoom_ylim=c(.3,.43),ratio=.4,          zoom_xaxis = NULL, zoom_gp = gpar(cex=3))

Return the First or Last Part of an Object

Description

Returns the upper left or lower right part of aspam object.

Usage

## S4 method for signature 'spam'head(x, n = 6L, m = n, ...)## S4 method for signature 'spam'tail(x, n = 6L, m = n, addrownums = TRUE, ...)

Arguments

Details

For matrices, 2-dim tables and data frames,head() (tail()) returnsthe first (last)n rows andm columns whenn > 0 or all but thelast (first)n rows whenn < 0 (with similar behaviorform).

tail() will add row and column names ofthe form"[n,]" and"[,n]" to the result, so that it looks similar to thelast lines and columns ofx when printed. Settingaddrownums = FALSE suppresses this behaviour.

A method forspam.chol.NgPeyton objects isexported as well.

Value

An regular matrix.

Author(s)

Reinhard Furrer

Examples

head( precmat.RW2( 10))tail( precmat.season(n=10, season=3), n=4, m=10)

Display a Sparse Matrix as Color Image

Description

The function creates a grid of colored rectangles withcolors corresponding to the values of thespam matrix.

Usage

## S4 method for signature 'spam'image(x, cex = NULL, ...)

Arguments

Details

getOption("spam.imagesize") determines if the sparse matrix iscoerced into a matrix and the plotted similarly toimage.default or ifthe matrix is simply represented as a scatterplot withpch=".". The points are scaled according tocex*getOption("spam.cex")/(nrow+ncol).For some devices or for non-square matrices,cex needs probably some adjustment.
The a zero matrix inspam format has as (1,1) entry the valuezero and only missing entries are interpreted asNA in thescatter plot.

Author(s)

Reinhard Furrer

See Also

display andspam.options.
The code is based onimage ofgraphics.

Examples

set.seed(13)smat <- spam_random(8)par(mfcol=c(1,2),pty='s')options(spam.imagesize=1000)image(smat) # or use better color schemesoptions(spam.imagesize=10)image(smat, cex=.25)smat <- spam_random(2^14, distribution=rnorm, density=1e-5, verbose=TRUE)par(mfcol=c(1,1), mai=c(.4,.4,.1,.1), pty='s')image(smat)

Read External Matrix Formats

Description

Read matrices stored in the Harwell-Boeing or MatrixMarket formats.

Usage

read.HB(file)read.MM(file)

Arguments

Details

The names of files storing matrices in theHarwell-Boeing format usually end in".rua" or".rsa".Those storing matrices in the MatrixMarket format usually end in".mtx".

Currently, only real assembled Harwell-Boeing can be read withread.HB. Reading MatrixMarket formats is more flexible.However, as entries ofspam matrices are of modedouble,integers matrices are coerced to doubles, patterns lead to matricescontaining ones and complex are coerced to the real part thereof. Inthese aforementioned cases, a warning is issued.

MatrixMarket also defines an array format, in which case a (possibly)densespam object is return (retaining only elements which arelarger thanoptions('spam.eps'). A warning is issued.

Value

A sparse matrix of classspam.

Note

The functions are based onreadHB andreadMM fromthe libraryMatrix to build the connection and read the rawdata.At present,read.MM is more flexible thanreadMM.

Author(s)

Reinhard Furrer based onMatrix functions byDouglas Batesbates@stat.wisc.edu and Martin Maechlermaechler@stat.math.ethz.ch

References

https://math.nist.gov/MatrixMarket/

https://sparse.tamu.edu/

Examples

## Not run: image(read.MM(gzcon(url(  "ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/bcspwr/bcspwr01.mtx.gz"))))## End(Not run)## Not run: ## Datasets supplied within Matrixstr(read.MM(system.file("external/pores_1.mtx",package = "Matrix")))str(read.HB(system.file("external/utm300.rua", package = "Matrix")))str(read.MM(system.file("external/lund_a.mtx", package = "Matrix")))str(read.HB(system.file("external/lund_a.rsa", package = "Matrix")))## End(Not run)

Test if a Sparse Matrix is Symmetric

Description

Efficient function to test if 'object' is symmetric or not.

Usage

# isSymmetric.spam(object, ...)## S3 method for class 'spam'isSymmetric(object, tol = 100 * .Machine$double.eps, ...)

Arguments

Details

symmetry is assessed by comparing the sparsity structure ofobject andt(object) via the functionall.equal.spam. If a difference is detected, the matrix iscleaned withcleanup and compared again.

Value

logical indicating ifobject is symmetric or not.

Author(s)

Reinhard Furrer

See Also

all.equal.spam,cleanup.

Examples

obj <- diag.spam(2)isSymmetric(obj)obj[1,2] <- .Machine$double.epsisSymmetric(obj)all.equal(obj, t(obj))

Kronecker Products on Sparse Matrices

Description

Computes the generalised kronecker product of two arrays,X andY.

Usage

kronecker.spam(X, Y, FUN = "*", make.dimnames = FALSE, ...)

Arguments

Details

The sparsity structure is determined by the ordinary%x%. Hence, the result ofkronecker(X, Y, FUN = "+") isdifferent depending on the input.

Value

An arrayA with dimensionsdim(X) * dim(Y).

Author(s)

Reinhard Furrer

Examples

# Starting with non-spam objects, we get a spam matrixkronecker.spam( diag(2), array(1:4, c(2, 2)))kronecker( diag.spam(2), array(1:4, c(2, 2)))# Notice the preservation of sparsity structure:kronecker( diag.spam(2), array(1:4, c(2, 2)), FUN="+")

Large 64-bit matrices require the R packagespam64

Description

The R packagespam can handle sparse matrices with upto 2^31-1 non-zero elements. For matrices with more non-zero elementsit is necessary to load thespam64 package in addition.

Details

With the help of the R packagedotCall64 spam interfaceseither the compiled code with 32-bit integers provided inspam or the compiled code with 64-bit integers provided inspam64.
To mimick 64-bit behavior, setoptions(spam.force64 = TRUE). The subsequent matrix indices are then stored in double format.

Author(s)

Reinhard Furrer, Florian Gerber, Kaspar Moesinger, Daniel Gerber

References

F. Gerber, K. Moesinger, R. Furrer (2017),Extending R packages to support 64-bit compiled code: An illustration with spam64 and GIMMS NDVI3g data,Computer & Geoscience 104, 109-119, https://doi.org/10.1016/j.cageo.2016.11.015.

See Also

spam64-package,dotCall64.

Examples

## Not run: ## the following matrices are very large, and hence,## require much memory and cpu time.library("spam64")s1 <- spam(1, ncol=2^30)        # 32-bit matrixs1s2 <- cbind(s1, s1)             # 64-bit matrixs2s3 <- spam(1, ncol=2^31)        # 64-bit matrixs3## End(Not run)

Lower and Upper Triangular Part of a Sparse Matrix

Description

Returns the lower or upper triangular structure orentries of a sparse matrix.

Usage

lower.tri(x, diag = FALSE)upper.tri(x, diag = FALSE)

Arguments

Details

Often not only the structure of the matrix is required but theentries as well. For compatibility, the default is only a structureconsisting of ones (representingTRUEs). Setting the flaggetOption( "spam.trivalues") toTRUE,the function returns the actualvalues.

See Also

spam.options anddiag

Examples

smat <- spam( c( 1,2,0,3,0,0,0,4,5),3)upper.tri( smat)upper.tri( smat, diag=TRUE)options(spam.trivalues=TRUE)upper.tri( smat)

Create Precision Matrices

Description

Creates precision matrices for gridded GMRF.

Usage

precmat.GMRFreglat(n,m, par, model = "m1p1",  eps = getOption("spam.eps"))

Arguments

Details

The function should be illustrative on how to create precisionmatrices for gridded GMRF. Hence, no testing (positive definiteness isdone).

The model specification"m" determines the complexity and"p" the number of parameters.

Please see the examples on the meaning of the different models.

Value

Aspam matrix of dimensionprod(dims)xprod(dims).

Author(s)

Reinhard Furrer

See Also

precmat,toeplitz.spam,kronecker.spam

Examples

as.matrix(precmat.GMRFreglat(4, 3, c(.4),         'm1p1'))as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3),      'm1p2'))as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3,.2),   'm2p3'))as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3,.2,.1),'m2p4'))# up to the diagonal, the following are equivalent:cleanup( precmat.IGMRFreglat(3,4) -             precmat.GMRFreglat(3,4,1, 'm1p1'))

Administrative districts of Germany

Description

Displaying dataover the administrative districts of Germany

Usage

map.landkreis(data, col=NULL, zlim=range(data), add=FALSE,              legendpos=c( 0.88,0.9,0.05,0.4))

Arguments

Details

The functiongermany.plot super-seedsmap.landkreis (it is several factors faster).

The perfect position of the legend is an art per se and depends onvariouspar parameters. See also the source code of the functionimage.plot offields.

Author(s)

Reinhard Furrer

References

The code ofmap.landkreis is very similar togermany.map from the packageINLA.

See Also

germany.plot super-seedingmap.landkreis.

Examples

## Not run: data( Oral)par( mfcol=c(1,2))germany.plot( log( Oral$Y), legend=TRUE)map.landkreis( log( Oral$Y))## End(Not run)

Maximum likelihood estimates

Description

Maximum likelihood estimates of a simple spatial model

Usage

neg2loglikelihood.spam(y, X, distmat, Covariance,                 beta, theta, Rstruct = NULL, cov.args = NULL, ...)neg2loglikelihood(y, X, distmat, Covariance,                 beta, theta, cov.args = NULL, ...)neg2loglikelihood.nomean(y, distmat, Covariance,                 theta, cov.args = NULL, ...)mle.spam(y, X, distmat, Covariance,     beta0, theta0, thetalower, thetaupper, optim.control=NULL,     Rstruct = NULL, hessian = FALSE, cov.args = NULL, ...)mle(y, X, distmat, Covariance,     beta0, theta0, thetalower, thetaupper, optim.control=NULL,     hessian = FALSE, cov.args = NULL, ...)mle.nomean.spam(y, distmat, Covariance,     theta0, thetalower, thetaupper, optim.control=NULL,     Rstruct = NULL, hessian = FALSE, cov.args = NULL, ...) mle.nomean(y, distmat, Covariance,     theta0, thetalower, thetaupper, optim.control=NULL,     hessian = FALSE, cov.args = NULL, ...)

Arguments

Details

We provide functions to calculate thenegative-2-log-likelihood and maximum likelihood estimates for the model

y ~ N_n( X beta, Sigma(h;theta) )

in the case of a sparse or ordinary covariance matrices.

In the case of the*.spam versions, the covariance function hasto return aspam object. In the other case, the methods arecorrectly overloaded and work either way, slightly slower than the*.spam counterparts though.

When working on the sphere, the distance matrix has to be transformed by

h -> R / 2 sin(h/2)

where R is the radius of the sphere.

The covariance function requires that the first argument is the distancematrix and the second the parameters. One can image cases in which thecovariance function does not take the entire distance matrix but onlysome partial information thereof. (An example is the use of a kroneckertype covariance structure.) In case of a sparse covariance constructionwhere the argumentRstruct is not given, the first parameterelement needs to be the range parameter. (This results from the fact,that a sparse structure is constructed that is independent of theparameter values to exploit the fast Choleski decomposition.)

In the zero-mean case, theneg2loglikelihood is calculated by settingthe parametersX orbeta to zero.

Value

The negative-2-loglikelihood or the output from the functionoptim.

Author(s)

Reinhard Furrer

See Also

covmat,rmvnorm.spam

Examples

# True parameter values:truebeta <- c(1,2,.2)    # beta = (intercept, linear in x, linear in y)truetheta <- c(.5,2,.02) # theta = (range, sill, nugget)# We now define a grid, distance matrix, and a sample:x <- seq(0,1,l=5)locs <- expand.grid( x, x)X <- as.matrix( cbind(1,locs))  # design matrixdistmat <- nearest.dist( locs, upper=NULL) # distance matrixSigma <- cov.sph( distmat, truetheta)    # true covariance matrixset.seed(15)y <- c(rmvnorm.spam(1,X %*% truebeta,Sigma)) # construct sample# Here is the negative 2 log likelihood:neg2loglikelihood.spam( y, X, distmat, cov.sph,                       truebeta, truetheta)# We pass now to the mle:res <- mle.spam(y, X, distmat, cov.sph,         truebeta, truetheta,thetalower=c(0,0,0),thetaupper=c(1,Inf,Inf))# Similar parameter estimates here, of course:mle.nomean.spam(y-X%*%res$par[1:3], distmat, cov.sph,         truetheta, thetalower=c(0,0,0), thetaupper=c(1,Inf,Inf))

Distance Matrix Computation

Description

This function computes and returns specific elements of distancematrix computed byusing the specified distance measure.

Usage

nearest.dist( x, y=NULL, method = "euclidean",             delta = 1, upper = if (is.null(y)) FALSE else NULL,             p = 2, miles = TRUE, R = NULL, fortran = FALSE)

Arguments

Details

For great circle distance, the matricesx andycontain the degrees longitudes in the first and the degrees latitudesin the second column.delta is indegrees. Hence to restrict to distances smaller thandelta.km,one has to specifydelta=delta.km*360/(6378.388*2*pi).

The distances are calculated based on spherical law of cosines.Care is needed for ‘zero’ distances due to the final acosin:acos(1-1e-16), especially with an actual radius.

Default value of Earth's radius is 3963.34miles (6378.388km).

There are many other packages providing distance functions. Especiallyfor great circle distances there are considerable differences betweenthe implementations. For high precision results,sp::spDists isa good candidate and distances of large amount of locations can beprocessed in parallel with theparallelDist package.

The formerly depreciated argumentseps anddiag are noweliminated.

x andy can be any object with an existingas.matrix method.

The Fortran code is based on a idea of Doug Nychka.

Value

Aspam object containing the distances spanned betweenzero anddelta. The sparse matrix may contain many zeros(e.g., collocated data). However, to calculate covariances, these zerosare essential.

Author(s)

Annina Cincera (C++ code), Reinhard Furrer

See Also

spam_rdist

Examples

# Note that upper=T and using t(X)+X is quicker than upper=NULL;#     upper=T marginally slower than upper=F.# To compare nearest.dist with dist, use as.dist(...)nx <- 4x <- expand.grid(as.double(1:nx),as.double(1:nx))sum( ( as.dist(nearest.dist( x, delta=nx*2))-          dist(x)                            )^2)# Create nearest neighbor structures:par(mfcol=c(1,2))x <- expand.grid(1:nx,1:(2*nx))display( nearest.dist( x, delta=1))x <- expand.grid(1:(2*nx),1:nx)display( nearest.dist( x, delta=1))

Options Settings

Description

Allow the user to set and examine a variety ofoptionswhich affect the way in whichR computes and displays sparsematrix results.

Details

Invokingoptions() with no arguments returns a list with thecurrent values of the options. To access the value of a single option, one shouldusegetOption("spam.eps"), e.g., rather thanoptions("spam.eps") which is alist of length one.

Of course, printing is still subordinate togetOption("max.print") or similar options.

Value

ForgetOption, the current value set for optionx, orNULL if the option is unset.

Foroptions(), a list of all set options sorted by category. Foroptions(name), a list of length one containing the set value,orNULL if it is unset. For uses setting one or more options,a list with the previous values of the options changed (returnedinvisibly).

Options used for the packagespam

A short description with the default values follows.

spam.eps=.Machine$double.eps:

values smaller than this areconsidered as zero. This is only used when creating spam objects.

spam.drop=FALSE:

default parameter fordrop when subsetting

spam.printsize=100:

the max number of elements of a matrix which wedisplay as regular matrix.

spam.imagesize=10000:

the max number of elements of a matrix we displayas regular matrix withimage ordisplay. Larger matrices are representedas dots only.

spam.cex=1200:

default dot size forimage ordisplay.

spam.structurebased=FALSE:

should operations be carried out onthe nonzero entries (the structure) or including the zeros.

spam.inefficiencywarning=1e6:

issue a warning when inefficientoperations are performed and the matrix exceeds the specified size.Valid value is a postive integer or a logical.TRUE correspondsto 1 (always),FALSE toInf.

spam.trivalues=FALSE:

a flag whether to return the structure(FALSE) or the values themselves (TRUE) when returning theupper and lower triangular part of a matrix.

spam.listmethod="PE":

algorithm forspam.list. Defaultis suggestion by Paul Eilers (thanks). Any other specification uses abubble sort algorithm which is only slightly faster for very sparse matrices.

spam.dopivoting=TRUE:

default parameter for "solve" routines.FALSEwould solve the system without using the permutation.

spam.NAOK=FALSE:

logical determines ifNA,NaN andInf areallowed to Fortan. Setting toTRUE allows to work with these butfull functionality has not been tested.

spam.safemodevalidity=TRUE:

logical determines if sanity checkis peformed when constructing sparse matrices.Default is safer but somewhat slower.

spam.cholsymmetrycheck=TRUE:

for the Cholesky factorization,verify if the matrix is symmetric.

spam.cholpivotcheck=TRUE:

for the Cholesky factorization,when passing a permutation, should a minimum set of checks be performed?

spam.cholupdatesingular="warning":

for a Cholesky update, whathappens if the matrix is singular:"warning" only andreturning the not updated factor,"error" or return simply"NULL".

spam.cholincreasefactor=c(1.25,1.25):

If not enought memorycould be allocated, these are the steps to increase it.

spam.nnznearestdistnnz=c(400^2,400):

Memory allocationparameters fornearest.dist.

spam.nearestdistincreasefactor=1.25:

If not enought memorycould be allocated, this is the step to increase it.

See Also

Functions influenced by these options include:print.spam,display.spam,image.spam,upper.tri.spam,chol.spam,nearest.dist, etc.
powerboost

Examples

smat <- diag.spam( 1:8)smatoptions(spam.printsize=49)smat# List all spam options:options()[grep("spam",names(options()))]# Reset to default values:options(spam.eps=.Machine$double.eps,        spam.drop=FALSE,        spam.printsize=100,        spam.imagesize=10000,        spam.cex=1200,        spam.structurebased=FALSE,        spam.inefficiencywarning=1e6,        spam.trivalues=FALSE,        spam.listmethod="PE",        spam.NAOK=FALSE,        spam.safemodevalidity=TRUE,        spam.dopivoting=TRUE,        spam.cholsymmetrycheck=TRUE,        spam.cholpivotcheck=TRUE,        spam.cholupdatesingular="warning",        spam.cholincreasefactor=c(1.25,1.25),        spam.nearestdistincreasefactor=1.25,        spam.nearestdistnnz=c(400^2,400))

Extract the permutation

Description

Extract the (inverse) permutation used by the Choleskydecomposition

Usage

ordering( x, inv=FALSE)

Arguments

Details

Recall that calculating a Cholesky factor from a sparse matrixconsists of finding a permutation first, then calculating the factorsof the permuted matrix. The ordering is important when working withthe factors themselves.

The ordering from a full/regular matrix is1:n.

Note that there exists many different algorithms to findorderings.

See the examples, they speak more than 10 lines.

Author(s)

Reinhard Furrer

See Also

chol.spam,solve.spam.

Examples

# Construct a pd matrix S to work with (size n)n <- 100    # dimensionS <- .25^abs(outer(1:n,1:n,"-"))S <- as.spam( S, eps=1e-4)I <- diag(n)  # Identity matrixcholS <- chol( S)ord <- ordering(cholS)iord <- ordering(cholS, inv=TRUE)R <- as.spam( cholS ) # R'R = P S P', with P=I[ord,],  # a permutation matrix (rows permuted).RtR <- t(R) %*% R# the following are equivalent:as.spam( RtR -            S[ord,ord],    eps=1e-15)as.spam( RtR[iord,iord] - S,             eps=1e-15)as.spam( t(R[,iord]) %*% R[,iord] - S, eps=1e-15)# we use 'eps' to avoid issues close to machine precision# trivially:as.spam( t(I[iord,]) - I[ord,])  # (P^-1)' = P  as.spam( t(I[ord,]) - I[,ord])  # as.spam( I[iord,] - I[,ord])as.spam( I[ord,]%*%S%*%I[,ord] - S[ord,ord] )   # pre and post multiplication with P and P' is ordering

Padding a Sparse Matrix

Description

Resets the dimension of aspam matrix by truncation or zero padding.

Usage

pad(x) <- value

Arguments

Details

It is important to notice the different behavior of the replacementmethod for ordinary arrays andspam objects (see‘Examples’). Here, the elements are not simply rearranged butan entirely new matrix is constructed. If the new column dimension issmaller than the original, the matrix is also cleaned (withoption("spam.eps") as filter).

Value

A (spam) matrix of dimensionvalue where trunction orpadding hasbeen used.

Author(s)

Reinhard Furrer

See Also

dim.spam.

Examples

x <- diag(4)dim(x)<-c(2,8)xs <- diag.spam(4)pad(s) <- c(7,3)  # any positive value can be useds <- diag.spam(4)pad(s) <- c(2,8)  # result is different than x

Permute a Matrix

Description

Row and/or column permutes a (sparse) matrix.

Usage

permutation.spam(A, P=NULL, Q=NULL, ind=FALSE, check=TRUE)

Arguments

Details

If P and Q are permutation matrices, the result isPAQ. However, it is also possible to specify the indicesand to perform in a very efficient wayA[rowind, colind], see examples.

A row permutation is much faster than a colum permutation.For very large matrices, a double transpose might be faster.

The spam optionspam.checkpivot determines if the permutationis verified.

Value

A permuted matrix.

Author(s)

Reinhard Furrer

See Also

ordering,spam.options.

Examples

A <- spam(1:12,3)P <- c(3,1,2)Q <- c(2,3,1,4)permutation(A,P,Q)-A[order(P),order(Q)]permutation(A,P,Q,ind=TRUE)-A[P,Q]

Specific Options Setting

Description

Sets several options for speed-up.

Usage

powerboost(flag)

Arguments

Details

The options turn checking off ("safemode","cholsymmetrycheck" and"cholpivotcheck") and switch tosingle precision for"eps".

Value

NULL in any case.

Author(s)

Reinhard Furrer, after receiving too much C.mc.st adds.

See Also

spam.options.

IGMRF Precision Matrices

Description

Fast ways to create sparse precision matrices for various IGMRF.

Usage

precmat(n, season=12, m=n, A=NULL, order=1, ... , type="RW1")precmat.RW1(n)precmat.RW2(n)precmat.RWn(n, order=3)precmat.season(n, season=12)precmat.IGMRFreglat(n, m, order=1, anisotropy=1)precmat.IGMRFirreglat(A, eps=getOption("spam.eps"))

Arguments

Details

precmat is a wrapper that calls the other functionsaccording to the argumenttype.
Implements many of the precision matrices discussed in Chapter3 of Rue and Held (2005). For example,precmat.RW1,precmat.RW2 andprecmat.season are given inequations (3.22), (3.40) and (3.59);precmat.IGMRFreglat onpage 107. Note that for the latter we reverse the order of the dimension here!
If adjacency matrix is a regular matrix, it is coerced to aspam object. Only the structure is used. Make sure, that thediagonal is empty.

Value

A sparse precision matrix.

Author(s)

Reinhard Furrer

References

Rue and Held (2005).

See Also

precmat.GMRFreglat,rmvnorm.prec,adjacency.landkreis.

Examples

n <- 10Q <- precmat.RW2( n)# rmvnorm.prec(1, Q=Q) # does not work, because the matrix is singular.Q%*%cbind(1,1:n)

Printing and Summarizing Sparse Matrices

Description

Printing (non-zero elements) of sparse matricesand summarizing the sparsity structure thereof.

Usage

  ## S4 method for signature 'spam'print(x, ...)## S4 method for signature 'spam'summary(object, ...)

Arguments

Details

getOption('spam.printsize') determines if the sparse matrix iscoerced into a matrix and the printed as an array or ifonly the non-zero elements of the matrix are given.

Value

NULL forprint, because the information is printed withcat thereis no real need to pass any object back.
A list containing the non-zero elements and the density forsummaryfor classspam.
A list containing the non-zero elements of the factor, the density andthe fill-in forsummary for classspam.chol.NgPeyton.

Author(s)

Reinhard Furrer

See Also

display orimage for a graphical visualization;spam.options

Examples

set.seed(13)smat <- spam_random(8)par(mfcol=c(1,2),pty='s')options(spam.printsize=1000)print(smat)options(spam.printsize=10)print(smat)summary(smat)summary(smat)$nnzsmat@entries[1:5] <- 0print(smat, minimal = FALSE)print(smat, minimal = FALSE, rowpointer = TRUE)smat@rowpointersprint_nnzpos(smat)

Create Random Sparse Matrices

Description

Creates random spam matrix given the dimension and other parameters.

Usage

spam_random(nrow = 1L, ncol = nrow, density = 0.5, distribution = NULL, digits = NULL,            sym = FALSE, spd = FALSE, verbose = FALSE, ...)

Arguments

Details

To create a random spam64 matrix, setoptions(spam.force64 = TRUE).

Value

A random matrix inspam format.

Author(s)

Florian Gerber, Roman Flury, Reinhard Furrer

See Also

spam-class anddisplay.spam

Examples

set.seed(42)rspam <- spam_random(500, digits = 2, distribution = rnorm, sd = 2, mean = 10, density = .01)display.spam(rspam, cex = 2)

Draw From a Gaussian Random Field

Description

Fast and intuitive ways to draw from a Gaussian random field.

Usage

rgrf( n,   locs, nx, ny=nx, xlim=c(0,1), ylim=c(0,1), tau=0,   Covariance, theta, beta=0, X,    method=c('chol'),  method.args=list(sparse=FALSE),    eps = getOption("spam.eps"), drop=TRUE, attributes=TRUE, ...)

Arguments

Details

If no locations are given, the function constructs theseaccording a regular or a regular perturbed grid. The perturbation isdetermined bytau, which has to be greater than zero (noperturbation) and strictly smaller than 1/2 (max perturbation).

The regular grid has spacing (here for x)dx=diff(xlim)/nx and runsfromxlim[1]+dx/2 toxlim[2]-dx/2.The locations are at least(1/nx-2*tau*dx) separated.

Currently, the only method implemented is a Cholesky factorizationroutine, (much as inrmvnorm).

Therdist() from thefields package is awefullyfast. Unless one has very sparse covariance matrices, a sparseapproach is not bringing a lot of improvements.

The methods may use different covariance construction approaches andthus the nesting ofcov.args inmethod.args.

Author(s)

Reinhard Furrer

See Also

rgrf,chol andordering.

Examples

require(fields)# Regular grid with constant mean:nx <- 10field <- rgrf(1, nx=nx,  Covariance="cov.wend2", theta=c(.5, 1), beta=5)quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=nx)points(attr(field,"locs"))# Irregluar grid:field <- rgrf(1, nx=10, tau=0.3, Covariance="cov.mat", theta=c(.2, 1, 1.5))fields::quilt.plot(attr(field,"locs"), field)

Draw Multivariate Normals

Description

Fast ways to draw multivariate normals when the variance or precision matrixis sparse.

Usage

rmvnorm(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, ..., mean, sigma)rmvnorm.spam(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, Rstruct=NULL, ..., mean, sigma)rmvnorm.prec(n, mu=rep.int(0, dim(Q)[1]), Q, Rstruct=NULL, ...)rmvnorm.canonical(n, b, Q, Rstruct=NULL, ...)

Arguments

Details

All functions rely on a Cholesky factorization of thecovariance or precision matrix.
The functionsrmvnorm.prec andrmvnorm.canonicaldo not require sparse precision matricesDepending on the the covariance matrixSigma,rmvnormorrmvnorm.spam is used. If wrongly specified, dispatching tothe other function is done.
Default mean is zero. Side note: mean is added viasweep() andno gain is accieved by distinguishing this case.
Often (e.g., in a Gibbs sampler setting), the sparsity structure ofthe covariance/precision does not change. In such setting, theCholesky factor can be passed viaRstruct in which only updatesare performed (i.e.,update.spam.chol.NgPeyton instead of afullchol).

Author(s)

Reinhard Furrer

References

See references inchol.

See Also

rgrf,chol andordering.

Examples

# Generate multivariate from a covariance inverse:# (usefull for GRMF)set.seed(13)n <- 25    # dimensionN <- 1000  # sample sizeSigmainv <- .25^abs(outer(1:n,1:n,"-"))Sigmainv <- as.spam( Sigmainv, eps=1e-4)Sigma <- solve( Sigmainv)  # for verificationiidsample <- array(rnorm(N*n),c(n,N))mvsample <- backsolve( chol(Sigmainv), iidsample)norm( var(t(mvsample)) - Sigma, type="m")# compare with:mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample, n)   #### ,n as patchnorm( var(t(mvsample)) - Sigma, type="m")# 'solve' step by step:b <- rnorm( n)R <- chol(Sigmainv)norm( backsolve( R, forwardsolve( R, b))-      solve( Sigmainv, b) )norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma )

Draw Conditional Multivariate Normals

Description

Fast way to draw conditional multivariate normals when the covariancematrix is sparse.

Usage

rmvnorm.conditional(n, y, mu = rep.int(0, dim(SigmaXX)[1]+dim(SigmaYY)[1]),                    SigmaXX, SigmaYY, SigmaXY, noise, RstructYY = NULL, ...)

Arguments

Details

Quite often, we want to draw condional observationsX|yfrom the modelY=X+e, whereX has covariance matrixSigmaXX ande has white noise.

Covariance ofY can be specified bySigmaYY orSigmaXX+diag(noise,). IfY andX do not have thesame dimensions,SigmaXY needs to be specified.

The function also implmements a general multivariate model, where thewe only observe part of the vector. The components are firstX thenY.

The functionrmvnorm.cond() is a wrapper tormvnorm.conditional() and included to increase similaritieswith other packages.

Author(s)

Reinhard Furrer

See Also

rmvnorm.spam.

Examples

set.seed(12)N <- 300y <- c(5, -5, -5, 5)SigmaXX <- as.spam(.95^abs(outer(1:N, 1:N, "-")), eps=1e-4)sel <- c(10, 100, 120, 300)        # where we observe ySigmaXY <- SigmaXX[, sel]SigmaYY <- SigmaXX[sel,sel] + diag.spam(.01, length(y)) # some noisex <- rmvnorm.conditional(3, y, SigmaXX=SigmaXX, SigmaXY=SigmaXY,                         SigmaYY=SigmaYY)# unconditional sample:ux <- rmvnorm(1, Sigma=SigmaXX)matplot(t(rbind(x, ux)), type='l', lty=1)points(sel, y, pch=19)

Draw Constrainted Multivariate Normals

Description

Fast ways to draw multivariate normals with linear constrains when the variance or precision matrixis sparse.

Usage

rmvnorm.const(n, mu = rep.int(0, dim(Sigma)[1]), Sigma, Rstruct = NULL,              A = array(1, c(1,dim(Sigma)[1])), a=0, U=NULL,  ...)rmvnorm.prec.const(n, mu = rep.int(0, dim(Q)[1]), Q, Rstruct = NULL,              A = array(1, c(1,dim(Q)[1])), a=0, U=NULL,  ...)rmvnorm.canonical.const(n, b, Q, Rstruct = NULL,              A = array(1, c(1,dim(Q)[1])), a=0, U=NULL, ...)

Arguments

Details

The functionsrmvnorm.prec andrmvnorm.canonicaldo not requrie sparse precision matrices.Forrmvnorm.spam, the differences between regular and sparsecovariance matrices are too significant to be implemented here.
Often (e.g., in a Gibbs sampler setting), the sparsity structure ofthe covariance/precision does not change. In such setting, theCholesky factor can be passed viaRstruct in which only updatesare performed (i.e.,update.spam.chol.NgPeyton instead of afullchol).

Author(s)

Reinhard Furrer

References

See references inchol.

See Also

rmvnorm.spam.

Examples

# to be filled in

Draw From a Multivariate t-Distribution

Description

Fast ways to draw from a multivariate t-distribution the scale (covariance) matrixis sparse.

Usage

rmvt(n, Sigma, df = 1, delta = rep(0, nrow(Sigma)),    type = c("shifted", "Kshirsagar"), ..., sigma)rmvt.spam(n, Sigma, df = 1, delta = rep(0, nrow(Sigma)),    type = c("shifted", "Kshirsagar"), ..., sigma)

Arguments

Details

This function is very much likermvt() from the packagemvtnorm. We refer to the help of the afore mentioned.

Author(s)

Reinhard Furrer

References

See references inmvtnorm::rmvt().

See Also

rmvnorm.

Form Row and Column Sums and Means

Description

Form row and column sums and means for sparsespam matrices

Usage

rowSums(x, na.rm = FALSE, dims = 1, ...)colSums(x, na.rm = FALSE, dims = 1, ...)rowMeans(x, na.rm = FALSE, dims = 1, ...)colMeans(x, na.rm = FALSE, dims = 1, ...)

Arguments

Details

Depending on the flag .

Value

Vector of appropriate length.

Author(s)

Reinhard Furrer

See Also

apply.spam,spam.options.

Examples

x <- spam( rnorm(20), 5, 4)rowSums( x)c( x %*% rep(1,4))

Wappers for Sparse Matrices

Description

These functions are convenient wrappers forspam objectsto classical matrix operations.

Usage

var.spam(x, ...)## S3 method for class 'spam'var(x, ...)

Arguments

Details

There is probably no point in fully defining methodshere. Typically, these functions do not exploit sparsitystructures. Hence, for very large matrices, warnings may be posted.

Value

Depends on function...

Author(s)

Reinhard Furrer

See Also

Option"inefficiencywarning" inspam.options.

Sparse Matrix Class

Description

This group of functions evaluates and coerces changes in class structure.

Usage

spam(x, nrow = 1, ncol = 1, eps = getOption("spam.eps"))as.spam(x, eps = getOption("spam.eps"))is.spam(x)

Arguments

Details

The functionsspam andas.spam act likematrixandas.matrixto coerce an object to a sparse matrix object of classspam.

Ifx is a list, it should contain either two or three elements.In case of the former, the list should contain an by twomatrix of indicies (calledind) and the values.In case of the latter, the list should contain three vectorscontaining the row, column indices (calledi andj) and the values. In both cases partial matching is done.In case there are several triplets with the samei,j,the values are added.

eps should be at least as large as.Machine$double.eps.

IfgetOption("spam.force64") isTRUE, a 64-bitspam matrix is returned in any case. IfFALSE, a 32-bitmatrix is returned when possible.

Value

A validspam object.
is.spam returnsTRUE ifx is aspam object.

Note

The zero matrix has the element zero stored in (1,1).

The functions do not test the presence ofNA/NaN/Inf. Virtuallyall call a Fortran routine with theNAOK=NAOKargument, which defaults toFALSE resulting in an error.Hence, theNaN do not always properly propagate through (i.e.spam is not IEEE-754 compliant).

Author(s)

Reinhard Furrer

References

Reinhard Furrer, Stephan R. Sain (2010)."spam: A Sparse Matrix R Package with Emphasis on MCMCMethods for Gaussian Markov Random Fields.",Journal of Statistical Software, 36(10), 1-25,doi:10.18637/jss.v036.i10.

See Also

SPAM for a general overview of the package;spam_random to create matrices with a random sparsity pattern;cleanup to purge a sparse matrix;spam.options for details about thesafemode flag;read.MM andforeign to createspammatrices from MatrixMarketfiles and from certainMatrix orSparseM formats.

Examples

# old message, do not loop, when you create a large sparse matrixset.seed(13)nz <- 128ln <- nz^2smat <- spam(0,ln,ln)is <- sample(ln,nz)js <- sample(ln,nz)## IGNORE_RDIFF_BEGINsystem.time(for (i in 1:nz) smat[is[i], js[i]] <- i)system.time(smat[cbind(is,js)] <- 1:nz)## IGNORE_RDIFF_ENDgetClass("spam")options(spam.NAOK=TRUE)as.spam(c(1, NA))

Spam internal and auxiliary functions

Description

The functions listed below are auxiliary functions and are not exported by the NAMESPACE. The user should not require to call these directly.

Details

printSize(x, size=8, digits=2L) constructs a string forpretty printing the object sizex.

Methods for Sparse Matrices

Description

Methods without any additional parameters for sparse matrices.

Details

If a method forspam objects takes the same arguments,produces the intuitive output. We do not provide additional helppages.

However, such methods are usually linked to axzy.spamfunction, that could also be called directly.

Author(s)

Reinhard Furrer

See Also

Corresponding base help functions.

Class "spam"

Description

Thespam class is a representation of sparse matrices.

Objects from the Class

Objects can be created by calls of the formnew("spam", entries, colindices, rowpointes, dimension).The standard "old Yale sparse format" is used to store sparse matrices.
The matrixx is stored in row form. The first element of rowi isx@rowpointers[i]. The length of rowi is determined byx@rowpointers[i+1]-x@rowpointers[i]. The column indices ofx are stored inthex@colindices vector. The column index for elementx@entries[k] isx@colindices[k].

Slots

entries:

Object of class"numeric" contains thenonzero values.

colindices:

Object of class"integer" ordered indicesof the nonzero values.

rowpointers:

Object of class"integer" pointer to the beginningof each row in the arraysentries andcolindices.

dimension:

Object of class"integer" specifyingthe dimension of the matrix.

Methods

as.matrix

signature(x = "spam"):transforming a sparsematrix into a regular matrix.

as.vector

signature(x = "spam"):transforming a sparsematrix into a vector (dependings onstructurebased) seeas.vector.spam for details.

as.spam

signature(x = "spam"):cleaning of a sparse matrix.

[<-

signature(x = "spam", i,j, value):assigning asparse matrix. The negative vectors are not implemented yet.

[

signature(x = "spam", i, j):subsetting asparse matrix. The negative vectors are not implemented yet.

%*%

signature(x, y):matrix multiplication, all combinations of sparse with fullmatrices or vectors are implemented.

c

signature(x = "spam"):vectorizes the sparse matrix and takes account of the zeros. Hencethe lenght of the result isprod(dim(x)).

cbind

signature(x = "spam"): binds sparse matrices, seecbind for details.

chol

signature(x = "spam"):seechol for details.

diag

signature(x = "spam"):seediag for details.

dim<-

signature(x = "spam"): rearrangesthe matrix to reflect a new dimension.

dim

signature(x = "spam"): gives the dimension of thesparse matrix.

pad<-

signature(x = "spam"): truncates or augmentsthe matrix seedim for details.

image

signature(x = "spam"):seeimage for details.

display

signature(x = "spam"):seedisplay for details.

length<-

signature(x = "spam"): Is not implemented andcauses an error.

length

signature(x = "spam"): gives the number ofnon-zero elements.

lower.tri

signature(x = "spam"): seelower.tri for details.

Math

signature(x = "spam"): seeMath for details.

Math2

signature(x = "spam"): seeMath2 for details.

norm

signature(x = "spam"): calculates the norm of a matrix.

plot

signature(x = "spam", y): same functionality asthe ordinaryplot.

print

signature(x = "spam"): seeprint for details.

rbind

signature(x = "spam"): binds sparsematrices, seecbind for details.

solve

signature(a = "spam"): seesolve for details.

summary

signature(object = "spam"): small summarystatement of the sparse matrix.

Summary

signature(x = "spam"):All functions of theSummary class (likemin,max,range...) operate on the vectorx@entries and return theresult thereof. See Examples orSummary for details.

t

signature(x = "spam"): transpose of a sparse matrix.

upper.tri

signature(x = "spam"): seelower.tri for details.

Details

The compressed sparse row (CSR) format is often described with thevectorsa,ia,ja. To be a bit morecomprehensive, we have chosen longer slot names.

Note

The slotscolindices androwpointers aretested for proper integer assignments. This is not true forentries.

Author(s)

Reinhard Furrer, some of the Fortran code is based on A. George,J. Liu, E. S. Ng, B.W Peyton and Y. Saad (alphabetical)

Examples

showMethods("as.spam")smat <- diag.spam(runif(15))range(smat)cos(smat)

Defunct Objects in Packagespam

Description

The functions or variables listed here are defunct, i.e. thorw an error when used.

Usage

  validspamobject(...)

Arguments

See Also

Deprecated,Defunct

Basic Linear Algebra for Sparse Matrices

Description

Basic linear algebra operations for sparse matricesof classspam.

Details

Linear algebra operations for matrices of classspam are designed to behave exactly as forregular matrices. In particular, matrix multiplication, transpose, addition,subtraction and various logical operations should work as with the conventionaldense form of matrix storage, as does indexing, rbind, cbind, and diagonalassignment and extraction (see for examplediag).Further functions with identical behavior aredim and thusnrow,ncol.

The functionnorm calculates the (matrix-)norm of the argument.The argumenttype specifies thel1 norm,sup or maxnorm (default), or the Frobenius or Hilbert-Schmidt(frobenius/hs) norm. Partial matching can be used. For example,norm is used to check for symmetry in the functionchol bycomputing the norm of the difference between the matrix and itstranspose

The operator%d*% efficiently multiplies a diagonal matrix (invector form) and a sparse matrix and is used for compatibility with thepackage fields. More specifically, this method is used in the internalfunctions ofKrig to make the code more readable. It avoidshaving a branch in the source code to handle the diagonal or nondiagonalcases. Note that this operator is not symmetric: a vector inthe left argument is interpreted as a diagonal matrix and a vector inthe right argument is kept as a column vector.

The operator%d+% efficiently adds a diagonal matrix (in vectorform) and a sparse matrix, similarly to the operator%d+%.

References

Some Fortran functions are based onhttps://github.com/johannesgerer/jburkardt-f/blob/master/sparsekit/sparsekit.html

See Also

spam for coercion and other class relations involving thesparse matrix classes.

Examples

# create a weight matrix and scale it:## Not run: wij <- distmat# with distmat from a nearest.dist(..., upper=TRUE) calln <- dim(wij)[1]wij@entries <- kernel( wij@entries, h) # for some function kernelwij <- wij + t(wij) + diag.spam(n)     # adjust from diag=FALSE, upper=TRUEsumwij <- wij %*% rep(1,n)    # row scaling:    #   wij@entries <- wij@entries/sumwij[ wij@colindices]    # col scaling:wij@entries <- wij@entries/sumwij[ rep(1:n, diff(wij@rowpointers))]## End(Not run)

Linear Equation Solving for Sparse Matrices

Description

backsolve andforwardsolve solve a systemof linear equations where the coefficient matrixis upper or lower triangular.
solve solves a linear system or computes the inverseof a matrix if the right-hand-side is missing.

Usage

## S4 method for signature 'spam'solve(a, b, Rstruct=NULL, ...)## S4 method for signature 'spam'backsolve(r, x, ...)## S4 method for signature 'spam'forwardsolve(l, x, ...)## S4 method for signature 'spam'chol2inv(x, ...)

Arguments

Details

We can solveA %*% x = b by first computing the Cholesky decompositionA = t(R)%*%R), then solvingt(R)%*%y = b fory, andfinally solvingR%*%x = y forx.solve combineschol, a Cholesky decomposition of asymmetric positive definite sparse matrix, withforwardsolve andthenbacksolve.

In casea is from achol call, thensolve is anefficient way to calculatebacksolve(a, forwardsolve( t(a), b)).

However, fora.spam anda.mat from achol callwith a sparse and ordinary matrix, note thatforwardsolve( a.mat, b, transpose=T, upper.tri=T)is equivalent toforwardsolve( t(a.mat), b)andbacksolve(a.spam, forwardsolve(a.spam, b, transpose=T, upper.tri=T))yields the desired result. Butbacksolve(a.spam,forwardsolve(t(a.spam), resid)) iswrong becauset(a.spam) is aspam and not aspam.chol.NgPeyton object.

forwardsolve andbacksolve solve a system of linearequations where the coefficient matrix is lower (forwardsolve) orupper (backsolve) triangular. Usually, the triangular matrix isresult from achol call and it is not required to transpose itforforwardsolve. Note that arguments of the defaultmethodsk,upper.tri andtranspose do not have anyeffects here.

Notice that it is more efficient to solve successively the linearequations (both triangular solves) than to implement these in theFortran code.

If the right-hand-side insolve is missing it will computethe inverse of a matrix. For details about the specific Cholseskydecomposition, seechol.

Recall that the Cholesky factors are from ordered matrices.

chol2inv(x) is a faster way tosolve(x).

Note

There is intentionally noS3 distinction between the classesspam andspam.chol.method.

Author(s)

Reinhard Furrer, based on Ng and Peyton (1993) Fortran routines

References

See references inchol.

See Also

chol.spam andordering.

Examples

# Generate multivariate form a covariance inverse:# (usefull for GRMF)set.seed(13)n <- 25    # dimensionN <- 1000  # sample sizeSigmainv <- .25^abs(outer(1:n,1:n,"-"))Sigmainv <- as.spam( Sigmainv, eps=1e-4)Sigma <- solve( Sigmainv)  # for verificationiidsample <- array(rnorm(N*n),c(n,N))mvsample <- backsolve( chol(Sigmainv), iidsample)norm( var(t(mvsample)) - Sigma)# compare with:mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample, n)   #### ,n as patchnorm( var(t(mvsample)) - Sigma)# 'solve' step by step:b <- rnorm( n)R <- chol(Sigmainv)norm( backsolve( R, forwardsolve( R, b))-      solve( Sigmainv, b) )norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma )# 'update':R1 <- update( R, Sigmainv + diag.spam( n))

Class "spam.chol.NgPeyton"

Description

Result of a Cholesky decomposition with theNgPeytonmethod

Details

It is not possible to directly change the length, dimension and thediagonal entries of a"spam.chol.NgPeyton" object.

Objects from the Class

Objects are created by calls of the formchol(x,method="NgPeyton", ...)and should not be created directly with anew("spam.chol.NgPeyton", ...) call.
At present, no proper print method is defined. However, the factor canbe transformed into aspam object.

Methods

as.matrix

signature(x = "spam.chol.NgPeyton"): Transform the factorinto a regular matrix.

as.spam

signature(x = "spam.chol.NgPeyton"): Transform the factorinto aspam object.

backsolve

signature(r = "spam.chol.NgPeyton"): solvinga triangular system, seesolve.

forwardsolve

signature(l = "spam.chol.NgPeyton"): solvinga triangular system, seesolve.

c

signature(x = "spam.chol.NgPeyton"): Coerce the factor into a vector.

determinant

signature(x = "spam.chol.NgPeyton"):Calculates the determinant from the factor, see alsodet.

diag

signature(x = "spam.chol.NgPeyton"): Extracts thediagonal entries.

dim

signature(x = "spam.chol.NgPeyton"): Retrieve thedimension. Note that"dim<-" is not implemented.

display

signature(x = "spam.chol.NgPeyton"): Transformationto aspam object and display, see alsodisplay.

image

signature(x = "spam.chol.NgPeyton"): Transformationto aspam object and display, see alsoimage.

length

signature(x = "spam.chol.NgPeyton"): Retrieve thedimension. Note that"length<-" is not implemented.

ordering

signature(x = "spam.chol.NgPeyton"):Retrieves the ordering, inordering.

print

signature(x = "spam.chol.NgPeyton"): Short description.

show

signature(object = "spam.chol.NgPeyton"): Short description.

summary

signature(object = "spam.chol.NgPeyton"):Description of the factor, returns (as a list)nnzR,nnzcolindices,the density of the factordensity, and fill-in ratiofillin. For the use of the first two, see ‘Examples’inchol.

t

signature(x = "spam.chol.NgPeyton"): Transformationto aspam object and transposition.

chol

signature(x = "spam.chol.NgPeyton"): Returnsx unchanged.

Author(s)

Reinhard Furrer

References

Ng, E. G. and B. W. Peyton (1993), "Block sparse Cholesky algorithmson advanced uniprocessor computers",SIAM J. Sci. Comput.,14,pp. 1034-1056.

See Also

print.spamordering andchol

Examples

x <- spam( c(4,3,0,3,5,1,0,1,4),3)cf <- chol( x)cfas.spam( cf)# Modify at own risk...slotNames(cf)

Create Toeplitz Matrices

Description

Creates symmetric and asymmetric Toeplitz matrices.

Usage

toeplitz.spam(x, y = NULL, eps = getOption("spam.eps"))

Arguments

Details

The vectory has to be of the same length asxand its first element is discarded.

Value

The Toeplitz matrix inspam format.

Author(s)

Reinhard Furrer

See Also

toeplitz,circulant.spam

Examples

toeplitz.spam(c(1,.25,0,0,0))

Transform a "spam" Format to Triplets

Description

Returns a list containing the indices and elements of aspam object.

Usage

triplet(x, tri=FALSE)

Arguments

Details

The elements are row (column) first ifx is aspam object (matrix).

Value

A list with elements

Author(s)

Reinhard Furrer

See Also

spam.creation for the inverse operation andforeign for other transformations.

Examples

x <- diag.spam(1:4)x[2,3] <- 5triplet(x)all.equal( spam( triplet(x, tri=TRUE)), x)

Validate a Sparse Matrix

Description

Checks if the sparse matrix has the correct structure.

Usage

validate_spam(object)

Arguments

Value

ReturnsTRUE ifobject is a valid spam objects.

See Also

spam.creation

Examples

validate_spam(spam(1, 20))

Spam Version Information

Description

spam.version is a variable (list) holdingdetailed information about the version ofspam loaded.

spam.Version() provides detailed information about the version ofspamrunning.

Usage

spam.version

Value

spam.version is a list with character-string components

spam.version is a list of class"simple.list" whichhas aprint method.

Author(s)

Reinhard Furrer

See Also

See the R counterpartsR.version.

Examples

spam.version$version.string