10.3 Vetítés és legkisebb négyzetek

A legkisebb négyzetes illesztést gyakran csak képletként tanuljuk meg, pedig mögötte nagyon egyszerű geometria van: egy vektort merőlegesen vetítünk egy altérre. Ebből a képből közvetlenül következik a normálegyenlet, és az is láthatóvá válik, miért veszélyes a majdnem kollineáris adat.

10.3.1 Cél és konkrét helyzet

Adott három megfigyelt pont: $(x_1,y_1)=(0,1)$, $(x_2,y_2)=(1,3)$, $(x_3,y_3)=(2,4)$. A megfigyelési egység egy $(x,y)$ pár; a koordináták dimenzió nélküli tanulóadatok, így itt nincs mértékegység. Célunk olyan $\hat y_i = w_0 + w_1 x_i$ egyenest illeszteni, amely minimalizálja a

$$Q(w_0,w_1)=\sum_{i=1}^{3}\bigl(y_i-\hat y_i\bigr)^2 =\sum_{i=1}^{3}\bigl(y_i-w_0-w_1 x_i\bigr)^2$$

négyzetes hibát. Mielőtt tovább olvasnál, gyors becslésként rajzold fel a pontokat: az első szakaszon $2$, a másodikon $1$ az $y$-emelkedés, ezért a legkisebb négyzetes meredekség kompromisszumként $w_1=1{,}5$ körül várható, a tengelymetszet pedig valamivel $1$ fölé esik. A levezetés után ellenőrizzük ezt a sejtést.

10.3.2 Merőlegesség és vetítés

Írjuk a modellt mátrixalakba, ahogyan 10.2 a lineáris transzformációt bevezette:

$$\mathbf X= \begin{pmatrix} 1 & 0\\ 1 & 1\\ 1 & 2 \end{pmatrix},\qquad \mathbf y= \begin{pmatrix} 1\\ 3\\ 4 \end{pmatrix},\qquad \mathbf w= \begin{pmatrix} w_0\\ w_1 \end{pmatrix}.$$

A design mátrix $\mathbf X\in\mathbb R^{3\times 2}$: minden sora egy megfigyelés, az első oszlop az intercept csupa $1$-ese, a második az $x_i$ értékek oszlopa. A súlyvektor $\mathbf w\in\mathbb R^2$, a predikció pedig $\hat{\mathbf y}=\mathbf X\mathbf w\in\mathbb R^3$. A $\hat{\mathbf y}$ vektor tehát az $\mathbf X$ oszlopainak lineáris kombinációja, ezért eleme az oszloptérnek (column space):

$$C(\mathbf X)=\{\mathbf X\mathbf v:\mathbf v\in\mathbb R^2\}\subseteq\mathbb R^3.$$

Jelen példában ez egy kétdimenziós sík a háromdimenziós megfigyelési térben, amelyet a $\mathbf c_1=(1,1,1)^T$ és $\mathbf c_2=(0,1,2)^T$ oszlopvektorok feszítenek ki. A legkisebb négyzetes illesztés geometriai feladata: válasszuk $\mathbf C(\mathbf X)$-ből az $\mathbf y$-hoz legközelebbi pontot.

A reziduum (residual) a megfigyelés és a predikció eltérése:

$$\mathbf r=\mathbf y-\hat{\mathbf y}=\mathbf y-\mathbf X\mathbf w,$$

ahol $\mathbf r\in\mathbb R^3$, és koordinátánként $r_i=y_i-\hat y_i$.

A kulcsfelismerés: akkor lesz $\hat{\mathbf y}$ a legközelebbi pont, ha az $\mathbf r$ reziduum merőleges az oszloptérre. Ha ugyanis $\mathbf r$ merőleges $\mathbf C(\mathbf X)$ minden vektorára, akkor bármely más $\mathbf v\in C(\mathbf X)$ esetén Pitagorasz tétele szerint

$$\|\mathbf y-\mathbf v\|^2 =\|\mathbf r\|^2+\|\hat{\mathbf y}-\mathbf v\|^2 \ge\|\mathbf r\|^2,$$

tehát $\hat{\mathbf y}$ valóban minimális távolságú. A merőlegességet elég az altér bázisvektoraira, vagyis $\mathbf X$ oszlopaira megkövetelni:

$$\mathbf X^T\mathbf r=0.$$

Ezt normálegyenletnek (normal equations) nevezzük. Behelyettesítve:

$$\mathbf X^T(\mathbf y-\mathbf X\mathbf w)=0 \implies \mathbf X^T\mathbf X\mathbf w=\mathbf X^T\mathbf y.$$

A shape-ek: $\mathbf X^T\mathbf X$ mérete $2\times 2$, $\mathbf X^T\mathbf y$ mérete $2\times 1$, tehát a megoldandó rendszer két egyenlet két ismeretlenre.

Ha $\mathbf X$ oszlopai lineárisan függetlenek, akkor $\mathbf X^T\mathbf X$ invertálható, és a legkisebb négyzetes megoldás egyértelmű:

$$\hat{\mathbf w}=(\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y.$$

A predikció ekkor az $\mathbf y$ vektor merőleges vetítése (projection) az oszloptérre:

$$\hat{\mathbf y}=\mathbf X(\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y =\mathbf P_{\mathbf X}\mathbf y.$$

A $\mathbf P_{\mathbf X}$ mátrix a vetítőmátrix (projection matrix), mérete $3\times 3$; szimmetrikus és idempotens, vagyis $\mathbf P_{\mathbf X}^2=\mathbf P_{\mathbf X}$. A vetítés fogalmát kis vektoros ábrán is érdemes látni. Az alábbi Matplotlib-kód egy $\mathbb R^2$-beli altérre vetítést rajzol; a piros reziduum merőleges a zöld altéregyenesre.

import matplotlib.pyplot as plt
import numpy as np

a = np.array([3.0, 1.0])          # az egydimenziós altér bázisvektora
y = np.array([2.0, 3.0])
y_hat = (a @ y) / (a @ a) * a     # merőleges vetítés
r = y - y_hat

fig, ax = plt.subplots()
ax.quiver(0, 0, *y,      angles='xy', scale_units='xy', scale=1,
          color='black', label='y')
ax.quiver(0, 0, *a,      angles='xy', scale_units='xy', scale=1,
          color='green', label='altér')
ax.quiver(0, 0, *y_hat,  angles='xy', scale_units='xy', scale=1,
          color='blue', label='vetített y')
ax.quiver(*y_hat, *r,    angles='xy', scale_units='xy', scale=1,
          color='red', label='reziduum')
ax.set_xlim(-1, 4); ax.set_ylim(-1, 4)
ax.set_aspect('equal'); ax.legend()

Az ábra a futás után szemlélteti, hogy a reziduum derékszöget zár be az altérrel. A kód itt nem futtatott; a leírt geometriai viszony a vetítési képletből következik.

10.3.3 Rang, pszeudoinverz és kondíció

A rang (rank) az $\mathbf X$ oszlopvektorai között található lineárisan független vektorok maximális száma; NumPy-ban np.linalg.matrix_rank. Ha a rang egyenlő az oszlopok számával ($\operatorname{rank}(\mathbf X)=p$), akkor $\mathbf X$ teljes oszloprangú, és $\mathbf X^T\mathbf X$ invertálható. Ha $\operatorname{rank}(\mathbf X)<p$, akkor az egyenletrendszernek több $\mathbf w$ megoldása is ugyanazt a $\hat{\mathbf y}$ predikciót adja; a predikció egyedi, a súlyvektor nem.

A pszeudoinverz (pseudoinverse, Moore–Penrose-inverz) ilyenkor is előállít egy egyértelmű legkisebb négyzetes megoldást: ez a legkisebb normájú megoldás. NumPy-ban np.linalg.pinv(X). A pszeudoinverz szinguláris felbontásra épül, és ranghiányos mátrixra is működik. A teljes rangú esetben ugyanazt adja, mint a normálegyenlet-megoldás, de numerikusan általában stabilabb.

A kondíciószám (condition number) azt méri, mennyire érzékeny a megoldás a bemeneti perturbációkra:

$$\kappa_2(\mathbf X)=\frac{\sigma_{\max}(\mathbf X)}{\sigma_{\min}(\mathbf X)},$$

ahol $\sigma_{\max}$ és $\sigma_{\min}$ az $\mathbf X$ legnagyobb és legkisebb szinguláris értéke; NumPy-ban np.linalg.cond. Ha $\kappa$ közel van $1$-hez, a feladat stabil; ha $\sigma_{\min}$ nagyon kicsi, $\kappa$ nagy, és kis adatváltozás is nagy súlyváltozást okozhat. Fontos: a normálegyenlet mátrixának kondíciója $\kappa(\mathbf X^T\mathbf X)=\kappa(\mathbf X)^2$, tehát a normálegyenletben még erősebb az instabilitás, mint magában az $\mathbf X$-ben.

10.3.4 Teljes minimális munkapélda

Az alábbi kód a hárompontos feladatot oldja meg, és minden fontos invariánst ellenőriz. A kód a kezdőlabor alap numerikus csomagjain kívül mást nem igényel. Itt nem futtatjuk; a kimenethez tartozó várt értékeket kézi számításból adjuk meg.

import numpy as np

X = np.array([
    [1.0, 0.0],
    [1.0, 1.0],
    [1.0, 2.0],
])
y = np.array([1.0, 3.0, 4.0])

w_lstsq, residuals, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
y_hat = X @ w_lstsq
r = y - y_hat

# Kézi normálegyenlet, az elméleti képlet ellenőrzésére.
w_direct = np.linalg.solve(X.T @ X, X.T @ y)
w_pinv = np.linalg.pinv(X) @ y

print("w_lstsq =", w_lstsq)
print("y_hat =", y_hat)
print("r =", r)
print("X.T @ r =", X.T @ r)
print("rank =", rank)
print("singular_values =", singular_values)
print("cond(X) =", np.linalg.cond(X))
print("w_direct =", w_direct)
print("w_pinv =", w_pinv)

A kézi vezetésből elméletileg várt kimenet: $\hat{\mathbf w}=(7/6,\ 3/2)^T$, azaz w_lstsq = [1.16666667, 1.5]. A predikciók: $\hat{\mathbf y}=(7/6,\ 8/3,\ 25/6)^T$, azaz körülbelül [1.1667, 2.6667, 4.1667]. A reziduumok: $(-1/6,\ 1/3,\ -1/6)$. Ellenőrizzük az invariánsokat:

$$\mathbf X^T\mathbf r =\begin{pmatrix} 1&1&1\\ 0&1&2 \end{pmatrix} \begin{pmatrix} -1/6\\ 1/3\\ -1/6 \end{pmatrix} = \begin{pmatrix} 0\\ 0 \end{pmatrix},$$

tehát a reziduum valóban merőleges mindkét oszlopra. A rang $2$, mert a két oszlop nem arányos. A szinguláris értékek $\sqrt{4+\sqrt{10}}\approx 2{,}676$ és $\sqrt{4-\sqrt{10}}\approx 0{,}915$, ezért $\kappa(\mathbf X)\approx 2{,}92$, ami stabil kis feladat. A solve, a lstsq és a pinv által adott súlyok ezen a jól kondicionált példán gyakorlatilag megegyeznek.

10.3.5 Majdnem kollineáris változók és az instabilitás

Hozzunk létre olyan adatot, ahol a második magyarázó változó majdnem konstansszorosa az elsőnek. Ez a gyakorlatban akkor fordul elő, ha két feature erősen összefügg. Az alábbi kód a jól értelmezhető $y=2x_1+x_2+0{,}5$ kapcsolatot használja; a zaj miatt a pontos súlyok nem reprodukálhatók gépfüggetlenül, ezért az itt közölt kimenet szemléltető jellegű, nem futtatott mérési jegyzék.

rng = np.random.default_rng(42)
n = 40
x1 = np.linspace(-2.0, 2.0, n)
x2 = 3.0 * x1 + 0.005 * rng.normal(size=n)
y_true = 2.0 * x1 + 1.0 * x2 + 0.5
y = y_true + 0.1 * rng.normal(size=n)
X = np.column_stack([np.ones(n), x1, x2])

w, *_ = np.linalg.lstsq(X, y, rcond=None)

# Kis perturbáció ugyanarra a design mátrixra.
y_perturbed = y + 0.02 * rng.normal(size=n)
w2, *_ = np.linalg.lstsq(X, y_perturbed, rcond=None)

print("cond(X) =", np.linalg.cond(X))
print("w =", w)
print("w_perturbed =", w2)
print("||w - w_perturbed|| =", np.linalg.norm(w - w2))
print("||y_hat - y_hat_perturbed|| =",
      np.linalg.norm(X @ w - X @ w2))

A megadott konstrukcióban $x_2$ majdnem $3x_1$, tehát a második és harmadik oszlop majdnem kollineáris. A kondíciószám várhatóan több száz vagy ezer körüli. A látványos jelenség: a $\mathbf w$ és $\mathbf w_2$ súlyvektorok önmagukban jelentősen eltérhetnek, mert a két súly kompenzálja egymást, miközben a predikciók különbsége kicsi marad. Ez nem hiba a kódban: az adat nem azonosítja egyértelműen a két együtthatót, csak a közelítő összegüket. A stabil eszköz itt is lstsq; a predikció marad stabil, a súlyvektorok értelmezése viszont óvatosságot igényel.

10.3.6 Feladatok

1. Értelmezés. A hárompontos példában számold ki a $\|\mathbf y\|^2$, $\|\hat{\mathbf y}\|^2$ és $\|\mathbf r\|^2$ értékeket, és igazold a Pitagorasz-összefüggést: $\|\mathbf y\|^2=\|\hat{\mathbf y}\|^2+\|\mathbf r\|^2$. Magyarázd el, miért pontosan a vetítésből következik ez.

2. Megvalósítás. Illessz egyenest a $(-1,0)$, $(0,2)$, $(1,1)$ pontokra. Írj NumPy-kódot lstsq-val, ellenőrizd a normálegyenletet solve-val, és írd ki a reziduumvektorra a $\mathbf X^T\mathbf r=0$ feltételt. Add meg a kondíciószámot is.

3. Hibadiagnózis. Egy kolléga megoldásként w_bad = np.linalg.inv(X.T @ X) @ X.T @ y sort használ. A majdnem kollineáris példán a predikciója mégsem rossz. Milyen diagnosztikai kísérlettel tudnád kimutatni, hogy a módszer instabilabb a kelleténél? Javasolj legalább egy konkrét numerikus összehasonlítást.

4. Megváltoztatott adat. A három eredeti ponthoz adj hozzá egy negyediket, amelyik a $(1,3)$ pont pontos duplikátuma. Futtasd az illesztést, és válaszold meg: változik-e a design mátrix rangja, a becsült egyenes, és a reziduumok négyzetösszege? Hogyan kapcsolódik ez ahhoz a duplikációs szerződéshez, amelyről 6.3 szólt?

10.3.7 Megoldások és indoklás

1. A kiszámolt vektorokkal:

$$\|\mathbf y\|^2=1^2+3^2+4^2=26,$$

$$\|\hat{\mathbf y}\|^2=\left(\frac76\right)^2+\left(\frac83\right)^2+\left(\frac{25}{6}\right)^2 =\frac{49+256+625}{36} =\frac{930}{36} =\frac{155}{6},$$

$$\|\mathbf r\|^2=\left(-\frac16\right)^2+\left(\frac13\right)^2+\left(-\frac16\right)^2 =\frac{1+4+1}{36} =\frac16,$$

és $\frac{155}{6}+\frac16=\frac{156}{6}=26$, tehát a Pitagorasz-összefüggés teljesül. A magyarázat: $\hat{\mathbf y}$ az $\mathbf y$ merőleges vetülete az oszloptérre, ezért $\mathbf r$ merőleges $\hat{\mathbf y}$-ra; az általános Pitagorasz-tételből épp ez az invariáns következik.

2. A feladat mátrixai:

$$\mathbf X=\begin{pmatrix}1&-1\\1&0\\1&1\end{pmatrix},\qquad \mathbf y=\begin{pmatrix}0\\2\\1\end{pmatrix}.$$

A normálegyenlet:

$$\mathbf X^T\mathbf X=\begin{pmatrix}3&0\\0&2\end{pmatrix},\qquad \mathbf X^T\mathbf y=\begin{pmatrix}3\\1\end{pmatrix},$$

ezért $\hat{\mathbf w}=(1,\ 0{,}5)^T$. A predikciók $(-1;\ 0{,}5)$, $(0;\ 1)$, $(1;\ 1{,}5)$, a reziduumok $(-0{,}5;\ 1;\ -0{,}5)$, és $\mathbf X^T\mathbf r=(0,0)^T$. A szinguláris értékek $\sqrt3$ és $\sqrt2$, ezért $\kappa(\mathbf X)=\sqrt{3/2}\approx1{,}225$; ez nagyon stabil. A kód:

X = np.array([[1.0, -1.0], [1.0, 0.0], [1.0, 1.0]])
y = np.array([0.0, 2.0, 1.0])
w, *_ = np.linalg.lstsq(X, y, rcond=None)
r = y - X @ w
print("w =", w)
print("X.T @ r =", X.T @ r)
print("cond(X) =", np.linalg.cond(X))
print("solve check =", np.linalg.solve(X.T @ X, X.T @ y))

3. A diagnosztikai kísérlet: legyártható egy kis perturbált $\mathbf y^*=\mathbf y+\varepsilon$, ahol $\varepsilon$ nagyon kicsi, például a fenti majdnem kollineáris példában $0{,}02$-szeres zaj. Hasonlítsd össze a $\|\hat{\mathbf w}-\hat{\mathbf w}^{*}\|$ eltérést a lstsq megoldás és az inv(X.T @ X) megoldás között. Az explicit inverz a normálegyenlet $\kappa(\mathbf X)^2$-szeres kondíciójával dolgozik, ezért a súlyeltérés ott általában sokkal nagyobb. A javítás: np.linalg.lstsq(X, y, rcond=None) vagy legalább np.linalg.solve(X.T @ X, X.T @ y); utóbbi még mindig a normálegyenletet oldja, ezért a solve is lehet érzékeny, de elkerüli az explicit inverzképzés plusz műveleteit. A legjobb a QR- vagy SVD-alapú lstsq vagy pinv.

4. A kibővített adat:

$$\mathbf X_4=\begin{pmatrix}1&0\\1&1\\1&2\\1&1\end{pmatrix},\qquad \mathbf y_4=\begin{pmatrix}1\\3\\4\\3\end{pmatrix}.$$

A rang $2$ marad, mert a duplikált sor nem hoz be új irányt. A normálegyenlet new mátrixa:

$$\mathbf X_4^T\mathbf X_4=\begin{pmatrix}4&4\\4&6\end{pmatrix},\qquad \mathbf X_4^T\mathbf y_4=\begin{pmatrix}11\\14\end{pmatrix},$$

amelynek megoldása $\hat{\mathbf w}=(5/4,\ 3/2)^T$. A meredekség nem változott, az intercept azonban $7/6\approx1{,}167$-ről $1{,}25$-ra nőtt, mert a $(1,3)$ pont kétszeres súlyt kapott. Ez összhangban van 6.3 tanulságával: a duplikátumot explicit döntéssel kell kezelni, mert ha nem valódi ismételt mérés, akkor torzítja a súlyozást.

10.3.8 Gyakori tévedés és ellenpélda

A tankönyvi képletben szereplő $(\mathbf X^T\mathbf X)^{-1}$ miatt sokan azt gondolják, hogy a helyes megvalósítás az explicit inverz:

w_bad = np.linalg.inv(X.T @ X) @ X.T @ y

Ez a megoldás kis, jól kondicionált mátrixokra gyakorlatilag rendben van, de nem általánosan a legjobb numerikus választás. Az ellenpélda: ha $\kappa(\mathbf X)=10^4$, akkor $\kappa(\mathbf X^T\mathbf X)=10^8$, tehát a normálegyenlet már eleve nyolc nagyságrendnyi érzékenységet visz a rendszerbe, és az explicit inverz további kerekítést ad. A lstsq QR- vagy SVD-felbontáson alapulva elkerüli az $\mathbf X^T\mathbf X$ explicit képzését. Ez ugyanaz a véges pontossági tanulság, amelyet 9.3 a lépésköz-optimumnál bemutatott: a kerekítési hiba a rossz képletből előbb-utóbb elrontja a numerikus eredményt.

10.3.9 Önellenőrzés

• Tudom-e bizonyítani az $\mathbf X^T\mathbf r=0$ feltételből, hogy a vetített pont minimális távolságú?
• Meg tudom-e mondani egy adott design mátrixról, mikor egyértelmű a súlyvektor, és mikor csak a predikció egyedi?
• Le tudom-e futattatni a majdnem kollineáris kísérletet úgy, hogy a seed és a seed-függő pontatlanság jelölt, ne pedig hallgatólagos?
• El tudom-e magyarázni, miért nézhet stabilnak egy instabil megoldó a predikción, miközben a súlyvektor erősen mozog?

10.3.10 Rövid összefoglaló

A legkisebb négyzetes illesztés geometriailag az $\mathbf y$ vektor merőleges vetítése az $\mathbf X$ mátrix oszlopterére; a reziduum merőlegessége adja a $\mathbf X^T\mathbf X\mathbf w=\mathbf X^T\mathbf y$ normálegyenletet. Teljes oszloprang esetén a súlyvektor egyértelmű, ranghiány esetén a pszeudoinverz a legkisebb normájú megoldást választja. A kondíciószám a megoldás érzékenységét méri; a majdnem kollineáris változóknál a súlyok instabilak lehetnek, miközben a predikció stabil marad. Tényleges számításhoz az explicit mátrixinverz helyett np.linalg.lstsq vagy np.linalg.pinv ajánlott. A következő alfejezetben a szinguláris felbontás és a dimenziócsökkentés épít erre a vetítési képre.