similar to: Numerical accuracy of matrix multiplication

>>>>> peter dalgaard <pdalgd at gmail.com> >>>>> on Fri, 16 Sep 2016 13:33:11 +0200 writes: > On 16 Sep 2016, at 12:41 , Alexis Sarda <alexis.sarda at gmail.com> wrote: >> Hello, >> >> while testing the crossprod() function under Linux, I noticed the following: >> >> set.seed(883)

>>>>> Alexis Sarda <alexis.sarda at gmail.com> >>>>> on Tue, 20 Sep 2016 17:33:49 +0200 writes: > I just realized that I was actually using a different random number > generator, could that be a valid reason for the discrepancy? > The code should be: > RNGkind("L'Ecuyer") > set.seed(883) > x <-

I just realized that I was actually using a different random number generator, could that be a valid reason for the discrepancy? The code should be: RNGkind("L'Ecuyer") set.seed(883) x <- rnorm(100) x %*% x - sum(x^2) # equal to 1.421085e-14 Regards, Alexis Sarda. On Tue, Sep 20, 2016 at 5:27 PM, Martin Maechler <maechler at stat.math.ethz.ch > wrote: >

R1.7.1/Win2k: Apologies if this posts twice - the first message was not in plain text. I have looked in help.start() and tried typing "crossprod" and "%*%" into the RGui to get an idea for what R is using as internal algorithms for its matrix computations/manipulations... to no avail. Could someone point me in the direction of some documentation? All I get for

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 +

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 } f2 <- function(a,X){ ignore <-

Dear addresses, I need perform a batch of 10 000 simulations for each of 4 options considered. (The idea is to obtain the parameter estimates in a heteroskedastic linear regression model - with additive or mixed heteroskedasticity - via the Kenward-Roger small-sample adjusted covariance matrix of disturbances). For this purpose I wrote an R program which would capture all possible options (true

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 that crossprod(x) is slower than crossprod(x,x): R>

Full_Name: Bill Simpson Version: 1.2.3 OS: linux Submission from: (NULL) (193.62.250.209) I think I have found an accuracy problem with sum() > a<-1:1000 > b<-1:1000 > > crossprod(a,a)*crossprod(b,b) [,1] [1,] 1.114448e+17 > drop(.Last.value) [1] 1.114448e+17 > #this should give same answer as above > sum(a*a)*sum(b*b) [1] -652010736 #seems to be due to

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 (x,

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

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

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

Dear list members; The code given below corresponds to the PCA-NIPALS (principal component analysis) algorithm adapted from the nipals function in the package chemometrics. The reason for using NIPALS instead of SVD is the ability of this algorithm to handle missing values, but that's a different story. I've been trying to find a way to improve (if possible) the efficiency of the code,

I conducted many tests and I noticed that I lose the domain-join on SMB1 soon as I joined SMB2 in the domain. Step 1: SMB1 "net ads join -Uadministrator" -> OK Step 2: SMB1 "net ads testjoin" -> OK Step 3: SMB2 "net ads join -Uadministrator" -> OK Step 4: SMB2 "net ads testjoin" -> OK Step 5: SMB1 "net ads testjoin" ->

R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 Hello, eigen(symmetric=TRUE) behaves strangely when given complex matrices. The following two lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being

The solve function in r-patched has been changed so that it applies a tolerance when using Lapack routines to calculate the inverse of a matrix or to solve a system of linear equations. A tolerance has always been used with the Linpack routines but not with the Lapack routines in versions 1.7.x and 1.8.0. (You can use the optional argument tol = 0 to override this check for computational

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 <-

Hello !   I’m trying to estimate a system of equation (demand and supply) using the systemfit package.  My program is:   library(systemfit) demand <- tsyud ~ tsyud1 + tsucp + tspo + tssn supply <- tscn ~ tsyn + tsqn + tsksn + tsucp system <- list(demand=eqdemand, learning = eqsupply) labels <- list(demand="eqdemand", learning="eqsupply") inst <- ~ tsupp1 + tsupp2

A new version of pqR is now available at pqR-project.org, which fixes several bugs that are also present in the latest R Core patch release (r66002). A number of bugs found previously during pqR development are also unfixed in the latest R Core release. Here is the list of these bugs that are unfixed in r66002 (including documentation deficiencies), taken from the pqR bug fix and documentation