Package 'fastmatrix'

Description

Small set of functions designed to speed up the computation of certain matrix operations that are commonly used in statistics and econometrics. It provides efficient implementations for the computation of several structured matrices, matrix decompositions and statistical procedures, many of which have minimal memory overhead. Furthermore, the package provides interfaces to C code callable by another C code from other R packages.

Details

The fastmatrix package provides functions to the efficient construction of duplication, commutation and symmetrizer matrices with minimal storage requeriments. Common matrix decompositions (e.g. LU, LDL), rank-1 updates (e.g. Cholesky update), iterative solvers for linear systems and other linear algebra utilities, as well as basic matrix operations, such as Hadamard (elementwise) and Kronecker products, the Sherman-Morrison formula and the power method. The package also offers several statistical procedures, namely: the sweep operator, weighted mean and covariance (using online algorithms), ordinary least squares via multiple strategies (Cholesky, QR, SVD, sweep and conjugate gradients), ridge regression (with procedures for selecting the ridge parameter), omnibus tests for univariate normality, multivariate skewness and kurtosis measures, Mahalanobis distance (with positive-definiteness checks), Wilson-Hilferty transformation of gamma variables (useful, for instance, for goodness-of-fit of multivariate normal data), and some random number generators. Finally, the package provides interfaces for code written in C, enabling other R packages (or user-written C code) to access the C routines in the fastmatrix package.

Author(s)

Felipe Osorio [email protected], Alonso Ogueda [email protected]

Description

Multiplication of 3-dimensional arrays was first introduced by Bates and Watts (1980). More extensions and technical details can be found in Wei (1998).

Usage

array.mult(a, b, x)

Arguments

Details

Let $\bold{X} = (x_{tij})$ be a 3-dimensional $n\times p\times q$ where indices $t, i$ and $j$ indicate face, row and column, respectively. The product $\bold{Y} = \bold{AXB}$ is an $n\times r\times s$ array, with $\bold{A}$ and $\bold{B}$ are $r\times p$ and $q\times s$ matrices respectively. The elements of $\bold{Y}$ are defined as:

$y_{tkl} = \sum\limits_{i=1}^p\sum\limits_{j=1}^q a_{ki}x_{tij}b_{jl}$

Value

array.mult returns a 3-dimensional array of dimension $n\times r\times s$.

References

Bates, D.M., Watts, D.G. (1980). Relative curvature measures of nonlinearity. Journal of the Royal Statistical Society, Series B 42, 1-25.

Wei, B.C. (1998). Exponential Family Nonlinear Models. Springer, New York.

See Also

array, matrix, bracket.prod.

Examples

x <- array(0, dim = c(2,3,3)) # 2 x 3 x 3 array
x[,,1] <- c(1,2,2,4,3,6)
x[,,2] <- c(2,4,4,8,6,12)
x[,,3] <- c(3,6,6,12,9,18)

a <- matrix(1, nrow = 2, ncol = 3)
b <- matrix(1, nrow = 3, ncol = 2)

y <- array.mult(a, b, x) # a 2 x 2 x 2 array
y

Description

Force a square matrix $\bold{X}$ to be symmetric

Usage

asSymmetric(x, lower = TRUE)

Arguments

Value

a square symmetric matrix.

Examples

a <- matrix(1:16, ncol = 4)
isSymmetric(a) # FALSE
a <- asSymmetric(a) # copy lower triangle into upper triangle

Description

Computes the Bezier curve based on $n+1$ control points using the De Casteljau's method.

Usage

bezier(x, y, ngrid = 200)

Arguments

Details

Given $\bold{p}_0,\bold{p}_1,\dots,\bold{p}_n$ control points the Bezier curve is given by $B(t)$ defined as

$B(t) = \left({\begin{array}{c} x(t) \\ y(t) \end{array}}\right) = \sum\limits_{k=0}^n {n \choose k} t^k (1 - t)^k\bold{p}_k$

where $t\in[0,1]$. To evaluate the Bezier curve the De Casteljau's method is used.

Value

A list containing xgrid and ygrid elements used to plot the Bezier curve.

Examples

# a tiny example
x <- c(1.0, 0.25, 1.25, 2.5, 4.00, 5.0)
y <- c(0.5, 2.00, 3.75, 4.0, 3.25, 1.0)
plot(x, y, type = "o")
z <- bezier(x, y, ngrid = 50)
lines(z$xgrid, z$ygrid, lwd = 2, lty = 2, col = "red")

# other simple example
x <- c(4,6,4,5,6,7)
y <- 1:6
plot(x, y, type = "o")
z <- bezier(x, y, ngrid = 50)
lines(z$xgrid, z$ygrid, lwd = 2, lty = 2, col = "red")

Description

Bracket product of a matrix and a 3-dimensional array.

Usage

bracket.prod(a, x)

Arguments

Details

Let $\bold{X} = (x_{tij})$ be a 3-dimensional $n\times p\times q$ array and $\bold{A}$ an $m\times n$ matrix, then $\bold{Y} = [\bold{A}][\bold{X}]$ is called the bracket product of $\bold{A}$ and $\bold{X}$, that is an $m \times p\times q$ with elements

$y_{tij} = \sum\limits_{k=1}^n a_{tk}x_{kij}$

Value

bracket.prod returns a 3-dimensional array of dimension $m\times p\times q$.

References

Wei, B.C. (1998). Exponential Family Nonlinear Models. Springer, New York.

See Also

array, matrix, array.mult.

Examples

x <- array(0, dim = c(2,3,3)) # 2 x 3 x 3 array
x[,,1] <- c(1,2,2,4,3,6)
x[,,2] <- c(2,4,4,8,6,12)
x[,,3] <- c(3,6,6,12,9,18)

a <- matrix(1, nrow = 3, ncol = 2)

