Regresja i kombinacje liniowe || Matematyka ∩ Programowanie

Ostatnio pomagałem w kursie algebry liniowej organizowanym przez Tai-Danae Bradleya i Jacka Hidary’ego i jedno z pytań, które pojawiło się kilka razy, brzmiało: „Dlaczego programiści powinni przejmować się koncepcją kombinacji liniowej?”

Dla tych, którzy nie wiedzą, dane wektory $ v_1, \dots, v_n$, a kombinacja liniowa wektorów to wybór niektórych współczynników $ a_i$, za pomocą których można ważyć wektory w sumie $ v = \sum_{i=1}^n a_i v_i$.

Muszę przyznać, że książki matematyczne słabo radzą sobie z przedstawianiem tej koncepcji jako czegoś więcej niż teoretycznego — być może kombinacje liniowe są potrzebne tylko do celów dowodowych, podczas gdy prawdziwym sednem jest mnożenie macierzy i iloczyny krzyżowe. Ale nie, kombinacje liniowe naprawdę leżą u podstaw wielu praktycznych zastosowań.

W niektórych przypadkach cały cel algorytmu polega na znalezieniu „użytecznej” kombinacji liniowej zbioru wektorów. Wektory to elementy składowe (często podstawa przestrzeni wektorowej lub podprzestrzeni), a zbiór kombinacji liniowych to legalne sposoby łączenia bloków. Prostsze bloki obsługują łatwiejsze i wydajniejsze algorytmy, ale ich liniowe kombinacje są mniej wyraziste. Stąd kompromis.

Konkretnym przykładem jest regresja. Większość ludzi myśli o regresji w kategoriach regresji liniowej. Szukasz funkcji liniowej, takiej jak $ y = mx+b$, która dobrze przybliża niektóre dane. W przypadku wielu zmiennych mamy np. $ \mathbf{x} = (x_1, x_2, x_3)$ jako wektor zmiennych wejściowych i $ \mathbf{w} = (w_1, w_2, w_3)$ jako wektor wag, a funkcją jest $ y = \mathbf{w}^T \mathbf{x} + b$.

Aby uniknąć przesunięcia o $ b$ (co sprawia, że funkcja jest afiniczna, a nie czysto liniowa; formuły funkcji czysto liniowych są łatwiejsze w użyciu, ponieważ przesunięcie jest jak nieznośny przypadek specjalny, z którym trzeba stale się liczyć), autorzy często dodają fałszywą zmienną wejściową $ x_0$, która jest zawsze stała na 1, i ponownie oznaczają $ b$ jako $ w_0$, aby otrzymać $ y = \mathbf{w}^T \mathbf{x} = \sum_i w_i x_i$ jako ostateczna forma. Problem optymalizacji do rozwiązania wygląda następująco, gdzie przybliżony zestaw danych wynosi $ \{ \mathbf{x}_1, \dots, \mathbf{x}_k \}$.

$$\displaystyle \min_w \sum_{i=1}^k (y_i – \mathbf{w}^T \mathbf{x}_i)^2$$

W tym przypadku wyuczona funkcja – wynik regresji – nie wygląda na kombinację liniową. Technicznie rzecz biorąc, jest, ale nie w interesujący sposób.

Staje się to bardziej oczywiste w przypadku kombinacji liniowych, gdy próbuje się modelować nieliniowość. Pomysł polega na zdefiniowaniu klasy funkcji tzw funkcje podstawowe $ B = \{ f_1, \dots, f_m \mid f_i: \mathbb{R}^n \to \mathbb{R} \}$ i pozwól, aby twoim przybliżeniem była dowolna liniowa kombinacja funkcji w $ B$, tj. dowolna funkcja z zakresu B.

$$\displaystyle \hat{f}(\mathbf{x}) = \sum_{i=1}^m w_i f_i(\mathbf {x})$$

Ponownie, zamiast ważyć każdą współrzędną wektora wejściowego za pomocą $ w_i$, ważymy wkład każdej funkcji bazowej (jeśli dany jest cały wektor wejściowy) w wynik. Jeśli funkcje bazowe miałyby wyprowadzić pojedynczą współrzędną ($ f_i(\mathbf{x}) = x_i$), wrócilibyśmy do regresji liniowej.

Następnie problemem optymalizacyjnym jest dobór wag tak, aby zminimalizować błąd aproksymacji.

$$\displaystyle \min_w \sum_{j=1}^k (y_j – \hat{f}(\mathbf{x}_j))^2$$

Załóżmy na przykład, że chcemy przeprowadzić regresję na podstawie wielomianów kwadratowych. Nasza baza dla trzech zmiennych wejściowych może wyglądać tak

$$\displaystyle \{1, x_1, x_2, x_3, x_1x_2, x_1x_3, x_2x_3, x_1^2, x_2^2, x_3^2 \}$$

