'What is R's crossproduct function?
I feel stupid asking, but what is the intent of R's crossprod function with respect to vector inputs? I wanted to calculate the cross-product of two vectors in Euclidean space and mistakenly tried using crossprod .
One definition of the vector cross-product is N = |A|*|B|*sin(theta) where theta is the angle between the two vectors. (The direction of N is perpendicular to the A-B plane). Another way to calculate it is N = Ax*By - Ay*Bx .base::crossprod clearly does not do this calculation, and in fact produces the vector dot-product of the two inputs sum(Ax*Bx, Ay*By).
So, I can easily write my own vectorxprod(A,B) function, but I can't figure out what crossprod is doing in general.
Solution 1:[1]
Here is a short code snippet which works whenever the cross product makes sense: the 3D version returns a vector and the 2D version returns a scalar. If you just want simple code that gives the right answer without pulling in an external library, this is all you need.
# Compute the vector cross product between x and y, and return the components
# indexed by i.
CrossProduct3D <- function(x, y, i=1:3) {
# Project inputs into 3D, since the cross product only makes sense in 3D.
To3D <- function(x) head(c(x, rep(0, 3)), 3)
x <- To3D(x)
y <- To3D(y)
# Indices should be treated cyclically (i.e., index 4 is "really" index 1, and
# so on). Index3D() lets us do that using R's convention of 1-based (rather
# than 0-based) arrays.
Index3D <- function(i) (i - 1) %% 3 + 1
# The i'th component of the cross product is:
# (x[i + 1] * y[i + 2]) - (x[i + 2] * y[i + 1])
# as long as we treat the indices cyclically.
return (x[Index3D(i + 1)] * y[Index3D(i + 2)] -
x[Index3D(i + 2)] * y[Index3D(i + 1)])
}
CrossProduct2D <- function(x, y) CrossProduct3D(x, y, i=3)
Does it work?
Let's check a random example I found online:
> CrossProduct3D(c(3, -3, 1), c(4, 9, 2)) == c(-15, -2, 39)
[1] TRUE TRUE TRUE
Looks pretty good!
Why is this better than previous answers?
- It's 3D (Carl's was 2D-only).
- It's simple and idiomatic.
- Nicely commented and formatted; hence, easy to understand
The downside is that the number '3' is hardcoded several times. Actually, this isn't such a bad thing, since it highlights the fact that the vector cross product is purely a 3D construct. Personally, I'd recommend ditching cross products entirely and learning Geometric Algebra instead. :)
Solution 2:[2]
The help ?crossprod explains it quite clearly. Take linear regression for example, for a model y = XB + e you want to find X'X, the product of X transpose and X. To get that, a simple call will suffice: crossprod(X) is the same as crossprod(X,X) is the same as t(X) %*% X. Also, crossprod can be used to find the dot product of two vectors.
Solution 3:[3]
In response to @Bryan Hanson's request, here's some Q&D code to calculate a vector crossproduct for two vectors in the plane. It's a bit messier to calculate the general 3-space vector crossproduct, or to extend to N-space. If you need those, you'll have to go to Wikipedia :-) .
crossvec <- function(x,y){
if(length(x)!=2 |length(y)!=2) stop('bad vectors')
cv <- x[1]*y[2]-x[2]*y[1]
return(invisible(cv))
}
Solution 4:[4]
Here is a minimalistic implementation for 3D vectors:
vector.cross <- function(a, b) {
if(length(a)!=3 || length(b)!=3){
stop("Cross product is only defined for 3D vectors.");
}
i1 <- c(2,3,1)
i2 <- c(3,1,2)
return (a[i1]*b[i2] - a[i2]*b[i1])
}
If you want to get the scalar "cross product" of 2D vectors u and v, you can do
vector.cross(c(u,0),c(v,0))[3]
Solution 5:[5]
There is a useful math operations package named pracma (https://rdrr.io/rforge/pracma/api/ or CRAN https://cran.r-project.org/web/packages/pracma/index.html). Easy to use and quick. The cross product is literally given by pracma::cross(x, y) for any two vectors.
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 | |
| Solution 2 | James Pringle |
| Solution 3 | Carl Witthoft |
| Solution 4 | Museful |
| Solution 5 |
