3  Matriser og matrise-vektor operasjoner

3.1 Matriser

En matrise er en rektangulær tabell av tall, ordnet i rader og kolonner

\[ A = \begin{bmatrix} 0 & -1.2 & 3 & 0.5 \\ 2.7 & 0 & -0.9 & -2 \\ -3 & 1.1 & 0 & 3.4 \end{bmatrix} \]

Dette er en såkalt \(3\times 4\) matrise (3 rader, 4 kolonner). Elementene i matriser er tallene, og disse skrives som \(a_{ij}\) (evt. \(a_{i,j}\)). For matrisen over \(a_{11} = 0\), \(a_{21} = 2.7\) og \(a_{34} =3.4\). Her er matriselementene vist eksplisitt:

\[ A = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \end{bmatrix} \]

Boyd–Vandenberghe bruker \(A_{ij}\)​ for elementer. Vi bruker store bokstaver for matriser og små bokstaver med indekser for elementer for å unngå forveksling med vektorer.

I kode: Julia

A = [
    0.0   -1.2   3.0   0.5;
    2.7    0.0  -0.9  -2.0;
   -3.0    1.1   0.0   3.4
]

Python

import numpy as np

A = np.array([
    [ 0.0, -1.2,  3.0,  0.5],
    [ 2.7,  0.0, -0.9, -2.0],
    [-3.0,  1.1,  0.0,  3.4]
])

Merk at Julia bruker semikolon ; for å indikere nye rader, mens Python bruker lister av lister \([[..],[..]]\) for å gjøre det samme.

Akkurat som med vektorer starter man å telle på 1 i Julia, mens 0 i Python, slik at i Julia

A[1,1]    # output: 0.0

mens dette må skrives A[0,0] i Python.

A[0, 0]    # output: 0.0

Kvadratisk matrise: En \(m\times m\) matrise kalles en kvadratisk matrise.

Høy matrise: En høy matriser har flere rader enn kolonner \(m>n\).

Bred matrise: En bred matrise har flere kolonner enn rader. Eksempelt \(A\) over er en bred matrise siden \(m=4\) og \(n=3\).

Vektorer som matriser En (kolonne) n-vektor kan forstås som en \(n\times 1\) matrise, mens en radvektor av samme dimensjon kan forstås som en \(1\times n\) matrise.

3.2 Submatriser og blokkmatriser

Submatriser: Om vi har en matrise \(A\), så er et utsnitt av den en submatrise. For matrisen \(A\) definert over, så kan vi skrive dette som:

\[ A_{1:2,1:3} = \begin{bmatrix} 0 & -1.2 & 3 \\ 2.7 & 0 & -0.9\\ \end{bmatrix} \] I kode:

Julia:

A[1:2, 1:3]

Python:

A[0:2, 0:3]

I Python må man starte på 0 og sluttindeksen inkluderes ikke.

Matrise på blokkform En effektiv måte å skrive og analysere store matiser er å bruke blokkform:

For eksempel om vi har den kvadratiske \(4\times 4\) matrisen,

\[ M = \begin{bmatrix} 1 & 2 & 0 & -1 \\ 3 & 4 & -2 & 0 \\ 5 & 6 & 9 & 0 \\ 7 & 8 & 0 & 9 \end{bmatrix} \]

så kan den skrives på blokkform (eller partisjonert form) som føler:

\[ M = \begin{bmatrix} A & B \\ C & D \end{bmatrix}, \]

hvor da:

\[ A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}, \; B = \begin{bmatrix} 0 & -1 \\ -2 & 0 \end{bmatrix}, \; C = \begin{bmatrix} 5 & 6 \\ 7 & 8 \end{bmatrix}, \; D = \begin{bmatrix} 9 & 0 \\ 0 & 9 \end{bmatrix} \]

\(A\), \(B\) osv. kan forstås som submateriser eller undermatriser av den fulle matrisen \(M\).

I kode: Julia:

A = [1 2;
     3 4]
B = [ 0 -1;
     -2  0]
C = [5 6;
     7 8]
D = [9 0;
     0 9]

M = [A B;
     C D]

Python (viser kun siste steg) bruker man np.block

M = np.block([[A, B],
              [C, D]])

Om man starter med en stor matrise \(M\), så kan man selv velge hvordan man vil partisjonere den så lenge blokkene passer inn i den store matrisen.

Matriser på kolonneform: Ofte er det nyttig å skrive matriser på kolonneform (altså at de består av en rad av kolonnevektorer)

\[ M = \begin{bmatrix} a& b & c & \ldots \end{bmatrix} \]

Der a, b og c hver for seg er \(n\)-vektorer. For eksempel, om

\[ M = \begin{bmatrix} 1& 2 & 1\\ 2& 2 & 1 \end{bmatrix} \]

så kan det skrives på denne formen, med \[a = \begin{bmatrix} 1 \\ 2\end{bmatrix},\qquad b = \begin{bmatrix} 2 \\ 2\end{bmatrix},\qquad c = \begin{bmatrix} 1 \\ 1\end{bmatrix}\].

Matriser på radform: Det er på tilsvarende vis mulig å skrive en matrise på radform, her for en \(3\times m\) matrise hvor \(b_1\), \(b_2\) , og \(b_3\) er radvektorer,

\[ M = \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} \]

For å unngå forvirring, så kan det være lurt å uttrykke matriser på radform på følgende vis:

\[ M = \begin{bmatrix} a_1^\top \\ a_2^\top \\ a_3^\top \end{bmatrix} \]

Da kan man forstå vektoeren \(a_1\), \(a_2\), \(a_3\) etc. hver for seg som en (vanlig) kolonnevektor, mens transponeringen gjør dem til radvektorer.

3.3 Grunnleggende matriseoperasjoner

Akkurat som vektorer kan matriser skaleres med et tall. Hvis vi multipliserer en matrise med en skalar \(c\), ganges hvert element med \(c\). For en \(3\times 4\)-matrise \(A=[a_{ij}]\):

\[ cA = \begin{bmatrix} c a_{11} & c a_{12} & c a_{13} & c a_{14} \\ c a_{21} & c a_{22} & c a_{23} & c a_{24} \\ c a_{31} & c a_{32} & c a_{33} & c a_{34} \end{bmatrix}. \]

Summen \(A+B\) er definert når \(A,B\in\mathbb{R}^{m\times n}\) (samme dimensjon), og beregnes elementvis: \[ (a+b)_{ij} = a_{ij} + b_{ij}. \]

Eksempel: \[ A = \begin{bmatrix} -1 & 0 \\ 2 & 1 \end{bmatrix}, \quad B = \begin{bmatrix} 1 & -2 \\ 2 & 2 \end{bmatrix} \;\Rightarrow\; A+B = \begin{bmatrix} 0 & -2 \\ 4 & 3 \end{bmatrix}. \]

Transponat
Transponatet bytter rader og kolonner: \((A^\top)_{ij}=a_{ji}\).

Julia
AT = transpose(A)

ellers så kan man adjugere (som er ekvivalent )

AT = A'

som er det samme transponering for reelle matriser

Python

AT = A.T

Egenskaper: \((A+B)^\top=A^\top+B^\top\), \((cA)^\top=cA^\top\), \((A^\top)^\top=A\).

NB (syntaks og semantikk): I Julia betyr * matriseprodukt, .* betyr elementvis produkt. I Python/NumPy betyr @ matriseprodukt, mens * er elementvis produkt.
Julia: C = A .* B (elementvis), C = A * B (matriseprodukt)
Python: C = A * B (elementvis), C = A @ B (matriseprodukt)


3.4 Egenskaper for addisjon og skalarmultiplikasjon

La \(A\), \(B\) og \(C\) være matriser av samme dimensjon \(m \times n\), og la \(c\) være en skalar. Da gjelder følgende egenskaper:

  • Kommutativitet: \[ A + B = B + A \]

  • Assosiciavitet: \[ (A + B) + C = A + (B + C) \]

  • Nøytralt element (nullmatrise):
    Nullmatrisen \(0\) er den unike matrisen der alle elementer er null, og den oppfyller: \[ A + 0 = A \]

  • Additiv invers:
    For hver matrise \(A\) finnes en additiv invers \(-A\) slik at: \[ A + (-A) = 0 \]

  • Distribusjon av skalarmultiplikasjon over addisjon:
    \[ c(A + B) = cA + cB \]

  • Transponat av sum:
    Transponatet av en sum er summen av transponatene: \[ (A + B)^\top = A^\top + B^\top \]

3.5 Matrisenorm

En matrisenorm gir et tall som måler “størrelsen” til en matrise. Den mest brukte i dette kurset er Frobenius-normen: \[ \|A\|_F = \sqrt{\sum_{i=1}^m \sum_{j=1}^n |A_{ij}|^2 }. \] Dette tilsvarer Euklidsk norm av matrisen “flatet ut” til én lang vektor.

