Berechnung der Projektions- / Hutmatrix über QR-Faktorisierung, SVD (und Cholesky-Faktorisierung?)
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.
2 Antworten
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
- SVD ist die stabilste Methode, aber zu teuer;
- QR ist befriedigend stabil, mit moderaten Rechenkosten und wird in der Praxis häufig verwendet;
- Pivoted Cholesky ist schnell, mit akzeptabler Stabilität. Für große Matrix wird diese bevorzugt.
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