search for: crossprod

...unction(alpha) diag(as.vector(alpha %*% t(z)), nrow=nrow(z), ncol=nrow(z)) sigma_MIX <- function(alpha) diag((as.vector(alpha %*% t(z)))^2, nrow=nrow(z), ncol=nrow(z)) sigma_EXP <- function(alpha) diag(exp(as.vector(alpha %*% t(z))), nrow=nrow(z), ncol=nrow(z)) fi <- function(sigma) solve(crossprod(x, crossprod(solve(sigma),x))) xs <- function(sigma) crossprod(x,solve(sigma)) sx <- function(sigma) crossprod(solve(sigma),x) fixs <- function(sigma) crossprod(fi(sigma),xs(sigma)) sxfixs <- function(sigma) crossprod(xs(sigma),crossprod(fi(sigma),xs(sigma))) # zav?dza (ncol(z) x 1) v...

I found out the other day that crossprod() will take a single matrix argument; crossprod(x) notionally returns crossprod(x,x). The two forms do not return identical matrices: x <- matrix(rnorm(3000000),ncol=3) M1 <- crossprod(x) M2 <- crossprod(x,x) R> max(abs(M1-M2)) [1] 1.932494e-08 But what really surprised me is...

The following is at least as much out of intellectual curiosity as for practical reasons. On reviewing some code written by novices to R, I came across: crossprod(x, y)[1,1] I thought, "That isn't a very S way of saying that, I wonder what the penalty is for using 'crossprod'." To my surprise the penalty was substantially negative. Handily the client had S-PLUS as well -- there the sign of the penalty was as I had expected, but the...