Eksempel: \[ A = \begin{bmatrix} 2 & -1 \\ 0 & 3 \end{bmatrix} \;\Rightarrow\; \|A\|_F=\sqrt{2^2+(-1)^2+0^2+3^2}=\sqrt{14}. \]

Siden \(A\) og \(A^\top\) har samme elementer, gjelder \(\|A\|_F=\|A^\top\|_F\).

Kolonneform-tolkning
Skriv \(A=[a_1\; a_2\;\cdots\; a_n]\) med \(a_j\in\mathbb{R}^m\). Da \[ \|A\|_F=\sqrt{\|a_1\|_2^2+\|a_2\|_2^2+\cdots+\|a_n\|_2^2}. \] Tilsvarende gjelder for radvektorer ved å bruke \(A^\top\) (som har samme norm).

Frobenius-normen er ikke det samme som matrise-2-normen.

Kode (Frobenius-norm) Julia:

using LinearAlgebra  
A = [2 -1; 0 3]  
norm_A = norm(A)    # Frobenius-norm

Python:

import numpy as np  
A = np.array([[2, -1], [0, 3]])  
norm_A = np.linalg.norm(A, 'fro')   # Frobenius-norm

3.6 Matrise-vektor multiplikasjon

Matrise-vektor multiplikasjon er helt sentralt for alle former for modellering. Om vi har en \(m\times n\) matrise

\[ A = \begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1n} \\ \vdots & & & \vdots\\ a_{m1} & a_{m2} & \ldots & a_{mn} \end{bmatrix} \]

og en vektor

\[ x = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n} \end{bmatrix} \]

så er produktet \(y = Ax\) gitt ved en ny vektor i \(\mathbb{R}^m\), der hver komponent er en lineærkombinasjon av elementene i \(x\) og radene i \(A\):

\[ Ax = \begin{bmatrix} a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n \\ a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n \\ \vdots \\ a_{m1} x_1 + a_{m2} x_2 + \ldots + a_{mn} x_n \end{bmatrix} \]

Om vi skriver matrisen på kolonneform:

\[ A = \begin{bmatrix} a_1 & a_2 & \ldots & a_n \end{bmatrix} \]

hvor \(a_1\) etc. er kolonnene med lengde \(m\), så kan vektoren \(y\) skrives som en lineær kombinasjon av kolonnevektorene \(a_j\), akkurat som ved et indreprodukt mellom to vektorer:

\[ Ax = x_1 a_1 + x_2 a_2 + \ldots + x_n a_n \]

Om vi derimot skiver mastrisen på radform

\[ A = \begin{bmatrix} b_1^\top \\ b_2^\top \\ \vdots \\ b_m^\top \end{bmatrix} \] så kan matrise-vektor-produktet skrives som

\[ A x= \begin{bmatrix} b_1^\top \\ b_2^\top \\ \vdots \\ b_m^\top \end{bmatrix} x = \begin{bmatrix} b_1^\top x \\ b_2^\top x \\ \vdots \\ b_m^\top x \end{bmatrix} \]

Dette samsvarer med en vanlig måten vi bruker å visualisere matrise-vektor-produkt

3.7 Sentrale matriser og matrise-vektor-produkter

3.7.1 Indreprodukt som matriseoperasjon

Vi har allerede i praksis sett mange eksempler på matrise-vektor-produkt og transponering: Et indreprodukt mellom to vektorer kan også ses som et spesialtilfelle av matrise-vektor-produkt, der vi transponerer én av vektorene:

\[ x\cdot a =x^\top a = \underbrace{\begin{bmatrix} x_1 & x_2 & \ldots & x_n \end{bmatrix}}_{1\times n \,\, \text{matrix}} \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix} = x_1 a_1 + x_2 a_2 + \ldots + x_n a_n \]

Her ser vi at transponering brukes for å gjøre en kolonnevektor om til en radvektor, slik at multiplikasjonen \(x^\top a\) gir et tall.

3.7.2 Plukke ut kolonnevektor av matrise

Om vi skiver en matrise på kolonneform,

\[A = \begin{bmatrix} a_1 & a_2 & \ldots & a_n \end{bmatrix} \]

da kan bruke en standard enhetsvektor \(e_i \in \mathbb{R}^n\) til å plukke ut kolonnen, her illustrert for \(e_2\) :

\[ A e_2 = \begin{bmatrix} a_1 & a_2 & \ldots & a_n \end{bmatrix} e_2 = \begin{bmatrix} a_1 & a_2 & \ldots & a_n \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} = a_2 \]

Vi kan tilsvarende plukke ut en rad om vi transponerer, slik at \((A^\top e_i)^\top\) svarer til rad \(i\) av matrisen \(A\).

3.7.3 Summasjon av kolonnevektorer

Om en matrise blir multiplisert men en en vektor av bare 1-tall (\(\mathbf{1}_n\)) svarer det til en summasjon av alle kolonnevektorene i matrisen

\[ A \mathbf{1}_n = \begin{bmatrix} a_1 & a_2 & \ldots & a_n \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} = \sum_{i=1}^n \,a_i \]

3.7.4 Null og en-matrise

Null matrisen er \(m\times n\) matrise, hvor alle emenetene er 0. Den skrives noen ganger som \(0_{m\times n}\) . Den gjør \(0_{m\times n}x = 0_m\). Den gjør en \(n\)-vektor til en 0-vektor med dimensjon \(m\).

En en-matrise hvor alle elementene er 1 kan noen ganger skrives som \(1_{m\times n}\).

3.7.5 Identitetsmatrisen.

En veldig viktig matrise er en Identitetsmatrisen. Den skrives ofte slik \(I\). Identitetsmatrisen \(I \in \mathbb{R}^{n \times n}\) er definert ved:

\[ I_{ij} = \begin{cases} 1 & \text{hvis } i = j \\\\ 0 & \text{ellers} \end{cases} \]

Det betyr at alle elementer på diagonalen er 1, og resten er 0. For eksempel, når \(n = 4\):

\[ I = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

Multiplikasjon med identitetsmatrisen bevarer en vektor:

\[ I x = x \]

I fysikk bruker man ofte den såkalte Kronecker-deltaen, som skrives \(\delta_{ij}\). Denne er definert som:

\[ \delta_{ij} = \begin{cases} 1 & \text{hvis } i = j \\\\ 0 & \text{ellers} \end{cases} \quad \Rightarrow \quad I_{ij} = \delta_{ij} \]

Den er derfor det samme som en identitetsmatrise, men ofte er ikke størrelsen på \(n\) spesifisert.

Kolonnene i identitetsmatrisen svarer til standard enhetsvektorene

\[ I = \begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} \quad \]

slik at for eksempel \[ e_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} \]

På indeksform kan dette skrives

\[ (e_2)_i = I_{i2} \]

I Python kan man gjøre dette slik:

import numpy as np

# Lag en 4x4 identitetsmatrise
I = np.eye(4)

# Pakk ut kolonnene direkte
e1, e2, e3, e4 = I.T

print("e2 =")
print(e2)

I Julia

using LinearAlgebra
# Lag en 4x4 identitetsmatrise
Id = Matrix(I, 4, 4)

# Pakk ut kolonnene direkte som enhetsvektorer
e1, e2, e3, e4 = eachcol(Id)

println("e2 = ")
display(e2)

Et praktisk aspekt med Julia er at I alene, virker som en identietsmatrise.

I +  [1 2; 2 3]
# 2×2 Matrix{Int64}:
#   2  2
#   2  4

``

3.7.6 Diagonalmatrisen

En diagonalmatrise er en kvadratisk matrise der alle elementer utenfor hoveddiagonalen er null. Den kan skrives som:

\[ D_{ij} = \begin{cases} d_i & \text{hvis } i = j \\ 0 & \text{ellers} \end{cases} \]

Der \(d_1, d_2, \dots, d_n\) er elementene på diagonalen. Vi kan skrive dette kompakt som: \[ D = \mathrm{diag}(d_1, d_2, \dots, d_n) \]

For eksempel, når \(n = 4\): \[ D = \begin{bmatrix} 2 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 0.5 \end{bmatrix} \]

Multiplikasjon med en diagonalmatrise tilsvarer skalering av hver komponent i en vektor: \[ Dx = \begin{bmatrix} d_1 x_1 \\ d_2 x_2 \\ \vdots \\ d_n x_n \end{bmatrix} \]

Diagonalmatrisen er som idenditetsmatrisen et trivielt eksempel på en symmetrisk matrise, siden \(D = D^\top\) alltid gjelder.


De følgende eksemplene viser hvordan en diagonalmatrise kan initialisers og brukes I Julia og i Python

Julia:

using LinearAlgebra

# Lag en diagonalmatrise fra en vektor
d = [2.0, -1.0, 3.0, 0.5]
D = Diagonal(d)

x = [1.0, 2.0, 3.0, 4.0]
y = D * x

println("D =")
display(Matrix(D))
println("Dx = ", y)

Python:

import numpy as np

# Lag en diagonalmatrise fra en liste
d = [2.0, -1.0, 3.0, 0.5]
D = np.diag(d)

x = np.array([1.0, 2.0, 3.0, 4.0])
y = D @ x

print("D =")
print(D)
print("Dx =", y)

3.7.7 Permutasjonsmatrisen

En permutasjonsmatrise er en kvadratisk matrise som representerer en omstokking (permutasjon) av elementene i en vektor.

Den er konstruert fra identitetsmatrisen ved å bytte om på radene (eller kolonnene). Hver rad og hver kolonne inneholder nøyaktig én \(1\) og resten nuller.

Hvis \(P\) er en permutasjonsmatrise og \(x\) en vektor, vil \(Px\) være vektoren \(x\) med elementene omordnet.

For eksempel, hvis vi starter med \(I_3\) og bytter plass på rad \(1\) og rad \(3\) får vi:

\[ P = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix} \]

For en vektor \[ x = \begin{bmatrix} 5 \\ 6 \\ 7 \end{bmatrix} \] får vi \[ Px = \begin{bmatrix} 7 \\ 6 \\ 5 \end{bmatrix} \] #### I kode

Julia:

using LinearAlgebra

# Start med 3x3 identitetsmatrise
P = Matrix(I, 3, 3)

# Bytt rad 1 og rad 3
P[[1, 3], :] = P[[3, 1], :]

x = [5.0, 6.0, 7.0]
y = P * x

println("P =")
display(P)
println("Px =", y)

Python:

import numpy as np

# Start med 3x3 identitetsmatrise
P = np.eye(3)

# Bytt rad 1 og rad 3
P[[0, 2]] = P[[2, 0]]

x = np.array([5.0, 6.0, 7.0])
y = P @ x

print("P =")
print(P)
print("Px =", y)

Permutasjonsmatriser er nyttige i algoritmer som krever omstokking av elementer eller rader, og de bevarer lengden på vektorer siden de bare endrer rekkefølgen på komponentene, slik at

\[\| Px\| = \|x\|\]

3.8 Ytreprodukt

Et indreprodukt er ikke den eneste matriseoperasjonen som er mulig mellom vektorer. En annen viktig matriseoperasjon er ytreproduktet.

Gitt to vektorer \(x \in \mathbb{R}^n\) og \(a \in \mathbb{R}^m\), er ytreproduktet definert som:

\[ x a^\top = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} \underbrace{ \begin{bmatrix} a_1 & a_2 & \ldots & a_m \end{bmatrix} }_{1 \times m} = \begin{bmatrix} x_1 a_1 & x_1 a_2 & \ldots & x_1 a_m \\ x_2 a_1 & x_2 a_2 & \ldots & x_2 a_m \\ \vdots & \vdots & \ddots & \vdots \\ x_n a_1 & x_n a_2 & \ldots & x_n a_m \end{bmatrix} \]

Resultatet er alltid en \(n \times m\)-matrise, der hver rad er en skalert kopi av vektoren \(a^\top\).

Symbolet \(\otimes\) brukes ofte i mer teoretiske sammenhenger for ytreprodukt, men i dette kompendiet vil ytreproduktet skrives som et matriseprodukt \(x a^\top\).


Kodekesempler

Julia:

# To vektorer med ulik lengde
x = [1.0, 2.0, 3.0]   # n = 3
a = [4.0, 5.0]        # m = 2

# Ytreprodukt
A = x * a'

println("x = ", x)
println("a = ", a)
println("x * a' =")
display(A)

Python:

import numpy as np

# To vektorer med ulik lengde
x = np.array([1.0, 2.0, 3.0])  # n = 3
a = np.array([4.0, 5.0])       # m = 2

# Ytreprodukt
A = np.outer(x, a)

print("x =", x)
print("a =", a)
print("np.outer(x, a) =")
print(A)

Ytreproduktet er en effektiv måte å bygge matriser fra vektorer på en systematisk måte, og brukes ofte i numeriske metoder, statistikk og maskinlæring.

3.9 Eksempler på matriser

Matriser gir ofte en god måte å lagre data og dette gir et godt utgangspunkt for videre beregninger.

3.9.1 Materialdata

En matrise gir en kompakt måte å samle og organisere mange datasett på.
Vi kan for eksempel beskrive ett materiale med en egenskapsvektor

\[ p = (\rho, E, \nu, \alpha, \kappa) \]

der komponentene kan være tetthet, elastisitetsmodul, Poissons tall, termisk ekspansjonskoeffisient og varmeledningsevne.

For \(n\) materialer samler vi disse vektorene i en egenskapsmatrise \(\mathcal{X}\), der hver rad tilsvarer ett materiale og hver kolonne en egenskap:

\[ \mathcal{X} = \begin{bmatrix} p_1^\top \\ p_2^\top \\ \vdots \\ p_n^\top \end{bmatrix} \]
Dette formatet gjør det mulig å sammenligne materialer, gjøre statistiske analyser eller bruke dataene i maskinlæring.


3.9.2 Beslutningstøtte: materialvalg med matrise–vektor-multiplikasjon

I forrige kapittel så vi hvordan vi kunne bruke indreprodukt til å gi en totalscore for ett materiale basert på en vektvektor \(w\) som uttrykker hvor viktige de ulike egenskapene er.
Når vi har mange materialer samlet i en matrise \(\mathcal{X}\), kan vi beregne score for alle materialene samtidig ved å bruke matrise–vektor-produktet:

\[ s = \mathcal{X} w \]

der \(s\) er en vektor med én totalscore per materiale.


Eksempel:

\[ \mathcal{X} = \begin{bmatrix} 1.00 & 0.00 & 0.90 & 0.50 \\ 0.30 & 0.60 & 0.40 & 0.10 \\ 0.60 & 1.00 & 0.00 & 1.00 \end{bmatrix} \]

Rad 1 = stål, rad 2 = aluminium, rad 3 = CFRP.

Vektvektoren er valgt som før:

\[ w = \begin{bmatrix} 0.4 \\ 0.2 \\ 0.1 \\ 0.3 \end{bmatrix} \]
som representerer hvor mye vi vektlegger henholdsvis høy stivhet, lav masse, lav pris og lav termisk ekspansjon.

Matrise–vektor-produktet gir:

\[ s = \begin{bmatrix} 0.73 \\ 0.36 \\ 0.60 \end{bmatrix} \]

slik at stål får høyest totalscore i dette eksemplet.


Python

import numpy as np

w = np.array([0.4, 0.2, 0.1, 0.3])
X = np.array([
    [1.00, 0.00, 0.90, 0.50],  # Stål
    [0.30, 0.60, 0.40, 0.10],  # Aluminium
    [0.60, 1.00, 0.00, 1.00]   # CFRP
])

s = X @ w
print(s)  # [0.73 0.36 0.60]

Julia:

w = [0.4, 0.2, 0.1, 0.3]
X = [
    1.00  0.00  0.90  0.50;  # Stål
    0.30  0.60  0.40  0.10;  # Aluminium
    0.60  1.00  0.00  1.00   # CFRP
]

s = X * w
println(s)  # [0.73, 0.36, 0.60]

Dette er et eksempel på matrise–vektor-operasjonen for materialvalg, men det benyttes ofte også for andre typer ingeniørbeslutninger.

3.9.3 Tidsserier

Matriser er et naturlig format for å lagre tidsseriedata, når vi har målinger av flere variable på flere tidspunkter.

Anta at vi måler temperaturen på \(n\) ulike punkter, ved \(m\) ulike tidspunkter
\(t = (t_1, t_2, \dots, t_n)\). Vi kan da samle dataene i en matrise:

\[ \mathcal{T} = \begin{bmatrix} T_1(t_1) & T_1(t_2) & T_1(t_3) & \dots & T_1(t_n) \\ T_2(t_1) & T_2(t_2) & T_2(t_3) & \dots & T_2(t_n) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ T_m(t_1) & T_m(t_2) & T_m(t_3) & \dots & T_m(t_n) \end{bmatrix} \]

  • Hver rad viser utviklingen av én variabel (her: temperatur på ett punkt) over tid.
  • Hver kolonne viser målinger for alle punktene på ett gitt tidspunkt.

Denne strukturen gjør det enkelt å: - Finne gjennomsnittstemperaturen for hvert tidspunkt (gjennomsnitt langs rader) - Analysere trender for et gitt punkt (kolonnevis) - Bruke lineære algebraoperasjoner for filtrering eller prediksjon

Slike matriseformater brukes i mange tekniske sammenhenger, fra værdata og prosesskontroll til signalbehandling og maskinlæring.

3.9.3.1 Hvorfor skrive det som en matrise?

Når dataene ligger i matriseform, kan vi utføre mange analyser med enkle operasjoner:

  • Gjennomsnittstemperatur per målepunkt (over tid):
    \[ \bar{T} = \frac{1}{n} \, \mathcal{T} \, \mathbf{1} \]

    der \(\mathbf{1}\) er en kolonnevektorer med lenge \(n\) gir antall tider .

  • Finn maksimalverdi per målepunkt (nyttig for å se toppbelastning):

import numpy as np
T = np.array([...])  # n x m målematrise
maksverdier = np.max(T, axis=1)
  • Sammenligne to målepunkter ved hjelp av korrelasjon for å se om de har likt mønster.

Her kommer eksemepl for gjennomsnitt over tid i kode:

Python:

import numpy as np

# n = 3 sensorer, m = 4 tidspunkter
T = np.array([
    [20.1, 20.3, 20.5, 20.4],  # Sensor 1
    [19.8, 19.9, 20.1, 20.2],  # Sensor 2
    [21.0, 21.2, 21.3, 21.1]   # Sensor 3
])

n = T.shape[1]
gjennomsnitt = (1/n) * T @ np.ones(n)
print(gjennomsnitt)  # [20.325 20.0 21.15]

Julia-eksempel: gjennomsnitt over tid

T = [
    20.1  20.3  20.5  20.4;  # Sensor 1
    19.8  19.9  20.1  20.2;  # Sensor 2
    21.0  21.2  21.3  21.1   # Sensor 3
]

n = size(T, 2)
gjennomsnitt = (1/n) * T * ones(n)
println(gjennomsnitt)  # [20.325, 20.0, 21.15]

3.9.4 Matriser som bildeformat

Matriser kan brukes til å representere bilder. Et gråtonebilde (svart–hvitt) kan lagres som en matrise der hver inngang representerer lysstyrken til én piksel, ofte med verdier mellom 0 (svart) og 1 (hvitt), eller mellom 0 og 255 ved heltallsrepresentasjon.

For eksempel: \[ \text{Gråtonebilde} \;\longleftrightarrow\; P \in \mathbb{R}^{h \times w} \]

der (h) er høyde (antall piksler i vertikal retning), og (w) er bredde (horisontal retning). Hvert element (P_{ij}) angir intensiteten i piksel nummer ((i, j)).

For fargebilder brukes tre matriser — én for hver av de tre fargekanalene: rød (R), grønn (G) og blå (B): \[ \text{Fargebilde} \;\longleftrightarrow\; (R, G, B) \]

Disse tre kanalene kombineres for å bestemme fargen i hver piksel. Dette formatet kalles RGB-representasjon og er standard i digitale bilder.

Nedenfor ser vi et eksempel med et bilde av tåke. Til høyre er den røde kanalen justert ved en enkel matriseoperasjon: \[ R' \;=\; 1.4\,R \;-\; 60 \] etterfulgt av klipping til intervallet ([0, 255]) for å sikre gyldige pikselverdier.

Original og endret bilde av tåke

Du kan selv velge et bilde fra https://www.publicdomainpictures.net (f.eks. dette: https://www.publicdomainpictures.net/en/view-image.php?image=232812&picture=fog) og lagre det som fog.jpg.


Python

Python-eksempelet (Pillow) bruker heltall i området 0–255.

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt

# Load
img = Image.open("fog.jpg").convert("RGB")
data = np.array(img).astype(float)  # 0..255 as float

# Separate channels
R, G, B = data[:, :, 0], data[:, :, 1], data[:, :, 2]

# Modify: exaggerate red
R_mod = R * 1.4 - 60
R_mod = np.clip(R_mod, 0, 255)

# Recombine
modified = np.stack([R_mod, G, B], axis=2).astype(np.uint8)

# Plot
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
axes[0].imshow(img);       axes[0].axis("off")
axes[1].imshow(modified);  axes[1].axis("off")
plt.tight_layout()
plt.show()

Julia: Julia-eksempelet bruker standardområdet [0,1][0,1][0,1]. Derfor tilsvarer -60 i Python et bias på -60/255 i Julia.

using FileIO, Images, ColorTypes

# Load image as RGB (values in [0,1])
img = load("fog.jpg")          # Matrix{RGB{N0f8}}
ch  = channelview(img)         # 3 × H × W

# Separate channels (copies for clarity)
R = ch[1, :, :]
G = ch[2, :, :]
B = ch[3, :, :]

# Modify red channel: R * 1.4 - 60 (scale bias to [0,1])
R_mod = clamp.(R .* 1.4 .- (60/255), 0, 1)

# Recombine
ch_mod = copy(ch)
ch_mod[1, :, :] .= R_mod
modified = colorview(RGB, ch_mod)

# Save side-by-side comparison
panel = hcat(img, modified)
save("fog_red_adjustment.png", panel)

Merk (Julia): Bruk av punkt-operatorer (f.eks. .*, .-, clamp.(...), .=) betyr at operasjonen utføres elementvis over hele matrisen. Uten punktene forsøker Julia å tolke * som matriseprodukt og clamp som en skalaroperasjon.

Bilder kan brukes som “test-arena” for å utforske ulike transformasjoner.

3.9.5 Støyfjerning med båndmatrise

Data fra eksperimenter – og ofte også fra datamodeller og simuleringer – inneholder gjerne en god del støy. For eksperimentelle data er det ikke uvanlig å ta gjennomsnittet av mange målinger (som ofte gjøres uten at brukeren et det), men det er ikke alltid mulig, og selv etter dette så kan det være mye gjenværende støy. Slik søye kan skape numeriske problemer og gi opphav til feilaktige resultater, f.eks. så har støy svært store numeriske derivater. Et vanlig første steg i støyfjerning er å bruke et støyfilter som reduserer tilfeldige variasjoner, men bevarer hovedtrenden i dataene.

En av de enkleste metodene er å ta et løpende gjennomsnittet av nærmeste naboene:

\[ y_i = \tfrac{1}{3}\big(x_{i-1} + x_i + x_{i+1}\big). \]

Samler vi datapunktene i en vektor \(x \in \mathbb{R}^n\), kan den gladdete vektoren \(y\) skrives som

\[ y = Sx, \]
der \(S\) er en \(n \times n\) båndmatrise av formen:

\[ S = \tfrac{1}{3} \begin{bmatrix} 1 & 1 & 0 & \cdots & 0 \\ 1 & 1 & 1 & & 0 \\ 0 & 1 & 1 & 1 & 0 \\ \vdots & & \ddots & \ddots & \ddots \\ 0 & \cdots & 0 & 1 & 1 \end{bmatrix}. \]

Julia-eksempel:

using LinearAlgebra, Random, Plots

# Discretize function f(x) = x * exp(-x) on interval [0, 5]
n = 200 
xs = range(0, 5, length=n)
f_true = xs .* exp.(-xs)

# Add some random noise
Random.seed!(0)
noise = 0.02 * randn(n)
f_noisy = f_true + noise

# Build smoothing matrix S (111 filter, with simple boundary handling)
S = zeros(n, n)
for i in 2:n-1
    S[i, i-1:i+1] .= 1/3 
end
# Boundary rows: just average available points
S[1, 1:2] .= 1/2 
S[n, n-1:n] .= 1/2 

# Apply filter
f_smooth = S * f_noisy

# Plot
plot(xs, f_true, lw=2, label="True function")
plot!(xs, f_noisy, lw=1.4, alpha=1, label="Noisy data")
plot!(xs, f_smooth, lw=1.4, label="Smoothed (111 filter)")
xlabel!("x"); ylabel!("f(x)")
title!("Smoothing with banded averaging matrix")

Dette gir følgende resultat

Som vi kan se så har det verste toppene forsvunnet, men det er fortsatt mye gjenværende støy.

Vi kan derfor prøve å fram med større bredder, her har vi brukt 15 og 31. Vi kan se at det reduserer støyen ytterligere, men når vi bruker så mange som 31 datapunkter mister vi mye av informasjonen spesielt rundt funksjonens maksimum.

![[denoise2.png]] Kode:

using LinearAlgebra, Random, Plots

# Funksjon og støy
n = 200
xs = range(0, 5, length=n)
f_true = xs .* exp.(-xs)

Random.seed!(1)
f_noisy = f_true .+ 0.02*randn(n)

# Hjelpefunksjon for smoothing-matriser
function smoothing_matrix(weights, n)
    k = length(weights)
    r = div(k,2)
    S = zeros(n,n)
    normed = weights ./ sum(weights)
    for i in 1:n
        for j in max(1,i-r):min(n,i+r)
            S[i,j] = normed[j-i+r+1]
        end
    end
    return S
end

# Tre uniforme filtre
S3  = smoothing_matrix(ones(3), n)
S15 = smoothing_matrix(ones(15), n)
S31 = smoothing_matrix(ones(31), n)

# Påfør
f_3  = S3 * f_noisy
f_15 = S15 * f_noisy
f_31 = S31 * f_noisy


# Plot
palette(:prism)

plot(xs, f_noisy, alpha=0.4, label="Støyete data", lw=1, color=:red)
plot!(xs, f_true, lw=2, label="f(x) = x e^{-x}", color=:black)
plot!(xs, f_15, lw=2, label="15-punkts gj.snitt", color=:blue, alpha=0.7)
plot!(xs, f_31, lw=2, label="31-punkts gj.snitt", alpha=0.7, color=:brown)
xlabel!("x"); ylabel!("f(x)")
title!("Uniforme filtre med ulik bredde")

Du kan selv prøve egne bredder. Eller kanskje varierer vekten, f.eks. vekt midten mer enn kantene.

3.9.6 Binning

Et alternativ til løpende gjennomsnitt, er å redusere antall datapunkter med bruk av “binning”. Det vi tar gjennomsnittet over et område og kollapser verdiene til et punkt.

Binning grupperer vi punktene i større bolker, og erstatter hele blokken med én verdi.
Da reduseres dimensjonen — matrisen blir rektangulær, f.eks. \(m \times n\) der \(m < n\).

Et enkelt eksempel: anta at vi har \(n = 12\) datapunkter, og vi vil samle dem i \(m = 4\) bokser à 3 punkter hver.
Da blir binning-matrisen

\[ B = \tfrac{1}{3} \begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \end{bmatrix}. \]

Da får vi \(y = Bx\) med bare 4 datapunkter igjen, hvor hvert punkt er gjennomsnittet av 3 originaler.

Nedenfor, ser vi resultater av det samme eksemplet, men “bin”-størrelser på 10.

using LinearAlgebra, Random, Plots

# Funksjon og støy
n = 200
xs = range(0, 5, length=n)
f_true = xs .* exp.(-xs)

Random.seed!(1)
f_noisy = f_true .+ 0.02*randn(n)

function binning_matrix(n, bin_size)
    m = div(n, bin_size)
    B = zeros(m, n)
    for i in 1:m
        B[i, (i-1)*bin_size+1:i*bin_size] .= 1/bin_size
    end
    return B
end

# Eksempel: 200 punkter, bin-størrelse 10
n = 200
bin_size = 10
B = binning_matrix(n, bin_size)

# Bruk
f_binned = B * f_noisy
xs_binned = range(0, 5, length=size(B,1))

plot(xs, f_noisy, alpha=0.4, label="Støyete data", lw=1, color=:red)
plot!(xs_binned, f_binned, lw=2, markershape=:circle, label="Binned (10 punkter)")
xlabel!("x"); ylabel!("f(x)")
title!("Uniforme filtre med ulik bredde")

Dette gir følgende figur:

Binning

Vi kan nå se at vi får en mye mer behagelig kurve som er lettere å visualisere.
Om vi i tillegg velger å interpolere disse datapunktene (f.eks. med splines) kan vi få virkelig pene kurver. Kurvetilpassing skal vi diskutere mer senere i kompendiet.

I kommersielle programmer, f.eks. for prosessering av optiske spektra, brukes ofte binning eller automatisk glatting.
Dette kan gjøre dataene enklere å håndtere og raskere å prosessere, men det reduserer også oppløsningen.

Binning og spesielt interpolering basert på kurvetilpassing kan være risikable strategier:
skarpe topper eller små skuldrer kan forsvinne helt, eller nye trekk kan se ut til å dukke opp.
Derfor er det noen ganger nødvendig å skru av slike forhåndsvalg og heller arbeide med rådata, for deretter å bruke mer kontrollerte metoder (som løpende gjennomsnitt eller mer avanserte filtre) til å redusere støy.

Det finnes mer avanserte metoder for støyreduksjon, som vi også kommer tilbake til senere i kompendiet.

3.9.7 Mange partikler:

Om vi skal uttrykke posisjonen til mange partikler, kan det gjøres men en matrise:

\[ \mathcal{R} = \begin{bmatrix} r_1 & r_2 & \ldots & r_n \end{bmatrix} \] eller mer eksplisitt

\[ \mathcal{R} = \begin{bmatrix} x_1 & x_2 & \cdots & x_n \\ y_1 & y_2 & \cdots & y_n \\ z_1 & z_2 & \cdots & z_n \\ \end{bmatrix} \]

hvor hver kolonne representerer én partikkel, og hver rad representerer én koordinatretning \((x, y, z)\), hvor \(\mathcal{R} \in \mathbb{R}^{3 \times n}\).

Men noen ganger er det nyttig å uttrykke dette som en lang \(3n\)-vektor,

\[ r = (x_1, y_1, z_1, x_2, y_2,z_2 \ldots, x_n, y_n, z_n ) \]

for eksempel hvis man ønsker å lagre posisjonene til alle partiklene ved ulike tidspunkter som rader i en matrise, der hver rad er én tilstand. Dette er vanlig i simuleringer og datastrukturer hvor hvert tidspunkt eller scenario lagres som én vektor.

For å stable tallene på denne måten, kan man bruke “flatten” i Python

import numpy as np
# Lage en 3 x 4 matrise for 4 partikler
R = np.array([
    [1.0, 2.0, 3.0, 4.0],   # x-koordinater
    [0.5, 0.0, -0.5, -1.0], # y-koordinater
    [9.0, 8.0, 7.0, 6.0]    # z-koordinater
])

print("R =")
print(R)

# Gjør om til 3n-vektor med rekkefølge: x1, y1, z1, x2, y2, z2, ...
r = R.T.flatten()

print("\nr =")
print(r)
# Lage en 3 x 4 matrise for 4 partikler
R = [
    1.0  2.0  3.0  4.0;   # x-koordinater
    0.5  0.0 -0.5 -1.0;   # y-koordinater
    9.0  8.0  7.0  6.0    # z-koordinater
]

println("R =")
display(R)

# Gjør om til 3n-vektor: x1, y1, z1, x2, y2, z2, ...
r = vec(R')

println("\nr =")
display(r)

I begge tilfellene måtte vi transponere matrisene før vi flatet dem ut.
I Python gjorde vi dette med R.T, mens vi i Julia skrev R'.

Ved en transponering bytter vi rader og kolonner i matrisen, slik at elementet på plass \((i, j)\) i \(A\) flyttes til plass \((j, i)\) i \(A^\top\) , slik at: \[ (A^\top)_{ij} = A_{ji} \]

Gitt for eksempel matrisen:

\[ A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix} \]

Så blir den transponerte matrisen:

\[ A^\top = \begin{bmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{bmatrix} \]

Transponering var nødvendig her når vi skulle flate ut matrisen med flatten eller vec pga. konvensjonen som disse metodene brukte. >Når man skal reorganisere dimensjonene til en matrise går det ofte galt og det er lurt å gjøre litt ekstra testing for å se om alt ble riktig.

3.9.8 Fjærkrefter og stivhetsmatrise: to partikler, tre fjærer

Koblede fjærer er viktige modellsystemer for å forstå spenninger og vibrasjoner i teknisk design, som f.eks. i FEM-analyse. De er også et eksempel på hvordan matrise-vektor-multiplikasjon oppstår naturlig i fysikk og ingeniørfag.


Eksempelproblem

  • To partikler kan bevege seg langs en linje.
  • De er koblet med tre identiske fjærer med fjærkonstant \(k\).
  • Systemet er i likevekt når begge partikler står i ro med fjærene i sin naturlige lengde.
  • Vi ser kun på små forskyvninger fra likevekt: \(x_1\) og \(x_2\).
  • En positiv verdi betyr at partikkelen er flyttet mot høyre.

Springs

Krefter på partiklene

Fjærene følger Hookes lov: kraft = \(k \cdot\) forlengelse.

Partikkel 1: - Venstre fjær gir kraft \(-k x_1\) - Midtre fjær gir kraft \(k(x_2 - x_1)\)

Totalt: \[ f_1 = -k x_1 + k(x_2 - x_1) = -2k x_1 + k x_2 \]

Partikkel 2: - Høyre fjær gir kraft \(-k x_2\) - Midtre fjær gir kraft \(-k(x_2 - x_1)\)

Totalt: \[ f_2 = -k(x_2 - x_1) - k x_2 = k x_1 - 2k x_2 \]


I matriseform kan dette skrives kompakt som

\[ \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} = k \begin{bmatrix} -2 & 1 \\ 1 & -2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \]

Dette gir et godt utgangspunkt for å løse dynamikken til slike ligningssystemer på en effektiv måte.

3.9.9 Eksempel: Diskretisering av deriverte

I forrige kapittel så vi hvordan en funksjon \(f(x)\) kan diskretiseres og gjøres til en vektor. Om vi velger ut (en vektor) av posisjoner \[ x_\text{vec} = (x_1, x_2, \ldots, x_n), \] med avstand \(h\) mellom punktene, så får vi vektoren \[ f_\text{vec} = (f(x_1), f(x_2), \ldots, f(x_n)). \]

Nå kan vi utrykke derivater av en funksjon med matrise-vektor produkt.


3.9.9.1 Førstederiverte

En enkel og mye brukt approksimasjon for den første deriverte er sentraldifferansen:

\[ f'(x_i) \approx \frac{f(x_{i+1}) - f(x_{i-1})}{2h}. \]

Dette kan vi samle i en derivert-matrise \(D\) som virker på:

\[ f'_\text{vec} \approx D f_\text{vec}. \]

Om vi tar for oss en periodske funksjon \(f(x) =f(x+L)\), så gjeldet et at \(f(0) =f(L)\), og begge disse punktetene blir gitt av elementet \(f_{\text{vec}, 1} = f(0) =f(L)\) , mens \(f_{\text{vec}, n} = f(L-h)\)

Da blir derivert-matrisen gitt ved:

\[ D = \tfrac{1}{2h} \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 & -1 \\ -1 & 0 & 1 & & & 0 \\ 0 & -1 & 0 & 1 & & \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ 0 & & & -1 & 0 & 1 \\ 1 & 0 & \cdots & 0 & -1 & 0 \end{bmatrix}. \]

TEMP: \[ y = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1}\\ y_n \\ y_{n+1} \\ \vdots\\ y_N\end{bmatrix}. \]

Punktene nededer til venstre og øverst til høyre svarer da til deriverte som trenger første eller siste indeks for å evalueres.


using LinearAlgebra, Plots

# Antall punkter og steglengde
n = 20
L = 2*pi     # intervall-lengde
h = L / n

# Gridpunkter og funksjonsverdier
x_vec = range(0, L - h, length=n)   # 0, h, 2h, ..., L-h
f_vec = sin.(x_vec)

# Bygg derivert-matrise (periodisk, sentral differanse)
D = zeros(n, n)
for i in 2:n-1
    D[i, i-1] = -1/(2*h)
    D[i, i+1] =  1/(2*h)
end
# Periodiske hjørner
D[1, 2]   =  1/(2*h)
D[1, n]   = -1/(2*h)
D[n, 1]   =  1/(2*h)
D[n, n-1] = -1/(2*h)

# Numerisk første-deriverte
f1_num = D * f_vec

# Analytisk første-deriverte
f1_true = cos.(x_vec)

# Plot sammenligning
plot(x_vec, f1_true, label="cos(x)", lw=2, color=:blue)
plot!(x_vec, f1_num, label="D * sin(x)", lw=2, ls=:dash, color=:red)
xlabel!("x")
ylabel!("derivert")
title!("Første-derivert: numerisk vs analytisk")

3.9.9.2 Andrederiverte

På tilsvarende måte kan vi bruke en 3-punkts formel for den andre deriverte:

\[ f^{\prime\prime}(x_i) \approx \frac{f(x_{i+1}) - 2f(x_i) + f(x_{i-1})}{h^2}. \]

Dette gir en Laplacian-matrise \(L\):

\[ L = \tfrac{1}{h^2} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 1 \\ 1 & -2 & 1 & & & 0 \\ 0 & 1 & -2 & 1 & & \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ 0 & & & 1 & -2 & 1 \\ 1 & 0 & \cdots & 0 & 1 & -2 \end{bmatrix}. \]

using LinearAlgebra, Plots

# Antall punkter og steglengde
n = 10
L = 2*pi          # intervall-lengde
h = L / n

# Gridpunkter og funksjonsverdier
x_vec = range(0, L - h, length=n)   # 0, h, 2h, ..., L-h
f_vec = sin.(x_vec)

# Bygg andre-derivert matrise (periodisk, 3-punkts formel)
Lmat = zeros(n, n)
for i in 2:n-1
    Lmat[i, i-1] =  1/h^2
    Lmat[i, i]   = -2/h^2
    Lmat[i, i+1] =  1/h^2
end
# Periodiske hjørner
Lmat[1, 1] = -2/h^2;   Lmat[1, 2] = 1/h^2;    Lmat[1, n] = 1/h^2
Lmat[n, n] = -2/h^2;   Lmat[n, n-1] = 1/h^2;  Lmat[n, 1] = 1/h^2

# Numerisk andre-deriverte
f2_num = Lmat * f_vec

# Analytisk fasit: f''(x) = -sin(x)
f2_true = -sin.(x_vec)

# Plot sammenligning
plot(x_vec, f2_true, label="-sin(x)", lw=2, color=:blue)
plot!(x_vec, f2_num, label="L * sin(x)", lw=2, ls=:dash, color=:red)
xlabel!("x")
ylabel!("andre derivert")
title!("Andre-derivert: numerisk vs analytisk")

I eksemplene her har vi brukt enkle differanseformler (såkalte stensiler).
Det finnes også mer avanserte derivat-matriser som gir høyere nøyaktighet.

Vi har også valgt periodiske randbetingelser, som passer godt for funksjoner
som repeterer seg. Andre randbetingelser er mulig, for eksempel faste ender
(Dirichlet) eller frie ender (Neumann), og det gir litt andre matriser.

3.10 Nettverksmatriser

Matriser er også veldig nyttige for å beskrive relasjoner mellom ting.
Studiet av slike relasjoner kalles grafteori, hvor noder representerer objekter og kanter representerer relasjoner. Lineær algebra er et viktig verktøy i slike studier.

Eksempel: Om vi har fire nettsider \(A,B,C,D\). Hvis en side linker til en annen, gir vi matrise-elemenet verdien \(1\), ellers får den verien \(0\).

Dette kan uttrykkes som en nettverksmatrise (adjacency-matrise)

\[ M = \begin{bmatrix} 0 & 1 & 1 & 0 \\ 0 & 0 & 1 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}. \]

Her betyr f.eks. \(M_{12}=1\) at side \(A\) linker til side \(B\), mens \(M_{21}=0\) viser at \(B\) ikke linker tilbake til \(A\).
Vi ser at matrisen ikke er symmetrisk, fordi lenker ofte går bare én vei.

Nettverksmatriser brukes i mange sammenhenger:

  • Søkemotorer: Googles opprinnelige PageRank-algoritme bygger på en slik linker for å rangere nettsider.
  • Sosiale nettverk: Hvem som følger hvem på Twitter eller Instagram.
  • Transport og logistikk: I veinettverk eller flyruter.
  • Infrastruktur og energi: Strømnett eller vannledninger kan modelleres som grafer med tilhørende matriser for å analysere stabilitet.

Et annet eksempel finner vi i kjemi. Molekyler kan beskrives som grafer.

For vannmolekylet H₂O kan vi lage en netteverksmatrise. Vi kan representere atomene som tre noder: ett O og to H. Hvis vi velger rekkefølgen (O, H₁, H₂), hvor subindeksen angir hvilket H-atom vi refererer til, får vi følgende nettverksmatrise:

\[ M = \begin{bmatrix} 0 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{bmatrix}. \]

Her indikerer \(M_{12} = 1\) at O er bundet til H₁, mens \(M_{23} = 0\) viser at H₁ og H₂ ikke er direkte bundet.

Slike matriser brukes bl.a. i maskinlæring for materialer og molekyler.

3.11 Funksjoner med matrise-vektorprodukt

3.11.1 Kvadratiske former

I mange sammenhenger møter vi funksjoner av typen

\[ f(x) = x^\top A x, \]
hvor \(A^\top = A\) er symmetrisk. Slike uttrykk kalles kvadratiske former.

Eksempel:
\[ f(x_1,x_2) = x_1^2 + 4x_1x_2 + x_2^2. \]

Dette kan skrives med matrisen \[ A = \begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix}, \]
slik at \(f(x) = x^\top A x\).