y <- bracket.prod(a, x) # a 3 x 3 x 3 array
y

Description

Calculates Lin's concordance correlation coefficient for evaluating the degree of agreement between measurements generated by two different methods.

Usage

ccc(x, data, method = "z-transform", level = 0.95, equal.means = FALSE, 
    ustat = TRUE, subset, na.action)

Arguments

Value

A list with class 'ccc' containing the following named components:

References

Bland, J., Altman, D. (1986). Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet 327, 307-310.

King, T.S., Chinchilli, V.M. (2001). A generalized concordance correlation coefficient for continuous and categorical data. Statistics in Medicine 20, 2131-2147.

King, T.S., Chinchilli, V.M. (2001). Robust estimators of the concordance correlation coefficient. Journal of Biopharmaceutical Statistics 11, 83-105.

Lin, L. (1989). A concordance correlation coefficient to evaluate reproducibility. Biometrics 45, 255-268.

Lin, L. (2000). A note on the concordance correlation coefficient. Biometrics 56, 324-325.

Vallejos, R., Osorio, F., Ferrer, C. (2025). A new coefficient to measure agreement between two continuous variables. doi: 10.48550/arXiv.2507.07913.

Examples

## data in Fig.1 from Bland and Altman (1986).
x <- list(Large = c(494,395,516,434,476,557,413,442,650,433,
          417,656,267,478,178,423,427),
          Mini  = c(512,430,520,428,500,600,364,380,658,445,
          432,626,260,477,259,350,451))
x <- as.data.frame(x)

plot(Mini ~ Large, data = x, xlim = c(100,800), ylim = c(100,800),
     xlab = "PERF by Large meter", ylab = "PERF by Mini meter")
abline(c(0,1), col = "gray", lwd = 2)

## estimating CCC
z <- ccc(~ Mini + Large, data = x, method = "asymp")
z
## output:
# Call:
# ccc(x = ~ Mini + Large, data = x, method = "asymp")
#
# Coefficients:
#  estimate  variance  accuracy precision 
#   0.9427    0.0008    0.9994    0.9433 
#
# Asymptotic 95% confidence interval:
#     CCC     SE  lower  upper 
#  0.9427 0.0286 0.8867 0.9988

Description

Conjugate gradients (CG) method is an iterative algorithm for solving linear systems with positive definite coefficient matrices.

Usage

cg(a, b, maxiter = 200, tol = 1e-7)

Arguments

Value

a vector with the approximate solution, the iterations performed are returned as the attribute 'iterations'.

Warning

The underlying C code does not check for symmetry nor positive definitiveness.

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

Hestenes, M.R., Stiefel, E. (1952). Methods of conjugate gradients for solving linear equations. Journal of Research of the National Bureau of Standards 49, 409-436.

See Also

jacobi, seidel, solve

Examples

a <- matrix(c(4,3,0,3,4,-1,0,-1,4), ncol = 3)
b <- c(24,30,-24)
z <- cg(a, b)
z # converged in 3 iterations

Description

Density, distribution function, quantile function and random generation for the chi distribution with df degrees of freedom.

Usage

dchi(x, df = 1, log = FALSE)
pchi(q, df = 1, lower.tail = TRUE, log.p = FALSE)
qchi(p, df = 1, lower.tail = TRUE, log.p = FALSE)
rchi(n, df = 1)

Arguments

Details

If df is not specified, they assume the default value of 1.

The chi distribution with $\nu$ degrees of freedom (also known as shape parameter) has density

$f(x) = \frac{1}{2^{\nu/2-1}\Gamma(\nu/2)} x^{\nu-1}\exp(-x^2/2),$

where $y > 0$, and $\nu > 0$.

rchi implements the routine proposed by Monahan (1987) for generating chi-distributed random variables.

Value

dchi, pchi, and qchi are respectively the density, distribution function and quantile function of the chi distribution. rchi generates random deviates drawn from the chi distribution, the length of the result is determined by n.

Author(s)

Felipe Osorio

References

Forbes, C., Evans, M., Hastings, N., Peacock, B. (2010). Statistical Distributions, 4th Ed. Wiley, New York.

Johnson, N.L., Kotz, S., Balakrishnan, N. (1994). Continuous Univariate Distributions, Volume 1, 2nd Ed. Wiley, New York.

Monahan, J.F. (1987). An algorithm for generating chi random variables. ACM Transactions on Mathematical Software 13, 168-172.

See Also

Distributions for other standard distributions.

Examples

x <- rchi(1000, df = 2)
## QQ-plot for chi data against true theoretical distribution:
qqplot(qchi(ppoints(1000), df = 2), x, main = "chi QQ-plot",
  xlab = "Theoretical quantiles", ylab = "Sample quantiles")
abline(c(0,1), lwd = 2, lty = 2, col = "red")

Description

function cholupdate returns the upper triangular Cholesky factor of $\bold{A} + \bold{xx}^T$, where $\bold{R}$ is the original Cholesky factor of $\bold{A} = \bold{R}^T\bold{R}$, with $\bold{x}$ a column vector of appropriate dimension.

Usage

cholupdate(r, x)

Arguments

References

Golub, G.H., Van Loan, C.F. (2013). Matrix Computations, 4th Edition. John Hopkins University Press.

See Also

chol

Examples

a <- matrix(c(1,1,1,1,2,3,1,3,6), ncol = 3)
r <- chol(a)
x <- c(0,0,1)
b <- a + outer(x,x)
r1 <- cholupdate(r, x)
r1
all(r1 == chol(b)) # TRUE

Description

Forms a symmetric circulant matrix using a backwards shift of its first column

Usage

circulant(x)

Arguments

Value

A symmetric circulant matrix.

Examples

x <- c(2,3,5,7,11,13)
circulant(x)