Hi the manpage says that crossprod(x,y) is formally equivalent to, but faster than, the call 't(x) %*% y'. I have a vector 'a' and a matrix 'A', and need to evaluate 't(a) %*% A %*% a' many many times, and performance is becoming crucial. With f1 <- function(a,X){ ignore <- t(a) %*% X %*% a...

Hi everyone, I do a lot of work with large variance matrices, and I like to use "crossprod" for speed and to keep everything symmetric, i.e. I often compute "crossprod(Q %*% t(A))" for "A %*% Sigma %*% t(A)", where "Sigma" decomposes as "t(Q) %*% Q". I notice in the code that "crossprod", current definition > crossprod function...

Dear R-users, I found a strange problem working with products of two matrices, say: a <- A[i, ] ; crossprod(a) where i is a set of integers selecting rows. When i is empty the result is in a sense random. After some trials the right answer (a matrix of zeros) appears. --------------- Illustration -------------------- R : Copyright 2003, The R Development Core Team Version 1.8.0 (2003-10-08) > A...

Hi I have a square matrix Ainv of size N-by-N where N ~ 1000 I have a rectangular matrix H of size N by n where n ~ 4. I have a vector d of length N. I need X = solve(t(H) %*% Ainv %*% H) %*% t(H) %*% Ainv %*% d and H %*% X. It is possible to rewrite X in the recommended crossprod way: X <- solve(quad.form(Ainv, H), crossprod(crossprod(Ainv, H), d)) where quad.form() is a little function that evaluates a quadratic form: quad.form <- function (M, x){ jj <- crossprod(M, x) return(drop(crossprod(jj, jj))) } QUESTION: how to calculate H %*%...

...copy the bug report (after finishing it) to # your favorite email program and send it to # # r-bugs@r-project.org # ###################################################### # The last line of following code produces a segmentation fault: x <- 1:10 f <- gl(5,2) mns <- tapply(x,f,mean) crossprod(mns) #to get sum of squares of mns. # Of course sum(mns^2) is more straightforward. --please do not edit the information below-- Version: platform = i686-pc-linux-gnu arch = i686 os = linux-gnu system = i686, linux-gnu status = major = 1 minor = 7.1 year = 2003 month = 06 day = 16 lan...

Dear Rdevelopers The background for this email is that I was helping a PhD student to improve the speed of her R code. I suggested to replace calls like t(AA)%*% BB by crossprod(AA,BB) since I expected this to be faster. The surprising result to me was that this change actually made her code slower. > ## Examples : > > AA <- matrix(rnorm(3000*1000),3000,1000) > BB <- matrix(rnorm(3000^2),3000,3000) > system.time(crossprod(AA,BB),gcFirst=TRUE) u...

ABCD are four matrix. A * Inverse((Transpose(A)*Tranpose(B)*B*A+C)) * Transpose(A) * Transpose(B) * D how to write in R in an efficient way? --------------------------------- [[alternative HTML version deleted]]

...<- NULL # matrix of scores and loadings for (h in 1:ncomp) { th <- Xh[, 1] ende <- FALSE while (!ende) { nr <- nr + 1 ph <- t(crossprod(th, Xh) * as.vector(1 / crossprod(th))) ph <- ph * as.vector(1/sqrt(crossprod(ph))) thnew <- t(tcrossprod(t(ph), Xh) * as.vector(1/(crossprod(ph)))) prec <- crossprod(th-thnew) th &l...

hi, the help page for crossprod states that crossprod(A,B) is faster than t(A) %*% B; experimentation certainly bears this out. more alarming is the evidence that crossprod(t(A), B) is faster than A %*% B: on a PII laptop, 128MB memory, win98, R-1.5.0.-patched precompiled (no ATLAS): > A <- matrix(rnorm(250000),500,500) &...

Hi all, I have been searching all sorts of documentation, reference cards, cheat sheets but can't find why R's crossprod(A, B) which is identical to A%*%B does not produce the same as Matlabs cross(A, B) Supposedly both calculate the cross product, and say so, or where do I go wrong? R is only doing sums in the crossprod however, as indicated by (z <- crossprod(1:4)) # = sum(1 + 2^2 + 3^2 + 4^2) in http:/...

> From: Douglas Bates > > How does the Opteron perform on floating point? Can you try something > like > > > mm = matrix(rnorm(1e6), nc = 1e3) > > system.time(crossprod(mm)) > [1] 0.51 0.02 0.53 0.00 0.00 > > system.time(crossprod(mm)) > [1] 0.37 0.03 0.40 0.00 0.00 > > system.time(crossprod(mm)) > [1] 0.38 0.02 0.40 0.00 0.00 > > system.time(crossprod(mm)) > [1] 0.38 0.02 0.40 0.00 0.00 > > (That was with R compiled to use Got...

Hello, Consider this function for generalized ridge regression: gre <- function (X,y,D){ n <- dim(X)[1] p <- dim(X)[2] intercept <- rep(1, n) X <- cbind(intercept, X) X2D <- crossprod(X,X)+ D Xy <- crossprod(X,y) bth <- qr.solve(X2D, Xy) } # suppose X is an (nxp) design matrix and y is an (nx1) response vector p <- dim(x)[2] D<- diag(rep(1.5,p)) bt <- gre(X,y,D) I am getting following error: Error in crossprod(X, X) + D : non-conformable arrays But when i...

...for to solve (X' V^{-1} X) ^{-1} X' V^{-1} y for the point estimates and (X' V^{-1} X) ^{-1} for the standard errors. score<-long$score X.mat<-model.matrix(score~time, long, row.names=F) var.mat<-getVarCov(fm1) I<-diag(sample.size) V <- kronecker(I,var.mat) pe<-solve(crossprod(X.mat,solve(V,X.mat)))%*%crossprod(X.mat,solve(V,sco re)) varcov<-solve(crossprod(X.mat,solve(V,X.mat))) This perfectly replicates the gls point estimates, but does not replicate the standard errors. The reason I am concerned is that this does not occur when I do this for a mixed model or for a...

...;- matrix(y,ncol=1) Design.data <- cbind(X, Y) for (j in 1:nrow(X)) { result <- vector("list", ) for( i in 1: 3100) { data <- Design.data[sample(50,50,replace=TRUE),] ##### and unique(Design.data.....) dataX <- data[,1:5] dataY <- data[,6] B.cap.simulation <- solve(crossprod(dataX)) %*% crossprod(dataX, dataY) P.simulation <- dataX %*% solve(crossprod(dataX)) %*% t(dataX) Y.cap.simulation <- P.simulation %*% dataY e.simulation <- dataY - Y.cap.simulation dX.simulation <- nrow(dataX) - ncol(dataX) var.cap.simulation <- crossprod(e.simulation) / (dX.simula...

...m/en-us/mkl ) With the last version of the MKL libraries Intel 2018, we are facing to an issue with *all INTERNAL command* that are executing in R. The R console is freezing executing a process at 100% and never stop!!! It?s really an issue for us. As example, we can reproduce the error with *crossprod. Crossprod *which is a wrapper of BLAS GEMM (optimized with MKL libraries), in this function it seems that variables are not protected ( PROTECT(); UNPROTECT() ), see the screenshot below, which is a recommendation for external commands: Picture1 *RECOMMANDATION* *Picture2* *Code of CROSSPRO...

...ternative ways to 'build' this function are as in f and g below: coeff<-sqrt(1:4) zeta<-function(i) {force(i); function(u){ zeta<-sqrt(2)+sin(pi*i*u) }} f<-function(u){ f<-0 for (j in 1:4){ f<-f+coeff[j]*(zeta(j)(u)) } f<-f } g<-function(u){ g<-crossprod(coeff,zeta(1:4)(u)) } Obviously, f and g are equivalent, but in the actual code I am writing, g is much faster. However, I can't seem to integrate (neither to plot) g. In fact, these are the error messages I get: > print(f(.1)) [1] 4.443642 > print(g(.1)) [,1] [1,] 4.443642 Ev...

I cleaned up my old benchmarking code and added checks for missing data to compare various ways of finding OLS regression coefficients. I thought I would share this for others. the long and short of it is that I would recommend ols.crossprod = function (y, x) { x <- as.matrix(x) ok <- (!is.na(y))&(!is.na(rowSums(x))) y <- y[ok]; x <- subset(x, ok) x <- cbind( 1, x) XtX <- crossprod(x) Xty <- crossprod(x, y) solve(...