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

\[A = U\,\Sigma\, V^T\]

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:

\[\begin{split} \begin{bmatrix} 1 & -2 \\ -2 & -2 \end{bmatrix} = \begin{bmatrix} -\frac{4\sqrt{5}}{15} & \frac{3\sqrt{5}}{10} \\ \frac{2\sqrt{5}}{15} & \frac{3\sqrt{5}}{5} \end{bmatrix} \begin{bmatrix} 3 & 0 \\ 0 & 2 \end{bmatrix} \begin{bmatrix} -\frac{2\sqrt{5}}{5} & \frac{\sqrt{5}}{5} \\ -\frac{\sqrt{5}}{5} & -\frac{2\sqrt{5}}{5} \end{bmatrix} \end{split}\]

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

\[\begin{split} \boldsymbol{x_1(t)} = \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix} \quad \text{og} \quad \boldsymbol{x_2(t)} = \begin{bmatrix} -\sin(t) \\ \cos(t) \end{bmatrix}\,. \end{split}\]
  1. Gør rede for at \(\,\boldsymbol{x_1(t)}\,\) og \(\,\boldsymbol{x_2(t)}\,\) er ortogonale for alle \(t.\)

  2. 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.

  3. 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)\,.\)

  4. 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_ellipse aktiveres.

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:

  1. Find med slideren et \(\,t_0,\,\) således at \(\,\boldsymbol{y_1(t_0)}\,\) og \(\,\boldsymbol{y_2(t_0)}\,\) er ortogonale.

  2. Sæt \(\,\boldsymbol{v_1}=\boldsymbol{x_1(t_0)}\,\) og \(\,\boldsymbol{v_2}=\boldsymbol{x_2(t_0)}\,\)

  3. 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\|\,\)

  4. 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)}\,\)

  5. 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.

  6. 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\,\)

  7. 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:

  1. Symmetrisk matrix, dvs. typen \(\,\left[\begin{array}{c}a & b\\b & c\end{array}\right]\)

  2. Skævsymmetrisk matrix af typen \(\,\left[\begin{array}{c}a & b\\-b & a\end{array}\right]\)

  3. 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:

LEGO-skib

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_colors erstatte white med 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 x med 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)
# 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 funktionen svd_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?