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

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

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

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

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

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

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