Es erschien in linearer Regression mit mehreren Variablen eines bestimmten Online-Kurses für maschinelles Lernen. Über die normale Gleichung.
Dr. Andrew Ng (im Folgenden als Dr. Ang abgekürzt) sagte: "Es ist ein Schmerz abzuleiten (freie Übersetzung)" und nur das Ergebnis wurde gezeigt, also habe ich etwas tiefer gegraben.
Unter ihnen habe ich einige Fragen, also werde ich sie teilen. Ich selbst verstehe es vielleicht noch nicht, also begrüße ich Tsukkomi.
[24.07.2015 23:10] Der Bestätigungscode wurde hinzugefügt und erheblich überarbeitet.
Überprüfen Sie zunächst.
$ X $ ist die $ m \ times (n + 1) $ -Matrix, die das gesamte Trainingsdatenmerkmal darstellt ($ m $ (Anzahl der Zeilen) ist die Anzahl der Daten, $ n $ ist die Anzahl der Merkmale (Merkmale)). $ y $ ist ein $ m $ -Dimensionsvektor, der die "korrekten Werte" der Trainingsdaten auflistet. $ \ theta $ repräsentiert $ (n + 1) $ repräsentiert die Parameter der Hypothesenfunktion ($ h_ {\ theta} (x) = \ sum_i \ theta_i x_i $, $ 0 \ le i \ le n $) Dimensionsvektor. Außerdem ist $ J (\ theta) = \ frac {1} {2m} \ sum_ {i = 1} ^ {m} (h_ {\ theta} (x ^ {(i)}) - y ^ {(i)} ) ^ 2 $ heißt die Kostenfunktion [^ 1](der linearen Regression), und der Zweck der linearen Regression besteht darin, sie zu minimieren.
[^ 1]: Diese Formel ist die sogenannte quadratische Fehlerfunktion. Mit anderen Worten, es ist die Methode des minimalen Quadrats.
Als Minimierungsmethode wird das letzte $ J (\ theta) $ teilweise durch $ \ theta_0, \ dots, \ theta_n $ unterschieden, und $ \ theta $ wird so berechnet, dass sie alle 0 sind. Wenn Sie beispielsweise teilweise mit einem bestimmten $ \ theta_j $ differenzieren:
Wenn Sie dies in $ m $ Zeilen umwandeln und in Form einer Matrix [^ 2] umschreiben:
[^ 2]: Sumimasen, ich habe es hier gefaltet. Es ist ärgerlich, eins nach dem anderen zu schreiben (^ - ^;
Dies wird in $ = {\ bf 0} $ übersetzt und für $ \ theta $ gelöst:
Dies ist die "normale Gleichung", die Dr. Ang sagt. Durch Berechnung der rechten Seite der letzten Gleichung kann der optimale Parameter $ \ theta $ (der den quadratischen Fehler minimiert) gefunden werden.
Es ist nicht immer so, dass $ (X ^ TX) ^ {-1} $ für $ X $ immer vorhanden ist. Wieder sagte Professor Ang in seinem Vortrag: "Es ist sehr selten, also mach dir keine Sorgen." Er erklärt auch, dass selbst wenn ein solcher Fall beispielsweise in der Oktave auftreten sollte, "pinv (X" * X) * X "* y" verwendet werden sollte. "X" ist eine transponierte Matrix der Matrix "X" ($ X ^ T $ für $ X $).
Octaves pinv ist" Pseudoinverse "(pseudoinverse Matrix, Generalisierung) Wird auch als inverse Matrix bezeichnet. Im Folgenden "[Generalisierte inverse Matrix](https://ja.wikipedia.org/wiki/%E6%93%AC%E4%BC%BC%E9%80%86%E8%A1%". 8C% E5% 88% 97) "wird verwendet).
Wenn die Argumentmatrix reversibel ist, wird die inverse Matrix (die übereinstimmende Matrix) zurückgegeben, andernfalls wird die "verallgemeinerte inverse Matrix" zurückgegeben, die bestimmte Eigenschaften erfüllt.
Die Details der verallgemeinerten inversen Matrix werden hier weggelassen, aber in einfachen Worten handelt es sich um eine "Matrix mit Eigenschaften, die der inversen Matrix ähnlich sind". Dies kann nicht nur für nicht reguläre (= nicht reversible) quadratische Matrizen definiert werden, sondern auch für nicht quadratische (= $ m \ mal n $ Matrizen wie $ m \ ne n $) Matrizen.
Auf diese Weise können Sie $ (X ^ TX) ^ {-1} $ zurückgeben und berechnen, ob $ X ^ TX $ reversibel ist, andernfalls durch eine ordnungsgemäß berechnete "verallgemeinerte Inverse". Es ist ein Mechanismus, der fehlerfrei berechnet und angemessene Berechnungsergebnisse erzielt.
Wie ich zuvor geschrieben habe, ist die verallgemeinerte inverse Matrix auch in der Matrix $ m \ times n $ ($ m \ ne n $) definiert. Insbesondere wenn $ m> n $ (= mehr Zeilen als Spalten = Hochformatmatrix), kann es tatsächlich wie folgt geschrieben werden:
Ich denke an lineare Regression (multivariate). Die Anzahl der Daten ist im Vergleich zur Anzahl der Merkmale (im Allgemeinen) tendenziell sehr groß. Mit anderen Worten sollte die Matrix eine vertikal lange Matrix mit mehr Zeilen sein.
Schauen wir uns nun noch einmal die Normalgleichung an.
Der andere Teil als $ y $ auf der rechten Seite. Es ist die gleiche Form, oder? Mit anderen Worten. Ist es nicht in Ordnung, das zu schreiben? Wann.
Im Oktavcode "pinv (X) * y". Ist das nicht in Ordnung? Wann.
des Weiteren. Betrachten Sie die folgende Gleichung für einen Moment.
Wenn $ X $ eine reguläre Matrix ist (reversibel mit einer quadratischen Matrix), kann sie mit $ \ theta = X ^ {-1} y $ gelöst werden. Wenn nicht, gibt es keine solche gewöhnliche Methode (es gibt überhaupt keine zufriedenstellende Lösung oder die Lösung ist nicht eindeutig bestimmt). Wenn wir jedoch $ \ theta = X ^ -y $ betrachten, das die inverse Matrix durch die verallgemeinerte inverse Matrix ersetzt, können wir die Lösung finden, die den Fehler (oder die Norm) minimiert. Mit anderen Worten bedeutet diese Formel, die minimale Fehlerlösung für $ X \ theta = y $ zu finden.
Und Octave bietet eine andere Möglichkeit, $ X \ theta = y $ so wie es ist zu lösen. Es wird mit dem Backslash-Operator als "X \ y" geschrieben. Dies bedeutet $ X ^ {-1} y $, wenn $ X $ eine reguläre Matrix ist, andernfalls scheint es die gleiche Berechnung wie $ X ^ -y $ durchzuführen.
Erinnern wir uns nun an die ursprüngliche Form der Normalgleichung.
Dies ist auch in Form von $ \ bigcirc \ theta = \ bigcirc $, nicht wahr? Es liegt in einer Form vor, die mit ○ \ ○ und dem Backslash-Operator gelöst werden kann. Bei Anwendung "(X" * X) \ (X "* y)".
Mit anderen Worten. Der Octave-Code zum Lösen der Normalgleichung bedeutet, dass es vier mögliche Typen gibt:
pinv(X'*X)*X'*ypinv(X)*yX \ y(X'*X) \ (X'*y)Von diesen wird in der Vorlesung nur das frühere "pinv (X" * X) * X "* y" vorgestellt. Warum werden der 2., 3. und 4. nicht eingeführt?
Also habe ich den eigentlichen Code geschrieben und überprüft.
solve_X_y.m
% solve_X_y.m
X = rand(10, 4)
# X =
#
# 0.033536 0.816107 0.996677 0.958327
# 0.683542 0.116498 0.614316 0.884338
# 0.734337 0.769245 0.696212 0.245270
# 0.216938 0.013297 0.885327 0.906086
# 0.630620 0.733668 0.820551 0.784664
# 0.138834 0.838178 0.216751 0.638286
# 0.100739 0.893597 0.891867 0.239482
# 0.362333 0.404999 0.018274 0.922847
# 0.102606 0.442110 0.744582 0.452299
# 0.590709 0.274452 0.459526 0.656588
y = rand(10, 1)
# y =
#
# 0.48518
# 0.13242
# 0.60525
# 0.31265
# 0.59250
# 0.47161
# 0.95971
# 0.44011
# 0.60115
# 0.75571
# calcuration
# [1]
pinv(X' * X) * X' * y
# ans =
#
# 0.1861915
# 0.5484641
# 0.2473279
# -0.0031611
# [2]
pinv(X) * y
# ans =
#
# 0.1861915
# 0.5484641
# 0.2473279
# -0.0031611
# [3]
X \ y
# ans =
#
# 0.1861915
# 0.5484641
# 0.2473279
# -0.0031611
# [4]
(X'*X) \ (X'*y)
# ans =
#
# 0.1861915
# 0.5484641
# 0.2473279
# -0.0031611
# time measurement (n = 10)
# [1]
tic();
for k=1:10000;
X = rand(40, 10);
y = rand(40, 1);
pinv(X' * X) * X' * y;
end;
toc()
# Elapsed time is 1.26513 seconds.
# [2]
tic();
for k=1:10000;
X = rand(40, 10);
y = rand(40, 1);
pinv(X) * y;
end;
toc()
# Elapsed time is 1.16283 seconds.
# [3]
tic();
for k=1:10000;
X = rand(40, 10);
y = rand(40, 1);
X \ y;
end;
toc()
# Elapsed time is 0.902037 seconds.
# [4]
tic();
for k=1:10000;
X = rand(40, 10);
y = rand(40, 1);
(X'*X) \ (X'*y);
end;
toc()
# Elapsed time is 0.689348 seconds.
# time measurement (n = 30)
# [1]
tic();
for k=1:10000;
X = rand(100, 30);
y = rand(100, 1);
pinv(X' * X) * X' * y;
end;
toc()
# Elapsed time is 5.79588 seconds.
# [2]
tic();
for k=1:10000;
X = rand(100, 30);
y = rand(100, 1);
pinv(X) * y;
end;
toc()
# Elapsed time is 7.11547 seconds.
# [3]
tic();
for k=1:10000;
X = rand(100, 30);
y = rand(100, 1);
X \ y;
end;
toc()
# Elapsed time is 3.64188 seconds.
# [4]
tic();
for k=1:10000;
X = rand(100, 30);
y = rand(100, 1);
(X'*X) \ (X'*y);
end;
toc()
# Elapsed time is 1.37039 seconds.
Was die Berechnungsergebnisse betrifft, berechnen zunächst alle [1], [2], [3] und [4] mit Sicherheit den gleichen Wert [^ 3].
[^ 3]: Bei näherer Betrachtung scheinen tatsächlich unterschiedliche Werte mit einem Fehler von etwa 1e-14 berechnet zu werden. Gleitkomma-Berechnungsfehler aufgrund unterschiedlicher Berechnungs- und Berechnungsmethoden, nicht wahr?
Bei der Messung der Zeit ist es natürlich, dass je größer der Wert von $ n $ (= Anzahl der Spalten in der Matrix = Anzahl der Features) ist, desto länger dauert es und zu diesem Zeitpunkt [1] pinv (X) '* X) * X' * y hat eine kürzere Ausführungszeit als [2] pinv (X) * y. Mit anderen Worten, die Berechnungsformel sieht komplizierter aus, aber der Rechenaufwand ist bei ersteren geringer.
Dies liegt daran, dass wenn $ X $ eine $ m \ mal n $ Matrix ist ($ m> n $), dann ist $ X ^ TX $ eine quadratische Matrix von $ n \ mal n $, die fast reversibel ist. Verglichen mit der Tatsache, dass die inverse Matrix mit hoher Geschwindigkeit berechnet wird (die Kosten für die Multiplikation von $ X ^ T $ von $ n \ mal m $ danach sind relativ gering), ist die verallgemeinerte inverse Matrix von $ X $, die in erster Linie keine quadratische Matrix ist, ( Außerdem ist es $ m \ mal n> n \ mal n $) Ich denke, dass die Berechnungskosten hoch sind.
Die Ergebnisse zeigen jedoch, dass die Berechnungszeit von "X \ y" in [3] schneller ist.
Außerdem ist das (X '* X) \ (X' * y) in [4] viel schneller, nicht wahr?
[4] ist aus dem gleichen Grund wie zuvor schneller als [3].
Da sich "rand ()" verwendet, ändern sich die Werte von "X" und "y" natürlich jedes Mal, wenn sie ausgeführt werden, aber das Berechnungsergebnis und die Ausführungszeit sind fast gleich. Bei dieser Rate scheint "(X" * X) \ (X "* y)" das Beste zu sein, aber ...
Früher habe ich versucht, mit Code unter Verwendung einer Zufallsmatrix zu überprüfen. Versuchen wir es jetzt mit einer Matrix aus absichtlich angeordneten Zahlen.
rank_deficient.m
X = [1 1 2 3;
1 2 3 5;
1 3 5 8;
1 5 8 13;
1 8 13 21;
1 13 21 34]
y = [8; 13; 21; 34; 55; 89]
# check rank of matrix
rank(X)
# => 3
# calcuration
# [1]
pinv(X'*X) * X'*y
# ans =
#
# 3.1974e-13
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# [2]
pinv(X) * y
# ans =
#
# 0.00000
# 0.33333
# 1.33333
# 1.66667
# [3]
X \ y
# ans =
#
# -1.3628e-14
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# [4]
(X'*X) \ (X'*y)
# warning: matrix singular to machine precision, rcond = 4.97057e-18
# ans =
#
# -1.8859e-13
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# Square Matrix
X = X(1:4, 1:4)
# X =
#
# 1 1 2 3
# 1 2 3 5
# 1 3 5 8
# 1 5 8 13
y = y(1:4)
# y =
#
# 8
# 13
# 21
# 34
# calcuration
# [1]
pinv(X'*X) * X'*y
# ans =
#
# 1.8119e-13
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# [2]
pinv(X) * y
# ans =
#
# -7.1054e-15
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# [3]
X \ y
# warning: matrix singular to machine precision, rcond = 0
# ans =
#
# -7.3807e-15
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
# [4]
(X'*X) \ (X'*y)
# warning: matrix singular to machine precision, rcond = 1.26207e-17
# ans =
#
# 1.5742e-14
# 3.3333e-01
# 1.3333e+00
# 1.6667e+00
Die zweite bis vierte Spalte der $ X $ -Matrix sind aufeinanderfolgende Fibonacci-Zahlen. Das heißt, $ x_2 + x_3 = x_4 $ (wenn Sie die Spalte $ i $ als $ x_i $ schreiben). Da die Matrix $ m \ times 4 $ ($ m \ ge 4 $) nur drei unabhängige Spalten hat (da die vierte Spalte durch eine lineare Verknüpfung der anderen Spalten dargestellt wird), ist der Rang der Matrix Es ist $ 3 $ [^ 4].
[^ 4]: In einer Matrix X mit m Zeilen und n Spalten wird X als "voller Rang" bezeichnet, wenn Rang (X) = min (m, n) ist. Wenn Rang (X) <min (m, n) ist, spricht man von "Rangmangel".
Zu diesem Zeitpunkt ist $ X ^ TX $, eine $ n \ times n $ -Matrix, keine reguläre Matrix (sie hat keine inverse Matrix). Wenn $ m = n $ ist, ist $ X $ selbst eine quadratische Matrix, aber auch keine reguläre Matrix.
Infolgedessen ist [4]((X '* X) \ (X' * y)), wenn $ m> n $, und [3], 4, wenn $ m = n $. Bei der Berechnung von X \ y, (X '* X) \ (X' * y)`) wurde eine Warnung ausgegeben.
Die Lösung durch den Operator "" scheint ein solches Problem zu haben.
Andererseits scheint es im Fall der Lösung durch die verallgemeinerte inverse Matrix so zu sein, dass die Berechnung ohne Vorwarnung durchgeführt wird. Das liegt daran, dass es von Anfang an "vernünftige Berechnungen durchführen soll, auch wenn es keine reguläre Matrix ist (auch wenn es keine quadratische Matrix ist)".
Dr. Ang führte "pinv (X" * X) * X "* y" ein, das dem abgeleiteten Ausdruck (aus der ursprünglichen Ausdruckstransformation) am nächsten kommt und besser ist als "pinv (X) * y". Ich denke, das liegt daran, dass es mit hoher Geschwindigkeit berechnet werden kann.
Außerdem habe ich in diesen Berechnungen nicht "X \ y" oder "(X" * X) \ (X "* y)" erwähnt, da die Matrix links vom Operator eine nicht reguläre quadratische Matrix (oder " Es scheint, dass das Problem im Fall der "Rangabfall" -Matrix auftritt (da kein Grund zur Sorge besteht, wenn "pinv ()" verwendet wird) [^ 5].
Dies kann jedoch als der Fall bezeichnet werden, wenn bei der Auswahl von Daten (Feature) überhaupt ein Problem auftritt. Wenn Sie also unabhängig von anderen Features nur unabhängige Features auswählen, ist dies einfacher und schneller (X * X). Ich denke nicht, dass es schlecht ist, \ (X '* y) `zu verwenden.
Konzentrieren Sie sich danach auf die Datenanalyse zu diesem Zweck (zahlen Sie die Kosten für die Arbeit) und beschleunigen Sie die Berechnung ((X '* X) \ (X' * y)) oder speziell für die Daten Führen Sie eine relativ schnelle und fehlerfreie Berechnung ohne Einfügen durch (pinv (X '* X) * X' * y)? Es wird dieses Schnäppchen sein.
[^ 5]: Bevor ich den Artikel änderte, fragte ich: "Wenn es keine Probleme gibt, ist X \ y besser? Was ist mit?", Aber dies ist wahrscheinlich das größte Problem. Ich habe es selbst gelöst, aber was ist mit?
Bei solchen Vorlesungen (bei denen die Sprache und das Verarbeitungssystem angegeben sind) und beim Lernen in Lehrbüchern versuche ich oft, in anderen Sprachen zu schreiben, um mein Verständnis zu vertiefen. Daher habe ich den diesmal überprüften Code auch mit Julia (v0.3.x / 0.4.0-dev) und Python (v2.7.x / 3.x) + NumPy überprüft.
In Python gibt es keinen entsprechenden Operator für "X \ y", sondern eine Funktion namens "numpy.linalg.solve (X, y)", aber es scheint, dass es nicht funktioniert, wenn "X" keine quadratische Matrix ist. ..
Julia v0.3.x / 0.4.0 Funktioniert mit beiden [^ 6]:
solve_X_y.jl
X = rand(10, 4)
# 10x4 Array{Float64,2}:
# 0.71148 0.968352 0.0952939 0.796324
# 0.915128 0.128326 0.630086 0.0635579
# 0.351199 0.131409 0.934867 0.501701
# 0.165645 0.874088 0.173725 0.976326
# 0.765261 0.790716 0.760362 0.204496
# 0.544099 0.156464 0.041718 0.507071
# 0.764964 0.852837 0.230312 0.134783
# 0.0738597 0.75529 0.693856 0.0107293
# 0.621861 0.56881 0.66972 0.163911
# 0.9471 0.453841 0.466836 0.10744
y = rand(10, 1)
# 10x1 Array{Float64,2}:
# 0.389321
# 0.436261
# 0.308189
# 0.734617
# 0.410237
# 0.4969
# 0.0708882
# 0.0840005
# 0.944711
# 0.14718
# calcuration
# [1]
pinv(X' * X) * X' * y
# 4x1 Array{Float64,2}:
# 0.169937
# -0.0365938
# 0.273122
# 0.55004
# [2]
pinv(X) * y
# 4x1 Array{Float64,2}:
# 0.169937
# -0.0365938
# 0.273122
# 0.55004
# [3]
X \ y
# 4x1 Array{Float64,2}:
# 0.169937
# -0.0365938
# 0.273122
# 0.55004
# [4]
(X'*X) \ (X'*y)
# 4x1 Array{Float64,2}:
# 0.169937
# -0.0365938
# 0.273122
# 0.55004
# time measurement (n = 10)
# [1]
@time for k=1:10000
X = rand(40, 10)
y = rand(40, 1)
pinv(X' * X) * X' * y
end
# elapsed time: 1.087745051 seconds (283600016 bytes allocated, 17.28% gc time)
# [2]
@time for k=1:10000
X = rand(40, 10)
y = rand(40, 1)
pinv(X) * y
end
# elapsed time: 1.278193773 seconds (334800016 bytes allocated, 17.29% gc time)
# [3]
@time for k=1:10000
X = rand(40, 10)
y = rand(40, 1)
X \ y
end
# elapsed time: 1.014968744 seconds (324320000 bytes allocated, 20.29% gc time)
# [4]
@time for k=1:10000
X = rand(100, 30)
y = rand(100, 1)
(X'*X) \ (X'*y)
end
# elapsed time: 0.163586767 seconds (62720032 bytes allocated, 41.51% gc time)
# time measurement (n = 30)
# [1]
@time for k=1:10000
X = rand(100, 30)
y = rand(100, 1)
pinv(X' * X) * X' * y
end
# elapsed time: 5.820615493 seconds (1557840000 bytes allocated, 19.02% gc time)
# [2]
@time for k=1:10000
X = rand(100, 30)
y = rand(100, 1)
pinv(X) * y
end
# elapsed time: 7.518744844 seconds (1914480016 bytes allocated, 16.51% gc time)
# [3]
@time for k=1:10000
X = rand(100, 30)
y = rand(100, 1)
X \ y
end
# elapsed time: 3.455976006 seconds (1292000000 bytes allocated, 22.67% gc time)
# [4]
@time for k=1:10000
X = rand(100, 30)
y = rand(100, 1)
(X'*X) \ (X'*y)
end
# elapsed time: 0.777771618 seconds (407840016 bytes allocated, 32.71% gc time)
[^ 6]: Julias Art, Matrizen / Vektoren zu handhaben und zu formatieren, wurde MATLAB / Octave entnommen, so dass es fast ohne Umschreiben funktionierte.
Ich habe den Eindruck, dass es nicht so schnell ist, wie ich erwartet hatte, aber ich denke, es liegt daran, dass ich "für" direkt auf der obersten Ebene geschrieben habe.
Es wird sicherlich schneller sein, wenn Sie es in eine Funktion oder einen anderen Schreibstil ändern, der auf Leistung achtet.
Aber (X '* X) \ (X' * y) ist schnell genug, um andere zu überwältigen! Aber ...
rank_deficient.jl
X = [1 1 2 3;
1 2 3 5;
1 3 5 8;
1 5 8 13;
1 8 13 21;
1 13 21 34]
y = [8; 13; 21; 34; 55; 89]
# check rank of matrix
rank(X)
# => 3
# calcuration
# [1]
pinv(X'*X) * X'*y
# 4-element Array{Float64,1}:
# -7.10543e-15
# 0.333333
# 1.33333
# 1.66667
# [2]
pinv(X) * y
# 4-element Array{Float64,1}:
# 3.55271e-15
# 0.333333
# 1.33333
# 1.66667
# [3]
X \ y
# 4-element Array{Float64,1}:
# -4.35117e-15
# 2.0
# 3.0
# 0.0
# [4]
(X'*X) \ (X'*y)
# 4-element Array{Float64,1}:
# 3.22113e-13
# -1.50024
# -0.500244
# 3.50024
# Square Matrix
X = X[1:4, 1:4]
# 4x4 Array{Int64,2}:
# 1 1 2 3
# 1 2 3 5
# 1 3 5 8
# 1 5 8 13
y = y[1:4]
# 4-element Array{Int64,1}:
# 8
# 13
# 21
# 34
# calcuration
# [1]
pinv(X'*X) * X'*y
# 4-element Array{Float64,1}:
# 8.52651e-14
# 0.333333
# 1.33333
# 1.66667
# [2]
pinv(X) * y
# 4-element Array{Float64,1}:
# 3.55271e-15
# 0.333333
# 1.33333
# 1.66667
# x[3]
# X \ y
# @> SingularException(4)
# x[4]
# (X'*X) \ (X'*y)
# @> SingularException(4)
Im Falle eines Rangabfalls. Im Fall einer vertikal langen Matrix ergeben sich die Berechnungsergebnisse, bei denen [3] und [4]("X \ y", "(X" * X) \ (X "* y)") sich erheblich von denen unterscheiden, die "pinv ()" verwenden Es ist geworden. Außerdem erscheint im Fall einer quadratischen Matrix ein Fehler mit der Aufschrift "keine reguläre Matrix" und eine Berechnung ist nicht möglich. N (> _ <) (Übrigens möchte ich kurz erwähnen, dass die Ergebnisse zwischen Julia v0.3.x und v0.4.0 leicht unterschiedlich waren.)
Funktioniert auch mit Python v2.7.x / 3.x, die Ergebnisse sind ähnlich:
solve_X_y.py
import numpy as np
X = np.random.rand(10, 4)
# array([[ 0.61009055, 0.71722947, 0.48465025, 0.15660522],
# [ 0.02424431, 0.49947237, 0.60493258, 0.8988653 ],
# [ 0.65048106, 0.69667863, 0.52860957, 0.65003537],
# [ 0.56541266, 0.25463788, 0.74047536, 0.64691215],
# [ 0.03052439, 0.47651739, 0.01667898, 0.7613639 ],
# [ 0.87725831, 0.47684888, 0.44039111, 0.39706053],
# [ 0.58302851, 0.20919564, 0.97598994, 0.19268083],
# [ 0.35987338, 0.98331404, 0.06299533, 0.76193058],
# [ 0.625453 , 0.70985323, 0.62948802, 0.627458 ],
# [ 0.64201569, 0.22264827, 0.71333221, 0.53305839]])
y = np.random.rand(10, 1)
# array([[ 0.99674247],
# [ 0.66282312],
# [ 0.68295932],
# [ 0.14330449],
# [ 0.17467666],
# [ 0.90896029],
# [ 0.65385071],
# [ 0.00748736],
# [ 0.93824979],
# [ 0.91696375]])
# calcuration
# [1]
np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
# array([[ 0.32591078],
# [ 0.46479763],
# [ 0.6684976 ],
# [-0.26695783]])
# [2]
np.linalg.pinv(X).dot(y)
# array([[ 0.32591078],
# [ 0.46479763],
# [ 0.6684976 ],
# [-0.26695783]])
# x[3]
# np.linalg.solve(X, y)
# @> LinAlgError
# [4]
np.linalg.solve(X.T.dot(X), X.T.dot(y))
# array([[ 0.32591078],
# [ 0.46479763],
# [ 0.6684976 ],
# [-0.26695783]])
# time measurement (n = 10)
from timeit import timeit
# [1]
def test_a():
X = np.random.rand(40, 10)
y = np.random.rand(40, 1)
np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
timeit("test_a()", setup="from __main__ import test_a", number=10000)
# 1.1948060989379883
# [2]
def test_b():
X = np.random.rand(40, 10)
y = np.random.rand(40, 1)
np.linalg.pinv(X).dot(y)
timeit("test_b()", setup="from __main__ import test_b", number=10000)
# 1.2698009014129639
# [4]
def test_c():
X = np.random.rand(40, 10)
y = np.random.rand(40, 1)
np.linalg.solve(X.T.dot(X), X.T.dot(y))
timeit("test_c()", setup="from __main__ import test_c", number=10000)
# 0.4645709991455078
# time measurement (n = 30)
# [1]
def test_d():
X = np.random.rand(100, 30)
y = np.random.rand(100, 1)
np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
timeit("test_d()", setup="from __main__ import test_d", number=10000)
# 4.615994930267334
# [2]
def test_e():
X = np.random.rand(100, 30)
y = np.random.rand(100, 1)
np.linalg.pinv(X).dot(y)
timeit("test_e()", setup="from __main__ import test_e", number=10000)
# 5.413921117782593
# [4]
def test_f():
X = np.random.rand(100, 30)
y = np.random.rand(100, 1)
np.linalg.solve(X.T.dot(X), X.T.dot(y))
timeit("test_f()", setup="from __main__ import test_f", number=10000)
# 0.9642360210418701
In NumPys ndarray ist "*" das Produkt jedes Elements, und das Produkt der Matrix muss die Methode ".dot ()" verwenden. Selbst komplizierte Ausdrücke werden noch komplizierter ... aber schnell! lösen () ist auch schnell! Aber ...
rank_deficient.jl
import numpy as np
X = np.array([
[1, 1, 2, 3],
[1, 2, 3, 5],
[1, 3, 5, 8],
[1, 5, 8, 13],
[1, 8, 13, 21],
[1, 13, 21, 34]])
y = np.array([[8], [13], [21], [34], [55], [89]])
# check rank of matrix
np.linalg.matrix_rank(X)
# => 3
# calcuration
# [1]
np.linalg.pinv(X.T.dot(X)).dot(X.T.dot(y))
# array([[ 2.27373675e-13],
# [ 3.33333333e-01],
# [ 1.33333333e+00],
# [ 1.66666667e+00]])
# [2]
np.linalg.pinv(X).dot(y)
# array([[ 3.55271368e-15],
# [ 3.33333333e-01],
# [ 1.33333333e+00],
# [ 1.66666667e+00]])
# [4]
np.linalg.solve(X.T.dot(X), X.T.dot(y))
# array([[ -8.12048841e-14],
# [ -5.00000000e+00],
# [ -4.00000000e+00],
# [ 7.00000000e+00]])
# Square Matrix
X = X[0:4]
# array([[ 1, 1, 2, 3],
# [ 1, 2, 3, 5],
# [ 1, 3, 5, 8],
# [ 1, 5, 8, 13]])
y = y[0:4]
# array([[ 8],
# [13],
# [21],
# [34]])
# calcuration
# [1]
np.linalg.pinv(X.T.dot(X)).dot(X.T.dot(y))
# array([[ -1.13686838e-13],
# [ 3.33333333e-01],
# [ 1.33333333e+00],
# [ 1.66666667e+00]])
# [2]
np.linalg.pinv(X).dot(y)
# array([[ 4.44089210e-15],
# [ 3.33333333e-01],
# [ 1.33333333e+00],
# [ 1.66666667e+00]])
# [4]
np.linalg.solve(X.T.dot(X), X.T.dot(y))
# array([[ -1.47008842e-14],
# [ 1.00000000e+00],
# [ 2.00000000e+00],
# [ 1.00000000e+00]])
Im Falle eines Rangabfalls. Im Fall von NumPy gab es keinen Fehler bei der Verwendung von "lösen ()", aber das Berechnungsergebnis unterschied sich immer noch stark von dem bei "pinv ()".
Recommended Posts