Description

This function provides the minimum information required to create the commutation matrix.

The commutation matrix is a square matrix of order $mn$ that, for an $m\times n$ matrix $\bold{A}$, transform vec$(\bold{A}$) to vec$(\bold{A}^T)$.

Usage

comm.info(m = 1, n = m, condensed = TRUE)

Arguments

Details

This function returns a list containing two vectors that represent an element of the commutation matrix and is accesed by the indexes in vectors row and col. This information is used by function comm.prod to do some operations involving the commutation matrix without forming it. This information also can be obtained using function commutation.

Value

A list containing the following elements:

References

Magnus, J.R., Neudecker, H. (1979). The commutation matrix: some properties and applications. The Annals of Statistics 7, 381-394.

See Also

commutation, comm.prod

Examples

z <- comm.info(m = 3, n = 2, condensed = FALSE)
z # where are the ones in commutation matrix of order '3,2'?

K32 <- commutation(m = 3, n = 2, matrix = TRUE)
K32 # only recommended if m and n are very small

Description

Given the row and column dimensions of a commutation matrix $\bold{K}$ of order $mn$ and a conformable matrix $\bold{x}$, performs one of the matrix-matrix operations:

• $\bold{Y} = \bold{KX}$, if side = "left" and transposed = FALSE, or

• $\bold{Y} = \bold{K}^T\bold{X}$, if side = "left" and transposed = TRUE, or

• $\bold{Y} = \bold{XK}$, if side = "right" and transposed = FALSE, or

• $\bold{Y} = \bold{XK}^T$, if side = "right" and transposed = TRUE.

The main aim of comm.prod is to do this matrix multiplication without forming the commutation matrix.

Usage

comm.prod(m = 1, n = m, x = NULL, transposed = FALSE, side = "left")

Arguments

Details

Underlying Fortran code only uses information provided by comm.info to performs the matrix multiplication. The commutation matrix is never created.

See Also

commutation

Examples

K42 <- commutation(m = 4, n = 2, matrix = TRUE)
x <- matrix(1:24, ncol = 3)
y <- K42 %*% x

z <- comm.prod(m = 4, n = 2, x) # K42 is not stored
all(z == y) # matrices y and z are equal!

Description

This function returns the commutation matrix of order $mn$ which transforms, for an $m\times n$ matrix $\bold{A}$, vec$(\bold{A})$ to vec$(\bold{A}^T)$.

Usage

commutation(m = 1, n = m, matrix = FALSE, condensed = FALSE)

Arguments

Details

This function is a wrapper function for the function comm.info. This function provides the minimum information required to create the commutation matrix. If option matrix = FALSE the commutation matrix is stored in two vectors containing the coordinate list of indexes for rows and columns. Option condensed = TRUE only returns vector of indexes for the rows of commutation matrix.

Warning: matrix = TRUE is not recommended, unless the order m and n be small. This matrix can require a huge amount of storage.

Value

Returns an $mn$ by $mn$ matrix (if requested).

References

Magnus, J.R., Neudecker, H. (1979). The commutation matrix: some properties and applications. The Annals of Statistics 7, 381-394.

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

comm.info

Examples

z <- commutation(m = 100, condensed = TRUE)
object.size(z) # 40.6 Kb of storage

z <- commutation(m = 100, condensed = FALSE)
object.size(z) # 80.7 Kb of storage

K100 <- commutation(m = 100, matrix = TRUE) # time: < 2 secs
object.size(K100) # 400 Mb of storage, do not request this matrix!

# a small example
K32 <- commutation(m = 3, n = 2, matrix = TRUE)
a <- matrix(1:6, ncol = 2)
v <- K32 %*% vec(a)
all(vec(t(a)) == as.vector(v)) # vectors are equal!

Description

This function is a constructor for the corAR1 correlation matrix representing an autocorrelation structure of order 1.

Usage

corAR1(rho, p = 2)

Arguments

Value

Returns a $p$ by $p$ matrix.

Examples

R <- corAR1(rho = 0.8, p = 5)

Description

This function is a constructor for the corCS correlation matrix representing a compound symmetry structure corresponding to uniform correlation.

Usage

corCS(rho, p = 2)

Arguments

Value

Returns a $p$ by $p$ matrix.

Examples

R <- corCS(rho = 0.8, p = 5)

Description

Returns a list containing the mean and covariance matrix of the data.

Usage

cov.MSSD(x)

Arguments

Details

This procedure uses the Holmes-Mergen method using the difference between each successive pairs of observations also known as Mean Square Successive Method (MSSD) to estimate the covariance matrix, which is given by

$\bold{S}_{HD} = \frac{1}{2(n-1)} \sum\limits_{i=2}^n (\bold{x}_i - \bold{x}_{i-1})(\bold{x}_i - \bold{x}_{i-1})^T.$

Value

A list containing the following named components:

References

Holmes, D.S., Mergen, A.E. (1993). Improving the performance of the $T^2$ control chart. Quality Engineering 5, 619-625.

See Also

cov and var.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z0 <- cov(x)
z0
z1 <- cov.MSSD(x)
z1

Description

Returns a list containing estimates of the weighted mean and covariance matrix of the data.

Usage

cov.weighted(x, weights = rep(1, nrow(x)))

Arguments

Details

The covariance matrix is divided by the number of observations, which arise for instance, when we use the class of elliptical contoured distributions. Thus,

$W_n = \sum\limits_{i=1}^n w_i, \qquad \overline{\bold{x}}_n = \frac{1}{W_n} \sum\limits_{i=1}^n w_i\bold{x}_i \qquad \bold{S}_n = \frac{1}{n} \sum\limits_{i=1}^n w_i (\bold{x}_i - \overline{\bold{x}}_n)(\bold{x}_i - \overline{\bold{x}}_n)^T.$

