Ich versuche in R eine Projektionsmatrix P
einer beliebigen N x J Matrix S
zu berechnen:
Ich habe versucht, dies mit der folgenden Funktion auszuführen:
%Vor%Aber wenn ich das benutze, bekomme ich Fehler, die so aussehen:
%Vor%Ich denke, dass dies ein Ergebnis von numerischem Unterlauf und / oder Instabilität ist, wie an vielen Stellen diskutiert wurde, wie r-help und hier , aber Ich bin nicht erfahren genug, SVD oder QR-Zerlegung zu verwenden, um das Problem zu beheben, oder den vorhandenen Code in die Tat umzusetzen. Ich habe auch den vorgeschlagenen Code ausprobiert, der das Schreiben als ein System lösen soll:
%Vor%Aber es geht immer noch nicht. Irgendwelche Vorschläge würden geschätzt.
Ich bin ziemlich sicher, dass meine Matrix invertierbar sein sollte und keine Ko-Linearitäten hat, schon allein deshalb, weil ich versucht habe, dies mit einer Matrix von orthogonalen Dummy-Variablen zu testen, und sie funktioniert immer noch nicht.
>Ich möchte das auch auf ziemlich große Matrizen anwenden, also suche ich nach einer einfachen allgemeinen Lösung.
Obwohl OP seit mehr als einem Jahr nicht mehr aktiv ist, entscheide ich mich immer noch für eine Antwort. Ich würde X
anstelle von S
verwenden, da wir in der Statistik oft Projektionsmatrix im linearen Regressionskontext wollen, wobei X
die Modellmatrix ist, y
der Antwortvektor ist, während H = X(X'X)^{-1}X'
ist hat / Projektionsmatrix, so dass Hy
prädiktive Werte ergibt.
Diese Antwort nimmt den Kontext der gewöhnlichen kleinsten Quadrate an. Für gewichtete kleinste Quadrate siehe Holen Sie sich die Hut-Matrix aus der QR-Dekomposition für die gewichtete Least-Square-Regression .
Eine Übersicht
solve
basiert auf der LU-Faktorisierung einer allgemeinen quadratischen Matrix. Für X'X
(sollte von crossprod(X)
anstatt t(X) %*% X
in R berechnet werden, lesen Sie ?crossprod
für mehr), das symmetrisch ist, können wir chol2inv
verwenden, das auf der Choleksy-Faktorisierung basiert.
Allerdings ist die Dreiecksfaktorisierung weniger stabil als QR
Faktorisierung. Das ist nicht schwer zu verstehen. Wenn X
die bedingte Anzahl kappa
hat, wird X'X
die bedingte Anzahl kappa ^ 2
haben. Dies kann große numerische Schwierigkeiten verursachen. Die Fehlermeldung erhalten Sie:
sagt dies nur. kappa ^ 2
ist ungefähr e-28
, viel viel kleiner als die Maschinengenauigkeit bei e-16
. Mit der Toleranz tol = .Machine$double.eps
wird X'X
als rank deficient angesehen, daher wird die LU- und Cholesky-Faktorisierung zusammenbrechen.
Im Allgemeinen wechseln wir in dieser Situation zu SVD oder QR, aber die
Im Folgenden werde ich alle drei Methoden erklären.
QR-Faktorisierung verwenden
Beachten Sie, dass die Projektionsmatrix permutationsunabhängig ist, d. h. es spielt keine Rolle, ob wir eine QR-Faktorisierung mit oder ohne Pivotierung durchführen.
In R kann qr.default
die LINPACK-Routine DQRDC
für die nicht-geschwenkte QR-Faktorisierung und die LAPACK-Routine DGEQP3
für die blockierte QR-Faktorisierung aufrufen. Lassen Sie uns eine Spielzeugmatrix erstellen und beide Optionen testen:
Beachten Sie, dass $pivot
für zwei Objekte unterschiedlich ist.
Nun definieren wir eine Wrapper-Funktion, um QQ'
zu berechnen:
Wir werden sehen, dass qr_linpack
und qr_lapack
dieselbe Projektionsmatrix haben:
Verwenden der Singulärwertzerlegung
In% berechnet svd
die Singulärwertzerlegung. Wir verwenden immer noch die obige Beispielmatrix X
:
Auch hier erhalten wir die gleiche Projektionsmatrix.
Verwenden der Pivoted Cholesky-Faktorisierung
Zur Demonstration verwenden wir immer noch das obige Beispiel X
.
Auch hier erhalten wir die gleiche Projektionsmatrix.
Wie wäre es mit chol
auf S'S
, dann auf chol2inv
? Sollte stabiler sein:
Tags und Links r regression svd projection-matrix qr-decomposition