4.3 Multiple lineare Regression & Normalengleichung

Mehrere Eingabegrößen

Meist hängt die Zielgröße von mehreren Merkmalen ab. Beispiel auf 📄 04_linear_methods_ws26.pdf, Folie 16: Der Advertising-Datensatz mit 200 Zeilen, Werbeausgaben für TV, Radio und Zeitung (Newspaper) sowie die erzielten Verkäufe (Sales). Das Modell:

$$\widehat{\text{sales}} = 1\cdot w_0 + w_1\,\text{TV} + w_2\,\text{Radio} + w_3\,\text{Newspaper}$$

Das ist multiple (auch multivariate) lineare Regression (multiple linear regression).

Matrixschreibweise und Designmatrix

📄 04_linear_methods_ws26.pdf, Folie 17 Man schreibt die Daten in eine Matrix, eine Zeile pro Beobachtung, eine Spalte pro Merkmal, und ergänzt vorne eine Spalte aus Einsen (Intercept). Diese Matrix heißt Designmatrix (design matrix) $X$:

$$\underbrace{\begin{pmatrix} 1 & 230{,}1 & 37{,}8 & 69{,}2\\ 1 & 44{,}5 & 39{,}3 & 45{,}1\\ 1 & 17{,}2 & 45{,}9 & 69{,}3\\ \vdots & \vdots & \vdots & \vdots \end{pmatrix}}_{X \;(200\times 4)} \underbrace{\begin{pmatrix} w_0\\ w_1\\ w_2\\ w_3\end{pmatrix}}_{\boldsymbol w} = \underbrace{\begin{pmatrix} 20{,}24\\ 12{,}43\\ 12{,}31\\ \vdots\end{pmatrix}}_{\hat{\boldsymbol y}} \qquad \text{vgl. } \boldsymbol y = \begin{pmatrix} 22{,}1\\ 10{,}4\\ 9{,}3\\ \vdots\end{pmatrix}$$

Kompakt: $X\cdot\boldsymbol w = \hat{\boldsymbol y}$. Eine einzige Matrix-Vektor-Multiplikation liefert alle 200 Vorhersagen.

Warum die Einserspalte? Für Zeile $i$ ergibt die Multiplikation $1\cdot w_0 + x_{i1}w_1 + \dots$, der Achsenabschnitt wird also automatisch addiert. So muss man $b$ nicht gesondert behandeln.

Das Modell allgemein

📄 04_linear_methods_ws26.pdf, Folie 18 Mit $d$ Merkmalen:

$$\hat y = f_\theta(x) = 1\cdot w_0 + w_1x_1 + w_2x_2 + \dots + w_dx_d = \sum_{j=1}^d w_jx_j + w_0 = \boldsymbol w^T\boldsymbol x + w_0 = \langle\boldsymbol w, \boldsymbol x\rangle + w_0$$

• Als neuronales Netz: ein Neuron mit Eingängen $1, x_1, \dots, x_d$, Gewichten $w_0, \dots, w_d$ und Summe $\Sigma$, keine Aktivierung.
• Geometrisch: Bei einem Merkmal ist $f_\theta$ eine Gerade, bei zwei Merkmalen eine Ebene im 3D-Raum, bei $d$ Merkmalen eine Hyperebene im $(d+1)$-dimensionalen Raum. Die Residuen sind die senkrechten (in $y$-Richtung gemessenen) Abstände der Punkte zur Ebene.

Notation

Die Literatur ist hier uneinheitlich 📄 04_linear_methods_ws26.pdf, Folie 19:

• $w_0$ heißt oft $b$ (Bias).
• Oft steckt der Bias mit im Gewichtsvektor, $\boldsymbol w^T = (w_0, w_1, \dots, w_d)$, und $\boldsymbol x$ bekommt eine führende 1 („erweiterte Datenmatrix“). Dann gilt einfach $\hat y = \boldsymbol w\cdot\boldsymbol x$.
• In der Statistikliteratur (Tibshirani u. a.) heißt der Gewichtsvektor $\boldsymbol\beta$.
• Häufig werden die Vektorpfeile/Fettdruck weggelassen.

