Bedste linje og bedste docking

Bedste linje og bedste docking#

For en given punktmængde i planen ønsker vi at finde den linje, som approksimerer punktmængden bedst muligt. Godheden af approksimationen udtrykkes ved en fejlfunktion, som simpelthen angiver summen af kvadraterne på de vinkelrette afstande fra de enkelte punkter til de linjer, der kandiderer til at være den bedste i den forstand. Approksimationen fås altså IKKE ved den (måske velkendte) lineære regression, hvor den lodrette afstand fra de enkelte punkter til linjen benyttes. Lodret betyder her: parallel med \(y\)-aksen, så lineær regression er meget sensitiv overfor valg af koordinatsystem. Det er approksimationen med de vinkelrette afstande til linjen ikke. Den approksimation er derfor meget mere geometrisk.

Hvorfor afhænger den vinkelrette afstand fra et punkt til en linje ikke af det koordiantsystem, vi benytter til at beregne den afstand?

Opsætning#

Vi starter med at importere de nødvendige kodepakker.

# importer standardbiblioteker
import numpy as np 
import matplotlib.pyplot as plt

# print instillinger for numpy matricer
np.set_printoptions(
    precision=3,       # Number of decimal places
    suppress=True,     # Suppress scientific notation
    linewidth=100,     # Width of output line
    threshold=20,      # Max number of array elements before truncation
    edgeitems=5        # Number of array elements at beginning/end when truncated
)

# importer ekstra funktioner
from ekstra_funktioner import plot_Q_Le

Punktmængder og linjer#

Vi begynder med opgaven fra invitationen til workshoppen: I et \((x, y)\)-koordinatsystem i planen har vi givet et punkt \(p=(x,y)\) og en linje \(L_e\) igennem \((0, 0)\) med enhedsretningsvektor \(\boldsymbol e=(\cos(v), \sin(v))\).

Linjen langs den korteste strækning fra \(p\) til linjen \(L_e\) danner en ret vinkel med \(L_e\). Dermed udspænder skæringspunktet, Origo og \(p\) tilsammen en retvinklet trekant.

Hvad er afstanden \(r\) fra \(p\) til linjen \(L_e\) udtrykt ved \(p\)’s koordinater \((x,y)\) og vinklen \(v\)?

Kan du udtrykke kvadratet \(r^{2}\) som et matrix produkt på følgende måde (se modul \(1\) i Pyodide om matrix-produkter):

(1)#\[\begin{equation} r^{2} = \left[\begin{array}{cc} e_{1}&e_{2} \end{array}\right] \cdot \left[\begin{array}{cc} y^{2} & -xy\\ -yx&x^{2} \end{array}\right] \cdot \left[\begin{array}{c} e_{1}\\ e_{2} \end{array}\right] \end{equation} \]

Vi ønsker at definere en funktion, der kan beregne denne værdi for ethvert givet punktet \(p\) og givet vinkel \(v\). I Python bruger vi def til at skrive en funktion. Her er et eksempel på en funktion gennemsnit, der regner gennemsnittet af tre tal.

Kør cellen for at definere funktionen.

def gennemsnit(a,b,c):
    sum = a + b + c
    resultat = sum / 3
    return resultat

I følgende celle tester vi funktionen.

Kør cellen. Ser du, hvad du forventer?

gns = gennemsnit(-2,6,11)
print("Resultatet er ", gns)

I følgende kodecelle skal vi definere en funktion afstand2, der netop beregner den ønskede kvadrerede afstand \(r^{2}\) fra punktet \(p\) til linjen \(L_e\). Som input har vi \(p\)’s koordinater \(x\) og \(y\) samt vinklen \(v\).

Brug dit resultat fra tidligere til at færdiggøre funktionen.

Vink: Du får måske brug for udtrykkene

\(z^2=\)z**2,

\(\mathrm{cos}(z)=\)np.cos(z),

\(\mathrm{sin}(z)=\)np.sin(z).

def afstand2(x,y,v):
    '''
    Beregner den kvadrerede afstand fra punktet (x,y) til linjen givet 
    ved enhedsvektoren e = (cos(v),sin(v)).
    '''
    resultat = "INDSÆT KODE HER"
    return resultat

I følgende kodecelle kan du teste din funktion. Du kan f.eks. teste det simple tilfælde, hvor \(v=0\), dvs. \(L_e\) ligger på \(x\)-aksen.

Hvad forventer du afstanden fra et punkt \((x,y)\) er i det tilfælde?

Kan du eftervise det i koden?