Dowolny wielomian kwadratowy trzech zmiennych można zapisać jako kombinację liniową tych funkcji bazowych. Zauważ również, że jeśli potraktujemy to jako podstawę przestrzeni wektorowej, wówczas wektor będzie krotką złożoną z 10 liczb — dziesięciu współczynników wielomianu. Jest to to samo, co $ \mathbb{R}^{10}$, tylko z inną interpretacją znaczenia wpisów wektora. Dzięki temu możemy zobaczyć, jak obliczalibyśmy iloczyny skalarne, rzuty i inne przydatne rzeczy, chociaż mogą one nie mieć tej samej wrażliwości geometrycznej.

Nie są to zwykłe funkcje bazowe używane w praktyce w regresji wielomianowej (patrz uwaga na końcu tego artykułu), ale już pisząc algorytmy regresji możemy wyrządzić pewne szkody.

Proste stochastyczne zejście gradientowe

Chociaż istnieje rozwiązanie w formie zamkniętej wielu problemów regresyjnych (w tym problem regresji kwadratowej, choć z niewielkim skrętem), opadanie gradientowe jest wystarczająco prostym rozwiązaniem, aby pokazać, jak moduł optymalizacyjny może znaleźć użyteczną kombinację liniową. Ten kod zostanie napisany w Pythonie 3.9. Jest na Githubie.

Najpierw zaczniemy od kilku przydatnych aliasów typów

from typing import Callable, Tuple, List

Input = Tuple(float, float, float)
Coefficients = List(float)
Gradient = List(float)
Hypothesis = Callable((Input), float)
Dataset = List(Tuple(Input, float))

Następnie zdefiniuj prostą klasę opakowania dla naszych funkcji podstawowych

class QuadraticBasisPolynomials:
    def __init__(self):
        self.basis_functions = (
            lambda x: 1,
            lambda x: x(0),
            lambda x: x(1),
            lambda x: x(2),
            lambda x: x(0) * x(1),
            lambda x: x(0) * x(2),
            lambda x: x(1) * x(2),
            lambda x: x(0) * x(0),
            lambda x: x(1) * x(1),
            lambda x: x(2) * x(2),
        )

    def __getitem__(self, index):
        return self.basis_functions(index)

    def __len__(self):
        return len(self.basis_functions)

    def linear_combination(self, weights: Coefficients) -> Hypothesis:
        def combined_function(x: Input) -> float:
            return sum(
                w * f(x)
                for (w, f) in zip(weights, self.basis_functions)
            )

        return combined_function

basis = QuadraticBasisPolynomials()

The linear_combination funkcja zwraca funkcję, która oblicza sumę ważoną funkcji bazowych. Teraz możemy zdefiniować błąd na zbiorze danych, a także dla pojedynczego punktu

def total_error(weights: Coefficients, data: Dataset) -> float:
    hypothesis = basis.linear_combination(weights)
    return sum(
        (actual_output - hypothesis(example)) ** 2
        for (example, actual_output) in data
    )


def single_point_error(
        weights: Coefficients, point: Tuple(Input, float)) -> float:
    return point(1) - basis.linear_combination(weights)(point(0))

Następnie możemy zdefiniować gradient funkcji błędu w odniesieniu do wag i pojedynczego punktu danych. Przypomnijmy, że funkcja błędu jest zdefiniowana jako

$$\displaystyle E(\mathbf{w}) = \sum_{j=1}^k (y_j – \hat{f}(\mathbf{x}_j})^2$$

gdzie $ \hat{f}$ jest liniową kombinacją funkcji bazowych

$$\hat{f}(\mathbf{x}_j) = \sum_{s=1}^n w_s f_s(\mathbf{x}_j)$$

Ponieważ będziemy wykonywać stochastyczne opadanie w gradiencie, wzór na błąd jest nieco prostszy. Obliczamy to nie dla całego zbioru danych, ale tylko dla pojedynczego losowego punktu na raz. Więc błąd jest

$$\displaystyle E(\mathbf {w}) = (y_j – \kapelusz {f}(\mathbf {x}_j))^2$$

Następnie obliczamy gradient w odniesieniu do poszczególnych wpisów $ \mathbf{w}$, korzystając z reguły łańcuchowej i zauważając, że jedynym wyrazem kombinacji liniowej, który ma niezerowy wkład w gradient dla $ \frac{\partial E}{\partial w_i}$ jest człon zawierający $ w_i$. Jest to jedna z głównych zalet stosowania kombinacji liniowych: obliczenia gradientu są łatwe.

$ $ \ Displaystyle \ frac {\ częściowe E} {\ częściowe w_i} = -2 (y_j – \ kapelusz {f} (\ mathbf {x} _j)) \ frac {\ częściowe \ kapelusz {f}} {\ częściowe w_i} (\ mathbf {x} _j) = –2 (\ mathbf {x} _y) = –2 \hat{f}(\mathbf{x}_j)) f_i(\mathbf{x}_j)$$

