Dlaczego nie skorzystać z „równań normalnych”, aby znaleźć proste współczynniki najmniejszych kwadratów?

17

Widziałem tę listę tutaj i nie mogłem uwierzyć, że istnieje tak wiele sposobów rozwiązania najmniejszych kwadratów. „Normalne równania” na Wikipedii wydawał się być dość prosty sposób do

\begin{aligned} \hat{α} & = \bar{y} - \hat{β} \bar{x}, \\ \hat{β} & = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y})}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}} \end{aligned}

${\begin{aligned}{\hat {\alpha }}&={\bar {y}}-{\hat {\beta }}\,{\bar {x}},\\{\hat {\beta }}&={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}\end{aligned}}$

Dlaczego więc ich nie użyć? Zakładam, że musi istnieć problem obliczeniowy lub precyzyjny, biorąc pod uwagę, że w pierwszym linku powyżej Mark L. Stone wspomina, że SVD lub QR są popularnymi metodami w oprogramowaniu statystycznym i że równania normalne są „OKAZISTE z punktu widzenia niezawodności i dokładności numerycznej”. Jednak w poniższym kodzie równania normalne dają mi dokładność do ~ 12 miejsc po przecinku w porównaniu z trzema popularnymi funkcjami pytona: polyfit numpy ; regres linowy Scipy'ego ; i scikit-learn'S regresja liniowa .

Co ciekawsze, metoda równania normalnego jest najszybsza, gdy n = 100000000. Czasy obliczeniowe dla mnie wynoszą: 2,5 s dla regresji liniowej; 12,9s dla polyfit; 4.2s dla regresji liniowej; i 1,8 s dla równania normalnego.

Kod:

import numpy as np
from sklearn.linear_model import LinearRegression
from scipy.stats import linregress
import timeit

b0 = 0
b1 = 1
n = 100000000
x = np.linspace(-5, 5, n)
np.random.seed(42)
e = np.random.randn(n)
y = b0 + b1*x + e

# scipy                                                                                                                                     
start = timeit.default_timer()
print(str.format('{0:.30f}', linregress(x, y)[0]))
stop = timeit.default_timer()
print(stop - start)

# numpy                                                                                                                                      
start = timeit.default_timer()
print(str.format('{0:.30f}', np.polyfit(x, y, 1)[0]))
stop = timeit.default_timer()
print(stop - start)

# sklearn                                                                                                                                    
clf = LinearRegression()
start = timeit.default_timer()
clf.fit(x.reshape(-1, 1), y.reshape(-1, 1))
stop = timeit.default_timer()
print(str.format('{0:.30f}', clf.coef_[0, 0]))
print(stop - start)

# normal equation                                                                                                                            
start = timeit.default_timer()
slope = np.sum((x-x.mean())*(y-y.mean()))/np.sum((x-x.mean())**2)
stop = timeit.default_timer()
print(str.format('{0:.30f}', slope))
print(stop - start)

regression least-squares scikit-learn

— Oliver Angelil
źródło

Odpowiedzi są dość przesadzone. To nie jest takie straszne, jeśli po prostu unikniesz jawnego obliczenia odwrotności.

— matematyk

3

Kilka uwag na temat prędkości: patrzysz tylko na jedną zmienną towarzyszącą, więc koszt odwrócenia macierzy wynosi zasadniczo 0. Jeśli spojrzysz na kilka tysięcy zmiennych towarzyszących, to się zmieni. Po drugie, ponieważ masz tylko jedną zmienną towarzyszącą, munging danych zajmuje tak naprawdę dużo czasu u zapakowanych konkurentów (ale to powinno skalować się tylko liniowo, więc nie jest to wielka sprawa). Normalne rozwiązanie równań nie przetwarza danych, więc jest szybsze, ale nie ma żadnych dzwonków ani gwizdków związanych z wynikami.

— Cliff AB

23

$Ax \approx b$ $A$ $A^TA$ $log_{10}(cond)$ $A^TA$ $A^TAx = A^Tb$ cyfr dokładności. Tzn., Utworzenie równania normalnego podwoiło liczbę utraconych cyfr dokładności, od samego początku. $log_{10}(cond(A^TA)) = 2*log_{10}(cond(A))$

$10^8$ $10^{16}$

Czasami unika się normalnych równań, a czasem nie.

— Mark L. Stone
źródło

2

Prostszym sposobem, aby to zobaczyć (jeśli nie znasz / nie obchodzi Cię liczba warunków), jest to (zasadniczo) mnożenie czegoś samodzielnie („kwadratura”), co oznacza, że możesz spodziewać się utraty około połowy swoich bitów precyzja. (Powinno to być bardziej oczywiste, jeśli A jest skalarem, i powinno być łatwo zauważyć, że utworzenie macierzy A tak naprawdę nie zmienia podstawowego problemu.)

— user541686 27.04.2018

Czy oprócz różnic w dokładności istnieje duża różnica prędkości między równaniami QR i normalnymi? ponieważ w tym drugim przypadku możesz rozwiązać (X'X) -1 * X'Y, co jest wolne z powodu odwrotności? Pytam, bo nie jestem pewien, jak działa QR, więc może jest coś tak powolnego jak odwracanie matrycy. A może jedynym punktem do rozważenia jest utrata dokładności?

— Simon

4

A^{T} A

$A^TA$

A^{T} b

$A^Tb$

8

Jeśli musisz rozwiązać tylko ten jeden zmienny problem, skorzystaj z formuły. Nie ma w tym nic złego. Widziałem na przykład, jak piszesz kilka wierszy kodu w ASM dla urządzenia osadzonego. W rzeczywistości użyłem tego rodzaju rozwiązania w niektórych sytuacjach. Oczywiście nie trzeba przeciągać dużych bibliotek statystycznych, aby rozwiązać ten jeden mały problem.

Niestabilność numeryczna i wydajność są problemami większych problemów i ogólnych ustawień. Jeśli rozwiążesz wielowymiarowe najmniejsze kwadraty itp. W przypadku ogólnego problemu oczywiście tego nie użyłbyś.

— Aksakal
źródło

0

Żaden nowoczesny pakiet statystyczny nie rozwiązałby regresji liniowej za pomocą równań normalnych. Równania normalne istnieją tylko w księgach statystycznych.

Równań normalnych nie należy używać, ponieważ obliczanie odwrotności macierzy jest bardzo problematyczne.

Po co używać opadania gradientu do regresji liniowej, gdy dostępne jest rozwiązanie matematyczne w formie zamkniętej?

... chociaż dostępne jest bezpośrednie równanie normalne. Zauważ, że w równaniu normalnym należy odwrócić macierz. Teraz odwracanie macierzy kosztuje O (N3) do obliczeń, gdzie N jest liczbą wierszy w macierzy X, tj. Obserwacjami. Co więcej, jeśli X jest źle warunkowany, spowoduje to błędy obliczeniowe w oszacowaniu ...

— SmallChess
źródło