Singulární rozklad – aplikace v image deblurring M. Plešinger, Z. Strakoš TUL, Fakulta mechatroniky, Liberec AV ČR, Ústav informatiky, Praha
1
Úvod
Uvažujme obecnou reálnou matici A ∈ Rn×m ,
n m,
r = rank(A) ≤ min{n, m} .
Pak existuje rozklad A = U ΣV T,
kde U ∈ Rn×n ,
Σ ∈ Rn×m ,
V ∈ Rm×m ,
(1)
takový, že matice U , V jsou ortogonální, a matice Σ má jediný nenulový diagonální blok řádu r Σr 0 , Σr = diag(σ1 , . . . , σr ) , přičemž platí σ1 ≥ σ2 ≥ . . . ≥ σr > 0 . (2) Σ = 0 0 Důkaz nalezneme v dnes již klasické literatuře, například [1], [2]. Rozklad (1, 2) se nazývá singulární rozklad (SVD) matice A. Sloupce ui matice U jsou levé singulární vektory a sloupce vi matice V jsou pravé singulární vektory. Čísla σi se nazývají singulární čísla. Singulární rozklad (1, 2) lze přepsat do dyadického tvaru A = U ΣV T =
r
ui σi viT .
(3)
i=1
Užitím (3) lze snadno řešit obecné soustavy lineárních algebraických rovnic ve smyslu nejmenších čtverců r T ui b † (4) vi , Ax ≈ b, x = A b = σi i=1 kde matice A† = ri=1 vi σi−1 uTi je tzv. Moore-Penroseova pseudoinverze a vektor x je řešení ve smyslu nejmenších čtverců minimální v normě (je-li A čtvercová nesingulární, platí A† = A−1 a x je řešením v klasickém smyslu), viz [2, Chap. 4.3].
2
Image deblurring
Obrázek vysoký n a široký m pixelů (pro jednoduchost uvažujme černobílý) lze interpretovat jako matici dimenze n×m. Tento bude hledaný vzor, příslušnou matici označíme X. Rozmazaný, nebo jinak zdeformovaný obraz tohoto vzoru, v praxi například v důsledku průchodu optickou soustavou, pak označíme B. Je-li proces rozmazání lineární a nezávislý po sloupcích a řádcích, lze zapsat jako maticová rovnice (viz [3]) AC X ATR = B ,
kde AC ∈ Rn×n , 1
AR ∈ Rm×m ,
X, B ∈ Rn×m .
(5)
Maticová rovnice (5) představuje soustavu lineárních algebraických rovnic, lze tedy přepsat v obvyklém tvaru (6) A x = b , kde A = AR ⊗ AC ∈ Rnm×nm je čtvercová, za předpokladu rank(AC ) = n ,
rank(AR ) = m ,
(7)
regulární matice a x = [ xT1 , . . . , xTm ]T ∈ Rnm ,
b = [ bT1 , . . . , bTm ]T ∈ Rnm .
(8)
Symbolem ⊗ značíme Kroneckerův součin a xi , bi jsou i-té sloupce matic X, B. right−hand side B
50
100
150
200
250 50
100
150
200
250
300
Obrázek 1: Rozmazaný obrázek – pravá strana B maticové rovnice (5). Je jasné, že pro nalezení X musíme znát matice AC , AR , respektive A. U reálných úloh vždy vycházíme z matematického popisu fyzikální soustavy, která původní vzor X deformovala na obraz B, a na základě tohoto popisu konstruujeme příslušné matice. U reálných úloh tak matici A nikdy neznáme přesně, máme k dispozici pouze její odhad. Pozorovaný obrázek navíc může být, a často bývá, zatížen šumem, který není v matematickém modelu zahrnut, tedy platí b = b EXACT + b NOISE . V případě pravé strany na obrázku 1 známe matice AC , AR popisující model přesně, podmínka (7) je splněna. Užitím vztahu (4) se pokusíme nalézt řešení soustavy (6, 8), tedy nalézt vzor X. Toto řešení označíme jako „naivní nm T ui b −1 (9) vi = x + A−1 b NOISE . x NAIVE = A b = σi i=1
Z obrázku 2, je patrné, že tento postup selhává, to co jsme získali je pouze šum, který pravá strana obsahovala. 2
naive solution
50
100
150
200
250 50
100
150
200
250
300
Obrázek 2: Naivní řešení x = A−1 b nemá se skutečným řešením nic společného.
3
Diskrétní Picardova podmínka
Když budeme zvětšovat rozlišení obrázku B, nejmenší singulární čísla matice A se budou blížit k nule, tedy lim i → ∞ σi = 0. Budeme-li považovat matici A za přesný model dané fyzikální soustavy a budeme-li pracovat s přesnou pravou stranu (bez šumu), musí limita nm uT b EXACT i vi lim n, m → ∞ σi i=1
reprezentovat původní spojitý obrázek, musí tedy být omezena. Picardova podmínka říká, že projekce uTi b EXACT přesné pravé strany do levých singulárních podprostorů matice A, musí s rostoucím indexem klesat k nule rychleji než singulární čísla σi matice A. Musí klesat tak rychle, aby platilo T ui b EXACT vi = 0 lim i→∞ σi pro rozlišení zvětšující se nade všechny meze. Diskrétní Picardova podmínka je její diskrétní aproximací (tedy, vágně formulováno, projekce musí klesat „rozumně rychle, viz obrázek 3). Picardova podmínka vychází z nutnosti existence zobrazovaného reálného vzoru X. Šum b NOISE zatěžující pravou stranu Picardovu podmínku splňovat nemusí, a často také nesplňuje. Protože singulární čísla klesají k nule, šum je zesílen a v naivním řešením (9) dominuje nad hledaným x, platí následující nerovnost x2 = A−1 b EXACT 2 A−1 b NOISE 2 .
3
singular values of A and projections uTi b
4
10
right−hand side projections on left singular subspaces uTi b singular values σi
2
10
noise level
0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
0
1
2
3
4
5
6
7
8 4
x 10
Obrázek 3: Plnou čárou jsou vykreslena singulární čísla matice, body jsou projekce pravé strany do odpovídajících levých singulárních podprostorů. Čárkovanou čárou je naznačena hladina šumu, cca 10−3 , projekce uTi b od indexu i ∼ 0.3 × 104 díky přítomnosti šumu neklesají.
4
Regularizace řešení
Nalezení vzoru X přesným řešením soustavy (6)-(8) není možné, díky šumu, který je obsažen v pravé straně. Regularizací řešení chápeme snahu potlačit vliv šumu maximálním možným způsobem. Klasickou regularizační technikou je truncated SVD (TSVD). TSVD řešení soustavy (6)-(8) spočívá v nahrazení matice A maticí Ak s nižší hodností k. Přičemž tato matice je (v normě) nejlepší aproximací A pro dané k, platí (viz [2]) Ak =
k i=1
ui σi viT .
Parametr k volíme tak, abychom co nejvíce omezili vliv šumu. Modifikovanou soustavu Ak x ≈ b řešíme ve smyslu nejmenších čtverců, hledáme řešení minimální v normě xk =
A†k
k T ui b b = vi . σi i=1
4
(10)
TSVD solution, k = 2983
50
100
150
200
250 50
100
150
200
250
300
Obrázek 4: Řešení (10) regularizované užitím TSVD, parametr k = 2983.
5
Závěr
Důležitou otázkou v oblasti regularizačních metod je problém zastavení, v případě TSVD problém volby optimálního k. Protože ne vždy máme představu o tom, jak má řešení vypadat, ne vždy lze spočítat celý singulární rozklad (tak jako na obrázku 3), abychom nalezli hladinu šumu. Obvykle se pro zastavení používají L-curve kritérium nebo generalized cross-validation kritérium (viz [3]). Acknowledgement: This work was supported by the National Program of Research “Information Society” under project 1ET400300415.
Reference [1] Gene H. Golub, Charles F. Van Loan: Matrix Computations, 3rd Edtition. The Johns Hopkins University Press, Baltimore, 1996. [2] David S. Watkins: Fundamentals of Matrix Computations, 2nd Edtition. Wiley-Interscience, New York, 2002. [3] Per Christian Hansen, James G. Nagy, Dianne P. O’Leary: Deblurring Images – Matrices, Spectra and Filtering. SIAM, 2006.
5
