import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
Modul 2 Elever#
Introduktion til SVD#
Singular Value Decomposition (SVD) af reelle matricer baserer sig på følgende sætning:
Enhver reel \(m \times n\) matrix \(A\) kan skrives på SVD-formen
hvor \(V\) er en ortogonal \(n\times n\) matrix, \(\Sigma\) er en \(m\times n\) diagonal matrix, og \(U\) er en ortogonal \(m\times m\) matrix. Elementerne i \(\Sigma\) kaldes singulærværdierne, de er alle \(\geq 0\) og står i ikke-stigende rækkefølge.
Vi undersøger først grundigt \(2\times 2\) matricer med fuld rang. Bemærk af der ved SVD af en \(2\times 2\) matrix kun indgår \(2\times 2\) matricer. Her er et eksempel på SVD af en \(2\times 2\) matrix:
Opsætning#
Vi kalder de python pakker og procedurer vi får brug for nedenfor:
# Biblioteker vi bruger
import numpy as np
from numpy.linalg import svd, matrix_rank
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.patches import Circle
from matplotlib.image import imread
from PIL import Image
from ipywidgets import interact, FloatSlider, Layout
from IPython.display import clear_output
from scipy.optimize import brentq
# Visninger
np.set_printoptions(precision=3, suppress=True)
plt.rcParams.update({'figure.max_open_warning': 0})
Eksistens af løsning.#
Lad \(A\) være en reel 2x2 matrix med fuld rang. Vi skal vise et vigtigt resultat på vejen mod SVD: Der findes en ortonormal basis \(\,(\boldsymbol{v_1},\boldsymbol{v_2})\,\) for \(\Bbb R^2\) som opfylder at basisvektorernes billeder \(\,A\boldsymbol{v_1}\,\) og \(\,A\boldsymbol{v_2}\,\) er ortogonale.
For vilkårligt \(\,t\in\Bbb R\,\) betragter vi vektorerne
Gør rede for at \(\,\boldsymbol{x_1(t)}\,\) og \(\,\boldsymbol{x_2(t)}\,\) er ortogonale for alle \(t.\)
Vis at der findes et \(\,t_0\in \left[0,\frac {\pi}2\right]\) der gør billedvektorerne \(\, \boldsymbol{y_1(t_0)} = A\,\boldsymbol{x_1(t_0)}\,\) og \(\,\boldsymbol{y_2(t_0)} = A \,\boldsymbol{x_2(t_0)}\,\) ortogonale. Vink: Gør rede for at funktionen \(\,s(t)=\boldsymbol{y_1(t)}\cdot \boldsymbol{y_2(t)}\,\) er kontinuert, overvej dens værdi i intervallets endepunkter og brug Bolzanos sætning.
Med tallet \(t_0\) og funktionen \(\,s(t)\,\) fra forrige spørgsmål: Vis at \(\,s(t_0+ \frac{\pi}{2})=0\,.\)
Vink: Udnyt at der for alle \(t\) gælder \(\,\cos(t+\frac{\pi}{2})=-\sin(t)\,\) og \(\,\sin(t+\frac{\pi}{2})=\cos(t)\,.\)Forklar at vi hermed har løst opgaven.
Geometrisk løsning af SVD for 2x2 matrix \(\,A\,\) med fuld rang#
Vi opstiller her en procedure “enhedscirkel_ellipse” som visualiserer transformationen af enhedscirklen ved en given 2x2 matrix \(\,A\,:\)
def enhedscirkel_ellipse(A, slider_width='80%', figsize=(14, 7)):
"""
Opretter en interaktiv visualisering af enhedscirklen og dens billede under matrix A.
Parametre:
- A: 2x2 numpy array, transformationsmatrix
- slider_width: bredde på slider (f.eks. '80%')
- figsize: tuple med figurstørrelse (bredde, højde)
"""
# Basis- og billedfunktioner
def x1(t): return np.array([np.cos(t), np.sin(t)])
def x2(t): return np.array([-np.sin(t), np.cos(t)])
def y1(t): return A @ x1(t)
def y2(t): return A @ x2(t)
# Forudberegn ellipsepunkter
ts = np.linspace(0, 2*np.pi, 361)
ellipse_pts = np.column_stack([A @ x1(tt) for tt in ts]).T
# Plot-funktion
def plot_vectors(t):
clear_output(wait=True)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize)
for ax in (ax1, ax2):
ax.set_aspect('equal')
ax.grid(True, linestyle='--', linewidth=0.4, alpha=0.5)
# Fælles aksebegrænsninger
all_pts = np.vstack([ellipse_pts, x1(t), x2(t), y1(t), y2(t)])
lim = 1.15 * np.max(np.abs(all_pts))
for ax in (ax1, ax2):
ax.set_xlim(-lim, lim)
ax.set_ylim(-lim, lim)
# Venstre: enhedscirkel
circle = Circle((0, 0), radius=1, edgecolor='lightgrey', facecolor='none', linewidth=1.5, zorder=0)
ax1.add_patch(circle)
ax1.plot(np.cos(ts), np.sin(ts), color='lightgrey', linewidth=1.5, zorder=0)
# Højre: ellipse
ax2.plot(ellipse_pts[:, 0], ellipse_pts[:, 1], color='lightgrey', linewidth=1.5, zorder=0)
# Funktion til pile og labels
def draw_arrow(ax, v, color, label, z):
ax.quiver(0, 0, v[0], v[1], angles='xy', scale_units='xy', scale=1,
color=color, width=0.012, zorder=z)
ax.text(v[0], v[1], f' {label}{tuple(np.round(v, 2))}',
va='bottom', ha='left', fontsize=9, color=color)
# Tegn vektorer
draw_arrow(ax1, x1(t), 'tab:blue', 'x₁', 3)
draw_arrow(ax1, x2(t), 'tab:blue', 'x₂', 3)
draw_arrow(ax2, y1(t), 'tab:red', 'y₁', 4)
draw_arrow(ax2, y2(t), 'tab:red', 'y₂', 4)
ax1.set_title('Enhedscirklen S')
ax2.set_title('Ellipsen A.S')
fig.suptitle(rf"$t={t:.2f}$ (træk slideren)", fontsize=14)
plt.show()
# Slider
slider = FloatSlider(
value=0,
min=0,
max=2*np.pi,
step=0.01,
description='t',
continuous_update=True,
readout_format='.2f',
layout=Layout(width=slider_width)
)
interact(plot_vectors, t=slider)
I den følgende celle skal du indskrive den matrix du vil afprøve.
Prøv \(\,A=\left[\begin{array}{c}1 & 2\\-1 & 1\end{array}\right]\,\) for grafisk løsning med slider.
Kør cellen, så proceduren
enhedscirkel_ellipseaktiveres.
A = np.array([[1,2],
[-1,1]])
enhedscirkel_ellipse(A)
Vi tager udgangspunkt i situation og betegnelserne som i øvelse 1. Gå således frem:
Find med slideren et \(\,t_0,\,\) således at \(\,\boldsymbol{y_1(t_0)}\,\) og \(\,\boldsymbol{y_2(t_0)}\,\) er ortogonale.
Sæt \(\,\boldsymbol{v_1}=\boldsymbol{x_1(t_0)}\,\) og \(\,\boldsymbol{v_2}=\boldsymbol{x_2(t_0)}\,\)
Lad \(\sigma_1\) og \(\sigma_2\) være længderne \(\,\sigma_1=\left\|\boldsymbol{y_1(t_0)}\right\|\,\) og \(\,\sigma_2=\left\|\boldsymbol{y_2(t_0)}\right\|\,\)
Sæt \(\,\boldsymbol{u_1}=\frac{1}{\sigma_1}\,\boldsymbol{y_1(t_0)}\,\) og \(\,\boldsymbol{u_2}=\frac{1}{\sigma_2}\,\boldsymbol{y_2(t_0)}\,\)
Gør rede for at \(\,V=\left[\boldsymbol{v_1}\,\,\boldsymbol{v_2}\right]\,\) og \(\,U=\left[\boldsymbol{u_1}\,\,\boldsymbol{u_2}\right]\,\) er ortogonale matricer.
Lad \(\Sigma\) være diagonalmatricen \(\,\Sigma=\begin{bmatrix}\sigma_1&0\\0&\sigma_2\end{bmatrix},\) og gør rede for at \(\,A\,V=U\,\Sigma\,\)
Gør rede for og verificér i det givne eksempel, at \(\,A=U\,\Sigma\,V^T\,\)
Konvention: Det er fast konvention i SVD, at man ønsker singulærværdierne i faldende orden, det vil her sige, at \(\,\sigma_1 \geq \sigma_2\,.\)
Meget bekvemt fandt vi i første eksempel et \(\,t_0\in \left[0,\frac {\pi}2\right],\,\) hvor \(\,\left\|\boldsymbol{y_1(t_0)}\right\|>\left\|\boldsymbol{y_2(t_0)}\right\|\,.\)
Hvis det modsatte gælder, trækker man blot slideren videre til næste ortogonalitetspunkt. Det opstår ved \(\,t\)-vinklen \(\,t_0+\frac{\pi}2\,\) hvor længderne byttes om.
Prøv for eksempel \(\,A=\left[\begin{array}{c}1 & -1\\0 & 2\end{array}\right]\,\) eller \(\,A=\left[\begin{array}{c}1 & -1\\0 & -2\end{array}\right]\,.\)
Håndregningsversion hvor du har mulighed for at regne med eksakte tal:
Udfør SVD på \(\,A=\left[\begin{array}{c}5 & 3\\0 & 4\end{array}\right]\,.\)
Vi opstiller nu en procedure SVD_2x2_matrix som samler punkterne 1-7, og som du kan bruge til at teste dine 2x2 matricer:
Bemærk specielt at vi sørger for at \(\,\|\boldsymbol{y_1(t_0)}\|\,\) er større end eller lig med \(\,\|\boldsymbol{y_2(t_0)}\|\,.\)
def SVD_2x2_matrix(A, print_output = False):
# Basis vektorer som funktion af t
def x1(t): return np.array([np.cos(t), np.sin(t)])
def x2(t): return np.array([-np.sin(t), np.cos(t)])
def y1(t): return A.dot(x1(t))
def y2(t): return A.dot(x2(t))
# Dotprodukt-funktion
f = lambda t: np.dot(y1(t), y2(t))
# Find rod i (0, π/2)
t0 = brentq(f, 1e-6, np.pi/2 - 1e-6)
# Sørg for ||y1|| >= ||y2||
if np.linalg.norm(y1(t0)) < np.linalg.norm(y2(t0)):
t0 += np.pi/2
# Så kan vi sætte værdierne direkte:
sigma1 = np.linalg.norm(y1(t0))
sigma2 = np.linalg.norm(y2(t0))
Sigma = np.diag([sigma1, sigma2])
V = np.column_stack([x1(t0), x2(t0)])
U = np.column_stack([y1(t0)/sigma1, y2(t0)/sigma2])
# printer outputtet, hvis tilvalgt
if print_output:
print("t0 =", t0)
print("Σ =\n", Sigma)
print("V =\n", V)
print("U =\n", U)
return t0, U, Sigma, V
I den følgende celle afprøves proceduren SVD_2x2_matrix. Som output får du i denne rækkefølge: \(\,t_0,\, U, \,\Sigma\,\) og \(\,V\,.\)
Indtast din test matrix \(A\) og tjek i en følgende celle, at vi opnår den ønskede dekomposition \(\,A = U\, \Sigma\,V^T\,.\)
A = np.array([[1,2],
[-1,1]])
t0, U, Sigma, V = SVD_2x2_matrix(A, print_output=True)
Afprøvning af særlige matricer#
Find på simple test-eksempler på særlige 2x2 matricer:
Symmetrisk matrix, dvs. typen \(\,\left[\begin{array}{c}a & b\\b & c\end{array}\right]\)
Skævsymmetrisk matrix af typen \(\,\left[\begin{array}{c}a & b\\-b & a\end{array}\right]\)
Ortogonal matrix (du kan f.eks. bruge matrix \(U\) lagret fra foregående øvelse)
Test matricerne skrives i følgende celle. Når du kører cellen aktiverer du begge procedurer enhedscirkel_ellipse og SVD_2x2_matrix.
Eksperimenter via slideren og studér det algebraiske output på matricer af type 1-3. Diskutér hvad der er det særlige!
A = np.array([[1,2],
[-1,1]])
enhedscirkel_ellipse(A)
_ = SVD_2x2_matrix(A, print_output=True)
SVD kommando for en \(m\times n\) matrix \(A\)#
En \(m\times n\) matrix med tilfældige heltal mellem eksempelvis \(1\) og \(5\) som indgange kan laves med np.random.randint().
Vælg \(m\) og \(n\) og dan en tilfældigt generet matrix
m = 'INDSÆT KODE HER'
n = 'INDSÆT KODE HER'
np.random.seed(0) # fastlåser "tilfældigheden"
A = np.random.randint(1,6,size=(m,n))
A
SVD-kommandoen np.linalg.svd(A) for en matrix \(A\) giver (\(m\times m\))-matricen \(U\), \(\,\) (\(n\times n\))-matricen \(V^T\) og en vektor \(\boldsymbol{s}\) med singulærværdier i ikke-stigende orden.
Benyt SVD-kommandoen til at fremkalde \(U\), \(V^T\) og vektoren \(\boldsymbol{s}\) for matricen \(A\).
# Beregn SVD
U, s, Vt = 'INDSÆT KODE HER'
print('\nSingulærværdier (s):', s)
print('\nU =\n', U)
print('\nVt =\n', Vt)
Dan \(\Sigma\) matricen ud fra singulærværdierne i \(\boldsymbol{s}\) ved at færdiggøre koden nedenfor.
# Vi opretter en nulmatrix med samme dimensioner som A (m x n).
Sigma = np.zeros((m, n))
# Indæst de singulære værdier (s) på diagonalen.
for i in range(min(m, n)):
Sigma[i, i] = 'INDSÆT KODE HER'
print('\nSigma =\n', Sigma)
Vis at SVD-kommadoen passer ved at danne produktet \(U\Sigma V^T\).
# Husk at bruge @ til at lave matrixmultiplikation
recon = 'INDSÆT KODE HER'
print('\nRekonstruktion (U @ Sigma @ Vt) =\n', np.round(recon,3))
print('\nA =\n', A)
LEGO-skibet#
Vi tager udgangspunkt i et billede af et LEGO-skib fundet på internet.
Undervejs kommer vi til at:
Undersøge SVD og rang-\(k\) approksimation.
Implementere og visualisere \(k\)-rangs-rekonstruktion af billeder i Python.
Bruge SVD til billedkompression og tolke singulærværdier.
Vi betragter følgende billede fra internettet:
Modellering af LEGO-enheder#
Vi kan definere en LEGO-enhed som et lille rektangel, der indeholder to cylindre, set fra oven.
Det er nemt at se, at der er \(11\) enheder lodret.
Man kan også tælle, at der er \(15\) enheder vandret.
Nu skal vi vælge farver til 2D-rekonstruktionen af LEGO-figuren. Vi kan bruge et udvalg af farver fra CSS Colors: Matplotlib CSS Colors
Når hver LEGO-enhed og baggrunden får tildelt et tal svarende til dens farve, opstår der en \(15\times11\) matrix \(A\), som repræsenterer LEGO-figuren i 2D.
Vi antager, at du får brug for en farve til baggrund og ca. syv forskellige farver til LEGO-figuren.
Gå frem således:
I den følgende celle skal du i listen
base_colorserstattewhitemed de farver du har valgt. Derved tildeles farverne et tal fra \(1\) til \(8\), svarende til deres plads i listen. Den første farve du vælger, svarer til tallet \(1\), det vil sige baggrundsfarven.Derefter skal du i matricen \(A\) erstatte hver
xmed det tal, som du har valgt til LEGO-enheden, så enheden får den farve, som du ønsker.Når du er færdig, skal du køre cellen og de følgende celler for at se din LEGO-figur i 2D.
base_colors = ['white','white','white','white','white','white','white','white']
A = np.array([
[1,1,1,1,1,1,x,x,1,1,1,1,1,1,1],
[1,1,1,1,1,x,x,x,1,1,1,1,1,1,1],
[1,1,1,1,x,x,x,x,1,1,1,1,1,1,1],
[1,1,1,x,x,x,x,x,1,1,1,1,1,1,1],
[1,1,x,x,x,x,x,x,1,1,1,1,x,1,1],
[1,x,x,x,x,x,x,x,1,1,x,x,x,x,1],
[x,x,x,x,x,x,x,x,1,1,1,x,x,x,1],
[1,1,1,1,1,1,x,x,1,1,1,x,x,x,1],
[1,1,1,x,x,x,x,x,x,x,x,x,x,x,x],
[1,1,1,1,x,x,x,x,x,x,x,x,x,x,1],
[1,1,1,1,1,x,x,x,x,x,x,x,x,1,1]
], dtype=int)
Hint
Forslag til farvelægning.
# Knyt en farve til tallene + fejl-farver i hver ende
base_colors = ['palegreen','white','red','forestgreen','blue','yellow','black','cyan']
# Hver Lego-enhed tildeles et tal
A = np.array([
[1,1,1,1,1,1,3,3,1,1,1,1,1,1,1],
[1,1,1,1,1,4,4,4,1,1,1,1,1,1,1],
[1,1,1,1,3,3,3,3,1,1,1,1,1,1,1],
[1,1,1,4,4,4,4,4,1,1,1,1,1,1,1],
[1,1,3,3,3,3,3,3,1,1,1,1,2,1,1],
[1,4,4,4,4,4,4,4,1,1,7,7,7,7,1],
[3,3,3,3,3,3,3,3,1,1,1,8,8,6,1],
[1,1,1,1,1,1,5,5,1,1,1,6,6,6,1],
[1,1,1,2,2,2,2,2,2,2,2,2,2,2,2],
[1,1,1,1,2,2,2,2,2,2,2,2,2,2,1],
[1,1,1,1,1,2,2,2,2,2,2,2,2,1,1]
], dtype=int)
# KØR DENNE CELLE FOR AT SE DIN LEGO-FIGUR
extended_colors = ['magenta'] + base_colors + ['magenta']
cmap = ListedColormap(extended_colors)
# Definér grænser så værdier <0.5 og >8.5 bliver vist som magenta
bounds = np.arange(-0.5, len(extended_colors) - 0.5 + 1, 1)
norm = BoundaryNorm(bounds, cmap.N)
# 2) Figur og akse
fig, ax = plt.subplots(figsize=(6,5))
# 3) Tegn billedet – først billedet, så grid (så grid’en ligger øverst)
im = ax.imshow(A, cmap=cmap, norm=norm,
extent=[-0.5, A.shape[1]-0.5, -0.5, A.shape[0]-0.5])
# 4) Lav grid med meshgrid + plot
# x går fra -0.5 til 14.5 (15 kolonner), y fra -0.5 til 10.5 (11 rækker)
x = np.arange(-0.5, A.shape[1] + 0.5, 1)
y = np.arange(-0.5, A.shape[0] + 0.5, 1)
X, Y = np.meshgrid(x, y)
ax.plot(X, Y, color='lightgrey', linewidth=0.7) # vertikale
ax.plot(X.T, Y.T, color='lightgrey', linewidth=0.7) # horisontale
# 5) Fjern ticks, sæt ensartet aspect, og giv en sort ramme
ax.set_xticks([]); ax.set_yticks([])
ax.set_aspect('equal')
# 6) Titel og show
ax.set_title("Originalt billede af LEGO skib")
plt.tight_layout()
plt.show()
Rang-\(k\) matrix approksimation#
For at komprimere vores LEGO-billede kan vi bruge SVD på matricen \(A\):
U, s, Vt = np.linalg.svd(A.astype(float))
Ved at plotte singulærværdierne får vi et klart indtryk af, hvor meget hvert enkelt led i summen i \((*)\) bidrager til den samlede matrix.
# Plot singulærværdier
plt.figure(figsize=(8,5))
plt.scatter(np.arange(1, len(s)+1), s, color='blue')
plt.title("Singulærværdiers størrelse for A")
plt.xlabel("Index")
plt.ylabel("Singulærværdi størrelse")
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()
De største singulærværdier repræsenterer de mest betydningsfulde mønstre i LEGO-billedet, mens de små bidrager med mindre detaljer. Plottet hjælper os derfor med at beslutte, hvor mange singulærværdier vi kan beholde, uden at miste for meget information, hvilket er centralt, når vi laver en \(k\)-rang-approksimation af billedet.
Vi ønsker at visualisere det originale LEGO-billede og rang-\(k\) SVD-approksimationen side om side for at se effekten af approksimationen. Funktion svd_approksimation(A,k) nedenfor laver en rang-\(k\) approksimation af billedet.
def svd_approksimation(A, k):
"""
Rank-k SVD-approximation
"""
U, s, Vt = np.linalg.svd(A)
A_k = np.zeros_like(A, dtype=float) # nulmatrix med samme dimensioner som A
for i in range(k):
A_k += s[i] * np.outer(U[:, i], Vt[i, :]) # sum af ydreprodukter givet antal af medtagne singulærværdier
return A_k
Færdiggør funktionen
plot_svd(k)ved at indsætte funktionensvd_approksimation().
def plot_svd(k):
# Beregn k-rangs approksimation af A (kaldes A_recon)
A_recon = 'INDSÆT KODE HER'
fig, axes = plt.subplots(1, 2, figsize=(12, 5), dpi=100)
# Lav koordinatarrays til pcolormesh - definerer kanterne på hver pixel
x_edges = np.arange(-0.5, A.shape[1] + 0.5, 1)
y_edges = np.arange(-0.5, A.shape[0] + 0.5, 1)
# Lav meshgrid for kanterne
X, Y = np.meshgrid(x_edges, y_edges)
# Approksimation
ax = axes[0]
# Brug pcolormesh som naturligt aligner med griddet
mesh = ax.pcolormesh(X, Y, A_recon, cmap=cmap, norm=norm, shading='flat')
# Tilføj gridlinjer præcist på kanterne
ax.vlines(x_edges, y_edges[0], y_edges[-1], colors='lightgrey', linewidth=0.7)
ax.hlines(y_edges, x_edges[0], x_edges[-1], colors='lightgrey', linewidth=0.7)
ax.set_xlim(x_edges[0], x_edges[-1])
ax.set_ylim(y_edges[0], y_edges[-1])
ax.invert_yaxis() # Retter op-ned vendt billede
ax.set_aspect('equal')
ax.set_xticks([])
ax.set_yticks([])
ax.set_title(f"Approksimation med k={k} singulærværdier")
# Originalt
ax = axes[1]
mesh = ax.pcolormesh(X, Y, A, cmap=cmap, norm=norm, shading='flat')
ax.vlines(x_edges, y_edges[0], y_edges[-1], colors='lightgrey', linewidth=0.7)
ax.hlines(y_edges, x_edges[0], x_edges[-1], colors='lightgrey', linewidth=0.7)
ax.set_xlim(x_edges[0], x_edges[-1])
ax.set_ylim(y_edges[0], y_edges[-1])
ax.invert_yaxis() # Retter op-ned vendt billede
ax.set_aspect('equal')
ax.set_xticks([])
ax.set_yticks([])
ax.set_title("Originalt billede")
plt.tight_layout()
plt.show()
Prøv forskellige værdier af \(k\) med funktionen
plot_svd(k)og se, hvordan billedet rekonstrueres. Hvilken effekt har det på billedets kvalitet?
# INDSÆT KODE HER
Prøv at rekonstruere matrixen ved at medtage alle singulærværdier. Hvorfor ser det præcis ud som det originale LEGO-billede?