x = "INDSÆT KODE HER"
y = "INDSÆT KODE HER"
v = "INDSÆT KODE HER"
r2 = afstand2(x,y,v)
print("Afstanden er", r2)

Vi ønsker nu naturligvis at udvide til at kigge på flere punkter på én gang. Vi vil vælge en punktmængde \(Q\) med \(7\) punkter \(q_i = (x_i, y_i), \, i = 1, \cdots, 7\), i planen. Målet er at bestemme summen af kvadrat-afstandene \(r^2\) fra \(q_i\) til \(L_e\), altså udtrykket

\[ S(Q, v) = \sum_{i=1}^{i=7} r^2(v)(q_i). \]

I Python er det praktisk at gemme punkterne i en \(2\times 7\) matrix

\[\begin{split} Q = \begin{bmatrix}x_1&x_2&\dots&x_7\\ y_1&y_2&\dots & y_7\end{bmatrix}, \end{split}\]

så den \(i\)’te søjle i \(Q\) indeholder koordinaterne til punktet \(q_i\). I følgende kodecelle defineres \(Q\) med tre punkter \((4,8)\), \((-3,5)\) og \((0,9)\).

Vælg \(7\) punkter med koordinater mellem \(-10\) og \(10\). Modificer følgende kodecelle så \(Q\) indeholder dine \(7\) punkter.

# Opsætning af data-matricen - INDSÆT EGNE PUNKTER
Q = np.array([[4,8], 
              [-3,5],
              [0,9]]).T

# printer data-matricen
print("Matrix Q:\n", Q)

I følgende kodecelle laver vi et plot af punktmængden i planen.

# Plotter punktmængden Q
plt.scatter(Q[0,:], Q[1,:], color='blue', s=50)

# Tilføjer origo til plottet
plt.scatter(0, 0, color='red', s=100, marker='x', label='Origo')  # Origo

# Opsætning og udseende af plottet
plt.xlabel('x')
plt.ylabel('y')
plt.title('Punkter i planen')
plt.legend()
plt.grid(True)
plt.axis('equal')  # equal scaling
plt.show()

Når vi arbejder med en datamatrix \(Q\) i Python (et numpy-array), tilgår vi elementer ved at skrive Q[i, j].

Her betyder

i rækken, dvs. hvilken vandret linje vi vælger.
j søjlen, dvs. hvilken lodret linje vi vælger.

Python tæller fra 0, så Q[0, 1] betyder “første række, anden søjle”. Vi kan også bruge et kolon : til at sige “tag det hele” i en given retning.

Q[:, j] betyder hele søjle \(j\).
Q[i, :] betyder hele række \(i\).

Det betyder altså, at udtrykket Q[:, i] giver hele punktet nummer i, fordi det henter alle rækker i søjle i. Følgende kodecelle tester disse udtryk.

Færdiggør kodecellen, så sidste linje printer det andet punkt \(q_2\).

# printer data-matricen
print("Matrix Q:\n", Q)

# Et tal i Q: 
print("y1 = ", Q[1,0])

# Alle x-koordinater:
print("x-koordinaterne =", Q[0,:])

# Andet punkt:
print("Punkt 2 =", "INDSÆT KODE HER")

For at kigge på alle punkterne i punktmængden, kan vi bruge en for-løkke. Kodecellen nedenfor printer alle punkterne i \(Q\) vha. en for-løkke. For-løkken går igennem tallene \(j=0,1,...,6\) og for hver printes den \(j\)’te søjle i \(Q\) altså \(q_j\).

# finder antal punkter i Q - vi har N = 7
N = Q.shape[1]  

# for-løkke til at printe alle punkter
for j in range(N): # løkke over tallene i = 0,1,...,N-1
    # finder punkt j
    q = Q[:,j]

    # printer punkt i
    print("Punkt", j+1, ":", q)

Vi vil nu definere en funktion sum_afstand2, der beregner udtrykket \(S(Q,v)\). Vi bruger en for-løkke til at tage ét punkt ad gangen. Vi kan bruge afstand2(x,y,v) til at beregne den kvadrerede afstand for det enkelte punkt og lægge det til resultatet.

Færdiggør funktionen.