This differs from the behaviour of function cov.wt.

Value

A list containing the following named components:

References

Clarke, M.R.B. (1971). Algorithm AS 41: Updating the sample mean and dispersion matrix. Applied Statistics 20, 206-209.

See Also

cov.wt, cov and var.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z0 <- cov.weighted(x) # all weights are 1
D2 <- Mahalanobis(x, center = z0$mean, cov = z0$cov)
p <- ncol(x)
wts <- (p + 1) / (1 + D2) # nice weights!
z1 <- cov.weighted(x, weights = wts)
z1

Description

Given the order of two duplication matrices and a conformable matrix $\bold{X}$, this function performs the operation: $\bold{Y} = \bold{D}_n^T\bold{X}\bold{D}_k$, where $\bold{D}_n$ and $\bold{D}_k$ are duplication matrices of order $n$ and $k$, respectively.

Usage

dupl.cross(n = 1, k = n, x = NULL)

Arguments

Details

This function calls dupl.prod to performs the matrix multiplications required but without forming any duplication matrices.

See Also

dupl.prod

Examples

D2 <- duplication(n = 2, matrix = TRUE)
D3 <- duplication(n = 3, matrix = TRUE)
x <- matrix(1, nrow = 9, ncol = 4)
y <- t(D3) %*% x %*% D2

z <- dupl.cross(n = 3, k = 2, x) # D2 and D3 are not stored
all(z == y) # matrices y and z are equal!

x <- matrix(1, nrow = 9, ncol = 9)
z <- dupl.cross(n = 3, x = x) # same matrix is used to pre- and post-multiplying x
z # print result

Description

This function provides the minimum information required to create the duplication matrix.

Usage

dupl.info(n = 1, condensed = TRUE)

Arguments

Details

This function returns a list containing two vectors that represent an element of the duplication matrix and is accesed by the indexes in vectors row and col. This information is used by function dupl.prod to do some operations involving the duplication matrix without forming it. This information also can be obtained using function duplication

Value

A list containing the following elements:

See Also

duplication, dupl.prod

Examples

z <- dupl.info(n = 3, condensed = FALSE)
z # where are the ones in duplication of order 3?

D3 <- duplication(n = 3, matrix = TRUE)
D3 # only recommended if n is very small

Description

Given the order of a duplication and a conformable matrix $\bold{X}$, performs one of the matrix-matrix operations:

• $\bold{Y} = \bold{DX}$, if side = "left" and transposed = FALSE, or

• $\bold{Y} = \bold{D}^T\bold{X}$, if side = "left" and transposed = TRUE, or

• $\bold{Y} = \bold{XD}$, if side = "right" and transposed = FALSE, or

• $\bold{Y} = \bold{XD}^T$, if side = "right" and transposed = TRUE,

where $\bold{D}$ is the duplication matrix of order $n$. The main aim of dupl.prod is to do this matrix multiplication without forming the duplication matrix.

Usage

dupl.prod(n = 1, x, transposed = FALSE, side = "left")

Arguments

Details

Underlying C code only uses information provided by dupl.info to performs the matrix multiplication. The duplication matrix is never created.

See Also

duplication

Examples

D4 <- duplication(n = 4, matrix = TRUE)
x <- matrix(1, nrow = 16, ncol = 2)
y <- crossprod(D4, x)

z <- dupl.prod(n = 4, x, transposed = TRUE) # D4 is not stored
all(z == y) # matrices y and z are equal!

Description

This function returns the duplication matrix of order $n$ which transforms, for a symmetric matrix $\bold{A}$, vech$(\bold{A})$ into vec$(\bold{A})$.

Usage

duplication(n = 1, matrix = FALSE, condensed = FALSE)

Arguments

Details

This function is a wrapper function for the function dupl.info. This function provides the minimum information required to create the duplication matrix. If option matrix = FALSE the duplication matrix is stored in two vectors containing the coordinate list of indexes for rows and columns. Option condensed = TRUE only returns vector of indexes for the columns of duplication matrix.

Warning: matrix = TRUE is not recommended, unless the order n be small. This matrix can require a huge amount of storage.

Value

Returns an $n^2$ by $n(n + 1)/2$ matrix (if requested).

References

Magnus, J.R., Neudecker, H. (1980). The elimination matrix, some lemmas and applications. SIAM Journal on Algebraic Discrete Methods 1, 422-449.

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

dupl.info

Examples

z <- duplication(n = 100, condensed = TRUE)
object.size(z) # 40.5 Kb of storage

z <- duplication(n = 100, condensed = FALSE)
object.size(z) # 80.6 Kb of storage

D100 <- duplication(n = 100, matrix = TRUE)
object.size(D100) # 202 Mb of storage, do not request this matrix!

# a small example
D3 <- duplication(n = 3, matrix = TRUE)
a <- matrix(c( 1, 2, 3,
               2, 3, 4,
               3, 4, 5), nrow = 3)
upper <- vech(a)
v <- D3 %*% upper
all(vec(a) == as.vector(v)) # vectors are equal!

Description

Equilibrate a rectangular or symmetric matrix using 2-norm.

Usage

equilibrate(x, scale = TRUE)

Arguments

Value

For scale = TRUE, the equilibrated matrix. The scalings and an approximation of the condition number, are returned as attributes "scales" and "condition". If $\bold{X}$ is a rectangular matrix, only the columns are equilibrated.

Examples

x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
z <- equilibrate(x)
apply(z, 2, function(x) sum(x^2)) # all 1

xx <- crossprod(x)
equilibrate(xx)

Description

The Floyd-Warshall algorithm finds all shortest paths (if exist) in a directed graph.

Usage

floyd(x)

Arguments

Value

Returns a list with final costs and the shortest path between two nodes.

