Mengapa tidak menggunakan "persamaan normal" untuk menemukan koefisien kuadrat terkecil sederhana?

17

Saya melihat daftar ini di sini dan tidak percaya ada begitu banyak cara untuk menyelesaikan kuadrat terkecil. "Persamaan normal" di Wikipedia tampaknya merupakan cara yang cukup lurus ke depan:

\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}}$

Jadi mengapa tidak menggunakannya saja? Saya berasumsi pasti ada masalah komputasi atau presisi mengingat bahwa dalam tautan pertama di atas Mark L. Stone menyebutkan bahwa SVD atau QR adalah metode populer dalam perangkat lunak statistik dan bahwa persamaan normal adalah "TERRIBLE dari keandalan dan sudut pandang akurasi numerik". Namun , dalam kode berikut, persamaan normal memberi saya akurasi hingga ~ 12 tempat desimal bila dibandingkan dengan tiga fungsi python populer: polyfit numpy ; linregress scipy ; dan LinearRegression scikit- learn .

Yang lebih menarik adalah bahwa metode persamaan normal adalah yang tercepat ketika n = 100000000. Waktu komputasi bagi saya adalah: 2,5 detik untuk linregress; 12.9 untuk polyfit; 4.2s untuk LinearRegression; dan 1,8 untuk persamaan normal.

Kode:

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
sumber

Jawabannya cukup dilebih-lebihkan. Tidak begitu mengerikan jika Anda hanya menghindari secara eksplisit menghitung kebalikannya.

mathreadler

3

Beberapa catatan tentang kecepatan: Anda hanya melihat kovariat tunggal, sehingga biaya inversi matriks pada dasarnya adalah 0. Jika Anda melihat beberapa ribu kovariat, itu akan berubah. Kedua, karena Anda hanya memiliki kovariat tunggal, data-munging adalah apa yang sebenarnya menghabiskan banyak waktu Anda di pesaing yang dikemas (tetapi ini seharusnya hanya skala secara linear, jadi bukan masalah besar). Solusi persamaan normal tidak melakukan munging data, jadi lebih cepat, tetapi tidak memiliki lonceng dan peluit yang melekat pada hasil itu.

Cliff AB

22

$Ax \approx b$ $A$ $A^TA$ $log_{10}(cond)$ $A^TA$ $A^TAx = A^Tb$ $log_{10}(cond(A^TA)) = 2*log_{10}(cond(A))$

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

Kadang-kadang Anda lolos dengan persamaan Normal, dan kadang-kadang tidak.

Mark L. Stone
sumber

2

Cara yang lebih sederhana untuk melihatnya (jika Anda tidak tahu / peduli tentang angka-angka kondisi) adalah bahwa Anda (pada dasarnya) mengalikan sesuatu dengan sendirinya ("mengkuadratkannya"), artinya Anda dapat kehilangan sekitar setengah dari presisi (Ini seharusnya lebih jelas jika A adalah skalar, dan harus mudah dilihat bahwa membuat matriks A tidak benar-benar mengubah masalah yang mendasarinya.)

user541686

Selain perbedaan dalam akurasi, apakah ada juga perbedaan kecepatan besar antara QR dan persamaan normal? karena dalam kasus terakhir Anda mungkin memecahkan (X'X) -1 * X'Y, yang lambat karena kebalikannya? Saya bertanya karena saya tidak yakin bagaimana QR bekerja, jadi mungkin ada sesuatu di sana yang sama lambatnya dengan membalikkan sebuah matriks. Atau satu-satunya titik pertimbangan hilangnya akurasi?

Simon

4

A^{T} A

$A^TA$

A^{T} b

$A^Tb$

8

Jika Anda hanya perlu menyelesaikan masalah variabel yang satu ini, maka silakan dan gunakan rumus. Tidak ada yang salah dengan itu. Saya bisa melihat Anda menulis beberapa baris kode di ASM untuk perangkat yang disematkan, misalnya. Bahkan, saya menggunakan solusi semacam ini dalam beberapa situasi. Anda tidak perlu menyeret pustaka statistik besar hanya untuk menyelesaikan satu masalah kecil ini, tentu saja.

Ketidakstabilan numerik dan kinerja adalah masalah masalah yang lebih besar dan pengaturan umum. Jika Anda menyelesaikan kuadrat multivariat, dll. Untuk masalah umum, tentu saja Anda tidak akan menggunakan ini.

Aksakal
sumber

0

Tidak ada paket statistik modern yang akan menyelesaikan regresi linier dengan persamaan normal. Persamaan normal hanya ada di buku statistik.

Persamaan normal tidak boleh digunakan sebagai komputasi invers dari matriks sangat bermasalah.

Mengapa menggunakan gradient descent untuk regresi linier, ketika solusi matematika bentuk-tertutup tersedia?

... meskipun persamaan normal langsung tersedia. Perhatikan bahwa dalam persamaan normal kita harus membalikkan sebuah matriks. Sekarang pembalikan biaya matriks O (N3) untuk perhitungan di mana N adalah jumlah baris dalam matriks X yaitu pengamatan. Selain itu, jika X tidak dikondisikan maka akan membuat kesalahan perhitungan dalam estimasi ...

Catur kecil
sumber

Mengapa tidak menggunakan "persamaan normal" untuk menemukan koefisien kuadrat terkecil sederhana?

Jawaban: