# Optymalizacja 2024/25
## Regresja liniowa

Zajmiemy się zastosowaniem optymalizacji liniowej w dziale statystyki nazywanym regresją liniową. Następnie wykorzystamy regresję liniową do oszacowania złożoności metody sympleks.

## Jak mierzyć średnie?

Zaczniemy od przykładu. Poniżej podane są wyniki egzaminu z Optymalizacji Liniowej (w procentach):

$$28, 62, 80, 84, 86, 86, 92, 95, 98$$

Niech $m$ oznacza liczbę ocen ($m = 9$) i niech $b_i$, $i = 1, 2, \ldots, m$ oznacza ocenę (w porządku rosnącym, jak wyżej). Najprostszą średnią jaką możemy policzyć jest *średnia arytmetyczna*

$$
\bar x = \frac 1m \sum_{i=1}^m b_i = 79.0
$$

Jest to przykład *statystyki*.

Czy to dobra średnia? Czy wynik $80$ jest rzeczywiście powyżej średniej?

Drugą (dobrze znaną z gimnazjum) średnią jest *mediana*:

$$
\hat x = b_{(m+1)/2} = 86.
$$

(Na potrzeby tych zajęć, dla uproszczenia, przyjmiemy, że $m$ jest zawsze nieparzyste.)

Rzeczywiście, $86$ wydaje się bardziej "średnie" dla tego zestawu danych.

## Średnia artmetyczna jako zagadnienie optymalizacyjne

Średnia arytmetyczna jest rozwiązaniem problemu minimalizacyjnego (nie jest to problem liniowy)

$$
\bar x = {\arg\min}_{x \in \mathbb{R}} \sum_{i=1}^m (x - b_i)^2
$$

Dlaczego?

$$
f(x) = \sum_{i=1}^m (x - b_i)^2
$$

$$
f'(\bar x) = \sum_{i = 1}^m 2(\bar x-b_i) = 0
$$

$$
m \bar x = \sum_i b_i
$$

$$
\bar x = \frac 1m (\sum_i b_i)
$$

## Mediana jako zagadnienie optymalizacyjne

Mediana jest rozwiązaniem problemu minimalizacyjnego

$$
\hat x = {\arg\min}_{x \in \mathbb{R}} \sum_{i = 1}^m | x - b_i |.
$$

Dlaczego?

$$
f'(x) = \sum_{i=1}^m sgn (x - b_i)
$$

Znak pochodnej zmienia się z ujemnego na dodatni dla środkowej wartości $b_i$ - czyli dla mediany.

<center>

</center>

## Regresja liniowa

Zakładamy, że na "losową" obserwację $b$ składają się nieznana wielkość $x$ i losowa fluktuacja $\varepsilon$:

$$
b = x + \varepsilon.
$$

Przeprowadzamy ciąg obserwacji i indeksujemy wyniki:

$$
b_i = x + \varepsilon_i, i = 1, 2, \ldots, m.
$$

Jak widzieliśmy wyżej, średnia arytmetyczna minimalizuje sumę kwadratów $\varepsilon_i$, mediana minimalizuje sumę modułów $\varepsilon_i$.

Załóżmy, że zebraliśmy dane o ciśnieniu pacjentów i chcemy zbadać, jak zależy ono od wieku, wagi, płci, etc. Tak więc z każdą obserwacją $b$ (ciśnienie) związane są wartości *zmiennych kontrolnych* - $a_1$ to wiek, $a_2$ to BMI, $a_3$ to płeć ($0$ lub $1$), etc. Niech $n$ oznacza liczbę zmiennych kontrolnych, które rozważamy. W regresji liniowej **zakładamy**, że obserwacja $b$ zależy liniowo od zmiennych kontrolnych. Zakładamy więc, że istnieją (nieznane) wartości $x_j$, $j= 1,2,\ldots, n$, dla których

$$
b = \sum_{j=1}^n a_j x_j + \varepsilon.
$$

Oczywiście, dane gromadzimy w tysiącach przypadków, mamy więc ciąg obserwacji, które indeksujemy

$$
b_i = \sum_{j=1}^n a_{ij} x_j + \varepsilon_i.
$$

W zapisie macierzowym to

$$
b = Ax + \varepsilon
$$

W regresji linowej celem jest znalezienie wektora $x$, który najlepiej "tłumaczy" obserwacje $b$. Tak jak dla średniej arytmetycznej lub mediany, możemy zrobić to minimalizując sumę kwadratów lub modułów $\varepsilon_i$. Są też inne możliwości, które przedyskutujemy poniżej.

## $L^2$-regresja

Rozważamy $p$-tą normę wektora $y \in \mathbb{R}^n$:

$$
\| y \|_p = (\sum_i |y_i|^p)^{\frac 1p}.
$$

Przypadek $p = 2$ to odległość euklidesowa. Drugim najważniejszym przypadkiem jest $p = 1$. (Trzecim $p = \infty$.)

Minimalizując $\varepsilon$ w $L^2$-normie otrzymujemy zagadnienie $L^2$-regresji:

$$
\bar x = {\arg\min}_{x} \| b - Ax \|^2_2.
$$

Tak jak dla średniej arytmetycznej możemy wyznaczyć jawny wzór na $\bar x$.

$$
f(x) = \| b - Ax \|^2_2 = \sum_i (b_i - \sum_j a_{ij} x_j)^2
$$

$$
\frac{\partial \bar x}{\partial x_k} = \sum_i 2 (b_i - \sum_j a_{ij} \bar x_j) (-a_{ik}) = 0
$$

Po uproszczeniu mamy dla $k= 1, 2, \ldots, n$:
$$
\sum_i a_{ik} b_i = \sum_i \sum_j a_{ik} a_{ij} \bar x_j
$$

$$
A^T b = A^T A \bar x
$$

$$
\bar x = (A^T A)^{-1} A^T b.
$$

### Przykład


(Jedna zmienna kontrolna $a$, modelem jest $b = ax_1 + x_2$.)

Obserwacje w zapisie macierzowym:

$$
\begin{bmatrix}1\\\frac 52\\ 3\end{bmatrix} =
\begin{bmatrix}0 & 1 \\ 2 & 1 \\ 4 & 1\end{bmatrix}
\begin{bmatrix}x_1 \\ x_2\end{bmatrix} +
\begin{bmatrix}\varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \end{bmatrix}
$$

## $L^1$-regresja

Problem minimalizacyjny

$$
\hat x = {\arg \min}_{x} \| b - Ax \|_1.
$$

Mimo modułów jest to problem liniowy, bo

$$
\min \sum_i | b_i - \sum_j a_{ij} x_j |
$$

można zapisać jako

$\sum_i t_i \to \min$
$t_i - |b_i - \sum_j a_{ij} x_j | = 0$

co jest równoważne problemowi minimalizacyjnemu

$\sum_i t_i \to \min$

$-t_i \leq b_i - \sum_j a_{ij} x_j \leq t_i$ dla $i = 1, 2, \ldots, m$.

### Przykład

(Ten sam przykład, tym razem w przypadku $L^1$.)

(...)

## Jak szybka jest metoda sympleks?

Kryterium optymalności są odpowiednie znaki w tabelce problemu: znaki funkcji celu mają być ujemne, a problem ma być dopuszczalny, tzn. znaki w kolumnie wyrazów wolnych mają być dodatnie.

Tych znaków jest $m + n$. Szacujemy zakładając, że wierzchołek nie jest zdegenerowany, tzn. żadna z tych liczb nie jest równa zero.

Spodziewamy się, że w losowym wierzchołku połowa znaków będzie prawidłowa: $\frac{m+n}2$. Spodziewamy się, że jeden krok metody sympleks poprawi jeden znak, a pozostałe zmieni "losowo", czyli średnio nic nie zmieni.

W tym momencie możemy argumentować, że metoda sympleks będzie średnio wykonywać $\frac{m+n}2$ kroków, by dojść do wierzchołka optymalnego.

Hipoteza:

$$
T = 2^\alpha (m+n)^\beta
$$

Logarytmując stronami i oznaczająć przez $\varepsilon$ odchylenie obserwacji od modelu otrzymujemy

$$
\log T = \alpha \log 2 + \beta \log (m+n) + \varepsilon
$$

Oczywiście wykonujemy ciąg eksperymentów

$$
\begin{bmatrix}
\log T_1 \\ \log T_2 \\ \vdots \\ \log T_k
\end{bmatrix} =
\begin{bmatrix}
\log 2 & \log (m_1 + n_1) \\
\log 2 & \log (m_2 + n_2) \\
\vdots & \vdots \\
\log 2 & \log (m_k + n_k)
\end{bmatrix}
\begin{bmatrix}
\alpha \\ \beta
\end{bmatrix} +
\begin{bmatrix}
\varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_k
\end{bmatrix}
$$

Przy odpowiednich oznaczeniach
$$
b = Ax + \varepsilon
$$

Eksperymentalne dane są następujące


(Z książki Vanderbei'a.)

Rozwiązując problem za pomocą $L^2$-regresji otrzymujemy

$$
\begin{bmatrix}\bar\alpha \\ \bar\beta\end{bmatrix} =
\begin{bmatrix}-1.03561 \\ 1.05152\end{bmatrix}.
$$

Inaczej mówiąc
$$
T \approx 0.488(m+n)^{1.052}.
$$

Dla $L^1$-regresji otrzymujemy
$$
\begin{bmatrix}\hat\alpha \\ \hat\beta\end{bmatrix} =
\begin{bmatrix}-0.9508 \\ 1.0491\end{bmatrix}.
$$

Czyli
$$
T \approx 0.517(m+n)^{1.049}.
$$

<center>


Wykres (skale logarytmiczne) $T$ w zależności od $m+n$ dla $L^1$ i $L^2$ regresji.
</center>