Institutionen för naturvetenskap och teknik
Att avsudda bilder
Filtrering och iterativa metoder
Örebro universitet
Institutionen för naturvetenskap och teknik
Självständigt arbete för kandidatexamen i matematik, 15 hp
Att avsudda bilder
Filtrering och iterativa metoder
Mattis Kienmayer Maj 2019
Handledare: Mårten Gulliksson Examinator: Andrii Dmytryshyn
Sammanfattning
Denna uppsats rör vid problemet med suddiga bilder och vad man kan göra för att försöka återställa den eftersökta skarpa bilden. Problemet tacklas med hjälp av linjär-algebraiska medel, både regularisering och iterativa metoder. Som resultat visar sig DFPM (dynamical functional particle method) jäm-förbar med både konjugerad gradient-metoden och LSQR-algoritmen (least squares QR), dock med användning av andra parametrar än de teoretiskt optimala. Utöver tidigare kända metoder presenteras och testas även ett al-ternativt minimeringsproblem.
Innehåll
1 Inledning 5 2 Modellering 6 2.1 En enkel modell . . . 6 2.2 Suddmatrisen . . . 6 2.3 Separabelt sudd . . . 73 Ett försök till avsuddning 9 3.1 Singulärvärdesfaktorisering . . . 9 3.2 Störningsanalys . . . 9 3.3 Regularisering . . . 10 4 Iterativa metoder 13 4.1 DFPM . . . 13 4.2 Symplektisk Euler . . . 15 4.3 Konvergenstest . . . 18 4.3.1 Parameteranalys - DFPM . . . 19 4.4 En alternativ formulering . . . 21 A Matlabkoder 23
Kapitel 1
Inledning
När man tar en bild och ljuset ska färdas hela vägen till kamerans sensor finns det olika störningsmoment som påverkar ljuset på vägen. Det mest uppenbara som påverkar ljusets färd är linsens utformning, men även saker som temperaturskillnader i luften. Om vi kallar den sanna scenen framför oss som vi försöker föreviga i form av en bild för skarp så är det i princip omöjligt att inte ta en suddig (eng. blurred) bild, så avsuddning (eng. deblurring) är viktigt om man vill ha skarpa bilder. Tyvärr finns det på grund av olika fel i processen ingen garanti att man kan återskapa den skarpa bilden, man får därav ofta nöja sig med något defekt. Att man beskriver verkligheten med ett ändligt antal färger bidrar bland annat till detta. I exempelvis [4] eller [10] kan man läsa om olika modeller.
För att tackla problemet med avsuddning används i den här uppsatsen förutom ett par kända filtreringsmetoder även DFPM (dynamical functio-nal particle method), en iterativ metod presenterad i bland annat [1] och [2]. För att testa DFPM jämförs metoden med två väletablerade iterativa metoder, LSQR (least squares QR) och konjugerad gradient, på några av-suddningsproblem. Alla bilder som behandlas i uppsatsen är data bifogade till [4].
Kapitel 2
Modellering
2.1
En enkel modell
En digital bild består av ett antal pixlar, alla med någon ljusintensitet. En svart-vit bild är därför enkel att representera i form av en matris där varje element representerar en pixels ljusintensitet på gråskalan, med 8-bitars färg kan dessa element anta heltalsvärden mellan 0 och 255. Vi betecknar den sanna bilden för X och den suddiga för B. För att kunna använda oss utav linjär-algebraiska metoder antar vi att suddprocessen är linjär, och så arran-gerar vi om matriserna X och B så de blir vektorer x och b. Låt X, B ∈ Rn×m och definiera x = vec(X) = x1 .. . xm ∈ RN, b = vec(B) = b1 .. . bm ∈ RN, där xj ∈ Rn och b
j ∈ Rn betecknar kolonnvektor j i X respektive B, och
N = nm. Eftersom vi antar att suddprocessen är linjär kan vi nu säga att det finns en matris A ∈ RN ×N sådan att Ax = b. Något som ofta förekommer vid fotografi är mätfel och brus, därför är det naturligt att utvidga modellen så att den tar hänsyn till detta. Om vi låter e beteckna felen som kan uppkomma får vi modellen
Ax = b + e.
2.2
Suddmatrisen
Frågan som uppstår näst är hur matrisen A ser ut. A bestäms av en så kallad ’point spread function’ (PSF) P ∈ Rn×m som beskriver hur ljuset som når sensorn sprider sig. För att få information om hur P ser ut kan man ta en bild B på en isolerad ljuskälla i storleksordning av en pixel och få B ≈ P . Ett par egenskaper som kännetecknar en PSF är att dess element summeras
till 1 samt att alla element pij ≥ 0. Antag att vi har en bild X och en PSF P , X = x11 x12 x13 x21 x22 x23 x31 x32 x33 , P = p11 p12 p13 p21 p22 p23 p31 p32 p33 .
Att ljus sprids åt något håll kan även ses som att ljus tas från motsatt håll. För att bestämma element bij i bilden B roteras därför P 180 grader, varefter
centrum av P , p22 i vårat fall, placeras över element xij. Slutligen bestäms
bij genom summation av parvis multiplicerade element. Vi får
bij = p33· xi−1,j−1+ p32· xi−1,j + p31· xi−1,j+1
+ p23· xi,j−1 + p22· xij + p21· xi,j+1
+ p13· xi+1,j−1 + p12· xi+1,j + p11· xi+1,j+1.
Värdet på de element xij som ligger utanför bilden, d.v.s. de xij där i /∈ [1, n] eller j /∈ [1, m], bestäms av ens randvillkor. Om vi antar noll randvärdesvill-kor, d.v.s. att allt utanför bilden är 0, implicerar detta att
p22 p12 p21 p11 p32 p22 p12 p31 p21 p11 p32 p22 p31 p21 p23 p13 p22 p12 p21 p11 p33 p23 p13 p32 p22 p12 p31 p21 p11 p33 p23 p32 p22 p31 p21 p23 p13 p22 p12 p33 p23 p13 p32 p22 p12 p33 p23 p32 p22 x11 x21 x31 x12 x22 x32 x13 x23 x33 = b11 b21 b31 b12 b22 b32 b13 b23 b33 , (2.1)
eller på enklare form
Avec(X) = vec(B) ⇔ Ax = b.
Matrisen A är av typen block Toeplitz med Toeplitz block (BTTB), och kännetecknas av konstanta blockdiagonaler (i vårat fall är blocken 3 × 3) som i sig har konstanta diagonaler.
Något som antogs när vi tog fram A är att suddprocessen inte är ett lokalt fenomen, d.v.s. att sättet ljuset sprider sig på inte beror på dess position i bilden. I syfte av att begränsa vad uppsatsen rör vid för problem kommer detta att antas framöver.
2.3
Separabelt sudd
Om vi antar att rader och kolumner suddas oberoende av varandra, även kallat separabelt sudd, behövs en annan tackling av problemet. I linjäralge-braiska termer kan separabelt sudd uttryckas som
där Ac ∈ Rn×n, A
r ∈ Rm×m. För att kunna använda vår modell Ax = b
nyttjar vi Kroneckerprodukten, definierad som
Ar⊗ Ac= a(r)1,1Ac · · · a(r)1,mAc .. . . .. ... a(r)m,1Ac · · · a(r)m,mAc ,
där a(r)ij betecknar element tillhörande matrisen Ar. Vi har att
AcXATr = B ⇐⇒ (Ar⊗ Ac)x = b.
Kapitel 3
Ett försök till avsuddning
3.1
Singulärvärdesfaktorisering
Ett mycket användbart verktyg vid avsuddning är singulärvärdesfaktorise-ringen, eller SVD (singular value decomposition). Singulärvärdesfaktorise-ringen uttrycker en matris A ∈ RN ×N som
A = U ΣVT, där UTU = I, VTV = I, Σ = diag(σ
1, . . . , σN), σ1 ≥ . . . ≥ σr > σr+1 =
· · · = σN = 0 och antalet nollskiljda singulärvärden σi är matrisens rank,
d.v.s. r = rank(A). Kolumnerna ui och vi, i = 1, . . . , N i matriserna U och V kallas för vänstra respektive högra singulärvektorer. En följd av denna faktorisering är att man kan skriva A som en summa av rank-1-matriser. Vi har att A = U ΣVT= r X i=1 σiuivTi .
Under antagande att A är inverterbar, d.v.s. r = N , fås
A−1= V Σ−1UT = N X i=1 1 σi viuTi .
Mer om singulärvärdesfaktoriseringen hittas i exempelvis [11].
3.2
Störningsanalys
Problemet med avsuddning är ett så kallat illa ställt problem. Vad som menas med detta är att lösningarna till två närliggande problem kan ligga väldigt långt ifrån varandra, detta leder till att små mätfel kan leda till stora fel i den eftersökta lösningen. Speciellt kännetecknas suddmatriserna av stora konditionstal.
20 40 60 80 100 10 20 30 40 50 60 70 80 90 100 20 40 60 80 100 10 20 30 40 50 60 70 80 90 100 20 40 60 80 100 10 20 30 40 50 60 70 80 90 100
Figur 3.1: Originalbilden till vänster, suddig och brusig bild i mitten, där bruset e är taget från en likformig fördelning med kekF = 0.01, och en naiv
lösning x = A−1(b + e) till höger.
Betrakta det störda problemet som är vår modell Ax = b + e,
som ger upphov till en störd lösning
x = ˆx + δx = A−1(b + e),
där ˆx är lösningen till Ax = b, vilket ger Aδx = e. Om vi nyttjar singulär-värdesfaktoriseringen kan vi uttrycka lösningen x som
x = A−1(b + e) = A−1b + A−1e = N X i=1 uT i b σi vi+ N X i=1 uT i e σi vi.
Eftersom det för alla i gäller att kvik = 1 medför detta att kxk ≤ N X i=1 1 σi |uTi b| + N X i=1 1 σi |uTi e|, där k·k betecknar den Euklidiska normen kxk2 = PN
i=1x2i. Vid inspektion
av den andra termen i högerled inses att små störningar e förstoras av små singulärvärden σi, något man försöker hantera med hjälp av olika regulari-seringsmetoder. I figur 3.1 visas hur en naiv lösning x = A−1(b + e) kan se ut.
3.3
Regularisering
Om vi inför ett filter φi, 0 ≤ φi ≤ 1, som styr hur mycket varje singulärvärde
σi bidrar till lösningen erhålls en så kallad spektral filtreringsmetod, och
lösningen till ett problem Ax = b ser då ut som xfilt= N X i=1 φi uTi (b + e) σi vi.
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 10-20 10-15 10-10 10-5 100 105
Figur 3.2: Picard plot för problemet i figur 3.1. Orangea linjen är singulär-värden σi, och den blå linjen med punkter är högerledskoefficienter |uTi b|.
20 40 60 80 100 10 20 30 40 50 60 70 80 90 100 Figur 3.3: TSVD lösning med k = 750 (N = 10201) till problemet i figur 3.1. Det finns många sätt att välja
filterfakto-rer φi på. En enkel metod är TSVD (trunkerad
SVD), då definerar vi filterfaktorerna som φi=
(
1, i = 1, . . . , k 0, i = k + 1, . . . , N.
I figur 3.3 demonstreras en TSVD lösning till samma problem som i figur 3.1. Valet av k motiveras genom inspektion av den så kallade Picard-ploten, figur 3.2, där singulärvärden σi
jämförs med absolutvärdet av högerledskoeffi-cienter uTi (b + e) av anledningen att de ’rätta’ koefficienterna uTi b tenderar att minska, se
ex-empelvis [3]. Efter ungefär i = 750 ser vi att koefficienterna planar av. Vi drar slutsatsen att koefficienter uTi (b + e), i > 750, domineras av brus och
att vi inte vill ha med dessa i vår lösning.
Tikhonovregularisering är en annan spektral filtreringsmetod, då define-ras filterfaktorerna som
φi= σ2 i σ2 i + α2 , i = 1, . . . , N,
där α > 0 är regulariseringsparametern. Hansen visar i [4] att dessa filter-faktorer satisfierar φi= 1 −α σi 2 + Oα σi 4 , σi α σi α 2 + O σi α 4 , α σi
samt att α bör ligga i intervallet [σN, σ1]. Detta medför att φi ≈ 1 för små
i, och φi ≈ 0 för stora i. TSVD och Tikhonovregularisering är därav relativt
lika. Enligt [10] erhålls Tikhonovlösningen med regulariseringsparameter α vid lösning av problemet
min
x kAx − bk 2
Kapitel 4
Iterativa metoder
Att ta fram en singulärvärdesfaktorisering kan ofta kräva alltför många be-räkningar. Man kan då använda iterativa metoder som i sitt stoppkriterie innefattar någon slags regularisering.
I det här kapitlet jämförs DFPM, en iterativ metod presenterad i [1], med de kändare konjugerad gradient metoden och LSQR ([7]) algoritmen på några avsuddningsproblem (läs mer om dessa metoder i exempelvis [8]). Vad som jämförs är hur lång tid och hur många iterationer metoderna kräver för att konvergera men utöver det även bildkvalitét. För att mäta bildkvalitét objektivt används måttet MSE (mean square error), definierat med den rätta bilden ˆx som referens enligt
MSE(ˆx, x) = 1
Nkˆx − xk
2 2.
4.1
DFPM
Tanken med DFPM (dynamical functional particle method, se [1]) är att lösa problemet Ax = b genom att introducera ett dynamiskt system av andra ordningens differentialekvationer. Betrakta ekvationen
¨
x + η ˙x = b − Ax, x(0) = x0, x(0) = ˙˙ x0, (4.1)
där η > 0 och ˙x betecknar tidsderivatan av x(t). Om A har positiva egen-värden följer att en lösning x(t) till 4.1 närmar sig lösningen till ekvationen Ax = b då t växer. Vi tittar närmre på detta i sats 4.1.1, men först behövs en lite striktare definition.
Definition 4.1.1. Vi kallar ett dynamiskt system ˙x = f (t, x) för asymp-totiskt stabilt i ˆx om f (t, ˆx) = 0 och om det finns något δ > 0 sådant att
kˆx − x(t0)k < δ ⇒ lim
Sats 4.1.1. Antag att A ∈ RN ×N är symmetrisk med alla egenvärden λ(A) > 0, då gäller att det dynamiska systemet 4.1 är asymptotiskt stabilt i lösningen till Ax = b.
Bevis. Antag att ˆx är lösningen till Ax = b och låt y = x − ˆx. Då kan vi, med hjälp av likheten −Ay = b − Ax, skriva om 4.1 som
¨
y + η ˙y = −Ay, η > 0.
Då matrisen A antas vara symmetrisk kan vi diagonalisera A med en ortogo-nal matris P , d.v.s. A = P DPT, där PTP = I samt D = diag(λ1, . . . , λN),
för att få ¨ y + η ˙y = −P DPTy ⇐⇒ PTy + ηP¨ Ty = −DP˙ Ty. Låt nu z = PTy, då erhålls ¨ z + η ˙z = −Dz,
som är ett system av oberoende linjära ODE vars enskilda ekvationer ges av ¨
zi+ η ˙zi= −λizi, i = 1, . . . , N, och som har lösningarna
zi(t) = C1(i)er
(i)
1 t+ C(i)
2 e
r(i)2 t, i = 1, . . . , N, (4.2)
där konstanterna C1(i) och C2(i) beror på begynnelsevärden x(0) samt ˙x(0) och r(i)1 , r(i)2 ges av den karakteristiska ekvationen
r2+ ηr + λi = 0 ⇐⇒ r(1,2)= − η 2 ± r η2 4 − λi.
Nu finns det två fall att studera, då r är reellt och då r är icke-reellt. Om r är reellt, d.v.s. η2/4 ≥ λi, medför detta att r(1,2) < 0, vilket implicerar
att zi(t) → 0 då t → ∞, och eftersom z = PTy = PT(x − ˆx) följer även
att x → ˆx. Antag nu att r är icke-reellt. Sätt r = a ± ib, där a = −η/2, b =pλi− η2/4, då kan vi skriva lösningen 4.2 som
zi(t) = D(i)1 eatcos bt + D2(i)eatsin bt,
där D1 = C1+ C2 och D2= i(C1− C2). Då a < 0 har vi att
lim
t→∞zi(t) = 0 ⇒ t→∞lim x(t) = ˆx.
Notera att vi i beviset inte specifierat något δ. Detta gäller alltså för alla lösningar x(t) till ekvation 4.1, oavsett kˆx − xk. Denna egenskap kallas ibland för global asymptotisk stabilitet. Systemet är globalt asymptotiskt stabilt i lösningen till Ax = b även för ickesymmetriska diagonaliserbara matriser A med positiva egenvärden, men beviset för detta täcks inte i den här uppsatsen, se [1].
4.2
Symplektisk Euler
Symplektisk Euler är en semi-implicit metod ämnad att lösa system av första ordningens differentialekvationer. Applicerad på ett dämpat system
˙ x = f (t, y) − ηx, x(t0) = x0 ˙ y = g(t, x), y(t0) = y0, (4.3)
ges iterationerna utav
xk+1 = xk+ f (tk, yk) − ηxk ∆t
yk+1= yk+ g(tk, xk+1)∆t,
där x, y ∈ RN, ∆t betecknar tidssteget och tk = t0 + k∆t, k = 0, 1, 2, . . .
Skillnaden mellan symplektisk Euler och klassisk Euler är att vi för att be-stämma yk+1 använder oss utav xk+1 istället för xk. Det är av den
anled-ningen metoden ibland kallas för semi-implicit. Se exempelvis [6] för mer om symplektisk Euler applicerat på dämpade dynamiska system.
För att applicera symplektisk Euler på ekvation 4.1 måste det först trans-formeras om till formen given av 4.3. Vi introducerar därför en term v = ˙x så att vi får systemet
(
˙v = b − Ax − ηv ˙
x = v. Iterationerna ser slutligen ut som
vk+1= vk+ (b − Axk− ηvk)∆t xk+1 = xk+ vk+1∆t. (4.4)
Optimala parametrar η och ∆t ges enligt [1] av ηopt = 2 √ λminλmax/ √ λmin+ √ λmax ∆topt = 2/ √ λmin+ √ λmax, (4.5)
där λmin och λmax är minsta respektive största egenvärdet tillhörande
ma-trisen A.
Definition 4.2.1. Vi kallar en iterativ metod xk+1= Cxk, där xk= Ckx0,
för stabil om det finns någon konstant K > 0 så att kxkk = kCkx
0k ≤ Kkxk0, k = 1, 2, 3, . . .
Lemma 4.2.1. Ett nödvändigt krav för stabilitet är att egenvärdena till ma-trisen C i definition 4.2.1 alla är mindre än eller lika med 1 i absolutvärde.
Kravet som beskrivs i lemma 4.2.1 brukar även kallas för Von Neumann stabilitet, och ett bevis för dess nödvändighet finns i exempelvis [5].
Sats 4.2.1. Antag att A ∈ RN ×N är symmetrisk med alla egenvärden λ(A) > 0, samt att parametrar η och ∆t ges av ekvation 4.5, då gäller att symplektisk Euler applicerat på problemet 4.1 är stabilt.
Bevis. Vi börjar med en transformation av systemet, precis som i beviset till sats 4.1.1, så att det ser ut som
¨
z + η ˙z = −Dz,
där D = diag(λ1, . . . , λN). Låt v = ˙z så att vi får systemet
(
˙v = −Dz − ηv ˙
z = v. (4.6)
Med symplektisk Euler applicerad på 4.6 ges iterationerna utav vk+1= vk+ ∆t(−Dzk− ηvk) zk+1= zk+ ∆tvk+1, eller på matrisform " vk+1 zk+1 # = " (1 − η∆t)I −D∆t (1 − η∆t)∆tI I − D(∆t)2 # " vk zk # .
Det är alltså egenvärdena tillhörande matrisen C ∈ R2N ×2N vi söker, där C = " (1 − η∆t)I −D∆t (1 − η∆t)∆tI I − D(∆t)2 # .
Speciellt ska vi visa att alla egenvärden |γi| ≤ 1, i = 1, . . . , 2N . Egenvärden
till C fås av lösningarna till ekvationen det(C − γI) = 0 ⇐⇒ det γ2I + (η∆t − 2)I + (∆t)2D γ + (1 − η∆t)I = 0. Eftersom B = γ2I +((η∆t−2)I +(∆t)2D)γ +(1−η∆t)I är en diagonalmatris följer att det(B) = 0 om och endast om något diagonalelement är 0. Vi har alltså N andragradsekvationer med två lösningar vardera, som ger oss 2N eftersökta egenvärden. Ekvationerna ges av
Vi har att Pi(γ) = 0 om och endast om γ(1,2) = 2 − η∆t − (∆t)2λi 2 ± r (2 − η∆t − (∆t)2λ i)2 4 − (1 − η∆t) ⇐⇒ γ(1,2)= 2 − η∆t − (∆t) 2λ i 2 ± ∆t 2 p (∆tλi+ η)2− 4λi.
Nu finns tre fall att undersöka, nämligen fallen då (∆tλi+ η)2 är större än, mindre än, eller lika med 4λi. Det första fallet, (∆tλi+η)2> 4λi, uppkommer
aldrig om vi använder parametrar η, ∆t givna av 4.5 eftersom
(∆tλi+η)2−4λi= 2λi+ 2 √ λminλmax √ λmin+ √ λmax !2 −4λi = 4(λi− λmin)(λi− λmax) (√λmin+ √ λmax)2 . Detta kan aldrig vara större än 0 eftersom alla egenvärden antas vara po-sitiva samt att λi ∈ [λmin, λmax] . Antag nu att (∆tλi + η)2 < 4λi, då
kan vi skriva γ(1,2) = x ± iy, där x = (2 − η∆t − (∆t)2λi)/2 och y =
(∆t/2)p4λi− (∆tλi+ η)2. Vi har att |γ| ≤ 1 ⇐⇒ |γ|2≤ 1, och
|γ(1,2)|2 = x2+ y2 = 2 − η∆t − (∆t) 2λ i 2 !2 +(∆t) 2 4λ i− (∆tλi+ η)2 4 = 1 − η∆t, som med insättning av parametrar η, ∆t givna av ekvation 4.5 ger
1 − η∆t = λmax+ λmin− 2 √ λmaxλmin λmax+ λmin+ 2 √ λmaxλmin < 1.
För fallet då (∆tλi+ η)2 = 4λi använder vi likheten ∆tλi+ η = 2
√ λi för att få γ(1,2)= 2 − η∆t − (∆t)2λi 2 = 2 − ∆t(∆tλi+ η) 2 = 1 − ∆t p λi.
Med insättning av ∆t given av ekvation 4.5 erhålls γ(1,2)= √ λmin+ √ λmax− 2 √ λi √ λmin+ √ λmax , och eftersom λi∈ [λmin, λmax] följer även här att |γ(1,2)| < 1.
Tabell 4.1: Konvergenstest för symmetrisk suddmatris A, tiden t är mätt i sekunder och tol = 10−4.
DFPM
Problem Iterationer Median(t) Varians(t) MSE
1 102 0.939 3.96 · 10−4 244.73 2 62 0.576 3.34 · 10−4 143.70 3 85 0.780 5.73 · 10−4 185.72 Konjugerad gradient 1 47 0.220 2.30 · 10−5 157.52 2 32 0.153 4.99 · 10−5 113.52 3 39 0.184 7.91 · 10−5 147.63 LSQR 1 181 2.292 2.05 · 10−3 278.82 2 131 1.670 1.09 · 10−3 157.59 3 151 1.900 1.14 · 10−3 205.01
Tabell 4.2: Konvergenstest för icke-symmetrisk suddmatris A, tiden t är mätt i sekunder och tol = 10−3.
DFPM anpassad för normalekvationerna
Problem Iterationer Median(t) Varians(t) MSE
1 112 0.645 4.79 · 10−3 84.35 2 98 0.557 1.46 · 10−3 79.68 3 102 0.584 2.69 · 10−4 110.27 BICG 1 56 0.324 2.15 · 10−4 4.23 2 43 0.250 6.66 · 10−5 14.15 3 54 0.314 1.34 · 10−4 14.89 LSQR 1 58 0.499 3.20 · 10−4 85.73 2 37 0.320 2.96 · 10−4 11.38 3 48 0.414 1.82 · 10−4 125.37
4.3
Konvergenstest
I tabell 4.1 samt tabell 4.2 finns resultat från konvergenstest utförda på 3 avsuddningsproblem med 2 olika suddmatriser, en symmetrisk positiv definit och en icke-symmetrisk men båda av samma storlek N × N , N = 2601. Den symmetriska suddmatrisen är framtagen från en Gaussisk PSF med noll rand-värdesvillkor, och den ickesymmetriska kommer från separabelt sudd. Bruset ligger på samma nivå kekF = 0.001 i båda testen och iterationerna stannar
när den relativa residualen kAxk− bk/kbk < tol. Den ickesymmetriska matri-sen visade sig ha negativa egenvärden, så av anledningen att både DFPM och konjugerad gradient metodens konvergens förutsätter positiv definit matris
(även symmetrisk i det senare fallet) kördes det ickesymmetriska problemet istället med en variant av konjugerad gradient kallad BICG (biconjugate gra-dient) och en anpassning av DFPM så att den löser ATAx = ATb istället för Ax = b. I figur 4.1 demonstreras hur en bild avsuddad med iterativ metod kan se ut. 10 20 30 40 50 5 10 15 20 25 30 35 40 45 50 10 20 30 40 50 5 10 15 20 25 30 35 40 45 50 10 20 30 40 50 5 10 15 20 25 30 35 40 45 50
Figur 4.1: Problem 2 i konvergenstest för symmetrisk suddmatris. Original-bilden till vänster, suddig och brusig bild i mitten, där bruset e är taget från en likformig fördelning med kekF = 0.001, och konjugerad gradient lösningen
till höger.
4.3.1 Parameteranalys - DFPM
Ett problem som stöttes på under konvergenstesten var att användning av teoretiskt optimala parametrar η och ∆t inte ledde till konvergens på de problem som testades, åtminstone inte inom rimlig tid. Experiment utfördes på en ekvidistant grid [0, 2] × [0, 2] med h = 0.05 och vad som visade sig fungera bäst i praktiken var η1 = 0.45, ∆t1 = 1.55 för den symmetriska
matrisen samt η2 = 0.3, ∆t2 = 0.45 för den icke-symmetriska matrisen.
Notera möjligheten att med en finare grid hitta parametrar som presterar bättre.
För att hitta vilka egenvärden de experimentellt framtagna parametrarna tillhör kan vi lösa ekvation 4.5 för λmin och λmax, då fås istället systemet
x + y = 4 − 2η∆t(∆t)2
xy = η2(∆t)2,
där vi bytt ut λmin, λmaxmed x, y. Notera att x och y är utbytbara, vi tolkar
alltså en senare lösning som att det minsta värdet tillhör λminoch det största
λmax. Vi har från den första ekvationen att y = (4 − 2η∆t)/(∆t)2− x, som
med insättning i den andra ekvationen ger en andragradsekvation i x given av
x
(4 − 2η∆t)/(∆t)2− x= η2/(∆t)2.
Vid insättning av experimentellt framtagna optimala parametrar för den symmetriska matrisen η = 0.45 och ∆t = 1.55 erhålls lösningarna λmax≈ 1
0 500 1000 1500 2000 2500 10-20 10-15 10-10 10-5 100
Figur 4.2: Semilog plot på egenvärden tillhörande den symmetriska suddma-trisen.
samt λmin ≈ 0.0843. Enligt matlab ges det största egenvärdet av λmax ≈
0.961 samt det minsta av λmin ≈ 4.76 · 10−20, en semilog plot av alla
egen-värden syns i figur 4.2.
Minns från beviset till sats 4.1.1 att z(t) → 0 ⇐⇒ x(t) → ˆx, där z = PT(x − ˆx), och att varje zi(t) = C1er1 + C2er2. Med begynnelsevärden
x(0) = 0, ˙x(0) = 0 fås att z(0) = −PTx, ˙ˆ z(0) = 0. Vid insättning i zi(t) =
C1er1+ C2er2 erhålls två ekvationer
(
C1+ C2 = −pTi xˆ
C1r1+ C2r2 = 0,
(4.7) där pi betecknar kolumn i i matrisen P . Lösningen till 4.7 ges av C1 = αr2,
C2 = −αr1, där α = pTi x/(rˆ 1−r2). Vi har alltså att zi(t) = αr2er1t−αr1er2t.
Minns även att r(1,2) = −η/2 ±pη2/4 − λ
i, och antag utan inskränkning
av allmängiltigheten att |r2| |r1| = . Då gäller att C1 C2 = α,
men C1er1t→ 0 långsamt. Detta kan vara en del av förklaringen till varför
att kunna säga något med säkerhet behöver detta tittas på närmre.
4.4
En alternativ formulering
Med kunskapen om att en naiv lösning till ett avsuddningsproblem, d.v.s. en numerisk lösning till minimeringsproblemet
min
x kAx − bk 2,
ofta lider av stora fel kˆx − xk, och med tanke på feltermen MSE(ˆx, x) som introducerades vid analysen av iterativa metoder kan vi formulera ett nytt minimeringsproblem min x 1 2kAx − bk 2 2+ µ 2kˆx − xk 2 2, µ > 0, (4.8)
där ˆx uppfyller Aˆx = b. Detta med motiveringen att vi vill minska residualen kAx − bk men samtidigt ha ett litet fel kˆx − xk.
Sats 4.4.1. Minimeringsproblemet 4.8 ett unikt minimum och detta antas i punkten x = (ATA + µI)−1(ATb + µˆx). Bevis. Låt Φ(x) = 1 2kAx − bk 2 2+ µ 2kˆx − xk 2 2 = 1 2(Ax − b) T(Ax − b) +µ 2(ˆx − x) T(ˆx − x),
vars gradient ges av
∇Φ(x) = ATAx − ATb + µ(x − ˆx).
Eftersom Φ(x) är en positiv kvadratisk funktion har den ett minimum då ∇Φ(x) = 0 ⇐⇒ (ATA + µI)x = ATb + µˆx. (4.9)
Då ATA är positiv semidefinit och µ > 0 följer det att ATA + µI är positiv definit, vilket medför att ekvation 4.9 har en unik lösning given av
x = (ATA + µI)−1(ATb + µˆx).
Notera att formuleringen 4.8 kräver kunskap om hur den sanna bilden ˆx ser ut, något som nödvändigtvis inte är ett problem. Vid exempelvis stjärn-foto kan det vara ett mindre problem då många pixlar torde vara svarta. Förslagsvis kan man begränsa feltermen MSE(ˆx, x) till att bara ta hänsyn
20 40 60 80 100 10 20 30 40 50 60 70 80 90 100 20 40 60 80 100 10 20 30 40 50 60 70 80 90 100 20 40 60 80 100 10 20 30 40 50 60 70 80 90 100
Figur 4.3: Test av alternativa formuleringen. Originalbilden till vänster, sud-dig och brusig bild i mitten, där bruset e är taget från en likformig fördelning med kekF= 0.001, och lösningen till höger (µ = 0.2).
till pixlar man med relativ säkerhet vet vilken färg de ska ha, alternativt kan man använda sig utav uppskattningar. Matrisen ATA + µI är dock sym-metrisk positiv definit vilket betyder att både DFPM och den kändare kon-jugerad gradient metoden går att nyttja. I figur 4.3 demonstreras ett test av den här formuleringen. I testet användes DFPM som lösare med parame-tern µ = 0.2, och som referens i feltermen användes x∗ = ˆx + δx, där δx är brus från en likformig fördelning med kδxkF = kˆxkF/10, detta med syfte att
Bilaga A
Matlabkoder
1 f u n c t i o n [ x , k ] = DFPMatatest (A, b , x0 , my, dt , t o l ) 2 %DFPM a l g o r i t h m m o d i f i e d t o h a n d l e non−p o s i t i v e 3 %d e f i n i t e m a t r i c e s u s i n g t h e normal e q u a t i o n s 4 5 %s t a r t i n g v a l u e s xdot0 and x0 6 x d o t o l d=z e r o s(s i z e( b ) ) ; 7 x o l d=x0 ; 8 d o l d=A∗ x o l d−b ; 9 10 r e l r e s =norm( d o l d ) /norm( b ) ; 11 12 k =0; 13 w h i l e r e l r e s > t o l
14 xdot=x d o t o l d −(A’ ∗ d o l d+my∗ x d o t o l d ) ∗ dt ; 15 16 x d o t d t=xdot ∗ dt ; 17 18 d=d o l d+A∗ x d o t d t ; 19 x=x o l d+x d o t d t ; 20 21 x o l d=x ; 22 x d o t o l d=xdot ; 23 d o l d=d ; 24 25 k=k +1; 26 r e l r e s =norm( d o l d ) /norm( b ) ; 27 end
1 f u n c t i o n [ x , k ] = DFPMtest (A, b , x0 , my, dt , t o l ) 2 %DFPM a l g o r i t h m 3 4 %s t a r t i n g v a l u e s xdot0 and x0 5 x d o t o l d=z e r o s(s i z e( b ) ) ; 6 x o l d=x0 ;
7 r e l r e s =norm(A∗ x0−b ) /norm( b ) ; 8
9 k =0;
10 w h i l e r e l r e s > t o l
11 xdot=x d o t o l d +(b−A∗ x o l d−my∗ x d o t o l d ) ∗ dt ; 12 x=x o l d+xdot ∗ dt ; 13 14 x o l d=x ; 15 x d o t o l d=xdot ; 16 17 k=k +1;
18 r e l r e s =norm(A∗ x o l d−b ) /norm( b ) ; 19 end 1 f u n c t i o n [ x , k ] = cg (A, b , x0 , t o l ) 2 %C o n j u g a t e g r a d i e n t method 3 x=x0 ; 4 r=b−A∗ x0 ; 5 k =0; 6 p=r ; 7 normb=norm( b ) ; 8 9 w h i l e norm( r ) /normb > t o l 10 x o l d=x ; 11 r t r=r ’ ∗ r ; 12 Ap=A∗p ; 13 a l f a =( r t r ) / ( p ’ ∗ Ap) ; 14 15 x=x o l d+a l f a ∗p ; 16 r=r−a l f a ∗Ap ; 17 b e t a=(r ’ ∗ r ) / ( r t r ) ; 18 p=r+b e t a∗p ; 19 20 k=k +1; 21 end
1 f u n c t i o n [ x , xhat , k ] = b i c o n c g (A, b , x0 , xhat0 , t o l ) 2 %B i c o n j u g a t e g r a d i e n t method
3 i f i s c o l u m n ( x0 ) ~= 1 | | i s c o l u m n ( xhat0 ) ~= 0 4 e r r o r( ’ make x0 column and xhat0 row ’) ; 5 end 6 x=x0 ; 7 xhat=xhat0 ; 8 9 r=b−A∗x ; 10 r h a t=b’− xhat ∗A ’ ; 11 12 k =0; 13 p=r ; 14 phat=r h a t ; 15 r t r n e w=r h a t ∗ r ; 16 normb=norm( b ) ; 17 18 w h i l e norm( r ) /normb > t o l 19 x o l d=x ; 20 x h a t o l d=xhat ; 21 22 r t r=r t r n e w ; 23 Ap=A∗p ; 24 25 a l f a =( r t r ) / ( phat ∗Ap) ; 26 27 x=x o l d+a l f a ∗p ; 28 xhat=x h a t o l d+c o n j( a l f a ) ∗ phat ; 29 30 r=r−a l f a ∗Ap ; 31 r h a t=r h a t −c o n j( a l f a ) ∗ phat ∗A; 32 r t r n e w=r h a t ∗ r ; 33 34 b e t a=( r t r n e w ) / ( r t r ) ; 35 p=r+b e t a∗p ; 36 phat=r h a t+c o n j(b e t a) ∗ phat ; 37 38 k=k +1; 39 end
1 f u n c t i o n [ x , k ] = l s q r 1 (A, b , t o l ) 2 %A l g o r i t h m LSQR 3 n=s i z e( b , 1 ) ; 4 5 i f s i z e(A, 1 ) ~= n 6 e r r o r( ’ make s u r e A i s s q u a r e ’) 7 end 8
9 u=b ; b e t a=norm( u ) ; u=u/b e t a; 10 v=A’ ∗ u ; a l f a=norm( v ) ; v=v/ a l f a ; 11 w=v ; x=z e r o s( n , 1 ) ; 12 p h i h a t=b e t a; r h o h a t=a l f a ; 13 k =0; 14 normb=norm( b ) ; 15
16 w h i l e norm(A∗x−b ) /normb > t o l
17 u=A∗v−a l f a ∗u ; b e t a=norm( u ) ; u=u/b e t a; 18 v=A’ ∗ u−b e t a∗v ; a l f a=norm( v ) ; v=v/ a l f a ; 19 20 rh o=s q r t( r h o h a t^2+b e t a^2) ; 21 c=r h o h a t / r ho ; 22 s=b e t a/ r ho ; 23 t h e t a=s ∗ a l f a ; 24 r h o h a t=−c ∗ a l f a ; 25 p h i=c ∗ p h i h a t ; 26 p h i h a t=s ∗ p h i h a t ; 27 28 x=x+( p h i / rho ) ∗w ; 29 w=v−( t h e t a / rh o ) ∗w ; 30 31 k=k +1; 32 end
Litteraturförteckning
[1] Sverker Edvardsson, Magnus Neuman, Per Edström, and Håkan Olin. Solving equations through particle dynamics. Computer Physics Com-munications, 197:169–181, 2015.
[2] Mårten Gulliksson, Magnus Ögren, Anna Oleynik, and Ye Zhang. Dam-ped Dynamical Systems for Solving Equations and Optimization Pro-blems, pages 1–44. Springer International Publishing, 2019.
[3] Per Christian Hansen. Discrete Inverse Problems: Insight and Algo-rithms. Society for Industrial and Applied Mathematics, 2010.
[4] Per Christian Hansen, James G. Nagy, and Dianne P. O’Leary. Deblur-ring Images: Matrices, Spectra, and FilteDeblur-ring (Fundamentals of Algo-rithms). Society for Industrial and Applied Mathematics, 2006.
[5] Peter D. Lax and Robert D. Richtmyer. Survey of the stability of li-near finite difference equations. Communications on Pure and Applied Mathematics, 9:267–293, 1956.
[6] Don Morgan and Sanzheng Qiao. Analysis of damped mass-spring systems for sound synthesis. EURASIP Journal on Audio, Speech, and Music Processing, 2009.
[7] Christopher C. Paige and Michael A. Saunders. Lsqr: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Softw., 8:43–71, 1982.
[8] Yousef Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2003. [9] Charles Van Loan. Computational Frameworks for the Fast Fourier
Transform. Society for Industrial and Applied Mathematics, 1992. [10] Curtis R. Vogel. Computational Methods for Inverse Problems. Society
for Industrial and Applied Mathematics, 2002.
[11] David S. Watkins. Fundamentals of Matrix Computations. John Wiley & Sons, Inc., 1991.
