\input{preamble}

\title{4 - 6 - OM DØDSLENGSEL}

\begin{document}
\pagenumbering{gobble}

%primtall

\maketitle

\begin{minipage}{125mm}
Hvis du skal bygge kai og har bestilt sekstoms bord,
får du såkalte \textit{råbord}.
De er bra trykkimpregnert, 
men ikke saga så nøye, 
så de fleste av dem er ikke nøyaktig 152.4 mm,
men kanskje et eller annet sted mellom 152mm og 157mm.
Jeg fiksa kaien til svigers en sommer,
og målte den faktiske bredden
på alle de tjueni sekstomsbordene som var kjøpt inn.
Gjennomsnittsbredden var $\hat \mu = 154.6$ mm,
og det empiriske standardavviket på $\hat \sigma = 1.1$ mm.
Til høyre er et plot av normalfordelingsfunksjonen
    oppå relative frekvenser,
    altså antall bord av hver lengde delt på totalt antall bord.
    Du finner dette som .csv-fil her:\\
\url{https://folk.ntnu.no/mortano/python/kai}.
\end{minipage}
\begin{minipage}{60mm}
\begin{center}
\includegraphics[scale=.45]{../figurer/kai-3}
\end{center}
\end{minipage}



\begin{minipage}{150mm}
    Nå skal vi komplisere det litt.
    Forestill deg at du har et måleapparat
    og at du måler en normalfordelt størrelse,
    men så vet du at den ene halen i normalfordelingen
    er under deteksjonsgrensen for måleapparatet.
    Hvis du ønsker å estimere $\mu$ og $\sigma$,
    kan du ikke ta gjennomsnitt og empirisk standardavvik,
    for gjennomsnittet til treffe for høyt
    og det empiriske standardavviket for lavt.
    I datasettet over bredden på kaibordene 
    har jeg tatt ut de laveste målingene
    (tre bord på 153 mm og ett på 152 mm).
    Da får man datasettet i margen.
    Hvilken normalfordelingskurve 
    \[
        f(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\left( x-\mu \right)^2/2\sigma^2}
    \]
    kommer dette fra?
    I denne økten skal vi prøve å finne det ut.
    
    I det venstre plottet under
    finner du histogrammet til kaibordene og 
    en normalfordelingskurve der $\mu$ og $\sigma$ er 
    estimert ved gjennomsnitt og empirisk standardavvik 
    på data uten 152- og 153-målingene.
    Det er kanskje ikke så lett å se med 
    det blotte øye at kurven er feil, 
    så jeg har for illustrasjonens skyld 
    tatt ut 154mm-målingene også og gjentatt eksperimentet;
    dette er plottet til høyre.
    I dette plottet er det ganske lett å se at noe er galt.\\[5mm]
    \begin{minipage}{70mm}
        \begin{center}
            \includegraphics[scale=.5]{../figurer/kai-7}
        \end{center}
    \end{minipage}
    \begin{minipage}{70mm}
        \begin{center}
            \includegraphics[scale=.5]{../figurer/kai-9}
        \end{center}
    \end{minipage}
    \end{minipage}
\begin{minipage}{30mm}
        \begin{center}
        \includegraphics[scale=.6]{../figurer/kai-4}
        \end{center}
\end{minipage}    

\clearpage

Problemet på forrige side er litt komplisert, 
så la oss begynne med noe mye enklere - 
Newtons metode i to dimensjoner.
La oss si at vi ønsker å løse det ikkelineære likningssettet
\begin{align*}
    \V f(\V x)= 
    \left(\; 
        \begin{matrix}
            f_1(\mathbf{\V x})\\[1mm]     
            f_2(\mathbf{\V x})
        \end{matrix}
     \;\right)
    =\V 0 
\end{align*}
der $\V f$ er en funksjon fra 
$\mathbb{R}^2$ til $\mathbb{R}^2$.
Dette er to likninger for nullnivåkurvene til komponentene til $\V f$, 
som begge er funksjoner fra 
$\mathbb{R}^2$ til $\mathbb{R}$.
En løsning $\V r$ av likningssettet må altså ligge på 
skjæringspunktet mellom disse nivåkurvene. 
Her er en figur som illustrerer dette med to sirkulære nivåkurver:
 \begin{center}
    \begin{tikzpicture}[thick]
            \draw[blue] (-2,0) circle (2cm) node[right,xshift=1cm,yshift=2cm]{$f_1(\V x)=0$};
            \draw[red] (0,-1) circle (1.5cm) node[right,xshift=1.5cm,yshift=0cm]{$f_2(\V x)=0$};
            \node[below right] at (-.3,.7){$\circ$};
            \node[below right] at (0,1.2){$\V r$};
            \draw[->] (-3.5,0) -- (2.5,0) node[right]{$x_1$};
            \draw[->] (-2.5,-1) -- (-2.5,1.3) node[above]{$x_2$};
      \end{tikzpicture}
  \end{center}
Numeriske metoder er stort sett iterative. 
Vi lar $\V x_n$ være iterasjon nummer $n$. 

\begin{oppgave}{1}
    Sett opp tangentplanene til $f_1$ og $f_2$ i iterasjon $\V x_n$, 
    og utled \textbf{Newtons metode}
    \begin{align*}
        \V x_{n+1}
        = 
        \V x_{n}
        - 
        \V f'(\V x_n)^{-1}
        \V f(\V x_n)
    \end{align*}
    (Hint: I bunn og grunn samme som den endimensjonale varianten vi hadde i TMA4101.)
\end{oppgave}

\begin{oppgave}{2}
    Test ved å finne skjæringspunktene mellom de to ellipsene 
    \begin{align*}
        3 x^2 + 2 x y + y^2 - 2x - y = 1
        \hspace*{10mm}
        \text{og}
        \hspace*{10mm}
        x^2 + 2x y + 6 y^2  - x  - 2y = 1
    \end{align*}
    Hint: Det går greit å bruke np.linalg.solve, 
    men for $2\times2$-matriser finnes det en formel:
    \begin{equation*}
        \begin{pmatrix} a & b \\ c & d \end{pmatrix}^{-1}=
        \frac{1}{ad-bc}\begin{pmatrix} d & -b \\ -c & a \end{pmatrix}.
    \end{equation*}        
\end{oppgave}

Apropos - \textbf{kofaktorekspansjon}.
Det finnes en generell formel for inversen til en inverterbar matrise:\\
\url{https://en.wikipedia.org/wiki/Invertible_matrix#Analytic_solution}\\
Det er en og annen gubbe (inkludert meg) som sporadisk synes dette er nyttig å vite om. 
\footnote{Jeg har hatt bruk for det én gang i livet - i den eneste artikkelen jeg har ført i pennen noen gang. 
Den var kjempebra og biveilederen min fikk den publisert i en obskur bulgarsk konferansejournal 
rett før han ble pensjonist.}




\begin{oppgave}{3}
    Forklar hvorfor det å løse oppgave 2 er det 
    samme som å lete etter kritiske punkter til funksjonen 
    \begin{align*}
        V(x,y)=x^3 + x^2y + xy^2 + 2y^3 - x^2 - y^2 - xy- x -y.
    \end{align*}
    (Hint: Se på gradienten til $V$.)
\end{oppgave}

\clearpage

En kjent og kjær teknikk for å lete etter kritiske punkter 
heter \textit{også} \textbf{Newtons metode}. 

\begin{oppgave}{4}
    Forklar hvorfor denne er gitt ved iterasjonen
    \begin{align*}
        \V x_{n+1}
        = 
        \V x_{n}
        - 
        V''(\V x_n)^{-1}
        V' (\V x_n)
    \end{align*}
    der $V: \mathbb{R}^2 \to \mathbb{R}$ er funksjonen vi leter etter kritiske punkter til. \\
    (Husk at $V'$ er gradienten og $V''$ er hessematrisen.)
\end{oppgave}


Antagelig oppdaget du at dødslengselmetoden konvergerer ganske sakte, 
til tross for at den prøver å gå rett utfor så bratt som mulig. 
Dette skyldes at den bratteste retningen ikke nødvendigvis
peker mot det kritiske punktet. 
La oss se litt på dette og samtidig sette lineæralgebraen i nytt lys. 

\begin{oppgave}{8}    
    Forklar at for alle lineære systemer $A\V {x} = \V b$ du har løst til nå,
    er $\V x$ minimumspunktet til 
    \begin{align*}
        V(\V x) =
        (A\V x - \V b)^T
        (A\V x - \V b)
        =
        \|A\V x - \V b\|^2
    \end{align*}
\end{oppgave}

I praktiske anvendelser må man ofte løse enorme lineære systemer
der de fleste komponentene i matrisen er null. 
Da er det sjelden effektivt å gausseliminere. 
Når du bruker np.linalg.solve, 
spør python bare LAPACK\footnote{\url{https://en.wikipedia.org/wiki/LAPACK}}
og der finnes det et tonn med forskjellige metoder.
Vektoren $A\V x - \V b$ kalles \textbf{residualet},
og noen av de mest kjente iterative teknikkene for lineære systemer 
(for eksempel GMRES\footnote{\url{https://en.wikipedia.org/wiki/Generalized_minimal_residual_method}})
baserer seg på å minimere denne i forskjellige underrom av $\mathbb R^n$. 

\begin{oppgave}{9}    
    I anvendelser er matrsien $A$ ofte symmetrisk og positivt definitt. 
    Da kan vi minimere en enklere funksjon. 
    Forklar at dersom $A \in \mathbb{R}^n$ er symmetrisk og positivt definitt,
    er det å løse $A\V x = \V b$ er det samme som å minimere funksjonen 
    \begin{align*}
        V(\V x) = \frac{1}{2}\V x^TA\V x-\V x^T \V b. 
    \end{align*}
\end{oppgave}

Kvadratiske polynomer er praktiske å jobbe med i optimering, 
for de er lette å analysere;
mange klassiske linjesøksmetoder 
baserer seg på å minimere andre ordens taylorpolynomer
til den funksjonen som egentlig skal minimeres. 
Hvis $V$ er kvadratisk slipper for eksempel backtracking.

\begin{oppgave}{10}
    Du kjører et steg med en linjesøksmetode
    på problemet over,
    og har valgt ut retningen $\V p_n$. 
    Finn minimumspunktet til 
    \begin{align*}
        V(\V x_n + \alpha_n \V p_n)
    \end{align*}
    som funksjon av $\alpha_n$.
\end{oppgave}


Nå kan vi se litt på hvorfor dødslengsel konvergerer sakte. 
Dersom $A$ er positivt definitt, 
er nivåkurvene til $f$ i oppgave 9 alltid ellipser.

\begin{oppgave}{11}
    Tegn opp noen litt lange og flate elliptiske nivåkurver og 
    tegn noen dødslengseliterasjoner 
    omtrentlig for hånd med $\alpha_n$ bestemt slik som som i oppgave 10. 
\end{oppgave}




\clearpage

En av de mest kjente metodene for dette kalles \textbf{konjugerte gradienter},
og denne er så bra at oppfinnerne Hestenes og Stiefel ble kjendiser. 
Den konvergerer som et uvær i mange tilfeller,
men som de fleste numeriske metoder vil den også feile spektakulært i andre. 
Hvis du går MTKJ og treffer på Svein Sunde kan du spørre ham om denne metoden, 
han har brukt den i atomkraftverk. 

\begin{oppgave}{13}
    Vis at 
    \begin{align*}
        (\V x,\V y)_A
        =
        \V x^T A \V y
    \end{align*}
    er et indreprodukt dersom $A$ er symmetrisk og positivt definitt. 
\end{oppgave}

Dersom $\V x$ og $\V y$ er ortogonale med hensyn på indreproduktet over,
sier vi at de er $A$-ortogonale. 
Et velkjent resultat fra abstrakt lineæralgebra er
at om du har et vektorrom og et indreprodukt,
finnes det alltid en ortogonal basis. 

\begin{oppgave}{14}
    Anta at du har en $A$-ortogonal basis $\{\V p_k\}$ for $\mathbb{R}^n$, 
    at $A\V x = \V b$, 
    og at 
    \begin{align*}
        \V x
        =
        \sum_k
        \alpha_k \V p_k.
    \end{align*}
    Utled at 
    \begin{align*}
        \alpha_k = 
        \frac{(\V b, \V p_k)}{(\V p_k,\V p_k)_A}
    \end{align*}
\end{oppgave}


Konjugerte gradienters metode består i å konstruere et sett 
med søksretninger som er $A$-ortogonale.
Det viser seg at man på magisk vis får dette til ved rekursjonen
\begin{align*}
    \V p_{k+1}
    =
    -\V r_{k+1}+
    \beta_{k+1} \V p_k
\end{align*}
der 
\begin{align*}
\V r_{k} = A\V x_k - \V b 
\end{align*}
er residualvektoren etter iterasjon $k$, 
og $\beta_k$ velges slik at $\V p_k$ og $\V p_{k+1}$ er 
$A$-ortogonale.

\begin{oppgave}{15}
    Vis at 
    \begin{align*}
        \beta_{k+1}
        = 
        \frac{(\V r_{k+1}, \V p_k)_A}{(\V p_k,\V p_k)_A}.
    \end{align*}
\end{oppgave}

Steglengden velges slik som i oppgave 10 for at det skal funke.
Du finner full algoritme på s 108 i Nocedal og Wright.

\begin{oppgave}{16}
    Vis at $\alpha_k$ over minimerer $V$ slik som i oppgave 10. 
\end{oppgave}



\begin{oppgave}{17}
    Prøv konjugerte gradienters metode på implisitt skjema for varmelikningen. 
\end{oppgave}



\clearpage

\section*{UKENS NØTTER}

\begin{oppgave}{1}
    Estimér $\mu$ og $\sigma$ i kaieksemplet på side 1. 
    Dette kan gjøres på mange måter. 
\end{oppgave}

\begin{oppgave}{2}    
    Har det noen hensikt å kjære Newtons metode på et lineært system?
\end{oppgave}

Dersom $A$ er symmetrisk og positivt definitt, 
finnes det en mye mer fornuftig retning å kjøre enn brattest utfor, 
nemlig langs med halvaksene til de elliptiske nivåkurvene. 

\begin{oppgave}{3}
    Forklar at med strategien fra oppgave 10 trenger du egentlig bare 
    to iterasjoner for å komme frem til minimumspunktet dersom 
    du velger $\V p_n$ slik. 
\end{oppgave}


\begin{comment}


varmelikningsmatrisen

conjugate gradients - definitte matriser 

fft

pseudoinvers - normalikningene

schurfaktorisering

qr - dette er kun gramschmidt

lu

hølder cauchy schwarz


\end{comment}



\end{document}