def sum_afstand2(Q,v):
    '''
    Beregner summen af de kvadrerede afstande fra punkterne i Q til linjen
    givet ved enhedsvektoren e = (cos(v),sin(v)).
    '''
    # finder antal punkter i Q
    N = Q.shape[1]  

    resultat = 0
    for j in range(N): # løkke over tallene j = 0,1,...,N-1
        # finder punkt j
        q = "INDSÆT KODE HER"

        # koordinaterne for punkt j
        x = q[0]
        y = q[1]

        # beregner den kvadrerede afstand fra punkt j til linjen
        r2 = "INDSÆT KODE HER"

        # lægger den kvadrerede afstand til resultatet
        resultat = resultat + r2

    # når vi er løbet gennem alle punkter kan det endelige resultat returneres
    return resultat

Målet er nu at finde den linje \(L_e\), der bedst approksimerer punktmængden. Dvs. vi ønsker at finde den vinkel \(v\), der giver den mindste værdi af \(S(Q,v)\). I følgende kodecelle genereres et plot over punktmængden \(Q\) og linjen \(L_e\), og \(S(Q,v)\) beregnes.

Eksperimenter med valget af vinkel \(v\) så du synes linjen bedst approksimerer punktmængden. Giver det også en lille værdi af \(S(Q,v)\)?

Husk: Vinklen er i radianer, så vælg \(v\in[0,2\pi]\). Alternativt kan man bruge funktionen v=np.deg2rad(w) til at omregne fra grader til radianer.

# Vælg en vinkel for linjen i radianer
v = np.pi/4 # vælg selv værdien her

# Genererer plot af punkter og linje
plot_Q_Le(Q, v)

print("Sum af kvadrerede afstande:", sum_afstand2(Q, v))

Den bedste linje#

Indtil videre har vi kun ledt efter linjer \(L_e\), der går igennem Origo. For at kunne finde den bedste af alle linjer i planen til at approksimere punktmængden er vi nødt til at se ud over den restriktion.

Det viser sig - ved Steiner’s parallakse sætning - at for parallelle linjer, så er fejlen ved linjeapproksimationen, altså størrelsen \(S(Q,v)\), mindst for linjen der går igennem punktmængdens massemidtpunkt. Af det kan vi konkludere, at den bedste linjeapproksimation nødvendigvis må gå igennem massemidtpunktet. I tilfælde med \(7\) punkter beregnes massemidtpunktet som

\[ CM(Q)=\frac{1}{7}\sum_{i=1}^7 q_i, \]

altså gennemsnittet af punkterne. For at lede igennem alle linjer, der går igennem \(CM(Q)\), kan vi f.eks. dreje linjer med en vinkel som før, men rykke omdrejningspunktet fra Origo til \(CM(Q)\). Alternativt, og helt ækvivalent, kan vi rykke hele punktmængden \(Q\) så massemidtpunktet ligger i Origo. Så kan vi bruge vores metode fra før. Vi kalder den rykkede punktmængde for centreret. Det svarer altså til at kigge på punkterne

\[ q_i-CM(Q),\quad i=1,...,7\;. \]

For at bestemme den centrerede punktmængde skal vi så først beregne massemidtpunktet. Vi kan bruge funktionen gennemsnitsfunktion np.mean til at beregne massemidtpunktet i Python.

# beregner massemidtpunktet CM(Q)
CM = np.mean(Q, axis=1, keepdims=True)

print("Massemidtpunktet CM =\n", CM)

Vi kan nu bruge massemidtpunktet til at beregne den centrerede punktmængde \(Q_c\). Det gøres let i Python ved at trække punktet \(CM(Q)\) fra alle punkt-elementerne i matricen \(Q\).

Qc = Q - CM
print("Centreret punktmængde Qc:\n", Qc)

Den centrerede punktmængde \(Q_c\) er altså dannet ved at flytte hele punktmængden \(Q\) med afstanden mellem massemidtpunktet og Origo. Lad os illustrere dette ved at plotte punktmængde \(Q\) sammen med \(Q_c\). Vi viser også både massemidtpunktet \(CM\) og Origo.

Ser resultatet korrekt ud?

# plotter punktmængden Q
plt.scatter(Q[0,:], Q[1,:], color='blue', s=50, label='Q')

# plotter den centrerede punktmængde Qc
plt.scatter(Qc[0,:], Qc[1,:], color='green', s=50, marker='x', label='Qc')

# plotter massemidtpunktet CM(Q)
plt.scatter(CM[0], CM[1], color='purple', s=100, marker='o', label='Massemidtpunkt') 

# plotter Origo
plt.scatter(0, 0, color='red', s=100, marker='x', label='Origo')  # Origo

# opsætning og udseende af plot
plt.xlabel('x')
plt.ylabel('y')
plt.title('Punkter i planen med massemidtpunkt')
plt.legend()
plt.grid(True)
plt.axis('equal') 
plt.show()