Kvadratiske former dukker opp i en rekke situasjoner:
- areal- og volum-beregninger (ellipser og ellipsoider),
- mekanikk (energiuttrykk, treghetsmomenter),
- statistikk og maskinlæring (kovariansmatriser, PCA).


Kvadratiske former oppstår f.eks. når vi vil uttrykke lengder og areal i ikke-kartesiske koordinater.

Tenk at vi har to vektorer \(u_1\) og \(u_2\) som spenner ut et plan.
For eksempel

\[ u_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}= e_1\,, \qquad u_2 = \begin{bmatrix} \tfrac{1}{2} \\ \tfrac{\sqrt{3}}{2} \end{bmatrix} = \tfrac{1}{2} e_1 + \tfrac{\sqrt{3}}{2} e_2 \,, \]

som tilsvarer en enhetlig lengde i \(x\)-retning og en vektor med 60° vinkel.

Hvis vi beskriver en vektor i dette koordinatsystemet som

\[ x = \alpha u_1 + \beta u_2, \]

som vi kan velge å skrive på vektorform (siden dette er en ordnet list med tall), i basisen \({u_1,u_2}\), ikke i standardbasis \({e_1,e_2}\).

\[ x= (\alpha, \beta) = \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \]

så er lengden kvadrert gitt ved den kvadratiske formen

