Estou tentando calcular em R uma matriz de projeçãoP de uma matriz N x J arbitráriaS:

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

Eu tenho tentado fazer isso com a seguinte função:

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

Mas quando eu uso isso, recebo erros parecidos com o seguinte:

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

Acho que isso é resultado de um fluxo insuficiente numérico e / ou instabilidade, conforme discutido em vários lugares como r-help eaqu, mas não tenho experiência suficiente usando decomposição SVD ou QR para corrigir o problema, ou coloque esse código existente em ação. Eu também tentei o código sugerido, que é escrever solucionar como um sistema:

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

Mas ainda não funciona. Todas as sugestões serão apreciadas.

Tenho certeza de que minha matriz deve ser invertível e não possui co-linearidades, mesmo porque tentei testá-lo com uma matriz de variáveis fictícias ortogonais e ainda não funcion

Além disso, eu gostaria de aplicar isso a matrizes bastante grandes, então estou procurando uma solução geral elegant