10

PRで任意の 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やhereなどの多くの場所で説明されているように、数値のアンダーフローおよび/または不安定性の結果であると思いますが、SVD または QR 分解を使用して問題を解決したり、この既存のコードをアクション。システムとしてソルブを書くことである提案されたコードも試しました：

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

しかし、それでもうまくいきません。任意の提案をいただければ幸いです。

直交ダミー変数の行列でこれをテストしようとしたが、まだ機能しないという理由だけで、行列が可逆であり、共線形性がないことはかなり確信しています。

また、これをかなり大きな行列に適用したいので、きちんとした一般的な解決策を探しています。

4

1 に答える 1

14

OP は 1 年以上アクティブではありませんが、回答を投稿することにしました。Xの代わりにを使用します。統計ではS、線形回帰のコンテキストで射影行列が必要になることがよくあります。XyH = X(X'X)^{-1}X'Hy

この回答は、通常の最小二乗法のコンテキストを前提としています。加重最小二乗法については、加重最小二乗回帰の QR 分解からハット行列を取得する を参照してください。

概要

solve一般的な正方行列の LU 分解に基づいています。対称である ( Rではなく でX'X計算する必要があります。詳しくは、こちらを参照してください) 。crossprod(X)t(X) %*% X?crossprodchol2inv

ただし、三角分解は因数分解よりも不安定ですQR。これを理解するのは難しくありません。X条件付き番号 がある場合kappa、条件付き番号X'Xがありますkappa ^ 2。これは、大きな数値上の困難を引き起こす可能性があります。表示されるエラー メッセージ:

# system is computationally singular: reciprocal condition number = 2.26005e-28

これを伝えているだけです。kappa ^ 2約e-28、機械精度よりかなり小さい程度e-16です。許容範囲tol = .Machine$double.epsでX'Xは、ランク不足と見なされるため、LU とコレスキー分解は失敗します。

通常、この状況では SVD または QR に切り替えますが、ピボットされたコレスキー分解も別の選択肢です。

  • SVD は最も安定した方法ですが、コストがかかりすぎます。
  • QR は、中程度の計算コストで十分に安定しており、実際に一般的に使用されています。
  • Pivoted Cholesky は高速で、許容できる安定性があります。大きな行列の場合、これが優先されます。

以下では、3 つの方法すべてについて説明します。

QR 因数分解の使用

QR因数分解

射影行列は順列に依存しないことに注意してください。つまり、ピボットを使用して QR 分解を実行するかどうかは問題ではありません。

R では、非ピボット QR 分解用のqr.defaultLINPACK ルーチンと、ブロック ピボット付き QR 分解用の LAPACK ルーチンを呼び出すことができます。おもちゃのマトリックスを生成して、両方のオプションをテストしましょう。DQRDCDGEQP3

set.seed(0); X <- matrix(rnorm(50), 10, 5)
qr_linpack <- qr.default(X)
qr_lapack <- qr.default(X, LAPACK = TRUE)

str(qr_linpack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0861 0.3509 0.3357 0.1094 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.37 1.03 1.01 1.15
# $ pivot: int [1:5] 1 2 3 4 5
# - attr(*, "class")= chr "qr"

str(qr_lapack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0646 0.2632 0.2518 0.0821 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.21 1.56 1.36 1.09
# $ pivot: int [1:5] 1 5 2 4 3
# - attr(*, "useLAPACK")= logi TRUE
# - attr(*, "class")= chr "qr"

$pivotは 2 つのオブジェクトで異なることに注意してください。

次に、計算するラッパー関数を定義しますQQ'。

f <- function (QR) {
  ## thin Q-factor
  Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
  ## QQ'
  tcrossprod(Q)
  }

それを見て、同じ射影行列qr_linpackを与えます。qr_lapack

H1 <- f(qr_linpack)
H2 <- f(qr_lapack)

mean(abs(H1 - H2))
# [1] 9.530571e-17

特異値分解の使用

SVD

R では、svd特異値分解を計算します。上記の例の行列を引き続き使用しますX。

SVD <- svd(X)

str(SVD)
# List of 3
# $ d: num [1:5] 4.321 3.667 2.158 1.904 0.876
# $ u: num [1:10, 1:5] -0.4108 -0.0646 -0.2643 -0.1734 0.1007 ...
# $ v: num [1:5, 1:5] -0.766 0.164 0.176 0.383 -0.457 ...

H3 <- tcrossprod(SVD$u)

mean(abs(H1 - H3))
# [1] 1.311668e-16

繰り返しますが、同じ射影行列が得られます。

ピボットコレスキー分解の使用

ピボットコレスキー分解

デモンストレーションのために、上記の例を引き続き使用しますX。

## pivoted Chol for `X'X`; we want lower triangular factor `L = R'`:
## we also suppress possible rank-deficient warnings (no harm at all!)
L <- t(suppressWarnings(chol(crossprod(X), pivot = TRUE)))

str(L)
# num [1:5, 1:5] 3.79 0.552 -0.82 -1.179 -0.182 ...
# - attr(*, "pivot")= int [1:5] 1 5 2 4 3
# - attr(*, "rank")= int 5

## compute `Q'`
r <- attr(L, "rank")
piv <- attr(L, "pivot")
Qt <- forwardsolve(L, t(X[, piv]), r)

## P = QQ'
H4 <- crossprod(Qt)

## compare
mean(abs(H1 - H4))
# [1] 6.983997e-17

繰り返しますが、同じ射影行列が得られます。

于 2016-09-02T17:38:51.157 に答える