The inverse of the cross product of a m x n matrix

Question

I would like to recreate the following calculation with a random matrix:

I started out with the following, which gives a result:

kmin1 <- cbind(1:10,1:10,6:15,1:10,1:10,6:15,1:10,1:10,6:15)
C <- cbind(1, kmin1)                # Column of 1s
diag(C) <- 1
Ccrosprod <-crossprod(C)            # C'C
Ctranspose <- t(C)                  # C'
CCtransposeinv <- solve(Ccrosprod)  # (C'C)^-1
W <- Ctranspose %*% CCtransposeinv  # W=(C'C)^-1*C'

My assumption is however that C should be able to be an m x n matrix, as there is no good reason to assume that factors equal observations.

EDIT: Based on the comment by Hong Ooi, I changed kmin1 <- matrix(rexp(200, rate=.1), ncol=20) into kmin1 <- matrix(rexp(200, rate=.1), nrow=20)

I checked the wikipedia and learned that an m x n might have a left or a right inverse. To put this into practice I attempted the following:

kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)                # Column of 1s
Ccrosprod <-crossprod(C)            # C'C
Ctranspose <- t(C)                  # C'
CCtransposeinv <- solve(Ccrosprod)  # (C'C)^-1
W <- Ctranspose %*% CCtransposeinv  # W=(C'C)^-1*C'

EDIT: Based on the comments below this questions everything works.

I would post this on stackexchange if I was sure this did not have anything to do with syntax, but as I am not experienced with matrices I am not sure.

Answer 1

If the columns of C are linearly independent then C'C is invertible and (C'C)^-1C' equals any of these:

set.seed(123)
kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)

r1 <- solve(crossprod(C), t(C))
r2 <- qr.solve(crossprod(C), t(C))

r3 <- chol2inv(chol(crossprod(C))) %*% t(C)

r4 <- with(svd(C), v %*% diag(1/d) %*% t(u))
r5 <- with(eigen(crossprod(C)), vectors %*% diag(1/values) %*% t(vectors)) %*% t(C)

r6 <- coef(lm.fit(C, diag(nrow(C))))

# check

all.equal(r1, r2)
## [1] TRUE

all.equal(r1, r3)
## [1] TRUE

all.equal(r1, r4)
## [1] TRUE

all.equal(r1, r5)
## [1] TRUE

dimnames(r6) <- NULL
all.equal(r1, r6)
## [1] TRUE

If C'C is not necessarily invertible then the answer is not necessarily unique (although if we were interested in C(C'C)^-C' then that would be unique even though the pseudoinverse of C'C may not be). At any rate we can form one pseudo inverse by taking the singular value decomposition (or the eigenvalue decomposition) and using the reciprocal of the singular values (or eigenvalues) and using 0 for those that are near 0. This is equivalent to using the Moore Penrose pseudo inverse. (The lm.fit approach shown above will work too but will generate some NAs in the result.)

set.seed(123)
kmin1 <- matrix(rexp(200, rate=.1), nrow=20)
C <- cbind(1, kmin1)
C[, 11] <- C[, 2] + C[, 3] # force singularity

eps <- 1.e-5 
s1 <- with(svd(C), v %*% diag(ifelse(abs(d) < eps, 0, 1/(d))) %*% t(u))
s2 <- with(eigen(crossprod(C)), 
  vectors %*% diag(ifelse(abs(values) < eps, 0, 1/values)) %*% t(vectors)) %*% t(C)

# check
all.equal(s1, s2)
## [1] TRUE