[R] Problems with crossprod

[Prof Brian Ripley]

On Fri, 17 Oct 2003, Liaw, Andy wrote:

> Somehow R creates `a' as a matrix with 0 rows and 5 columns.  I don't know
> how crossprod() or other linear algebra functions deals with such a
> degenerate matrix.
> 
> I'd suggest R Core to add checks for strictly positive dimensions in such
> functions.

Rather, we do try to ensure they give the right answers for 0 extents.
The code has lines like

    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
			x, &nrx, y, &nry, &zero, z, &ncx);
    }

and so does nothing for 0 extents!  The corresponding matprod code does, 
so this seems a simple oversight that will be fixed shortly.

> (Also, I find it strange that A[1,] is a vector, but A[numeric(0),] is a 0x5
> matrix...)

That's the point of drop = TRUE (the default): it drops extents of length 
1 (and not 0).


> 
> Andy
> 
> > From: Giovanni Marchetti [mailto:gmm at ds.unifi.it] 
> > 
> > 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  <-matrix(0, 5, 5)
> > > i <- c()
> > > a <- A[i, ] ; crossprod(a)
> >               [,1]          [,2]          [,3]          [,4] [,5]
> > [1,] 6.578187e-313           NaN           NaN           NaN  NaN
> > [2,]           NaN 1.273197e-313           NaN 1.485397e-313  NaN
> > [3,]           NaN 4.243992e-313 2.121996e-314           NaN  NaN
> > [4,]           NaN 1.697597e-313           NaN 4.880590e-313  NaN
> > [5,] 5.941588e-313           NaN           NaN 1.697597e-313  NaN
> > > a <- A[i, ] ; crossprod(a)
> >               [,1]          [,2]          [,3]          [,4]  
> >         [,5]
> > [1,] 2.121996e-314 5.729389e-313           NaN           NaN  
> >          NaN
> > [2,]           NaN           NaN           NaN           NaN 
> > 1.909796e-313
> > [3,] 2.970794e-313           NaN           NaN           NaN  
> >          NaN
> > [4,]           NaN           NaN           NaN 8.487983e-314  
> >          NaN
> > [5,]           NaN 6.365987e-313 2.546395e-313           NaN  
> >          NaN
> > > a <- A[i, ] ; crossprod(a)
> >               [,1]          [,2] [,3]          [,4]          [,5]
> > [1,]           NaN 1.485397e-313  NaN           NaN 2.970794e-313
> > [2,] 3.182994e-313           NaN  NaN 1.060998e-313           NaN
> > [3,]           NaN           NaN  NaN 1.697597e-313 2.737375e-312
> > [4,]           NaN           NaN  NaN           NaN  2.048394e+10
> > [5,]           NaN           NaN  NaN           NaN 2.970794e-313
> > > a <- A[i, ] ; crossprod(a)
> >               [,1]        [,2] [,3] [,4] [,5]
> > [1,] 1.591383e-266 20489834629    0    0    0
> > [2,] 5.031994e-266           0    0    0    0
> > [3,] 1.591205e-266           0    0    0    0
> > [4,] 1.264128e-267           0    0    0    0
> > [5,] 1.037656e-311           0    0    0    0
> > > a <- A[i, ] ; crossprod(a)
> >      [,1] [,2] [,3] [,4] [,5]
> > [1,]    0    0    0    0    0
> > [2,]    0    0    0    0    0
> > [3,]    0    0    0    0    0
> > [4,]    0    0    0    0    0
> > [5,]    0    0    0    0    0
> > --------------- End of illustration------------
> > 
> > The same problem does not appear using the matrix product:
> > 
> > > a <- A[i, ] ; t(a) %*% a
> >      [,1] [,2] [,3] [,4] [,5]
> > [1,]    0    0    0    0    0
> > [2,]    0    0    0    0    0
> > [3,]    0    0    0    0    0
> > [4,]    0    0    0    0    0
> > [5,]    0    0    0    0    0
> > 
> > Note that Splus 6 returns an error message:
> > 
> > > a <- A[i, ] ; crossprod(a)
> > 
> > Problem in .Fortran.ok.Internal(if(cmplx) "zcrossp1"..: 
> > subroutine dcrossp1: 
> > Argument 1 has zero length
> > 
> > 
> > Thank you, 
> > 
> > Giovanni
> > 
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list 
> > https://www.stat.math.ethz.ch/mailman/listinfo> /r-help
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595