Il est sorti en régression linéaire avec plusieurs variables d'un certain cours d'apprentissage automatique en ligne, À propos de l'équation normale.
Le Dr Andrew Ng (ci-après abrégé en Dr Ang) a dit: "C'est une douleur à dériver (traduction libre)" et seul le résultat a été montré, alors j'ai creusé un peu plus profondément.
Parmi eux, j'ai quelques questions, je vais donc les partager. Je ne le comprends peut-être pas encore moi-même, alors je souhaite la bienvenue à Tsukkomi.
[2015/07/24 23:10] Le code de vérification a été ajouté et a été considérablement révisé.
Tout d'abord, révisez.
$ X $ est la matrice $ m \ times (n + 1) $ qui représente la totalité de la fonctionnalité de données d'entraînement ($ m $ (lignes) est le nombre de données, $ n $ est le nombre de fonctionnalités). $ y $ est un vecteur de dimension $ m $ qui répertorie les «valeurs correctes» des données d'apprentissage. $ \ theta $ représente $ (n + 1) $ représentant les paramètres de la fonction d'hypothèse ($ h_ {\ theta} (x) = \ sum_i \ theta_i x_i $, $ 0 \ le i \ le n $) Vecteur de dimension. De plus, $ J (\ theta) = \ frac {1} {2m} \ sum_ {i = 1} ^ {m} (h_ {\ theta} (x ^ {(i)}) - y ^ {(i)} ) ^ 2 $ est appelée la fonction de coût [^ 1](de la régression linéaire), et le but de la régression linéaire est de la minimiser.
[^ 1]: Cette formule est la fonction dite d'erreur carrée. En d'autres termes, c'est la méthode du carré minimum.
Comme méthode de minimisation, le dernier $ J (\ theta) $ est partiellement différencié par $ \ theta_0, \ dots, \ theta_n $, et $ \ theta $ est calculé de sorte qu'ils soient tous à 0. Par exemple, si vous vous différenciez partiellement avec un certain $ \ theta_j $:
Si vous le transformez en $ m $ lignes et le réécrivez sous la forme d'une matrice [^ 2]:
[^ 2]: Sumimasen, je l'ai plié ici. C'est ennuyeux d'écrire un par un (^ - ^;
Ceci est traduit en $ = {\ bf 0} $ et résolu pour $ \ theta $:
C'est «l'équation normale» que dit le Dr Ang. En calculant le côté droit de la dernière équation, nous pouvons trouver le paramètre optimal $ \ theta $ (qui minimise l'erreur quadratique).
Ce n'est pas toujours le cas que $ (X ^ TX) ^ {-1} $ existe pour $ X
[Pinv] d'Octave (https://www.gnu.org/software/octave/doc/interpreter/Basic-Matrix-Functions.html#index-pinv) est" Pseudoinverse "(matrice pseudo-inverse, généralisation) Également appelée matrice inverse. Ci-après, "[matrice inverse généralisée](https://ja.wikipedia.org/wiki/%E6%93%AC%E4%BC%BC%E9%80%86%E8%A1%" 8C% E5% 88% 97) "est utilisé).
Si l'argument matrice est réversible, la matrice inverse (la matrice qui correspond) est renvoyée, sinon la "matrice inverse généralisée" qui satisfait certaines propriétés est renvoyée.
Les détails de la matrice inverse généralisée sont omis ici, mais en termes simples, il s'agit d'une "matrice avec des propriétés similaires à la matrice inverse". Cela peut être défini non seulement pour les matrices carrées non régulières (= non réversibles), mais aussi pour les matrices non carrées (= $ m \ fois n $ telles que $ m \ ne n $).
Cela vous permet de retourner $ (X ^ TX) ^ {-1} $ et de calculer si $ X ^ TX $ est réversible, sinon par un "inverse généralisé" correctement calculé. C'est un mécanisme qui peut être calculé sans erreur et peut obtenir des résultats de calcul raisonnables.
Comme je l'ai écrit précédemment, la matrice inverse généralisée est également définie dans la matrice $ m \ times n $ ($ m \ ne n $). Surtout si $ m> n $ (= s'il y a plus de lignes que de colonnes = s'il s'agit d'une matrice verticalement longue), cela peut en fait s'écrire:
Je pense à la régression linéaire (multivariée). Le nombre de données a tendance à être très important par rapport au nombre de fonctionnalités (généralement). En d'autres termes, la matrice doit être une matrice verticalement longue avec plus de lignes.
Examinons maintenant à nouveau l'équation normale.
La partie autre que $ y $ sur le côté droit. C'est la même forme, non? En d'autres termes. N'est-il pas correct d'écrire ceci? Quand.
Dans le code Octave, pinv (X) * y. N'est-ce pas bien? Quand.
plus loin. Considérez l'équation suivante pendant un moment.
Si $ X $ est une matrice régulière (réversible avec une matrice carrée), elle peut être résolue avec $ \ theta = X ^ {-1} y $. Sinon, il n'y a pas de méthode ordinaire (il n'y a pas de solution satisfaisante en premier lieu, ou la solution n'est pas déterminée de manière unique). Cependant, en considérant $ \ theta = X ^ -y $, qui remplace la matrice inverse par la matrice inverse généralisée, nous pouvons trouver la solution qui minimise l'erreur (ou la norme). En d'autres termes, cette formule signifie trouver la solution d'erreur minimale pour $ X \ theta = y $.
Et Octave fournit un autre moyen de résoudre $ X \ theta = y $ tel quel. Il s'écrit «X \ y» à l'aide de l'opérateur de barre oblique inverse. Cela signifie $ X ^ {-1} y $ si $ X $ est une matrice régulière, sinon il semble faire le même calcul que $ X ^ -y $.
Rappelons maintenant la forme originale de l'équation normale.
C'est aussi sous la forme $ \ bigcirc \ theta = \ bigcirc $, n'est-ce pas? C'est sous une forme qui peut être résolue avec «○ \ ○» et l'opérateur de barre oblique inverse. Lorsqu'il est appliqué, (X '* X) \ (X' * y).
En d'autres termes. Le code d'Octave pour résoudre l'équation normale signifie qu'il existe quatre types possibles:
pinv(X'*X)*X'*ypinv(X)*yX \ y(X'*X) \ (X'*y)Parmi ceux-ci, seul l'ancien pinv (X '* X) * X' * y est introduit dans la conférence.
Pourquoi les 2e, 3e et 4e ne sont-ils pas introduits?
Donc, j'ai écrit le code réel et l'ai vérifié.
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.
Tout d'abord, en ce qui concerne les résultats du calcul, [1], [2], [3] et [4] calculent certainement la même valeur [^ 3].
[^ 3]: Après un examen plus approfondi, il semble que différentes valeurs soient en fait calculées avec une erreur d'environ 1e-14. Erreur de calcul en virgule flottante due à un processus de calcul et à une méthode de calcul différents, n'est-ce pas?
De plus, lors de la mesure du temps, il est naturel que plus la valeur de $ n $ (= le nombre de colonnes dans la matrice = le nombre de caractéristiques) est élevée, plus cela prend du temps, et à ce moment [1] pinv (X) '* X) * X' * y a un temps d'exécution plus court que [2] pinv (X) * y. En d'autres termes, la formule de calcul semble plus compliquée, mais la quantité de calcul est plus petite dans la première.
En effet, si $ X $ est une matrice $ m \ fois n $ ($ m> n $), alors $ X ^ TX $ est une matrice carrée de $ n \ fois n $, ce qui est presque réversible. Par rapport au fait que la matrice inverse est calculée à grande vitesse (le coût de multiplication de $ X ^ T $ de $ n \ fois m $ après cela est relativement faible), la matrice inverse généralisée de $ X $, qui n'est pas une matrice carrée en premier lieu, est ( De plus, c'est $ m \ fois n> n \ fois n $) Je pense que le coût du calcul est élevé.
Cependant, les résultats montrent que le temps de calcul de «X \ y» dans [3] est plus rapide.
De plus, le (X '* X) \ (X' * y) dans [4] est beaucoup plus rapide, n'est-ce pas?
[4] est plus rapide que [3] pour la même raison qu'auparavant.
Au fait, bien sûr, puisque «rand ()» est utilisé, les valeurs de «X» et «y» changent à chaque fois qu'il est exécuté, mais le résultat du calcul et le temps d'exécution sont presque les mêmes.
À ce rythme, (X '* X) \ (X' * y) semble être le meilleur, mais ...
Auparavant, j'ai essayé de vérifier avec du code en utilisant une matrice aléatoire. Essayons maintenant avec une matrice de nombres intentionnellement arrangés.
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
Les deuxième à quatrième colonnes de la matrice $ X $ sont des nombres de Fibonacci consécutifs. Autrement dit, $ x_2 + x_3 = x_4 $ (si vous écrivez la colonne $ i $ comme $ x_i $). Puisque la matrice $ m \ times 4 $ ($ m \ ge 4 $) n'a que trois colonnes indépendantes (car la quatrième colonne est représentée par une jointure linéaire des autres colonnes), le rang de la matrice est C'est $ 3 $ [^ 4].
[^ 4]: Dans une matrice X avec m lignes et n colonnes, lorsque rang (X) = min (m, n), X est dit "rang complet". Si rang (X) <min (m, n), on dit qu'il est "rang déficient".
À ce stade, $ X ^ TX $, qui est une matrice $ n \ times n $, n'est pas une matrice régulière (elle n'a pas de matrice inverse). De plus, si $ m = n $, $ X $ est lui-même une matrice carrée, mais ce n'est pas non plus une matrice régulière.
En conséquence, [4]((X '* X) \ (X' * y)) quand $ m> n $, et [3], 4 quand $ m = n $. Un avertissement a été émis lors du calcul de X \ y, (X '* X) \ (X' * y)). La solution de l'opérateur ` semble avoir un tel problème.
En revanche, dans le cas de la solution par la matrice inverse généralisée, il semble que le calcul soit effectué sans avertissement. En effet, depuis le début, "il est conçu pour effectuer des calculs raisonnables même si ce n'est pas une matrice régulière (même si ce n'est pas une matrice carrée)".
Le Dr Ang a introduit pinv (X '* X) * X' * y, qui est le plus proche de l'expression dérivée (de la transformation de l'expression d'origine) et est meilleur que pinv (X) * y. Je pense que c'est parce qu'il peut être calculé à grande vitesse.
De plus, je n'ai pas mentionné X \ y ou (X '* X) \ (X' * y) dans ces calculs, car la matrice à gauche de l'opérateur est une matrice carrée non régulière (ou " Il semble qu'il y ait un problème dans le cas de la matrice "rank drop") (car il n'y a pas lieu de s'inquiéter si pinv () est utilisé) [^ 5].
Cependant, on peut dire que c'est le cas lorsqu'il y a un problème dans la sélection des données (fonctionnalité) en premier lieu, donc si vous ne sélectionnez que des fonctionnalités indépendantes indépendamment des autres fonctionnalités, ce sera plus simple et plus rapide (X '* X). Je ne pense pas qu'il soit mauvais d'utiliser \ (X '* y) .
Ensuite, concentrez-vous sur l'analyse des données à cet effet (payez le coût du travail) et accélérez le calcul ((X '* X) \ (X' * y)), ou spécialement pour les données Faites-vous un calcul raisonnablement rapide et sans erreur sans insérer (pinv (X '* X) * X' * y)? Ce sera cette affaire.
[^ 5]: Avant de modifier l'article, j'ai demandé: "S'il n'y a pas de problèmes, X \ y est mieux? Et?", Mais c'est probablement le plus gros problème. Je l'ai résolu moi-même, mais qu'en est-il?
Dans ce genre de cours (où la langue et le système de traitement sont spécifiés) et d'étudier dans des manuels, j'essaie souvent d'écrire dans d'autres langues afin d'approfondir ma compréhension. J'ai donc également vérifié le code vérifié cette fois avec Julia (v0.3.x / 0.4.0-dev) et Python (v2.7.x / 3.x) + NumPy.
En Python, il n'y a pas d'opérateur correspondant pour X \ y, mais à la place il y a une fonction appelée numpy.linalg.solve (X, y), mais il semble que cela ne fonctionnera pas à moins que X soit une matrice carrée. ..
Julia v0.3.x / 0.4.0 Fonctionne avec les deux [^ 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]: La manière de Julia de gérer et de formater les matrices / vecteurs a été tirée de MATLAB / Octave, donc cela a fonctionné comme il était avec presque aucune réécriture.
J'ai l'impression que ce n'est pas aussi rapide que ce à quoi je m'attendais, mais je pense que c'est parce que j'ai écrit «for» directement au plus haut niveau.
Ce sera sûrement plus rapide si vous le changez pour une fonction ou un autre style d'écriture qui prête attention à la performance.
Mais (X '* X) \ (X' * y) est assez rapide pour submerger les autres! Mais ...
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)
En cas de baisse de rang.
Dans le cas d'une matrice verticalement longue, [3], [4](X \ y,(X '* X) \ (X' * y)) sont significativement différents de ceux utilisant pinv (). Il est devenu.
De plus, dans le cas d'une matrice carrée, une erreur indiquant "pas une matrice régulière" apparaît et le calcul n'est pas possible. N (> _ <)
(En passant, je mentionnerai brièvement que les résultats étaient légèrement différents entre Julia v0.3.x et v0.4.0)
Fonctionne également avec Python v2.7.x / 3.x, les résultats sont similaires:
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
Dans le ndarray de NumPy, * est le produit de chaque élément, et le produit de la matrice doit utiliser la méthode .dot (). Même les expressions compliquées deviennent encore plus compliquées ... mais vite! résoudre () est aussi rapide! Mais ...
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]])
En cas de baisse de rang.
Dans le cas de NumPy, il n'y avait pas d'erreur lors de l'utilisation de résoudre (), mais le résultat du calcul était encore très différent de celui utilisant pinv ().
Recommended Posts