Kolejną zaletą liniowości jest to, że formuła ta jest niezależna od zawartości podstawowych funkcji bazowych. Będzie tak tak długo, jak wagi nie pojawią się we wzorze na funkcje bazowe. Jako ćwiczenie: spróbuj zmienić implementację, aby używać radialnych funkcji bazowych wokół każdego punktu danych. (zobacz uwagę na końcu, dlaczego byłoby to problematyczne w prawdziwym życiu)

def gradient(weights: Coefficients, data_point: Tuple(Input, float)) -> Gradient:
    error = single_point_error(weights, data_point)
    dE_dw = (0) * len(weights)

    for i, w in enumerate(weights):
        dE_dw(i) = -2 * error * basis(i)(data_point(0))

    return dE_dw

Wreszcie rdzeń gradientowy z pomocnikiem debugowania.

import random

def print_debug_info(step, grad_norm, error, progress):
    print(f"{step}, {progress:.4f}, {error:.4f}, {grad_norm:.4f}")


def gradient_descent(
        data: Dataset,
        learning_rate: float,
        tolerance: float,
        training_callback = None,
) -> Hypothesis:
    weights = (random.random() * 2 - 1 for i in range(len(basis)))
    last_error = total_error(weights, data)
    step = 0
    progress = tolerance * 2
    grad_norm = 1.0

    if training_callback:
        training_callback(step, 0.0, last_error, 0.0)

    while abs(progress) > tolerance or grad_norm > tolerance:
        grad = gradient(weights, random.choice(data))
        grad_norm = sum(x**2 for x in grad)

        for i in range(len(weights)):
            weights(i) -= learning_rate * grad(i)

        error = total_error(weights, data)
        progress = error - last_error
        last_error = error
        step += 1

        if training_callback:
            training_callback(step, grad_norm, error, progress)

    return basis.linear_combination(weights)

Następnie utwórz przykładowe dane i uruchom optymalizację

def example_quadratic_data(num_points: int):
    def fn(x, y, z):
        return 2 - 4*x*y + z + z**2

    data = ()
    for i in range(num_points):
        x, y, z = random.random(), random.random(), random.random()
        data.append(((x, y, z), fn(x, y, z)))

    return data


if __name__ == "__main__":
    data = example_quadratic_data(30)
    gradient_descent(
        data,
        learning_rate=0.01,
        tolerance=1e-06,
        training_callback=print_debug_info
    )

W zależności od losowości może to zająć kilka tysięcy kroków, ale zazwyczaj kończy się błędem

Jądra i regularyzacja

Zakończę wyjaśnieniami dotyczącymi powyższych nawiasów.

Prawdziwe jądro wielomianu. Wybraliśmy prosty zbiór funkcji wielomianowych. Jest to ściśle powiązane z koncepcją „jądra”, ale „prawdziwe” jądro wielomianowe wykorzystuje nieco inne funkcje bazowe. Skaluje niektóre funkcje bazowe o $ \sqrt{2}$. Jest to w porządku, ponieważ kombinacja liniowa może to kompensować za pomocą współczynników odpowiednio podzielonych przez $ \sqrt{2}$. Ale dlaczego ktoś miałby chcieć to zrobić? Odpowiedź sprowadza się do techniki wydajności obliczeniowej zwanej „sztuczką jądra”. Krótko mówiąc, umożliwia obliczenie iloczynu skalarnego między dwiema liniowymi kombinacjami wektorów w tej przestrzeni wektorowej bez bezpośredniego przedstawiania wektorów w przestrzeni. Jeśli algorytm regresji używa w swoim kodzie wyłącznie iloczynów skalarnych (co ma miejsce w przypadku rozwiązania regresji w formie zamkniętej), można uzyskać korzyści z modelowania cech nieliniowych bez kosztów bezpośredniego obliczania cech. Jest tu o wiele więcej teorii matematycznej do omówienia (por. Odtwarzanie przestrzeni jądra Hilberta), ale na razie muszę to zostawić.

Co jest nie tak z promieniowym ćwiczeniem funkcji bazowej? W tym ćwiczeniu należało utworzyć rodzinę funkcji bazowych, po jednej dla każdego punktu danych. Problem polega na tym, że posiadanie tak wielu funkcji bazowych sprawia, że ​​przestrzeń kombinacji liniowych jest zbyt wyrazista. Optymalizacja spowoduje nadmierne dopasowanie danych. Przypomina to tabelę przeglądową: każdemu punktowi danych przypisany jest jeden wpis. Nowe punkty danych, które nie są ujęte w szkoleniu, rzadko będą dobrze obsługiwane, ponieważ nie znajdują się w „tabeli przeglądowej” znalezionej przez algorytm optymalizacji. Aby obejść ten problem, w praktyce należałoby dodać dodatkowy człon do błędu odpowiadającego normie L1 lub L2 wektora wag. Pozwala to zapewnić, że całkowity rozmiar wag jest mały i w przypadku L1 zwykle odpowiada to większości wag wynoszących zero, a tylko kilku wagom (najważniejszym) jest niezerowych. Nazywa się proces karania „wielkości” kombinacji liniowej regularyzacja.