\[ \|x\|^2 = (\alpha u_1 + \beta u_2)^\top (\alpha u_1 + \beta u_2) = \begin{bmatrix} \alpha & \beta \end{bmatrix} \; G \; \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \]

der \(G\) er en Gram-matrise,

\[ G = \begin{bmatrix} u_1^\top u_1 & u_1^\top u_2 \\ u_2^\top u_1 & u_2^\top u_2 \end{bmatrix}. \]

For vårt valg av \(u_1\) og \(u_2\) får vi \[ G = \begin{bmatrix} 1 & \tfrac{1}{2} \\ \tfrac{1}{2} & 1 \end{bmatrix}. \]

Dette viser hvordan kvadratiske former kan uttrykker lengder og arealer i et koordinatsystem hvor aksene står i 60° vinkel, ikke 90°.

Gram matrisen kan altså brukes til å definere nye normer eller metrikker med indreprodukt av typen:

\[ \langle x, y \rangle_G = x^\top G y. \]

Med dette indreproduktet, kan vi også definere vinkelen mellom to vektorer \(x\) og \(y\) som

\[ \cos \theta = \frac{\langle x, y \rangle_G}{\sqrt{\langle x, x \rangle_G}\,\sqrt{\langle y, y \rangle_G}}. \]

Dette er helt analogt med hvordan vi definerer vinkler i det vanlige kartesiske koordinatsystemet.

I dette eksempelet startet vi med vektorer der 2-normen ga et naturlig avstandsmål, og dette er selvsagt passende i kartesiske koordinatsystemer.
Men i dette kompendiet bruker vi også indreprodukt og normer til å sammenligne likheter mellom datavektorer som representerer f.eks. bilder eller tekster. Er det gitt at en vanlig 2-norm alltid er den mest naturlige? På ingen måte!

Statistiske metoder som PCA og mange maskinlæringsmetoder kan nettopp forstås som forsøk på å oppdage et mer hensiktsmessig indreprodukt eller norm for å sammenligne likheter mellom vektorer på en måte som gjenspeiler strukturen i dataene. I praksis betyr dette ofte at en «skjult metrikk» blir brukt – enten gjennom en transformasjon (som i PCA eller nevralnett) eller en Gram-matrise (som i kernelmetoder).

For tekstdokumenter kan vi f.eks. starte med en enkel modell der hvert ord tilsvarer en uavhengig basisvektor \(e_1, e_2, \dots\) (som beskrevet i kapittel 1). Da blir to tekster sammenlignet via vinkelen mellom deres tekstvektorer.

Men vi kan ta ett steg videre: vi kan bruke en basis som ikke er ortogonal,
der enkelte ord allerede ligger «nærmere» hverandre. Slik kan ord som konge og dronning forstås som vektorer med liten vinkel mellom seg, mens ord som konge og bil har nesten 90° vinkel.

Moderne språkmodeller (LLMs) går enda lengre og behandler tekstfragmenter som vektorer,
som de så kan sammenligne ved å vurdere avstander og vinkler mellom dem.


3.12 Taylorrekke

I mange anvendelser er det nyttig å approksimere funksjoner som polynomer. Dette kalles en Taylorrekke (eller Taylor-ekspansjon) rundt et punkt \(x_0\).

3.12.0.1 Én variabel

For en funksjon \(f:\mathbb{R}\to\mathbb{R}\), så får vi til andre i orden i \((x-x_0)\),

\[ f(x) \approx f(x_0) + f'(x_0)(x-x_0) + \tfrac{1}{2} f''(x_0)(x-x_0)^2. \]

Eksempel:

\[ f(x) = (1+x)\,e^{-x}. \]

  • \(f(0)=1\)
  • \(f'(x) = -x e^{-x}\), så \(f'(0)=0\)
  • \(f''(x)= (x-1)e^{-x}\), så \(f''(0)=-1\)

Tilnærming rundt \(x_0=0\):

\[ f(x) \approx 1 - \tfrac{1}{2}x^2. \]


3.12.0.2 Flere variabler (vektor-notasjon)

For \(f:\mathbb{R}^n \to \mathbb{R}\) og et punkt \(x_0\), sett \(\Delta x = x - x_0\). Da har vi

\[ f(x) \approx f(x_0) + \nabla f(x_0)^\top \Delta x + \tfrac{1}{2}\,\Delta x^\top H(x_0)\, \Delta x, \]

der
- \(\nabla f(x_0)\) er gradienten, en vektor av første-deriverte:

\[ \nabla f(x_0) = \begin{bmatrix} \frac{\partial f}{\partial x_1}(x_0) \\ \frac{\partial f}{\partial x_2}(x_0) \\ \vdots \\ \frac{\partial f}{\partial x_n}(x_0) \end{bmatrix}, \]

  • \(H(x_0)\) er Hessianen, matrisen av andre-deriverte:

\[ H(x_0) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{bmatrix}_{x=x_0}. \]

3.12.0.3 Eksempel i n-dimensjoner

La

\[ f(x) = 1 + a^\top x \, e^{-b^\top x}, \quad x\in \mathbb{R}^2, \]

med gitte vektorer \(a,b\in\mathbb{R}^n\). og

  • Funksjonsverdi i \(x_0=0\):

\[ f(0) = 1. \]

  • For å regne ut gradienten, bruker vi både produktregelen og kjederegelen. Kjederegelen er gitt ved (\(f:\mathbb{R}\to\mathbb{R}\) og \(g:\mathbb{R}^n\to\mathbb{R}\), ) \[ h(x) = f(g(x)), \qquad x\in\mathbb{R}^n. \] da er gradienten \[ \nabla h(x) = f'(g(x)) \, \nabla g(x). \]
    Dette bruker vi når vi regner ut gradient av \(e^{-b^\top x}\), og vi har at \(\nabla (b^\top x)= b\) og tilsvarende er også \(\nabla (x^\top b)=b\), siden vi kan bytte rekkenfølgen på indreproduktet av to vektorerer.

Sammen med produktregelen får vi da

\[ \nabla f(x) = e^{-b^\top x}\big( a - (a^\top x)b \big), \]
og dermed

\[ \nabla f(0)=a. \]

  • Hessianen ved x=0 er gitt ved:

\[ H(x) = -e^{-b^\top x}\big( b a^\top + a b^\top \big) + e^{-b^\top x}\big( (a^\top x)\, b b^\top \big). \]

Ved \(x_0=0\):

\[ H(0) = -\big( b a^\top + a b^\top \big). \]

Taylor-tilnærmelse rundt \(x_0=0\):

\[ f(x) \approx 1 + a^\top x - \tfrac{1}{2} x^\top \big( b a^\top + a b^\top \big) x. \]


Denne typen utvidelser brukes for å bygge en lokal lineær eller kvadratisk modell av en mer komplisert funksjon.


I eksempelet nedenfor plotter vi den opprinnelige, linere, og kvardatiske funksjonen. Vi har GLMakie for å kunne lage dynamiske figurer i 3d som kan roteres.

using GLMakie, LinearAlgebra

# sample points
xs = range(0, 1, length=20)
ys = range(0, 1, length=20)

# function setup
a = [1.0, 1.0]
b = [0.3, 0.4]
f(x) = 1 + (a'x)[1] * exp(-(b'x)[1])
grad0 = a
H0 = -(b * a' + a * b')
f_lin(x)  = 1 + (grad0' x)[1]
f_quad(x) = 1 + (grad0' x)[1] + 0.5 * (x' * H0 * x)[1]

# evaluate functions on grid
Zf    = [f([x,y]) for x in xs, y in ys]
Zlin  = [f_lin([x,y]) for x in xs, y in ys]
Zquad = [f_quad([x,y]) for x in xs, y in ys]

# figure with three subplots
fig = Figure(resolution=(900,300))
ax1 = Axis3(fig[1,1], title="Original f(x)", aspect=:data, limits=(0,1, 0,1, 0,3))
surface!(ax1, xs, ys, Zf; colormap=:viridis)

ax2 = Axis3(fig[1,2], title="Linear approx", aspect=:data, limits=(0,1, 0,1, 0,3))
surface!(ax2, xs, ys, Zlin; colormap=:viridis)

ax3 = Axis3(fig[1,3], title="Quadratic approx", aspect=:data, limits=(0,1, 0,1, 0,3))
surface!(ax3, xs, ys, Zquad; colormap=:viridis)

fig

Taylor2d

I høyere dimensjoner (og vanskelig nok i 2d) er det ikke så lett å visuelt se hvor godt slike Taylor ekspansjoner beskriver hvor godt en funksjon beskriver funksjon. En måte å løse dette på er å plotte \(|x|\)

using LinearAlgebra, Plots

# parameters
a = [1.0, 1.0]
b = [0.3, 0.4]

f(x) = 1 + a'x * exp(-b'x)
grad0 = a
H0 = -(b * a' + a * b')
f_lin(x)  = 1 + grad0'x
f_quad(x) = 1 + grad0'x + 0.5 * x'(H0*x) 

# grid
N = 20
xs = range(0, 1, length=N)
ys = range(0, 1, length=N)

# evaluate functions on grid (matrices)
Zf    = [f([x,y])      for x in xs, y in ys]
Zlin  = [f_lin([x,y])  for x in xs, y in ys]
Zquad = [f_quad([x,y]) for x in xs, y in ys]

# distances and errors (flattened)
dists     = vec([norm([x,y]) for x in xs, y in ys])
errs_lin  = vec(abs.(Zf .- Zlin))
errs_quad = vec(abs.(Zf .- Zquad))

# plot error vs |x|
scatter(dists, errs_lin;  label="Linear",    alpha=0.6, ms=4)
scatter!(dists, errs_quad; label="Quadratic", alpha=0.6, ms=4)
xlabel!("‖x‖")
ylabel!("Error")
title!("Taylor error vs distance from 0")

Taylor vs distance

Om du ønsker å sammenligne i høyere dimensjoner, så sammenligner man istedet ved å velge tilfeldige tall, ellers får man en dimensjons-eksplosjon.

Automatic differentiation: Man kan ofte velge å behandle output (f.eks. krefter) fra en ekstern kode, som en FEM-simulering, en prosessmodell eller molekyldynamikk, som en funksjon av variable man gir til modellen (f.eks. trykk eller en last). I noen tilfeller ønsker man å lage enklere representasjoner av dette, for eksempel ved en Taylor-ekspansjon. Da finnes det ingen analytisk formel, og vi er nødt til å regne ut derivater numerisk (slik vi illustrerte tidligere i 1D). I høye dimensjoner blir dette ofte både tregt og unøyaktig, spesielt for Hessianen.

Dersom funksjonen derimot er skrevet direkte i Julia (eller i tilsvarende rammeverk i andre språk, som JAX i Python), finnes det et svært elegant alternativ: automatic differentiation (AD). Dette bygger på at alle beregninger på en datamaskin kan brytes ned i elementære operasjoner (addisjon, multiplikasjon, eksponentialer osv.) der vi kjenner derivatene. Ved systematisk å bruke produkt- og kjederegelen kan AD regne ut gradienter og Hessianer eksakt, uten at vi selv trenger å utføre algebraen.

AD er ikke bare nyttig i små eksempler som her. Hele den moderne utviklingen av maskinlæring med nevrale nettverk er faktisk mulig fordi rammeverk som PyTorch, TensorFlow og JAX kan bruke AD til å beregne gradienter med millioner av parametre.

using LinearAlgebra, Plots, ForwardDiff

# parameters
a = [1.0, 1.0]
b = [0.3, 0.4]

f(x) = 1 + (a'x)[1] * exp(-(b'x)[1])

