Вычислить матрицу проекции / шляпы с помощью QR-факторизации, SVD (и факторизации Холецкого?)

Я пытаюсь вычислить в R матрица проекции P произвольной матрицы N x J S :

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

Я пытался выполнить это с помощью следующей функции:

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

Но когда я использую это Я получаю следующие ошибки:

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

Я думаю, что это результат потери числового значения и / или нестабильности, о чем говорилось во многих местах, таких как r-help и здесь , но У меня недостаточно опыта в использовании SVD или QR-разложения, чтобы решить проблему или привести в действие существующий код. Я также пробовал предложенный код, который должен написать решение как систему:

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

Но все равно он не работает. Мы ценим любые предложения.

Я почти уверен, что моя матрица должна быть обратимой и не иметь коллинеарностей, хотя бы потому, что я пробовал проверить это с помощью матрицы ортогональных фиктивных переменных, но она все еще не работает.

Кроме того, я хотел бы применить это к довольно большим матрицам, поэтому я ищу аккуратное общее решение.