Vi vil nu finde den bedste linje til at approksimere den centrerede punktmængde \(Qc\). Vi kigger igen på linjer, der går igennem Origo (som er massemidtpunktet i det her tilfælde). Ligesom tidligere vælges en vinkel \(v\). Vi genererer så et plot med plot_Q_Le og beregner fejlen \(S(Q_c,v)\) med sum_afstand2.

Eksperimenter med valget af \(v\) for at finde den bedste linjeapproksimation. Kan du få en mindre værdi af fejlen \(S(Q_c,v)\) end for den ikke-centrerede punktmængde.

# Vælg en vinkel for linjen i radianer
v = np.pi/4 # vælg selv værdien her

# Genererer plot af punkter og linje
plot_Q_Le(Qc, v)

print("Sum af kvadrerede afstande:", sum_afstand2(Qc, v))

Den bedste linje med SVD#

Data-matricen \(Q_c\) kan (som alle andre matricer) SVD dekomponeres (her udtrykt for en \((2 \times 7)\)-data-matrix) - se modul \(2\) om SVD dekomponering:

\[ Q_c = U \cdot \Sigma \cdot V^\top , \tag{2.11} \]

hvor så \(U\) er en \((2 \times 2)\)-matrix, \(\Sigma\) er en \((2 \times 7)\)-‘diagonalmatrix’, og \(V^\top\) er en \((7 \times 7)\)-matrix med følgende udseende for passende værdier af \(\theta\) og \(\sigma_1 \geq \sigma_2 \geq 0\):

\[\begin{split} U = \begin{bmatrix} \vert&\vert\\ u_1 & u_2 \\ \vert&\vert \end{bmatrix} = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} \end{split}\]

\[\begin{split} V^\top = \begin{bmatrix} -\, v_1 \, - \\ -\, v_2 \, - \\ -\, v_3 \, - \\ -\, v_4 \, - \\ -\, v_5 \, - \\ -\, v_6 \, - \\ -\, v_7 \, - \end{bmatrix} \end{split}\]

\[\begin{split} \Sigma = \begin{bmatrix} \sigma_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \sigma_2 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \end{split}\]

Søjlevektorerne \(u_k\), \(k = 1, 2\), i \(U\) er udtrykt simpelt ved den ene vinkel \(\theta\). Søjlevektorerne er enhedsvektorer og de står vinkelret på hinanden (prikproduktet er 0). Rækkevektorerne \(v_\ell\), \(\ell = 1, \cdots, 7\), i matricen \(V^\top\) er også enhedsvektorer (i \(\mathbb{R}^7\) med 7 koordinater i hver vektor) som er parvis vinkelrette på hinanden (med hensyn til det naturlige prikprodukt i \(\mathbb{R}^7\)).

I Python kan vi beregne SVD af med matrix med funktionen np.linalg.svd.

U, sigma, VT = np.linalg.svd(Qc)
print("U-matrix:\n", U)
print("Singulær-værdier:\n", sigma)
print("V^T-matrix:\n", VT)

Den første søjle i \(U\) kan vi skrive som

\[\begin{split} u_1=\begin{bmatrix} \cos(\theta)\\\sin(\theta)\end{bmatrix} \end{split}\]

for en vinkel \(\theta\in[0,2\pi[\). Linjen \(L_e\) for \(e=u_1\) er den bedste linjeapproksimation for punktmængden \(Q_c\). Vi ønsker altså at finde den tilsvarende vinkel \(\theta\). Fra trigonometri har vi, at

\[ \tan(\theta) = \frac{\sin(\theta)}{\cos(\theta)} \]

og

\[ \theta = \tan^{-1}\left(\frac{\sin(\theta)}{\cos(\theta)}\right)=\tan^{-1}\left(\frac{y}{x}\right), \]

hvis \(x\) og \(y\) er koordinaterne for \(u_1\). I Python kan vi beregne \(\tan^{-1}\) med np.arctan().

Brug np.arctan() til at bestemme vinklen \(\theta\).

theta = np.arctan("INDSÆT KODE HER")

Vi kan nu prøve at bruge linjeapproksimationen givet ved vinklen \(\theta\).

Ser approksimationen god ud?

Får du en bedre værdi af fejlen end du fandt selv tidligere?

# Genererer plot af punkter og linje
plot_Q_Le(Qc, theta)

print("Sum af kvadrerede afstande:", sum_afstand2(Qc, theta))