References

Floyd, R.W. (1962). Algorithm 97: Shortest Path. Communications of the ACM 5 (6), 345.

Examples

x <- matrix(c(0,3,Inf,5,2,0,Inf,4,Inf,1,0,Inf,Inf,Inf,2,0), nrow = 4, 
            ncol = 4, byrow = TRUE)
z <- floyd(x)
z

Description

This function returns the Frank matrix of order $n$.

Usage

frank(n = 1)

Arguments

Details

A Frank matrix of order $n$ is a square matrix $\bold{F}_n = (f_{ij})$ defined as

$f_{ij} = \left\{ {\begin{array}{ll} n - j + 1, & i \le j, \\ n - j, & i = j + 1, \\ 0, & i \ge j + 2 \end{array}} \right.$

Value

Returns an $n$ by $n$ matrix.

References

Frank, W.L. (1958). Computing eigenvalues of complex matrices by determinant evaluation and by methods of Danilewski and Wielandt. Journal of the Society for Industrial and Applied Mathematics 6, 378-392.

Examples

F5 <- frank(n = 5)
det(F5) # equals 1

Description

It calculates the geometric mean using a Fused-Multiply-and-Add (FMA) compensated scheme for accurate computation of floating-point product.

Usage

geomean(x)

Arguments

Details

If x contains any non-positive values, geomean returns NA and a warning message is displayed.

The geometric mean is a measure of central tendency, which is defined as

$G = \sqrt[n]{x_1 x_2 \ldots x_n} = \Big(\prod_{i=1}^n x_i\Big)^{1/n}.$

This procedure calculates the product required in the geometric mean safely using a compensated scheme as proposed by Graillat (2009).

Value

The geometric mean of the sample, a non-negative number.

References

Graillat, S. (2009). Accurate floating-point product and exponentiation. IEEE Transactions on Computers 58, 994-1000.

Oguita, T., Rump, S.M., Oishi, S. (2005). Accurate sum and dot product. SIAM Journal on Scientific Computing 26, 1955-1988.

See Also

mean, median.

Examples

set.seed(149)
x <- rlnorm(1000)
mean(x)    # 1.68169
median(x)  # 0.99663
geomean(x) # 1.01688

Description

This function returns the Hadamard or element-wise product of two matrices $\bold{X}$ and $\bold{Y}$, that have the same dimensions.

Usage

hadamard(x, y = x)

Arguments

Value

A matrix with the same dimension of $\bold{X}$ (and $\bold{Y}$) which corresponds to the element-by-element product of the two matrices.

References

Styan, G.P.H. (1973). Hadamard products and multivariate statistical analysis, Linear Algebra and Its Applications 6, 217-240.

Examples

x <- matrix(rep(1:10, times = 5), ncol = 5)
y <- matrix(rep(1:5, each = 10), ncol = 5)
z <- hadamard(x, y)
z

Description

Forms a symmetric Hankel matrix of order $n$ from the values in vector $\bold{x}$ and optionally the vector $\bold{y}$.

Usage

hankel(x, y = NULL)

Arguments

Value

A symmetric Hankel matrix of order $n$.

Examples

x <- 1:4
y <- c(4,6,8,10)

# H4
hankel(x)

# H({1,2,3,4},{4,6,8,10})
hankel(x, y)

Description

Performs large-sample methods for testing equality of $p \ge 2$ correlated variables.

Usage

harris.test(x, test = "Wald")

Arguments

Value

A list of class 'harris.test' with the following elements:

References

Harris, P. (1985). Testing the variance homogeneity of correlated variables. Biometrika 72, 103-107.

Examples

x <- iris[,1:4]
z <- harris.test(x, test = "robust")
z

Description

This function returns the Helmert matrix of order $n$.

Usage

helmert(n = 1)

Arguments

Details

A Helmert matrix of order $n$ is a square matrix defined as

$\bold{H}_n = \left[ {\begin{array}{ccccc} 1/\sqrt{n} & 1/\sqrt{n} & 1/\sqrt{n} & \dots & 1/\sqrt{n} \\ 1/\sqrt{2} & -1/\sqrt{2} & 0 & \dots & 0 \\ 1/\sqrt{6} & 1/\sqrt{6} & -2/\sqrt{6} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \dots & -\frac{(n-1)}{\sqrt{n(n-1)}} \end{array}} \right].$

Helmert matrix is orthogonal and is frequently used in the analysis of variance (ANOVA).

Value

Returns an $n$ by $n$ matrix.

References

Lancaster, H.O. (1965). The Helmert matrices. The American Mathematical Monthly 72, 4-12.

Gentle, J.E. (2007). Matrix Algebra: Theory, Computations, and Applications in Statistics. Springer, New York.

Examples

n <- 1000
set.seed(149)
x <- rnorm(n)

H <- helmert(n)
object.size(H) # 7.63 Mb of storage
K <- H[2:n,]
z <- c(K %*% x)
sum(z^2) # 933.1736

# same that
(n - 1) * var(x)

Description

This function returns the Householder matrix (also called Householder reflection) or the Householder vector, which are constructed based on a nonzero vector $\bold{x}$.

Usage

house(x, matrix = FALSE)

Arguments

Details

A Householder transformation is a rank-1 modification of the identity matrix which can be used to zero out selected elements of a vector. A Householder matrix $\bold{P}$ adopts the form

$\bold{P} - \tau\bold{uu}^T,$

where $\bold{u}$ is called Householder vector, and $\tau = 2/(\bold{u}^T\bold{u}).$

Value

If matrix = TRUE an $n$ by $n$ matrix is returned, otherwise house function return a list containing the following components:

References

Golub, G.H., van Loan, C.F. (1996). Matrix Computations, 3rd Edition. The Johns Hopkins University Press, Baltimore.

Examples

x <- c(3,1,5,1)
z <- house(x)
z

mat <- house(x, matrix = TRUE)
v <- mat %*% x
v[1,1] == minkowski(x) # TRUE

Description

This function applies the Householder matrix defined by a nonzero vector $\bold{x}$ to a matrix $\bold{A}$. Thus, the output of this function is the matrix-matrix operation:

• $\bold{Y} = \bold{PA}$, if side = "left", or

• $\bold{Y} = \bold{AP}$, if side = "right".

Usage

house.prod(a, x, side = "left")

Arguments

Details

Underlying code never entail the explicit formation of the Householder matrix.

References

Golub, G.H., van Loan, C.F. (1996). Matrix Computations, 3rd Edition. The Johns Hopkins University Press, Baltimore.

Examples

x <- c(3,1,5,1)
u <- house(x)$u
house.prod(u, x)

y <- as.matrix(x)
house.prod(y, x)

Description

Returns TRUE if the given matrix is lower or upper triangular matrix.

Usage

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

Arguments

Value

Check if a matrix is lower or upper triangular. You can also include diagonal to the check.

See Also

lower.tri, upper.tri

Examples

x <- matrix(rnorm(10 * 3), ncol = 3)
  R <- chol(crossprod(x))

  is.lower.tri(R)
  is.upper.tri(R)

Description

Jacobi method is an iterative algorithm for solving a system of linear equations.

Usage

jacobi(a, b, start, maxiter = 200, tol = 1e-7)

Arguments

Details

Let $\bold{D}$, $\bold{L}$, and $\bold{U}$ denote the diagonal, lower triangular and upper triangular parts of a matrix $\bold{A}$. Jacobi's method solve the equation $\bold{Ax} = \bold{b}$, iteratively by rewriting $\bold{Dx} + (\bold{L} + \bold{U})\bold{x} = \bold{b}$. Assuming that $\bold{D}$ is nonsingular leads to the iteration formula

$\bold{x}^{(k+1)} = -\bold{D}^{-1}(\bold{L} + \bold{U})\bold{x}^{(k)} + \bold{D}^{-1}\bold{b}$

Value

a vector with the approximate solution, the iterations performed are returned as the attribute 'iterations'.

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

seidel

Examples

a <- matrix(c(5,-3,2,-2,9,-1,3,1,-7), ncol = 3)
b <- c(-1,2,3)
start <- c(1,1,1)
z <- jacobi(a, b, start)
z # converged in 15 iterations

Description

Performs an omnibus test for univariate normality.

Usage

JarqueBera.test(x, test = "DH")

Arguments

Value

A list of class 'JarqueBera.test' with the following elements:

References

Doornik, J.A., Hansen, H. (2008). An omnibus test for univariate and multivariate normality. Oxford Bulletin of Economics and Statistics 70, 927-939.

Gel, Y.R., Gastwirth, J.L. (2008). A robust modification of the Jarque-Bera test of normality. Economics Letters 99, 30-32.

Jarque, C.M., Bera, A.K. (1980). Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters 6, 255-259.

Urzua, C.M. (1996). On the correct use of omnibus tests for normality. Economics Letters 53, 247-251.

Examples

set.seed(149)
x <- rnorm(100)
z <- JarqueBera.test(x, test = "DH")
z

set.seed(173)
x <- runif(100)
z <- JarqueBera.test(x, test = "DH")
z

Description

Computes the kronecker product of two matrices, $\bold{X}$ and $\bold{Y}$.

Usage

kronecker.prod(x, y = x)

Arguments

Details

Let $\bold{X}$ be an $m\times n$ and $\bold{Y}$ a $p\times q$ matrix. The $mp\times nq$ matrix defined by

$\left[{\begin{array}{ccc} x_{11}\bold{Y} & \dots & x_{1n}\bold{Y} \\ \vdots & & \vdots \\ x_{m1}\bold{Y} & \dots & x_{mn}\bold{Y} \end{array}}\right],$

is called the Kronecker product of $\bold{X}$ and $\bold{Y}$.

Value

An array with dimensions dim(x) * dim(y).

References

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

kronecker function from base package is based on outer. Our C version is slightly faster.

Examples

# block diagonal matrix:
a <- diag(1:3)
b <- matrix(1:4, ncol = 2)
kronecker.prod(a, b)

# examples with vectors
ones <- rep(1, 4)
y <- 1:3
kronecker.prod(ones, y) # 12-dimensional vector
kronecker.prod(ones, t(y)) # 3 x 3 matrix

Description

Given $\bold{A}$ an $n$ by $n$ real matrix and an $n$-vector $\bold{b}$, this function constructs the Krylov matrix $\bold{K}$, where

$\bold{K} = [\bold{b},\bold{Ab},\dots,\bold{A}^{m-1}\bold{b}].$

Usage

krylov(a, b, m = ncol(a))

Arguments

Value

Returns an $n$ by $m$ matrix.

Examples

a <- matrix(c(1, 3, 2, -5, 1, 7, 1, 5, -4), ncol = 3, byrow = TRUE)
b <- c(1, 1, 1)
k <- krylov(a, b, m = 4)
k

Description

Functions to compute measures of multivariate skewness $(b_{1p})$ and kurtosis $(b_{2p})$ proposed by Mardia (1970),

$b_{1p} = \frac{1}{n^2}\sum\limits_{i=1}^n\sum\limits_{j=1}^n ((\bold{x}_i - \overline{\bold{x}})^T\bold{S}^{-1}(\bold{x}_j - \overline{\bold{x}}))^3,$

and

$b_{2p} = \frac{1}{n}\sum\limits_{i=1}^n ((\bold{x}_i - \overline{\bold{x}})^T \bold{S}^{-1}(\bold{x}_i - \overline{\bold{x}}))^2.$

Usage

kurtosis(x)

skewness(x)

Arguments

References

Mardia, K.V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika 57, 519-530.

Mardia, K.V., Zemroch, P.J. (1975). Algorithm AS 84: Measures of multivariate skewness and kurtosis. Applied Statistics 24, 262-265.

Examples

setosa <- iris[1:50,1:4]
kurtosis(setosa)
skewness(setosa)

Description

Compute the LDL decomposition of a real symmetric matrix.

Usage

ldl(x)

Arguments

Value

The factorization has the form $\bold{X} = \bold{LDL}^T$, where $\bold{D}$ is a diagonal matrix and $\bold{L}$ is unitary lower triangular.

The LDL decomposition of $\bold{x}$ is returned as a list with components:

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

chol

Examples

a <- matrix(c(2,-1,0,-1,2,-1,0,-1,1), ncol = 3)
z <- ldl(a)
z # information of LDL factorization

# computing det(a)
prod(z$d) # product of diagonal elements of D

# a non-positive-definite matrix
m <- matrix(c(5,-5,-5,3), ncol = 2)
try(chol(m)) # fails
ldl(m)

Description

lu computes the LU factorization of a matrix.

Usage

lu(x)
## Default S3 method:
lu(x)

## S3 method for class 'lu'
solve(a, b, ...)

is.lu(x)

Arguments

Details

The LU factorization plays an important role in many numerical procedures. In particular it is the basic method to solve the equation $\bold{Ax} = \bold{b}$ for given matrix $\bold{A}$, and vector $\bold{b}$.

solve.lu is the method for solve for lu objects.

is.lu returns TRUE if x is a list and inherits from "lu".

Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the Fortran code.

Value

The LU factorization of the matrix as computed by LAPACK. The components in the returned value correspond directly to the values returned by DGETRF.

Note

To compute the determinant of a matrix (do you really need it?), the LU factorization is much more efficient than using eigenvalues (eigen). See det.

LAPACK uses column pivoting and does not attempt to detect rank-deficient matrices.

References

Anderson. E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. Sorensen, D. (1999). LAPACK Users' Guide, 3rd Edition. SIAM.

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

extractL, extractU, constructX for reconstruction of the matrices, lu2inv

Examples

a <- matrix(c(3,2,6,17,4,18,10,-2,-12), ncol = 3)
z <- lu(a)
z # information of LU factorization

# computing det(a)
prod(diag(z$lu)) # product of diagonal elements of U

# solve linear equations
b <- matrix(1:6, ncol = 2)
solve(z, b)

Description

Returns the original matrix from which the object was constructed or the components of the factorization.

Usage

constructX(x)
extractL(x)
extractU(x)

Arguments

Value

constructX returns $\bold{X}$, the original matrix from which the lu object was constructed (because of the pivoting the $\bold{X}$ matrix is not exactly the product between $\bold{L}$ and $\bold{U}$).

extractL returns $\bold{L}$. This may be pivoted.

extractU returns $\bold{U}$.

See Also

lu.

Examples

a <- matrix(c(10,-3,5,-7,2,-1,0,6,5), ncol = 3)
z <- lu(a)
L <- extractL(z)
L
U <- extractU(z)
U
X <- constructX(z)
all(a == X)

Description

Invert a square matrix from its LU factorization.

Usage

lu2inv(x)

Arguments

Value

The inverse of the matrix whose LU factorization was given.

Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the Fortran code.

Source

This is an interface to the LAPACK routine DGETRI. LAPACK is from https://netlib.org/lapack/ and its guide is listed in the references.

References

Anderson. E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. Sorensen, D. (1999). LAPACK Users' Guide, 3rd Edition. SIAM.

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

lu, solve.

Examples

a <- matrix(c(3,2,6,17,4,18,10,-2,-12), ncol = 3)
z <- lu(a)
a %*% lu2inv(z)

Description

Returns the squared Mahalanobis distance of all rows in $\bold{x}$ and the vector $\bold{\mu}$ = center with respect to $\bold{\Sigma}$ = cov. This is (for vector $\bold{x}$) defined as

$D^2 = (\bold{x} - \bold{\mu})^T \bold{\Sigma}^{-1} (\bold{x} - \bold{\mu})$

Usage

Mahalanobis(x, center, cov, inverted = FALSE)

Arguments

Details

Unlike function mahalanobis, the covariance matrix is factorized using the Cholesky decomposition, which allows to assess if cov is positive definite. Unsuccessful results from the underlying LAPACK code will result in an error message.

See Also

cov, mahalanobis

Examples

x <- cbind(1:6, 1:3)
xbar <- colMeans(x)
S <- matrix(c(1,4,4,1), ncol = 2) # is negative definite
D2 <- mahalanobis(x, center = xbar, S)
all(D2 >= 0) # several distances are negative

## next command produces the following error:
## Covariance matrix is possibly not positive-definite
## Not run: D2 <- Mahalanobis(x, center = xbar, S)

Description

Performs Mardia's tests to assess multivariate normality based on the multivariate skewness and kurtosis coefficients.

Usage

mardia.test(x)

Arguments

Value

A list of class 'Mardia.test' with the following elements:

References

Mardia, K.V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika 57, 519-530.

Mardia, K.V. (1974). Applications of some measures of multivariate skewness and kurtosis in testing normality and robustness studies. Sankhya 36, 115-128.

Examples

setosa <- iris[1:50,1:4]
z <- mardia.test(setosa)
z

set.seed(149)
Sigma <- matrix(c(10,3,3,2), ncol = 2)
x <- rmnorm(n = 300, Sigma = Sigma)
z <- mardia.test(x)
z

Description

This function computes the matrix function $\bold{F} =f(\bold{A})$ where $\bold{A}$ is upper triangular by applying a Parlett recurrence.

Usage

matrix.fun(a, FUN = "log")

Arguments

Details

The used-defined function FUN is evaluated at the triangular matrix argument. This function can be used in conjunction with Schur decomposition to evaluate the function of a matrix.

References

Higham, N.J. (1986). Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia.

Examples

a <- matrix(c(1,2,3,0,3,4,0,0,5), ncol = 3, byrow = TRUE)
fnc <- function(x) (1 + x) / x
f <- matrix.fun(a, FUN = fnc)
f

a <- matrix(c(-49,24,-64,31), ncol = 2, byrow = TRUE)
z <- schur(a)
m <- z$m
u <- z$vectors
m <- matrix.fun(m, FUN = exp)
u %*% m %*% t(u) # exp(a)

Description

Computes the inner product between two rectangular matrices calling BLAS.

Usage

matrix.inner(x, y = x)

Arguments

Value

a real value, indicating the inner product between two matrices.

Examples

x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
y <- matrix(1, nrow = 6, ncol = 3)
matrix.inner(x, y)

# must be equal
matrix.norm(x, type = "Frobenius")^2
matrix.inner(x)

Description

Computes a matrix norm of x using LAPACK. The norm can be the one ("1") norm, the infinity ("inf") norm, the Frobenius norm, the maximum modulus ("maximum") among elements of a matrix, as determined by the value of type.

Usage

matrix.norm(x, type = "Frobenius")

Arguments

Details

As function norm in package base, method of matrix.norm calls the LAPACK function DLANGE.

Note that the 1-, Inf- and maximum norm is faster to calculate than the Frobenius one.

Value

The matrix norm, a non-negative number.

Examples

# a tiny example
x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
matrix.norm(x, type = "Frobenius")
matrix.norm(x, type = "1")
matrix.norm(x, type = "Inf")

# an example not that small
n <- 1000
x <- .5 * diag(n) + 0.5 * matrix(1, nrow = n, ncol = n)
matrix.norm(x, type = "Frobenius")
matrix.norm(x, type = "1")
matrix.norm(x, type = "Inf")
matrix.norm(x, type = "maximum") # equal to 1

Description

Given $c_0,c_1,\dots,c_n$ coefficients of the polynomial and $\bold{A}$ an $n$ by $n$ matrix. This function computes the matrix polynomial

$\bold{B} = \sum\limits_{k=0}^n c_k\bold{A}^k,$

using Horner's scheme, where $\bold{A}^0 = \bold{I}$ is the $n$ by $n$ identity matrix.

Usage

matrix.polynomial(a, coef = rep(1, power + 1), power = length(coef))

Arguments

Value

Returns an $n$ by $n$ matrix.

Examples

a <- matrix(c(1, 3, 2, -5, 1, 7, 1, 5, -4), ncol = 3, byrow = TRUE)
cf <- c(3, 1, 2)
b <- matrix.polynomial(a, cf)
b # 3 * diag(3) + a + 2 * a %*% a
b <- matrix.polynomial(a, power = 2)
b # diag(3) + a + a %*% a

Description

This function computes a square root of an $n\times n$ matrix $\bold{A}$.

Usage

matrix.sqrt(a, method = "DB", maxiter = 50, tol = 1e-8)

Arguments

Details

A square root of a square matrix $\bold{A}$ is obtained by solving the equation $\bold{X}^2 = \bold{A}$, considering the Newton iteration proposed by Denman and Beavers (1976), or alternatively is based on the Schur decomposition.

References

Denman, E.D., Beavers, A.N. (1976). The matrix sign function and computations in systems. Applied Mathematics and Computation 2, 63-94.

Higham, N.J. (1986). Newton's method for the matrix square root. Mathematics of Computation 46, 537-549.

Higham, N.J. (1986). Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia.

Examples

a <- matrix(c(35,17,3,17,46,11,3,11,12), ncol = 3)
root <- matrix.sqrt(a) # 8 iterations

# just checking
root %*% root

root <- matrix.sqrt(a, method = "schur")

Description

Compute the Cholesky factorization of a real symmetric but not necessarily positive definite matrix.

Usage

mchol(x)

Arguments

Value

The lower triangular factor of modified Cholesky decomposition, i.e., the matrix $\bold{L}$ such that $\bold{X} + \bold{E} = \bold{LL}^T$, where $\bold{E}$ is a nonnegative diagonal matrix that is zero if $\bold{X}$ es sufficiently positive definite.

References

Gill, P.E., Murray, W., Wright, M.H. (1981). Practical Optimization. Academic Press, London.

Nocedal, J., Wright, S.J. (1999). Numerical Optimization. Springer, New York.

See Also

chol, ldl

Examples

# a non-positive-definite matrix
a <- matrix(c(4,2,1,2,6,3,1,3,-.004), ncol = 3)
try(chol(a)) # fails
z <- mchol(a)
z # triangular factor

# modified 'a' matrix
tcrossprod(z)

Description

It calculates the mediancenter (or geometric median) of multivariate data.

Usage

mediancenter(x)

Arguments

Details

The mediancenter for a sample of multivariate observations is computed using a steepest descend method combined with bisection. The mediancenter invariant to rotations of axes and is useful as a multivariate generalization of the median of univariate sample.

Value

A list containing the following named components:

References

Gower, J.C. (1974). Algorithm AS 78: The mediancentre. Applied Statistics 23, 466-470.

See Also

cov.wt, median.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z <- mediancenter(x)$median # degenerate solution
xbar <- colMeans(x)
plot(x, xlab = "", ylab = "")
points(x = xbar[1], y = xbar[2], pch = 16, col = "red")
points(x = z[1], y = z[2], pch = 3, col = "blue", lwd = 2)