[R] proportional matrix rows

[Peter Dalgaard]

Gabor Grothendieck <ggrothendieck at myway.com> writes:

> Robin Hankin <r.hankin <at> soc.soton.ac.uk> writes:
> 
> : 
> : Hi
> : 
> : I have a two-column integer matrix like this:
> : 
> : R> jj
> : 
> :        [,1] [,2]
> :   [1,]   -1    1
> :   [2,]   -2    2
> :   [3,]   -7    6
> :   [4,]   -8    7
> :   [5,]   -6    5
> :   [6,]   -9    8
> :   [7,]   -5    4
> :   [8,]    3   -3
> :   [9,]  -10    9
> : [10,]   -4    3
> : 
> : I want a diagnostic that detects whether a row is a multiple of
> : the first row or not.  In this case, this would be rows 1,2, and 8.
> : 
> : How to do this nicely?  Sometimes the first row has a zero, so the
> : method would have to work on
> : 
> :       [,1] [,2]
> : [1,]    0    1
> : [2,]    1    1
> : [3,]    0    8
> : [4,]    0   -4
> : [5,]    0    0
> : [6,]    4    0
> :  >
> : 
> : in which case rows 1,3,4,5 are multiples.  It'd be nice to have
> : a solution that works for any number of columns.
> 
> 
> If rowi is a multiple of row1 then 
> crossprod(cbind(row1, rowi)) is singular so:
> 
> apply(mat, 1, function(x) det(crossprod(cbind(x, mat[1,])))) == 0
> 
> If your matrix only has two columns then crossprod(x) is singular
> iff x is so you can eliminate the crossprod.

..and in that case it might be easier just to calculate the dot
product with the orthogonal of mat[1,]:

> z <- jj[1,]
> jj %*% c(z[1], -z[2])
      [,1]
 [1,]    0
 [2,]    0
 [3,]    1
 [4,]    1
 [5,]    1
 [6,]    1
 [7,]    1
 [8,]    0
 [9,]    1
[10,]    1

(Notice that the value "1" just comes from Robin being singularly bad
at selecting "arbitrary" rows that are not proportional to the first
one!)

This generalizes to more than two columns as

jj %*% (diag(ncol(jj)) - z%*%solve(crossprod(z))%*%t(z)) 

having rows that are identically zero (the multiplier matrix is
obviously rank-deficient so you can leave out one column, but you need
to guard against leaving in a column of zeroes, so it is hardly worth
it).

In all cases you probably need some sort of guard against round-off.
I.e. something like

> rowSums(zapsmall(jj %*% (diag(2) - z%*%solve(crossprod(z),t(z))))^2) == 0
 [1]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE

or

> rowSums((jj %*% (diag(2) - z%*%solve(crossprod(z),t(z))))^2) < 1e-7
 [1]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907