Einsteinsche Summenkonvention

Eine Kurzschrift, die auf den nächsten Folien verwendet wird: Über doppelt auftretende Indizes wird summiert, das Summenzeichen lässt man weg.

$$\hat y_i = \sum_{j=0}^p X_{ij}w_j = X_{ij}w_j$$

In NumPy entspricht das np.einsum('ij,j->i', X, w), also schlicht X @ w.

Laut Feedback-Folie 36 ist die Summenkonvention „nicht soo wichtig“. Wichtig ist nur, die Herleitung unten nachvollziehen zu können. Im Folgenden werden die Summen deshalb ausgeschrieben.

Herleitung der Normalengleichung

📄 04_linear_methods_ws26.pdf, Folie 20 und 📄 04_linear_methods_ws26.pdf, Folie 21 leiten die geschlossene Lösung in vier Schritten her.

Schritt 1: Verlustfunktion

Summe der quadrierten Fehler (der Faktor $\frac1n$ ändert die Lage des Minimums nicht und wird weggelassen):

$$L = \sum_i (\hat y_i - y_i)^2 = \sum_i\Big(\sum_j X_{ij}w_j - y_i\Big)^2$$

Schritt 2: Ableiten nach einem Gewicht $w_k$

Kettenregel: äußere Ableitung von $(\cdot)^2$ ist $2(\cdot)$, dann mal innere Ableitung.

$$\frac{\partial L}{\partial w_k} = 2\sum_i\Big(\sum_j X_{ij}w_j - y_i\Big)\cdot\frac{\partial}{\partial w_k}\Big(\sum_j X_{ij}w_j - y_i\Big)$$

Die innere Ableitung: In $\sum_j X_{ij}w_j = X_{i0}w_0 + X_{i1}w_1 + \dots$ hängt nur ein Summand von $w_k$ ab, nämlich $X_{ik}w_k$. Alle anderen fallen weg, $y_i$ auch:

$$\frac{\partial}{\partial w_k}\Big(\sum_j X_{ij}w_j - y_i\Big) = \sum_j X_{ij}\,\delta_{jk} = X_{ik}$$

($\delta_{jk}$ ist das Kronecker-Delta: 1 für $j = k$, sonst 0. Die Folie schreibt $\delta_{kj}X_{ij}$ unter die Klammer.) Also:

$$\frac{\partial L}{\partial w_k} = 2\sum_i\Big(\sum_j X_{ij}w_j - y_i\Big)X_{ik}$$

Ergänzung: In Matrixform ist das die $k$-te Komponente von

$$\nabla_{\boldsymbol w}L = 2X^T(X\boldsymbol w - \boldsymbol y).$$

Das ist auch der Gradient, den man für Gradientenabstieg braucht. Intuition: $X\boldsymbol w - \boldsymbol y$ ist der Vektor der Residuen, $X^T$ gewichtet jedes Residuum mit dem zugehörigen Merkmalswert.

Die Folie bemerkt, dass die direkte Vektorschreibweise „scary and error-prone“ ist. Deshalb die Herleitung über Indizes.

Schritt 3: Gradient gleich null setzen

$$\frac{\partial L}{\partial w_k} = 0 \;\Longrightarrow\; \sum_i\sum_j X_{ik}X_{ij}w_j = \sum_i X_{ik}y_i \quad\text{für alle } k$$

Links steht $\sum_j (X^TX)_{kj}w_j$, denn $(X^TX)_{kj} = \sum_i X_{ik}X_{ij}$. Rechts steht $(X^T\boldsymbol y)_k$. Für alle $k$ gleichzeitig:

$$\boxed{X^TX\,\boldsymbol w = X^T\boldsymbol y}$$

Das ist die Normalengleichung (normal equation), ein lineares Gleichungssystem mit so vielen Gleichungen wie Gewichten.

Schritt 4: Auflösen

Falls $X^TX$ invertierbar ist:

$$\boxed{\boldsymbol w = (X^TX)^{-1}X^T\boldsymbol y}$$

Warum „Normalen“gleichung?

Ergänzung, nicht auf den Folien: Umgestellt lautet sie $X^T(\boldsymbol y - X\boldsymbol w) = \boldsymbol 0$. Der Residuenvektor steht also senkrecht (normal) auf allen Spalten von $X$. Geometrisch ist $\hat{\boldsymbol y} = X\boldsymbol w$ die Orthogonalprojektion von $\boldsymbol y$ auf den Raum, den die Spalten von $X$ aufspannen. Der kürzeste Abstand zu einem Unterraum ist immer der senkrechte.

Durchgerechnetes Beispiel

Ergänzung, nicht auf den Folien: Die drei Punkte $(1,2)$, $(2,3)$, $(3,5)$ aus 4.1:

$$X = \begin{pmatrix} 1 & 1\\ 1 & 2\\ 1 & 3\end{pmatrix},\quad \boldsymbol y = \begin{pmatrix}2\\3\\5\end{pmatrix},\quad X^TX = \begin{pmatrix} 3 & 6\\ 6 & 14\end{pmatrix},\quad X^T\boldsymbol y = \begin{pmatrix} 10\\ 23\end{pmatrix}$$

($X^TX$ enthält $n$, $\sum x_i$ und $\sum x_i^2$; $X^T\boldsymbol y$ enthält $\sum y_i$ und $\sum x_iy_i$.) Mit $\det(X^TX) = 42 - 36 = 6$:

$$(X^TX)^{-1} = \frac16\begin{pmatrix} 14 & -6\\ -6 & 3\end{pmatrix},\qquad \boldsymbol w = \frac16\begin{pmatrix} 140 - 138\\ -60 + 69\end{pmatrix} = \begin{pmatrix} 1/3\\ 3/2\end{pmatrix}$$

Also $w_0 = b = 1/3$ und $w_1 = a = 1{,}5$, wie in 4.1 und mit Gradientenabstieg in 4.2.

Wann ist $X^TX$ nicht invertierbar?

Ergänzung, nicht auf den Folien (wichtig zum Verständnis von Ridge in 4.4):

• Mehr Parameter als Datenpunkte ($p > n$): $X^TX$ ist $p\times p$, hat aber höchstens Rang $n$. Es gibt dann unendlich viele Lösungen, die die Daten perfekt treffen.
• Linear abhängige Spalten (perfekte Kollinearität), z. B. ein Merkmal doppelt, einmal in Metern und einmal in Zentimetern.
• Fast abhängige Spalten machen $X^TX$ schlecht konditioniert: Kleine Änderungen in den Daten verändern $\boldsymbol w$ stark.

In Python

Die Folien verweisen auf das Notebook linear_methods/LinearRegression.ipynb 📄 04_linear_methods_ws26.pdf, Folie 23.

Ergänzung: Typische Varianten:

import numpy as np
from sklearn.linear_model import LinearRegression

X1 = np.column_stack([np.ones(len(X)), X])        # Einserspalte anhängen

w = np.linalg.inv(X1.T @ X1) @ X1.T @ y           # wörtlich die Formel
w = np.linalg.solve(X1.T @ X1, X1.T @ y)          # besser: LGS lösen
w, *_ = np.linalg.lstsq(X1, y, rcond=None)        # am stabilsten

model = LinearRegression().fit(X, y)              # fügt den Intercept selbst hinzu
model.intercept_, model.coef_

Fallstrick: Inverse explizit berechnen

Die Formel $(X^TX)^{-1}X^T\boldsymbol y$ ist zum Herleiten gedacht. Numerisch ist es besser, das Gleichungssystem zu lösen (solve) oder direkt ein Kleinste-Quadrate-Verfahren (lstsq, QR-Zerlegung) zu verwenden, weil das Bilden von $X^TX$ die Kondition quadriert. Und bei LinearRegression aus scikit-learn keine eigene Einserspalte anhängen, der Intercept wird automatisch geschätzt.