Estoy tratando de calcular en R una matriz de proyecciónP de una matriz arbitraria N x JS:

P = S (S'S) ^ -1 S'

He estado intentando realizar esto con la siguiente función:

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

Pero cuando uso esto obtengo errores que se ven así:

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

Creo que esto es el resultado de un flujo inferior numérico y / o inestabilidad como se discutió en numerosos lugares como r-help yaqu, pero no tengo suficiente experiencia en el uso de SVD o descomposición QR para solucionar el problema, o de lo contrario poner este código existente en acción. También probé el código sugerido, que es escribir resolver como un sistema:

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

Pero aún así no funciona. Cualquier sugerencia sera apreciada

Estoy bastante seguro de que mi matriz debe ser invertible y no tiene ninguna linealidad, aunque solo sea porque he intentado probar esto con una matriz de variables ficticias ortogonales, y todavía no funciona.

Además, me gustaría aplicar esto a matrices bastante grandes, así que estoy buscando una solución general ordenada.