# expansion point
x0 = zeros(2)

# automatic differentiation
grad0 = ForwardDiff.gradient(f, x0)
H0    = ForwardDiff.hessian(f, x0)

# linear and quadratic approximations
f_lin(x)  = f(x0) + (grad0'(x - x0))[1]
f_quad(x) = f(x0) + (grad0' * (x - x0))[1] + 0.5 * ((x - x0)'H0(x - x0))[1]

# grid
N = 20
xs = range(0, 1, length=N)
ys = range(0, 1, length=N)

# evaluate functions on grid (matrices)
Zf    = [f([x,y])      for x in xs, y in ys]
Zlin  = [f_lin([x,y])  for x in xs, y in ys]
Zquad = [f_quad([x,y]) for x in xs, y in ys]

# distances and errors (flattened)
dists     = vec([norm([x,y]) for x in xs, y in ys])
errs_lin  = vec(abs.(Zf .- Zlin))
errs_quad = vec(abs.(Zf .- Zquad))

# plot error vs |x|
scatter(dists, errs_lin;  label="Linear",    alpha=0.6, ms=4)
scatter!(dists, errs_quad; label="Quadratic", alpha=0.6, ms=4)
xlabel!("‖x‖")
ylabel!("Error")
title!("Taylor error vs distance from 0")

# --- explicit check against manual formulas ---
println("Explicit check")
grad0_reg = a
H0_reg = -(b * a' + a * b')
f_lin_reg(x)  = 1 + (grad0_reg' * x)[1]
f_quad_reg(x) = 1 + (grad0_reg' * x)[1] + 0.5 * (x' * H0_reg * x)[1]

Zlin_reg  = [f_lin_reg([x,y])  for x in xs, y in ys]
Zquad_reg = [f_quad_reg([x,y]) for x in xs, y in ys]

println(norm(Zlin_reg - Zlin))
println(norm(Zquad_reg - Zquad))

Hvor mye lengre tid tar AD en analytisk differensiering (på datamaskinen)? Dette er et spørsmål du selv bør undersøke, man da bør du lage litt mer komplekse funksjoner enn eksemplet her.

3.13 Komplekse matriser

I mange anvendelser (signalbehandling, kvantemekanikk, bølgeproblemer) arbeider vi med komplekse matriser, der hvert element \(a_{ij}\) kan være et komplekst tall
\[ a_{ij} \in \mathbb{C}. \] ### Adjungert (konjugert transponert)

For en matrise \(A \in \mathbb{C}^{m\times n}\) defineres den adjungerte (eller konjugert transponerte) som
\[ A^\dagger = A^{* \top}, \] der \(A^*\) betyr at vi tar det komplekse konjugatet av hvert element, og \(\top\) betyr transponering.

Elementvis: \[ (A^\dagger)_{ij} = a_{ji}^*. \]

I Julia: A' gir adjungert (konjugert transponert), mens transpose(A) gir ren transponering uten konjugering.


3.13.1 Hermittisk matrise

En kvadratisk matrise \(A \in \mathbb{C}^{n\times n}\) er Hermittisk hvis
\[ A = A^\dagger. \] Dette betyr at \(a_{ii} \in \mathbb{R}\) (diagonalelementene er reelle) og at \(a_{ij} = a_{ji}^*\). Hermittiske matriser er komplekse analoger til reelle, symmetriske matriser, og har mange av de samme egenskapene.

3.13.2 Frobeniusnorm

Frobeniusnormen for komplekse matriser, ligner på den for reelle, \[ \|A\|_F = \sqrt{\sum_{i=1}^m \sum_{j=1}^n |a_{ij}|^2}. \] der \(|a_{ij}|^2\) er kvadratet av absoluttverdien til det komplekse tallet.

Eksempel i Julia I dette eksemplet sjekker vi om en matrise er Hermittisk og regner ut Frobenius normen.

using LinearAlgebra

# Lage en kompleks matrise
A = [1+2im   3-4im;
     3+4im   5+0im]

A_adj = A'             # adjungert (konjugert transponert)
is_hermitian = ishermittian(A) # sjekk Hermittisk
fro_norm = norm(A)

println("A =\n", A)
println("\nA_adj =\n", A_adj)
println("\nHermittisk? ", is_hermitian)
println("\nFrobeniusnorm = ", fro_norm)

Eksempel i Python

import numpy as np

# Lage en kompleks matrise
A = np.array([[1+2j,  3-4j],
              [3+4j,  5+0j]])

A_adj = A.conj().T          # adjungert (konjugert transponert)
is_hermitian = np.allclose(A, A_adj)
fro_norm = np.linalg.norm(A, 'fro')

print("A =\n", A)
print("\nA† =\n", A_adj)
print("\nHermittisk?", is_hermitian)
print("\nFrobeniusnorm =", fro_norm)