Eötvös Loránd Tudományegyetem Természettudományi Kar
A szinguláris érték felbontás alkalmazásai a jel- és képfeldolgozásban
Szakdolgozat
Készítette:
Réti Attila ELTE TTK Alkalmazott matematikus BSc
Témavezet®: Kovács Péter Tanársegéd ELTE IK Numerikus Analízis Tanszék
Budapest, 2016
Tartalomjegyzék
1. Bevezetés
4
2. Matematikai háttér
6
2.1. De níciók, problémafelvetés . . . . . . . . . . . . . . . 2.1.1. Hilbert térbeli approximáció . . . . . . . . . . . 2.1.2. PCA . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3. A kovarianciamátrix közelítése ML-becsléssel . . 2.2. Szinguláris érték szerinti felbontás . . . . . . . . . . . . 2.3. Numerikus módszerek . . . . . . . . . . . . . . . . . . 2.3.1. Givens forgatások . . . . . . . . . . . . . . . . . 2.3.2. Bidiagonalizálás Householder transzformációval 2.3.3. Golub-Reinsch algoritmus . . . . . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
3. Alkalmazások 3.1. 3.2. 3.3. 3.4.
Tömörítés . . . . . . . Adatsimítás . . . . . . Arcképek osztályozása Kitekintés . . . . . . .
6 6 8 10 11 21 21 23 25
30 . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
A. Algoritmusok
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
30 33 34 35
37
2
Köszönetnyilvánítás
Ezúton szeretném megköszönni témavezet®mnek, Kovács Péternek, hogy elvállalta a konzulensi teend®ket, bármikor a rendelkezésemre állt, mindig ösztönzött és mindenben segítette a munkámat. Nélküle ez a dolgozat nem jöhetett volna létre. Szeretnék továbbá köszönetet mondani barátn®mnek, szüleimnek és testvéremnek, akik mindig mellettem álltak és támogattak a munkám során.
3
1. fejezet
Bevezetés
A huszonegyedik századra az információáramlás olyannyira meghatározza a világunkat, hogy a mindennapi élet szerves részét képezik például a különböz® képek, zenék, videók és cikkek. A számítógépek és okostelefonok korában az emberek nagy többsége rendszeresen dolgozik nagy adathalmazokkal, digitális fotókkal, táblázatokkal. A hadiiparban már külön drón egységeket gyártanak, melyek feladata egy adott táj domborzatának és lakosságának feltérképezése az általuk rögzített jelek alapján. A modern orvostudományban is számos képet és jelet dolgoznak fel, pl. MRI felvételeket, EKG jeleket, röntgen vagy CT képeket. Ilyen magas igények mellett jogosan lehet kíváncsi az ember, hogy milyen eszközök, lehetséges eljárások és módszerek teszik lehet®vé ennek a nagyon magas színvonalú és folyamatosan fejl®d® rendszernek a m¶ködését A dolgozat célja a szinguláris felbontás és annak tulajdonságainak ismertetése, elméleti hátterének alapos és részletes áttekintése, illetve a szinguláris érték szerinti felbontás alkalmazásainak bemutatása a jel- és képfeldolgozásban. De mit is értünk digitális jelfeldolgozás alatt? Digitális jelfeldolgozásnak (angolul Digital Signal Processing, DSP) nevezzük azt a folyamatot, amikor egy zikai mennyiséget valamilyen módon olyan formára hozunk, hogy számítógéppel feldolgozhatóvá, alakíthatóvá váljon. Egy képet például tekinthetünk egy n × m-es mátrixnak (n, m ∈ N+ ), melynek elemei egész számok. A mátrixon belüli elhelyezkedés az adott pixel pozícióját jelenti, az érték pedig az adott színárnyalatot. Ily módon numerikusan reprezentáltunk egy képet és tudunk vele különböz® m¶veleteket végezni. Adódhat azonban a kérdés, hogy mi történik akkor, ha nagy mennyiség¶ adatot szeretnénk eltárolni vagy feldolgozni. Hogyan tudunk kevesebb adatot tárolni, mint amennyi valójában a rendelkezésünkre áll anélkül, hogy (lényeges) információt veszí4
tenénk? Hogyan tudjuk kevesebb m¶velettel ugyan azt a végeredményt elérni? Ezekre, és ehhez hasonló kérdésekre ad majd választ a dolgozat. Látni fogjuk, hogy a probléma valójában visszavezethet® különböz® matematikai problémák együttesére pl.: sajátérték feladat, szinguláris értékfelbontás (angolul Singular Value Decomposition, SVD), Hilbert térbeli approximáció, altérre való projekció, vagy éppen kovariancia mátrix közelítése. A szakdolgozat második fejezetében kezdetben a dolgozat megértéséhez szükséges matematikai hátteret ismertetem [1] és [2] alapján. Megismerkedünk az euklideszi terekkel, tetsz®leges ponthoz megszerkesztjük egy adott altér legjobban közelít® elemét, majd megadjuk, hogyan érdemes megválasztani az alteret, amire projekciót végzünk. Bepillantunk a f®komponens-analízis (angolul Principal Component Analysis, PCA) statisztikai eljárásába majd innen eljutunk a szinguláris értékfelbontásig. Találkozni fogunk a téma néhány f®bb tételével (pl. EckartYoung tétel) amiket részletesebben taglalok és bizonyítok is. A fejezet végén pedig a numerikus módszerekr®l fogok beszélni, a [3] irodalom alapján ismertetem a Gene Howard Golub és Christian Reinsch által kidolgozott SVD algoritmust, az ahhoz szükséges segédalgoritmusokat és lépéseket. Az utolsó, harmadik fejezetben látni fogjuk, mire is alkalmazható az el®z® fejezetben felhalmozott elmélet: az általam készített programmal néhány példán keresztül bemutatom, hogyan is m¶ködik az SVD a gyakorlatban, a bevezetésben felvetett képtömörítési problémát miként oldhatjuk meg. Ezen kívül szó esik még egydimenziós adatok simításáról illetve arcképek osztályozásáról. A dolgozatomat egy kitekintéssel zárom, melyben a szinguláris értékfelbontás egy további alkalmazását fogom szemléltetni. Megjegyzem, hogy a dolgozatban bemutatott numerikus módszereket, így például a szinguláris érték felbontás különböz® változatait
Matlab-ban implementáltam. A
kódok és a példaprogramok a csatolt CD mellékleten elérhet®ek.
5
2. fejezet
Matematikai háttér
2.1. De níciók, problémafelvetés Ez a fejezet a fontos de níciók, tételek és eljárások ismertetéséül szolgál. Az olvasóról feltételezzük a matematikai el®képzettséget, ezért a fejezetben használt alapfogalmakat ismertnek tekintjük. A veszteséges tömörítés az adattömörítési algoritmusok egy osztálya, ami a veszteségmentes tömörítéssel ellentétben nem teszi lehet®vé a tömörített adatból az eredeti adatok pontos rekonstrukcióját, ám egy elég jó rekonstrukciót igen. A mi esetünkben is ez történik, veszteségesen fogunk tömöríteni, viszont keressük a lehet® legjobb közelítését egy adott képnek. 2.1.1. Hilbert térbeli approximáció
2.1.1. De níció. Legyen (Z, ||.||) Banach-tér, Y ⊂ Z, továbbá legyen b ∈ Z. Ekkor azt mondjuk, hogy az yb ∈ Y a b elem legjobb közelítése Y -ban, ha ||b − yb || ≤ ||b − y|| ∀ y ∈ Y -ra. A mi esetünkben Hilbert téren belül gondolkodunk, azaz teljes euklideszi tér felett, hiszen a képünket egy Rm×n (m, n ∈ N+ ) feletti mátrixként reprezentáljuk, melyben minden pixel- elhelyezkedést®l függ®en- egy mátrixbeli elemnek felel meg. Az egyszer¶ség kedvéért tekintsünk szürkeárnyalatos képeket, melyek pixelein a 0-255 közötti egész számok intenzitás értékeket jelölnek. Így tehát valós számokat tartalmazó vektorokkal illetve mátrixokkal van dolgunk, Hilbert térben keressük majd egy kitüntetett altér legjobban közelít® elemét. A következ® állítás ennek létezésére és egyértelm¶ségére ad választ: 6
2.1.2. Tétel. Legyen (H, h.i) Hilbert tér, f ∈ H és H 0 ⊂ H zárt altér. Ekkor ∃! f 0 ∈ H 0 : ||f − f 0 || = inf h∈H 0 ||f − h|| és f − f 0 ⊥ H 0 (azaz hf − f 0 , hi = 0 ∀h 0 ∈ H 0 ). Mivel véges dimenziós vektorterekben dolgozunk, adódhat a kérdés, hogy milyen alakban tudjuk megkapni a kívánt közelítést. Természetesen a válasz ekvivalens azzal, hogy milyen távol van a közelítend® elemt®l a közelít® elemet tartalmazó altér. Gondoljunk például a 2 dimenziós esetre ahol értelemszer¶en egy pont és egy egyenes között a legrövidebb út a pontot az egyenessel összeköt® mer®leges szakasz. Az alábbi tétel ezt fogalmazza meg véges dimenziós alterekre:
2.1.3. Tétel. Legyen b az f és a gi -kb®l álló skalárszorzatokat tartalmazó vektor, azaz b = hf, gi i (i = 1 . . . n). Legyen d a legjobban közelít® elem távolsága. Ekkor az alábbi egyenl®ség teljesül: d2 := ||f − f 0 ||2 = ||f||2 − bT c. A H 0 véges dimenziós altér esetén H 0 = span(g1 , . . . , gn ), ahol n a dimenziószámot, gi (i ∈ 1, . . . , n) pedig a bázisvektorokat jelöli. Ezen jelölések mellett a legjobban P közelít® elemet:f 0 = ni=1 ci gi alakban keressük, ahol ci tetsz®leges valós skalár. Fon-
2.1. ábra. A v vektor legjobb közelítése a szürkével jelölt H altérre nem más, mint a vektor vetülete az altérre, azaz a ProjH v.1 tos észrevétel, hogy ha a bázisvektorok valamilyen kitüntetett tulajdonságú rendszert alkotnak, akkor a feladat leegyszer¶södik.
2.1.4. Megjegyzés. Ha g1 , . . . , gn ortogonális rendszer(OGR), akkor a legjobban közelít® elem: f0 = 1 Az
ábra
megtalálható
a
n X hf, gi i gi hg i , gi i j=1
http://tex.stackexchange.com/questions/117789/setting-camera-
projection-angle-in-tikz-to-enhance-the-depth-perception-when-dra weboldalon.
7
Ortonormált rendszer esetén, mivel hgi , gi i = 1, a fenti egyenl®ség egyszer¶södik: f0 =
n X hf, gi igi j=1
A kés®bbiekben látni fogjuk, hogy ezek a speciális alakok fontos szerepet játszanak az alkalmazás során. 2.1.2. PCA
A f®komponens-analízis, azaz a PCA egy többváltozós statisztikai eljárás, amely adatredukciós módszerként is használható. Lényege, hogy egy er®sen korreláló változókat tartalmazó nagy adathalmaz dimenzióját lecsökkenti a szórásirányok gyelembevételével. Ezt egy ortogonális transzformáció segítségével éri el: a lehetségesen korreláltatható változókat lineárisan korrelálatlanná alakítja át. Ezeket a változókat nevezzük f®komponenseknek, innen adódik az eljárás neve. A fejezet további részeit a [4] és [7] irodalom alapján fogom ismertetni. A PCA-nak számos vonzó el®nye van, ami miatt az egyik leggyakrabban használt dimenziócsökkent® eljárás. Többek között segítségével könnyebb megtalálni az adatokat legjobban jellemz® mintázatokat, így például az adatok nagyon nagy mértékben tömöríthet®k viszonylag kevés veszteséggel. Nem elhanyagolható, hogy a dimenziócsökkentés zajsz¶résre, simításra is alkalmas. Matematikailag gyakorlatilag nem más, mint egy ortogonális lineáris transzformáció, mely az adatot egy új koordináta rendszerbe transzformálja, azaz egy bázistranszformációt hajt végre. Most nézzük a feladat formális felírását. Legyen X = [x1 , x2 , . . . , xn ] egy n dimenziós valószín¶ségi vektorváltozó. Legyen ennek az X valószín¶ségi vektorváltozónak a várhatóértéke 0. Ez feltehet®, hiszen ha mondjuk egy tetsz®leges m konstans lenne, akkor a mintaátlagot a diszkrét adatokból levonva várhatóértéknek már zérust kapnánk. Azt szeretnénk, hogy a különböz® koordináták korrelációja 0 legyen, tehát:
0, ha i 6= j E(xi xj ) = σ2i , különben, 8
(2.1)
ahol σ2i az xi -hez tartozó szórásnégyzet és E(x) az x várható értéke. A feladat megoldásához szükségünk van az X kovariancia mátrixára:
Cx = E((X − m)(X − m)T ) = E(XXT ) A korábbiak alapján tehát szükségünk van egy lineáris ortogonális transzformációra úgy, hogy a keletkez® transzformált elem is teljesítse a (2.1) egyenl®ség feltételeit. Legyen φ : Rn → Rn lineáris ortogonális transzformáció, Y = φ(x) = φT x transzformált valószín¶ségi vektorváltozó, mely teljesíti a fent említett feltételt:
Cy = E(YY T ) = E(φT x(φT x)T ) = E(φT (xxT )φ) = φT Cx φ Balról φ-vel szorozva könnyen látszik, hogy:
Cx φ = φCy Figyelembe véve, hogy Cy a (2.1) egyenl®ség alapján egy diagonális mátrix, továbbá a sajátérték és sajátvektor de nícióját használva adódik, hogy:
Cx φi = λi φi Ebb®l következik, hogy a (2.1) egyenl®ség nem másra redukálódik, mint egy sajátérték feladatra, azon belül is Cx diagonalizálására. Feltehet®, hogy a diagonalizált elemek, tehát a sajátértékek fentr®l lefelé abszolút értékben csökken® sorrendben következnek. (Ellenkez® esetben ekvivalens mátrixm¶veletekkel olyan alakra hozzuk.) Ekkor, mivel euklideszi térben vagyunk, és a {φi | i = 1, . . . , n} rendszer ortogonális, az el®z® fejezet alapján X-et `2 normában az alábbi formában közelíthetjük legjobban: 0
X≈X =
k X
hX, φi iφi
j=1
Az is látszik, hogy ez nem más, mint az els® k legnagyobb sajátértékhez tartozó sajátaltérre vett projekció, ahol k-t növelve csökkenthetjük a hibát. A gyakorlatban a valódi kovariancia mátrixot nem ismerjük, ezt az n1 XXT formulával szokás közelíteni, ahol X oszlopait a valószín¶ségi vektorváltozó mintavételével kapjuk.
9
2.1.3. A kovarianciamátrix közelítése ML-becsléssel
Korábban azt állítottam, hogy a képek feldolgozása során a kovariancia mátrixot az 1 XXT -tal n
szokás közelíteni, most pedig a [5] alapján meg is magyarázom, hogy mi-
ért van így. Jelölje most µ a várható értéket, θ a paraméter vektort, ∇θ a gradiens operátort, Σ pedig a kovariancia mátrixot. Legyen X Gauss-normális eloszlású valószín¶ségi változó. Az egyszer¶ség kedvéért tekintsük el®ször azt az esetet, amikor a várható érték nem ismert, a kovariancia mátrix viszont igen. Ezen feltételek mellett tekintsük az Xk mintavételi pontot és legyen
1 1 log P(Xk |µ) = − log (2π)d |Σ| − (Xk − µ)T Σ−1 (Xk − µ) 2 2 és
∇θ log P(Xk |µ) = Σ−1 (Xk − µ). Azonosítsuk θ-t µ-vel. Az el®z® egyenletb®l látszik, hogy a µ maximum likelihood (ML) becslésének ki kell elégítenie a n X
Σ−1 (Xk − µ ^) = 0
k=1
egyenletet. Ezt alakítva az alábbi formát nyerjük.
1X Xk . µ ^= n k=1 n
Most pedig ezek alapján tekintsük azt az esetet, amikor sem µ, sem pedig Σ nem ismert. Így most ezek az elemek alkotják a θ paraméter vektor komponenseit. Tekintsük el®ször az egyváltozós esetet, ahol θ1 = µ és θ2 = σ2 . Ekkor egyetlen pont log-likelihoodja
1 1 log P(Xk |θ) = − log 2πθ2 − (Xk − θ1 )2 2 2θ2
és a deriváltja
" ∇θ l = ∇θ log P(Xk |θ) =
1 (Xk − θ1 ) θ2 −θ1 )2 1 + (Xk2θ 2 2θ2 2
# .
A likelihood egyenlet miatt a teljes log-likelihood függvényre teljesül, hogy n X 1 (Xk − θ^1 ) = 0 θ^2 k=1
10
és
−
n n X X 1 (Xk − θ^1 )2 = 0, + ^2 ^22 θ θ k=1 k=1
^1 és θ^2 a maximum likelihood becslése θ1 -nek és θ2 -nek. Cseréljük most ki ezeket ahol θ a becsléseket úgy, hogy µ ^ = θ^1 és σ ^ 2 = θ^2 legyen. Ekkor néhány alakítással adódik, hogy
és
1X µ ^= Xk n k=1 n
1X σ ^ = (Xk − µ ^ )2 . n k=1 n
2
Ebb®l következik, hogy a kovariancia mátrix közelítése:
X ^= 1 Σ (Xk − µ ^ )(Xk − µ ^ )T . n k=1 n
Ez pedig nem más, mint a kívánt egyenl®ség. Tehát a mintevételezéssel kapott X =
[X1 , . . . , Xn ] mátrixxal az n1 XXT formula várható értékben jól közelíthet® a valódi kovariancia mátrix. Megjegyzem, hogy a technikát számos területen alkalmazzák a gyakorlatban, például mintafelsimerés során, illetve klaszterezési, osztályozási feladatokban (b®vebben lásd könyv 3. fejezetét). Ezt a közelítést viszont igen nehéz kiszámítani, ha XXT nagyon nagy, ugyanis akár csak egy 640×480 pixelb®l álló kép esetén is hatalmas, több gigabyte-nyi memóriára lenne szükségünk a szorzat meghatározására. Erre lehet alkalmazni egy másik dimenziócsökkent® algoritmust: a szinguláris értékfelbontást. A PCA egy sajátértékprobléma, melyet az SVD-vel szinguláris érték problémává konvertálhatunk. Ez utóbbi megoldásához már léteznek olyan numerikus módszerek, melyekkel a feladat kezelhet® méret¶ válik. A következ® fejezetben az SVD algoritmusról és annak elméletéi hátterér®l fogok beszélni.
2.2. Szinguláris érték szerinti felbontás A szinguláris érték felbontást (angolul Singular Value Decomposition, SVD) eredetileg a dierenciál geometriában fedezték fel. Azt szerették volna megállapítani, hogy vajon egy valós bilineáris formát egyenl®vé lehet e tenni egy másikkal független orto11
gonális transzformációkat használva. Eugenio Beltrami és Camille Jordan egymástól függetlenül, már 1873-ban és 1874-ben felfedezte az SVD létezését. Harmadikként James Joseph Sylvester érkezett el a valós négyzetes mátrixok szinguláris érték szerinti felbontásához, nyilvánvalóan függetlenül az el®deihez. egy A mátrix szinguláris értékeit az A mátrix kanonikus szorzóinak nevezte. A negyedik matematikus Autonne volt 1915-ben, aki felfedezte a szinguláris felbontást. Érdekesség, hogy a polár felbontás során jutott el az SVD-ig. Az els® bizonyítás, tetsz®leges Cm×n -beli mátrixok szinguláris érték felbontására, Carl Eckart és Gale Young nevéhez f¶z®dik 1936-ból. Úgy látták, hogy ez a f®tengely transzformáció általánosítása hermitikus vagy Hermite-mátrixokra. A szinguláris felbontás, hasonlóan a PCA-hoz, egy valós, vagy komplex elemeket tartalmazó mátrix faktorizációja. A felbontást a tudomány számos területén alkalmazzák, például a statisztikában és a jelfeldolgozásban. Felhasználható többek közt a legkisebb négyzetek módszeréhez, kis rangú mátrixszal való közelítéshez, a Kabsch algoritmushoz, vagy az általánosított inverz számolására. El®ször is nézzük, mit is jelent egy mátrix szinguláris érték szerinti felbontása!
2.2.1. De níció (Szinguláris felbontás). Egy A ∈ Cn×m (n, m ∈ N+ ) mátrix szinguláris felbontásán egy A = UDV ∗
(2.2)
szorzattá bontást értünk, ahol U ∈ Cn×n , V ∈ Cm×m unitér, D ∈ Cn×m diagonális mátrixok.
2.2.2. De níció
(Szinguláris érték). A D mátrix f®átlójában lév® σ1 ≥ σ2 ≥ · · · ≥
σn ≥ 0 elemeket szinguláris értékeknek nevezzük. σi =
p λi (A∗ A)
(i = 1 . . . n),
ahol tetsz®leges B ∈ Cn×n mátrixra a λi (B) a B mátrix i. legnagyobb abszolút érték¶ sajátértékét jelenti.
2.2.3. De níció
(Szinguláris vektor). A felbontásból kapott U és V oszlopait bal-,
illetve jobboldali szinguláris vektoroknak nevezzük. Fentebb de niáltuk, mi is az SVD, most pedig a létezésére fogunk látni elégséges feltételt. 12
2.2. ábra. Az A mátrix szinguláris felbontása abban az esetben, ha n > m. Természetesen ellenkez® esetben is hasonló a felbontás.
2.2.4. Tétel
(Szinguláris érték felbontás létezése). Tetsz®leges A ∈ Cn×m (n, m ∈
N+ ) mátrixra létezik a szinguláris felbontás.
Bizonyítás. Legyen r = rang(A). Tudjuk, hogy A∗ A mátrix önadjungált, továbbá azt is, hogy pozitív szemide nit, hiszen
hA∗ Ax, xi = x∗ (A∗ Ax) = (Ax)∗ Ax = hAx, Axi ≥ 0. Az A∗ A mátrix nem nulla sajátértékeit jelöljük σ1 . . . σr -rel. Ekkor természetesen a maradék nemnegatív sajátértéke 0, melynek multiplicitása n − r. Legyen továbbá Σ = diag(σ1 . . . σr ) ∈ Cr×r . Mivel A∗ A önadjungált, ezért van diagonális Schur-felbontása, azaz létezik V ∈ Cn×n unitér mátrix, hogy A∗ A = VD2 V ∗ , ahol
" D2 = diag(λi (A∗ A)) =
Σ2 0 0 0
# .
Mivel A∗ AV = VD2 , ezért V oszlopai az A∗ A sajátvektorai. Bontsuk fel V -t két komponensre: legyen V = [V1 , V2 ], ahol V1 tartalmazza a nem nulla sajátértékekhez tartozó sajátvektorokat, V2 pedig a nullához tartozóakat. Vizsgáljuk most a nem nulla sajátértékekhez tartozó sajátvektorokat tartalmazó részt:
V1∗ A∗ AV1 = Σ2 Balról és jobbról a szigma Σ−1 -el beszorozva egy r × r nagyságú egységmátrixot ka13
punk.
Σ−1 V1∗ A∗ AV1 Σ−1 = Ir×r Legyen U∗1 = Σ−1 V1∗ A∗ ; U1 = AV1 Σ−1 . U1 -et kiegészítjük m − r darab ortonormált
Rm -beli vektorral, ezek alkotják majd az U2 -t. Legyen az U = [U1 , U2 ], melyre U∗ U = Im×m . Ekkor: " U∗ AV =
U∗1 U∗2
# A
h
V1 V2
i
" =
U∗1 AV1 U∗1 AV2 U∗2 AV1 U∗2 AV2
#
" =
Σ 0 0 0
# , m×n
ahol U1 megadásából következik, hogy AV1 = U1 Σ; U1 és U2 ortogonalitásából pedig az, hogy U∗2 (AV1 ) = U∗2 (U1 Σ) = 0, illetve AV2 = 0 ∗ V2 = 0 miatt a 2. oszlop elemei nullák.
2.2.5. Megjegyzés. A szinguláris felbontásból keletkez® U és V mátrixok nem egyértelm¶en meghatározottak, például a nulla szinguláris értékekhez tartozó szinguláris vektorok sorrendje nem egyértelm¶. Láthattuk tehát, hogy minden mátrixnak létezik a szinguláris érték felbontása, most pedig az SVD néhány kivételes tulajdonságáról lesz szó. Az alábbiakban látni fogjuk, hogyan lehet általánosított inverzet készíteni SVD-b®l; hogyan lehet a kovariancia mátrixot közelíteni SVD-vel, illetve hogy egy tetsz®leges mátrix legjobb közelítését kettes illetve Frobenius normában egy, a szinguláris érték felbontásból kapott mátrix adja meg. Mindezek el®tt de niálom a mátrixnorma fogalmát, illetve annak el®bb említett két speciális esetét.
2.2.6. De níció
(Mátrixnorma). Az f : Rm×n → R függvényt mátrixnormának ne-
vezzük, ha 1. f(x) ≥ 0 (∀x ∈ Rm×n ), 2. f(x) = 0 ⇔ x = 0, 3. f(λx) = |λ|f(x) (∀x ∈ Rm×n , ∀λ ∈ R), 4. f(x + y) ≤ f(x) + f(y) (∀x, y ∈ Rm×n ). A norma szokásos jelölése: f(x) = ||x||.
2.2.7. De níció (Spektrál, vagy 2-es norma). Egy tetsz®leges A mátrix spektrálnormáját az alábbi két alakban tudjuk felírni: n
||A||2 = max ||Ax||2 = max ||x||2 =1
i=1
14
p
λi (A∗ A),
ahol ||x||2 az x vektor euklideszi normája és λi (A∗ A) az A∗ A mátrix sajátértékeit jelöli.
2.2.8. De níció (Frobenius norma). v uX n u m X t ||A||F = |aij |2 i=1 j=1
2.2.9. De níció (Ak közelít®mátrix). Legyen T
A = UDV =
r X
σi ui vTi
i=1
az A mátrix szinguláris felbontása, ahol r a mátrix rangját jelöli. Ekkor: Ak =
k X
σi ui vTi ,
i=1
minden 0-nál nagyobb, r-nél kisebb egész k-ra. A konstrukcióból nyilvánvaló, hogy Ak rangja k. Most nézzünk egy állítást ezekre a kitüntetett tulajdonságokkal bíró Ak mátrixokra.
2.2.10. Lemma. Ak sorai az A sorainak projekciói a k legnagyobb szinguláris értékhez tartozó szinguláris vektorok által kifeszített Vk altérre. Bizonyítás. Legyen a az A egy független sorvektora. Mivel minden i-re a vi ortonormált, a projekcióját Vk -ra a
Pk
T i=1 ha, vi ivi
képlet adja. Tehát a mátrixot, aminek
a sorai az A matrix soraiból vett projekció Vk -ra, a következ® alakba írható föl: Pk T i=1 Avi vi . Az így kapott kifejezést egyszer¶síthetjük: k X i=1
Avi vTi
=
k X
σi ui vTi = Ak .
i=1
Az alábbi tétel Carl Eckart és Gale Young nevéhez f¶z®dik, melyet és annak bizonyítását a [6] irodalomban publikálták 1936-ban. A tétel lényege, hogy megadja tetsz®leges r rangú mátrix 2-es és Frobenius normában vett legjobb k rangú közelítését (k < r).
15
2.2.11. Tétel (Eckart-Young, 1936). Legyen A ∈ Rm×n , rang(A) = r, k ∈ {1, . . . , r}. Az A mátrix legjobb k rangú közelítése 2-es, illetve Frobenius normában az Ak =
k X
σi ui vTi ,
i=1
mátrix. A két mátrix távolsága kettes normában: ||A − Ak ||2 = σk+1 ,
Frobenius normában pedig: v u r uX σ2i . ||A − Ak ||F = t i=k+1
Bizonyítás. A tételt több részfeladatra szedve, segédállítások segítségével bizonyítom. El®ször nézzük a Frobenius normabeli legjobb közelítés tulajdonságot.
2.2.12. Tétel. Legyen A és Ak olyan, mint ahogy fentebb de niáltuk. Ekkor tetsz®leges B legfeljebb k rangú mátrixra: ||A − Ak ||F ≤ ||A − B||F
Bizonyítás. Legyen B az a mátrix, ami a k, vagy kevesebb rangú mátrixok között minimalizálja az ||A − B||2F kifejezést. Legyen L a B sorai által kifeszített tér. Mivel
B minimalizálja az ||A − B||2F normát, muszáj, hogy B minden egyes sora a megfelel® sor projekciója legyen A-ból L-be, különben B azon sorát kicserélve a hozzá tartozó sor projekciójára A-ból L-be az L nem változna, (ebb®l kifolyólag a B rangja sem,) viszont csökkentené ||A − B||2F értékét. Mivel B mindegyik sora az A megfelel® sorának projekciója, ebb®l következik, hogy ||A − B||2F nem más, mint a távolságok négyzetösszege A soraitól L-ig. Viszont Ak az, ami minimalizálja a távolságok négyzetösszegét A soraitól tetsz®leges k dimenziós altérig, ebb®l pedig következik, hogy ||A − Ak ||F ≤ ||A − B||F . Beláttuk, hogy Frobenius normában a legjobban közelít® k rangú mátrixa A-nak
Ak . Most pedig vizsgáljuk ugyan ezt kettes normában, kezdve a közelítés A-tól vett távolságával. 16
2.2.13. Lemma. ||A − Ak ||2F =
Pr i=k+1
σ2i
Bizonyítás. Legyen A = UDV T és legyen Ak =
Pk
σi ui vTi . Induljunk ki a két mátrix normabeli eltéréséb®l! Mivel ||A − Ak ||F = ||UDV − Ak ||F = ||D − UT Ak V||F , legyen N = UT Ak V egy m × n-es k rangú mátrix. N megadásából következik, hogy: ||D −
N||2F
=
X
|Di,j − Ni,j | = 2
i,j
r X
i=1 T
|σi − Ni,i |2 +
X
i=1
|Ni,i |2 +
i>r
X
|Ni,j |2 .
i6=j
Ez az összeg akkor minimális, ha az N nem diagonális elemei mind nullák, illetve ha i > r esetén a diagonális elemek szintén nullával egyeznek meg. Nyilvánvalóan a Pr 2 i=1 |σi − Ni,i | minimumát úgy kapjuk, hogy Ni,i = σi minden i = 1, . . . , , k számra és Ni,j = 0, bármely más i, j indexre. Ebb®l következik, hogy
||A −
Ak ||2F
= ||D −
N||2F
=
r X
σ2i .
i=k+1
2.2.14. Lemma. ||A − Ak ||22 = σ2k+1 Bizonyítás. Legyen A =
Pr
σi ui vTi az A mátrix szinguláris felbontása. Ekkor Ak = P r T T i=k+1 σi ui vi . Legyen x tetsz®leges vektor. Fejezzük ki i=1 σi ui vi , A − Ak pedig Pr az x-et a v1 , v2 , . . . , vr vektorok lineáris kombinációjaként: x = j=1 αi vi , ahol α1 , α2 , . . . , αr tetsz®leges skalárok. Ekkor Pk
i=1
r r r X X X ||(A − Ak )x||2 = σi ui vTi αi vi = αi σi ui vTi vi = i=k+1
i=1
i=k+1
v r r X u X u t = αi σi ui = α2i σ2i . i=k+1
i=k+1
Az x vektor úgy maximalizálja az el®bbi értéket a ||x||22 =
Pr
α2i = 1 megkötés mellett, hogy az αk+1 -et 1-nek, a többi αi -t 0-nak választjuk, ebb®l pedig következik az állítás. i=1
Végül az állításunk igazolásához lássuk be, hogy kettes norma esetén is az Ak a legjobb közelítés.
17
2.2.15. Tétel. Legyen A és Ak olyan, mint ahogy fentebb de niáltuk. Ekkor tetsz®leges B legfeljebb k rangú mátrixra: ||A − Ak ||2 ≤ ||A − B||2
Bizonyítás. Ha az A mátrix rangja kevesebb, mint k, akkor az állítás nyilvánvaló, hiszen ekkor ||A − Ak ||2 = 0. Legyen tehát rang(A) szigorúan nagyobb, mint k. Az el®z® lemma alapján ||A − Ak ||22 = σ2k+1 . Ezek után tegyük fel, hogy létezik olyan
B legfeljebb k rangú mátrix, ami kettes normában jobban közelíti A-t Ak -nál. Ez azt jelenti, hogy ||A − B||2 < σk+1 . Legyen null(B) a B nulltere, tehát azon vektorok halmaza, melyekre teljesül, hogy Bv = 0. Mivel a B rangja legfeljebb k, a nulltér dimenziója legalább m−k. Legyen v1 , v2 , . . . , vk+1 az els® k+1 legnagyobb szinguláris vektora A-nak. A dimenziók száma miatt létezik egy z 6= 0 eleme az alábbi térnek: null(B) ∩ span(v1 , v2 , . . . , vk+1 ). Skálázzuk z-t úgy, hogy |z| = 1. Most megmutatjuk, hogy a z-re ami az A els® k + 1 szinguláris vektora által kifeszített térben van teljesül, hogy (A−B)z ≥ σk+1 . Hiszen ebb®l az következik, hogy az A − B kettes normában legalább σk+1 , ami ellentmond annak a feltevésnek, hogy ||A − B||2 < σk+1 . El®ször is tudjuk, hogy:
||A − B||22 ≥ ||(A − B)z||22 . Mivel z ∈ null(B),
||A − B||22 ≥ ||Az||22 . Viszont úgy választottuk meg z-t, hogy benne legyen a v1 , v2 , . . . , vk+1 vektorok konvex burkában:
n 2 n k+1 k+1 X X X X 2 T 2 T 2 2 T 2 2 ||Az||2 = σi ui vi z = σi (vi z) = σi (vi z) ≥ σk+1 (vTi z)2 = σ2k+1 . i=1
i=1
i=1
i=1
Ez azt jelenti, hogy ||A − B||22 ≥ σ2k+1 , ez pedig ellentmondás, hiszen feltettük, hogy
||A − B||2 < σk+1 . Beláttuk tehát, hogy ||.||2 , illetve ||.||F normában Ak az A legjobb közelítése, illetve az eredeti és a közelít® mátrix távolságát is megállapítottuk. Így tehát a segédtételek segítségével a 2.2.11 tételt is beláttuk.
18
Mint korábban említettem, az SVD-nek számos alkalmazása van. A korábbiakban láttuk a felhasználását az alacsonyabb rangú mátrixszal való közelítésben, most pedig az általánosított inverz konstrukciójában mutatom be a szerepét. Mindenek el®tt de niáljuk az általánosított inverz fogalmát.
2.2.16. De níció
(Általánosított inverz). Az A ∈ Cm×n mátrix Moore-Penrose-féle
általánosított, vagy pszeudo inverze az A+ ∈ Cn×m mátrix, ha 1. AA+ önadjungált, 2. A+ A önadjungált, 3. AA+ A = A, 4. A+ AA+ = A+ .
2.2.17. Tétel. Az A ∈ Cm×n mátrix Moore-Penrose-féle általánosított inverze egyértelm¶.
2.2.18. Tétel. A D ∈ Cm×n diagonális mátrix általánosított inverze a D+ ∈ Cn×m diagonális mátrix, ahol d+ ii =
0, 1 , dii
ha dii = 0, különben.
Tekintsük most az Ax = b egyenletrendszert, illetve az A szinguláris felbontását. Ekkor:
Ax = b ⇔ Dy = c,
x := Vy,
c = U∗ b,
tehát akkor oldható meg a lineáris egyenletrendszer, ha
• σi = 0 ⇒ ci = 0,
∀i ∈ (1, 2, . . . , n), az ilyen indexekre az yi -k tetsz®legesek.
• i > n ⇒ ci = 0. Célszer¶ a minimális normájú x megoldást kiválasztani, ami egyben a minimális normájú y := V ∗ x kiválasztását is jelenti, mivel V ortogonális mátrix. Ezt a megoldást úgy kapjuk, hogy az yi -t nullának vesszük a fentebb említett esetben, így y = D+ c, amib®l pedig x-re az alábbi egyenl®séget kapjuk:
x = VD+ U∗ b.
19
2.2.19. Tétel (Általánosított inverz el®állítása). A szinguláris felbontás felhasználásával az A ∈ Cm×n mátrix általánosított inverze a következ® alakban írható fel: A+ = VD+ U∗
Bizonyítás. Legyen rang(A) = p. Ahhoz, hogy az állítást igazoljuk, elég csak leellen®rizni, hogy a konstrukciónk kielégíti e az általánosított inverz de níciójában lév® feltételeket,.
" 1. A+ A = (VD+ U∗ )(UDV ∗ ) = VD+ U∗ UDV ∗ = V
Ip×p 0 0 0
# V ∗.
"
# I 0 p×p 2. AA+ = (UDV ∗ )(VD+ U∗ ) = UDV ∗ VD+ U∗ = U U∗ . 0 0 " # " # I 0 I 0 p×p p×p V ∗ = UD V ∗ = A. 3. 1. alapján: AA+ A = UDV ∗ V 0 0 0 0 " # " # I 0 I 0 p×p p×p U∗ = A + . 4. 2. alapján: A+ AA+ = VD+ U∗ U U∗ = VD+ 0 0 0 0
Az SVD el®állítása nehéz/számításigényes feladat (lásd kés®bb a Golub-Reinsch algoritmsnál). Viszont speciális esetben, ha az A mátrix teljes rangú és túlhatározott, akkor a feladat leegyszer¶södik:
2.2.20. Tétel. Legyen A ∈ Cm×n , m > n és r = rang(A) = n. Ekkor A+ = (A∗ A)−1 A∗ . A probléma abban a speciális esetben is könnyebb lesz, ha az A mátrix teljes rangú és alulhatározott:
2.2.21. Tétel. Legyen A ∈ Cm|timesn , m < n és r = rang(A) = m. Ekkor A+ = A∗ (AA∗ )−1 . Ebben a fejezetben megismerkedhettünk a szinguláris felbontás fontos tételeivel és tulajdonságaival, viszont a numerikus el®állításáról nem esett szó. A következ® fejezet erre a kérdésre ad választ. 20
2.3. Numerikus módszerek Ebben a fejezetben az szinguláris érték felbontás és az ahhoz kell® eljárások numerikus megoldásával fogok foglalkozni a [3] alapján. Látni fogjuk, miként tudunk numerikus módszerekkel egy mátrixot bidiagonális alakra hozni, milyen haszna van az SVD meghatározásában a Givens vagy Jakobi forgató mátrixokkal való szorzásnak illetve a Householder transzformációnak, majd végül ismeretem a Golub-Reinsch féle SVD algoritmust és annak különböz® tulajdonságát. Az algoritmusok mindegyikét leprogramoztam és az alkalmazások fejezetben fel is használom azokat. Az eljárások pszeudokódja megtalálható az A fejezetben. 2.3.1. Givens forgatások
A Givens forgatásnak komoly szerepe van a numerikus analízis számos területén, például sajátértékek meghatározásában, vagy a legkisebb négyzetek módszerében. El®ször kezdjük a 2 × 2-es eset vizsgálatával.
2.3.1. De níció
(Forgatás). Egy 2 × 2-es Q ortogonális mátrixot forgatásnak neve-
zünk, ha megadható az alábbi alakban: "
cos(θ) sin(θ) −sin(θ) cos(θ)
# ,
ahol θ a forgatás szöge.
2.3.2. De níció (Tükrözés). Egy 2 × 2-es Q ortogonális mátrixot tükrözésnek nevezünk, ha megadható az alábbi alakban: "
cos(θ) sin(θ) sin(θ) −cos(θ)
# .
Legyen y := QT x = Qx, ahol x tetsz®leges vektortérbeli elem. Ekkor az y kiszámítható, mint az x vektor tükrözése az S egyenesre, ahol " S = span
cos sin
#
θ 2 θ 2
.
Láttuk kétszer kettes esetben, hogyan is néz ki egy forgató mátrix, most pedig nézzük általános esetben! 21
2.3.3. De níció (Givens mátrix). i
1
.. . 0 G(i, k, θ) = ... 0 .. . 0
k
···
0
···
0
···
0
···
c
···
s
···
0 i .. , . 0 k .. . 1
...
.. .
.. . . . .
.. . .. .
···
−s
···
c
···
···
0
···
0
···
.. .
.. . . . .
.. .
(2.3)
ahol c = cos(θ), s = sin(θ), θ pedig a forgatás szöge.
2.3.4. Megjegyzés. • A G(i, k, θ) tisztán ortogonális. • A G(i, k, θ)T -val történ® szorzás egy óramutató járásával ellentétes, θ szög¶ forgatásnak felel meg az (i, k) koordinátasíkon. Az el®z®ekb®l az következik, hogy ha x ∈ Rn és y = G(i, k, θ)T x, akkor az y elemeit az alábbi zárt alakban kapjuk:
cxi − sxk , ha j = i yj = sxi + cxk , ha j = k xj , különben. Ebb®l az explicit felírásból látszik, hogy az yk kinullázható, ha
xi c= p 2 xi + x2k
és
−xk s= p 2 . xi + x2k
Így tehát egyszer¶ számítással eliminálni tudunk egy bizonyos elemet egy vektorban a Givens forgatás segítségével. A kés®bbiekben az SVD során ki fogjuk használni a forgatás ezen el®nyös tulajdonságát. El®bb viszont nézzünk egy másik fontos eljárást, amely segítségével megkaphatjuk tetsz®leges mátrix bidiagonális alakját.
22
2.3.2. Bidiagonalizálás Householder transzformációval
2.3.5. De níció (Householder-transzformáció). Legyen A ∈ Rn×n reguláris mátrix és ||.||2 az euklideszi norma. Householder-transzformációnak hívunk egy olyan A → QA leképezést, ahol H szimmetrikus ér ortogonális (komplex esetben hermitikus és unitér)
mátrix, továbbá az alábbi speciális alakban írható fel: H = I − cuuT ,
ahol I := In×n az egységmátrixot jelöli, u tetsz®leges vektor, u 6= 0, és c :=
2 . ||u||2
Az el®z®ek alapján tehát Householder-mátrixnak nevezünk egy olyan Rn×n beli mátrixot, melyre teljesül, hogy H(v) = I − 2vvT , ahol v ∈ Rn és a vektor egységnyi hosszú. Mértanilag a H(v)x = x − 2(vT x)v leképezés azt jelenti, hogy meghatározzuk
x tükörképét egy olyan síkra, melyre mer®leges v.
2.3. ábra. A Householder leképezés nem más, mint az x vektor tükrözése egy v-re mer®leges síkra.2
Most pedig nézzük a Householder-mátrixok tulajdonságait:
2.3.6. Állítás. 1. HT = H, azaz szimmetrikus, 2. H2 = I, azaz ortogonális, 3. H(v)v = −v, 4. H(v)y = y,
∀y ⊥ v.
2 Az ábra megtalálható a [2] könyvben.
23
A Householder-mátrixszal való m¶veleteknek az is nagy el®nye, hogy nem kell el®állítani a H(v) mátrixot, anélkül alkalmazhatjuk vektorokra:
H(v)x = (I − 2vvT )x = x − 2vvT x. Ennek a meghatározásához elegend® O(n2 ) lépés, szemben a mátrixszal való szorzással, amihez O(n2 ) m¶velet szükséges. A hatékonyság mellett sok mindenre használható a Householder-transzformáció, például ki tudjuk számolni vele egy mátrix QR-felbontását, bidiagonalizálni tudunk vele, illetve a számos tulajdonsága miatt egyéb tükrözési eljárásokra is alkalmas. A következ® tétel azt mondja ki, hogy két tetsz®leges egyforma hosszú vektor egymásba vihet® Householder-tükrözéssel.
2.3.7. Tétel. Legyen a, b ∈ Rn , a 6= b és ||a||2 = ||b||2 6= 0. Ekkor, ha v=
a−b ||a − b||2
⇒
H(v)a = b.
Most pedig térjünk rá a Householder-féle bidiagonalizálásra. El®ször is de niáljuk a bidiagonális mátrix fogalmát!
2.3.8. De níció
(Bidiagonális mátrix). Egy bidiagonális B ∈ Rm×n , m, n ∈ N+
mátrix alatt egy olyan mátrixot értünk, aminek csak a f®átlójában és a f®átló feletti átlóban található nem nulla érték. Legyen A ∈ Rm×n és m ≥ n. Ekkor A-ból a B bidiagonális mátrixot az alábbi formula adja:
d1 f1 0 d2 .. . . . .
T UB AVB = 0 ... 0 ...
0 f2 .. .
0
...
0 0 .. .
.. . mxn dn−1 fn−1 ∈R . dn
(2.4)
Ebben a formulában mind UB -t, mind pedig VB -t Householder-mátrixok szorzataként kapjuk, azaz:
UB = U1 · · · U n
és
VB = V1 · · · Vn−2 .
Most nézzük meg a bidiagonalizálást szemléletesen, egy 4 × 4-es esetre. 24
× × × × ×
× × × × ×
× × × × ×
× × × × ×
× 0 0 0 0
× × × × ×
0 × × × ×
0 × × × ×
× × 0 0 0 × × 0 U3 0 0 × × → 0 0 × × 0 0 × ×
U1 →
U2 →
× 0 0 0 0
× × × × ×
× × × × ×
× × × × ×
V1 →
× × 0 0 0 × × × V2 0 0 × × → 0 0 × × 0 0 × ×
× × 0 0 0 × × 0 U4 0 0 × × → 0 0 0 × 0 0 0 ×
× × 0 0 0 × × 0 0 0 × × 0 0 0 × 0 0 0 0
(2.5)
Ebben az el®állításban az Ui és Vi mátrixokkal ha megszorozzuk az eredeti mátrixot (balról UTi -tal, jobbról Vi -vel), akkor a megfelel® sor illetve oszlop f®átló utáni második elemt®l kezdve illetve f®átló alól az elemek nullára változnak. Alapvet®en az Uk -k felel®sek a k. oszlopba bevezetett nullákért, a Vk -k pedig a k. sorba bevezetettekért. A bidiagonális alak eléréséhez 4mn2 − 4n3 /3 lépésre van szükség. A teljes algoritmus megtalálható az A fejezet 2 részében. 2.3.3. Golub-Reinsch algoritmus
Ebben a fejezetben elérkeztünk a Golub-Reinsch féle SVD algoritmushoz. Az algoritmus célja egy tetsz®leges A mátrix szinguláris érték szerinti felbontásának meghatá-
25
rozása. Az els® lépés, hogy írjuk fel A-t bidiagonális alakban:
" UTB AVB =
B 0
#
d1 f1 0 0 d2 f2 .. . . . . . . .
B= 0 ... 0 ...
... ..
0 0 .. .
.
dn−1 fn−1 dn
∈ Rmxn .
A fennmaradó probléma tehát kiszámolni a B mátrix szinguláris felbontását. Ehhez alkalmazzuk egy implicit-shift QR lépést a T = BT B tridiagonális mátrixra.
2.3.9. De níció (Implicit szimmetrikus QR lépés Wilkinson shifteléssel). Adott egy T ∈ Rn×n redukált szimmetrikus tridiagonális mátrix. A de niált lépés átalakítja T -t egy ZT TZ mátrixszá, ahol Z = G1 · · · Gn−1 a Givens forgatások eredménye azzal a tulajdonsággal, hogy ZT (T − µI) fels® háromszög mátrix, µ pedig a T mátrix f®átlójára illeszked® 2x2-es záró részmátrixának T (n, n) = tnn -hez közelebbi sajátértéke. Legyen tehát a BT B = T legalsó, f®átlóra illeszked® 2x2-es részmátrixának, azaz
" T (n − 1 : n, n − 1 : n) =
d2n−1 + f2n−2 dn−1 fn−1 dn−1 fn−1 d2n + f2n−1
# (2.6)
mátrixnak a d2n + f2n−1 -hez közelebbi sajátértéke µ. Most számítsuk ki a c1 = cos(θ1 ) és s1 = sin(θ1 ) értékeket aszerint, hogy
"
c1 s1 −s1 c1
#T "
d21 − µ d1 f1
#
" =
x 0
# ,
ahol x egy nem nulla elemet jelöl. Legyen továbbá G1 = G(1, 2, θ1 ). A további
G2 , . . . , Gn−1 forgatómátrixot aszerint számoljuk ki, hogy ha Q = G1 · · · Gn−1 , akkor QT TQ tridiagonális és Qe1 = G1 e1 . Szorozzuk most be a B mátrixot a G1 forgató mátrixszal! A keletkezett BG1 mátrixba bekerül a f®átló alá egy nem nulla elem.
26
Nézzük az alábbi egyszer¶, 6 × 6-os esetet.
B ← BG1 =
× × 0 0 0 0 + × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 ×
(2.7)
Ehhez pedig ki tudjuk számolni az U1 , V2 , U2 ,. . . , Vn−1 és Un−1 mátrixokat, hogy segítségükkel eltüntessük a nem kívánt elemet.
B ← UT1 B =
B ← BV2 =
T B ← U2 B =
× × + 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 ×
× × 0 0 0 0 0 × × 0 0 0 0 + × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 ×
× × 0 0 0 0 0 × × + 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 × × 0 0 0 0 0 ×
(2.8)
És így tovább. Az UTi mátrixokkal balról, míg a Vi -kkel jobbról szorozva egy új B
27
bidiagonális mátrixot kapunk, amit az alábbi módon lehet meghatározni: T BV. = (UT · · · UT )B(G1 V2 · · · Vn−1 ) = U B 1 n−1 Mivel minden Vi -re teljesül, hogy Vi = G(i, i + 1, θi ) (i = 2, . . . n − 1) ebb®l követ 1 = Qe1 . Az implicit Q tétel miatt feltehetjük, hogy V és Q teljesen kezik, hogy Ve
megegyezik. Így implicit át tudjuk alakítani T -t T -vé, ha direktbe csak a B mátrixszal dolgozunk. Persze ahhoz, hogy tartsuk ezeket az állításokat, szükségszer¶, hogy a tridiagonális mátrixunk redukált legyen. Ez azt jelenti, hogy meg kell néznünk, mi történik, ha a B valamely diagonális vagy szuperdiagonális eleme nulla. Ha a f®átló feletti k. elem 0 (fk , k < n), akkor a B felbontható
"
B1 0 0 B2
#
alakba, ahol B1 k × k-as, B2 pedig n − k × n − k-as bidiagonális mátrix, így az eredeti feladat két részre bomlik. Ha valamely digonális dk elem lenne 0, akkor viszont a mellette lév® diagonális elem kinullázható Givens forgatásokkal, így szintén két részfeladatra bomlik az eredeti probléma. A Golub-Reinsch algoritmus egyik legmeghatározóbb része a Golub-Kahan lépés. Egy bidiagonális B mátrixból indul ki, aminek sem a diagonális, sem a f®átló feletti átló nem tartalmaz nullát. Az algoritmus a (2.6) egyenl®séggel de niált µ sajátértékb®l indul ki, majd n − 1 lépés követi, a k. lépésben kiszámolja az aktuális G(k, k + 1, θ) forgatómátrixot, szorozza a jelenlegi B-t jobbról vele, majd újraszámolja a G-t és balról is megszorozza B-t a Givens mátrix transzponáltjával. Az Golub-Kahan lépés = UB V T mátrixot konstruál, melynek a f®átlón kívüli elemeinek a a B-b®l egy B négyzetösszege csökken a B-beli négyzetösszeghez képest.Hatásosan implementálva az eljárás 30n lépésben tárolja a B ritka mátrix elemeit, 6mn lépés szükséges U, 6n2 lépés pedig V meghatározásához. Tipikusan a fenti eljárás néhány egymásutáni alkalmazása során az fn−1 , azaz a legalsó szuperdiagonális elem elhanyagolhatóvá válik. Az elhanyagolhatóság határát az alábbi egyenl®tlenségek adják.
|fi | ≤ (|di | + |di+1 |) |di | ≤ ||B||, 28
(2.9)
ahol rögzített, alacsony érték¶ szám. A Golub-Reinsch SVD a fejezetben de niált algoritmusokat ötvözi egyben. A bemenetünk egy tetsz®leges A mátrix. Els® lépésben ezt bidiagonalizáljuk. Majd amíg nem lesz a kapott B mátrix diagonális, a következ®ket tesszük: 1. Végignézzük a szuperdiagonális elemeket és ha a (2.9) feltétel alapján elhanyagolhatónak bizonyul, akkor nullává alakítjuk. B11 0 0 2. Ezután B = 0 B22 0 alakban írható fel, ahol B33 q × q-s diagonális
0 0 B33 mátrix, B22 f®átló feletti átlójában pedig van nem nulla elem. 3. Ha B22 bármely diagonális eleme 0, akkor a mellette lev® szuperdiagonális elemet elimináljuk, ha nincs ilyen 0 érték¶ diagonális elem, akkor pedig alkalmazzuk
B22 -re a Golub-Kahan lépést. Végül megkapjuk az U és V ortogonális, illetve D diagonális mátrixok egy jó közelítését.
29
3. fejezet
Alkalmazások
Az el®z® fejezetben megismerkedtünk a szinguláris érték szerinti felbontással, annak tulajdonságaival és egy lehetséges el®állításával, ebben a fejezetben pedig néhány alternatív alkalmazását fogom szemléltetni.
3.1. Tömörítés Tömörítés során a célunk nem más, mint az adatok feldolgozása oly módon, hogy lehet®leg minél kevesebb tárhelyet kelljen használnunk. Ez általában akkor valósítható meg, ha az adat redundáns.
3.1. ábra. Az eredeti kép.
Mint említettem, egy képet meg lehet feleltetni egy Rm×n -es mátrixnak, ahol a mátrix i-edik sorának j-edik eleme reprezentálja a kép adott pixelének a színét/intenzitását.
30
Legyen A mátrix a 3.1 ábra reprezentációja:
A :=
a1,1 a1,2 · · · a1,n a2,1 a2,2 · · · a2,n .. .. . . .. . . . . an,1 an,2 · · · an,n
.
Az egyszer¶ség kedvéért szürkeárnyalatos képet vizsgálunk, azaz ai,j ∈ (0, 1, . . . , 255). Alakítsunk 1×2-es blokkokat az A mátrix szomszédos elemeib®l, majd pedig ábrázoljuk azokat koordináta rendszerben. Egy pont els® koordinátája tehát a páros baloldali értéke, a második koordináta pedig a páros jobboldali értéke. Általában egy képen az egymás melletti pixelek nagy valószín¶séggel azonos szín¶ek, így a pontok nagy része az identitás környékére fog esni. Mivel a szomszédos pixelek er®sen korrelálnak, 250
100 80
200
60 40
150
20 0
100
-20 -40
50
-60
0
-100
-80
0
50
100
150
200
250
0
(a) standard koordináta rendszerben
50
100
150
200
250
300
350
(b) transzformált bázissal
3.2. ábra. Az ábrázolt pontok különböz® koordináta rendszerekben. így célszer¶ nem a hagyományos, [0, 1], [1, 0] bázisú koordináta rendszerben ábrázolni a pontokat, hanem valamilyen másik transzformált koordináta rendszerben. Látszik, hogy a hagyományos kanonikus bázist használva mind a páros, mind pedig a páratlan pixelek kihasználják a teljes [0, . . . , 255] tartományt. Viszont ha az elforgatott rendszert nézzük, észrevehet®, hogy a páros index¶ pixelek lényegében a [−40, 40] intervallumra koncentrálódnak. A transzformáció nem más, mint a 2.3.1 fejezetben taglalt Givens forgató mátrix. Alapvet®en a transzformációval még nem érnénk el tömörítést, azonban, mivel a páros index¶ pontok a [−40, 40]-be esnek, így ezeket 6 biten tárolva olyan közelítést kapunk, melynek `2 normában hibája kicsi. Ezt a két dimenziós módszert ki tudjuk terjeszteni hasonlóan n dimenzióra is. Ekkor egy pont nem két változóból fog állni (egymás melletti pixelek), hanem maga a kép lesz. 31
Ezt úgy érjük el, hogy a reprezentáló mátrixot vektorizáljuk, azaz sorfolytonosan egy
nm × 1-es vektorba illesszük az értékeket. Legyen M a képek száma, N = 640 ∗ 480 a képek mérete. Ekkor MN adatot kéne eltárolnunk. Ehelyett alkalmazhatjuk a következ® lépéseket: válasszunk ki néhány képet, ezen képek halmazát nevezzük tanítóhalmaznak, melynek elemszáma k. Megjegyezzük, hogy akár az összes képet választhatnánk, ekkor a teljes adatbázisunk tanítóhalmaz lenne. A 2.1.2 jelöléseket használva legyen a tanítóhalmaz i. vektorizált képe xi ∈ RN (1 ≤ i ≤ k). Az ezekb®l készített RN×k mátrix pedig X = [x1 , x2 , . . . , xk ]. Alkalmazzuk a 2.1.2 fejezetben leírt f®komponens analízis technikáját az X mátrixra. Ebben az esetben a kovariancia mátrix közelítéséhez használt XXT ∈ RN×N mátrix olyan nagy, hogy a szorzást el sem tudjuk végezni. Ezért a PCA-t a 2.2 fejezetben bemutatott szinguláris érték szerinti felbontáson keresztül implementáljuk. Az így nyert SVD-ben vegyük a keletkez® l darab legnagyobb szinguláris értékéhez tartoP zó szinguláris vektort. A 2 fejezet alapján X ≈ li=1 hXi , Vi iVi , ahol V az SVD-b®l kapott nagyobbik ortogonális mátrix. Ekkor ez az összeg nem más, mint egy projekció span(V1 , V2 , . . . , Vl )-re. Nekünk a kovariancia mátrixra van szükségünk, viszont az el®z® fejezetben azt is beláttuk, hogy a kovariancia mátrixot jól tudjuk közelíteni az k1 XXT -tal. Így tehát összesen az NM adat helyett elegend® l darab szinguláris vektort, illetve Ml darab hXi , Vi i együtthatót tárolni. Nézzünk egy példát a [8] honlapról felhasznált Yale adatbázisból,amely 165 darab 64 × 64 pixeles képet tartalmaz. Ebb®l legyen mondjuk 30 elem¶ a tanítóhalmaz és nézzük az egyik kép közelítését
3.3. ábra. Az els® 20 saját arc.
32
úgy, hogy mindössze 20 szinguláris vektort használunk. A 3.3 ábrán láthatjuk az els® 20 saját arc ábráját (angolul: eigenfaces). A 3.4 ábrán láthatóak az eredeti kép közelítései 1,2,. . . ,20 szinguláris vektor esetén. Érdemes meg gyelni, hogy az els® 5 szinguláris vektor használatánál még igen nagy a hiba, szinte teljesen más alakot látunk, mint a kés®bbi képeken. Az els® 10-nél valamivel kisebb a hiba, már kezd
3.4. ábra. A kép közelítése 1. . . 20 szinguláris vektor használatával.
hasonlítani a közelítés az eredeti képhez. Viszont 15 sajátvektortól kezdve már szinte teljesen azonos képeket kapunk. Tehát ha tömöríteni szeretnénk, akkor elég csak a 20 legnagyobb szinguláris értékhez tartozó szinguláris vektort eltárolni, illetve minden képhez a skalárszorzatból keletkezett együtthatókat.
3.2. Adatsimítás A zaj az eredetit®l eltér® frekvenciájú és intenzitású jelek összessége, melyek gondot okozhatnak, hiszen általában távol vannak a valós értékekt®l. A zajok eltüntetésére több módszer is van, az egyik ilyen az adatsimítás. Ez egy eljárás, mely a kiugró értékeket a szomszédos adatok alapján módosítja. Nézzünk egy példát adatsimításra! Adott egy EKG felvétel, amely 60 szívütést tartalmaz, ezt szeretnénk zajmentesíteni. A tömörítés problémájához hasonlóan, használni tudjuk az SVD-t. Most egy dimenziós adatsorozatokkal dolgozunk, tehát nem kell vektorizálni. Az X mátrix oszlopai az EKG jelek szívütései. Tekintsünk két különböz® esetet egy-egy szívverésre! Az els® 33
10
10
8
8
6
6
4
4
2
2
0
0
-2
-2 0
50
100
150
200
250
300
350
400
0
(a) A jel közelítése 3 sajátvektorral.
50
100
150
200
250
300
350
400
(b) A jel közelítése 10 sajátvektorral.
3.5. ábra. EKG jelek simítása. esetben közelítsünk 3 sajátvektorral, a második esetben pedig 10 sajátvektorral. A 3.5 ábrán kékkel van jelölve az eredeti, zajos jel, pirossal pedig a közelítésünk. Látszik, hogy a kék görbe tele van apró zajokkal, "recékkel". Meg gyelhet®, hogy már 3 sajátvektort használva a jel kell®en sima, csak néhol van az eredeti gra kontól eltérés, 10 sajátértékre viszont már szinte teljesen pontos a közelítés. Tehát ebben az esetben simítás mellett tömörítést is kaptunk. A pontos alkalmazásáról és egyéb részletekr®l a [10] irodalomban olvashatunk.
3.3. Arcképek osztályozása Az arckép osztályozási feladat a következ®képpen írható fel: Adott egy arcképadatbázis és szeretnénk eldönteni, hogy egy adott kép benne van e az adatbázisban. Például adott egy kaszinó tagjainak az arcképeit tartalmazó adatbázis és meg szeretnénk határozni, hogy egy jelenleg megjelen® személynek van e tagsága a klubba, vagy sem. El®ször a pixelenkénti összehasonlítás jutna az eszünkben, viszont ez több okból sem járható út: el®ször is rettent®en m¶veletigényes, hiszen minden kép minden pixelével össze kéne hasonlítani, ebb®l kifolyólag lassú is, az esetleges kaszinótagunknak kiábrándítóan sokat kéne várnia. Másrészr®l pontosnak sem mondható, hiszen más fényviszonyok mellett teljesen eltér®ek lesznek a pixelek, illetve azt is gyelembe kell vennünk, hogy a képet az eredetit®l különböz® pozícióban is rögzíthetjük. A megoldást a szinguláris felbontás szolgáltatja: transzformáljuk másik térbe az adatokat, és vegyünk a kép szinguláris érték szerinti felbontását. A megkapott szinguláris felbontás legnagyobb k darab szinguláris értékéhez tartozó szinguláris vektorok által 34
meghatározott hxi , vi i skalárszorzatokat hasonlítjuk össze az adatbázisunk képeinek hasonló skalárszorzataival. Ha a skalárszorzatok eltérése, azaz a hiba kettes normában valamilyen rögzített értéknél kisebb, akkor elfogadjuk a képet, tehát a vendégünk nyugodt szívvel léphet a játékok mezejére, ha viszont nagyobb annál, akkor elutasítjuk azt, tehát sajnálatos módon távoznia kell. Az eljárást dimenziócsökkentésnek nevezzük, mivel az adatot alacsonyabb dimenzióba transzformáltuk.
3.4. Kitekintés Ebben a fejezetben az SVD egy alternatív el®állításáról lesz szó. Az úgynevezett 2DSVD egy kétdimenziós változata a szinguláris felbontásnak. Hasonlóan de niálhatjuk a kétdimenziós SVD-t 2D-s objektumokra, mint ahogyan azt tettük az SVDvel a 1D-s vektorokra. Legyenek X1 , . . . , Xn az eredeti képeket reprezentáló Rr×c -es mátrixok. De niáljuk az átlagolt sor-sor és oszlop-oszlop kovariancia mátrixokat a következ®képp:
n X T (Xi − X)(X F= i − X) , i=1
G=
n X
T (Xi − X), (Xi − X)
i=1
az Xi -k mintaátlaga. Ekkor F az XXT -nak, G pedig XT X-nek felel meg. Legyen Ahol X
Uk az F els® k legnagyobb sajátvektorát tartalmazó mátrix, Vs pedig tartalmazza a G els® s darab legnagyobb sajátvektort. Ekkor: F=
r X
λ` u` uT` ,
Uk ≡ (u1 , . . . , uk );
`=1
G=
c X
ξ` v` vT` ,
Vs ≡ (v1 , . . . , vs ).
`=1
Ezek alapján de niáljuk következ®képpen az SVD kiterjesztését kétdimenziós objektumokra:
e = Uk M i V T , X s
Mi = UTk Xi Vs ,
i = 1, . . . , n.
A standard SVD esetében az Uk biztosítja a közös alteret, amelyre végezzük az egydimenziós vektorok projekcióit. A 2DSVD esetében hasonlóan, az Uk és Vs mátrixok jelentik a projekciós alterét a 2D-s objektumoknak. Fontos észrevétel, hogy a kétdimenziós szinguláris felbontásban szerepl® Mi mátrixnak nem kell diagonálisnak 35
lennie, szemben az SVD-vel, ahol a Σk diagonális volt. Kísérletek igazolják, hogy az arcképek osztályozása esetén a 2DSVD hasonlóan hatékony tud lenni, mint az SVD. Viszont nagy el®nye, hogy sokkal kevesebb tárhelyet igényel, ugyanis míg a szinguláris felbontás esetén M darab m ∗ n pixel¶ képnél M × mn nagyságú volt a probléma, a 2DSVD-nél viszont ez a probléma redukálódik m × n nagyságúra. Megjegyzem, hogy a 2.2 fejezetben ismertetett Eckart-Young tételek analóg változatai beláthatóak a 2DSVD-re is. B®vebben lásd a [9] cikkben.
36
A. függelék
Algoritmusok
Algorithm 1 Householder 1: function Householder(x) 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:
n = length(x) σ = x(2 : n)T x(2 : n) v = [1, x(2 : n)]T if σ = 0 then β=0 elsey4
else
||v||2 = −σ/(||x|| + µ)
end if
β = 2||v||22 /(σ + ||v||22 ) v = v/||v||2
end if return v, β end function
37
Algorithm 2 Bidiagonalization 1: function Bidiag(A) 2: for j = 1 : n do 3: 4: 5: 6: 7: 8: 9: 10: 11: 12:
[v, β] = Householder(A(j : m, j)) A(j : m, j : n) = (Im−j+1 − βvvT )A(j : m, j : n) A(j + 1 : m, j) = v(2 : m − j + 1) if j ≤ n − 2 then [v, β] = Householder(A(j, j + 1 : n)T ) A(j : m, j + 1 : n) = A(j : m, j + 1 : n)(In−j − βvvT ) A(j, j + 2 : n) = v(2 : n − j)T
end if end for end function
Algorithm 3 Givens Rotation 1: function Givens(a, b) 2: if b = 0 then 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16:
c = 1; s = 0
else if |b| > |a| then r = −a/b √ s = 1/ 1 + r2 c = sr
else
r = −b/a √ c = 1/ 1 + r2 s = cr
end if end if return c, s end function
38
Algorithm 4 Golub-Kahan Step 1: function GKS(B) 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:
Let µ be the eigenvalue of the trailing 2-by-2 submatrix of T = BT B that is closer to tnn . y = t11 − µ z = t12 for k = 1 : n − 1 do Determine c = cos(θ) and s = sin(θ) from y and z B = BG(k, k + 1, θ) y = bkk z = bk+1,k Determine c = cos(θ) and s = sin(θ) from y and z B = G(k, k + 1, θ)T B if k < n − 1 then y = bk,k+1 z = bk,k+2
end if end for end function
Algorithm 5 Golub-Reinsch SVD 1: function SVD(A) 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20:
[U, B, V] = Bidiag(A) while q = n do for i = 1 : n do Set bi,i+1 to zero if |bi,i+1 | ≤ (|bii | + |bi+1,i+1 |)
end for
Findthe smallest p and the largest q such that if B11 0 0 0 B = 0 B22 0 0 B33 then B33 is diagonal and B22 has nonzero superdiagonal. if q < n then if any diagonal entry in B22 is zero then zero the superdiagonal entry in the same row
else
apply Golub-Kahan step to B22 , B = diag(Ip , U, Iq+m−n )T B diag(Ip , V, Iq )
end if end if end while return U, B, V end function
39
Irodalomjegyzék
[1] Gergó L., Numerikus módszerek, ELTE Eötvös Kiadó, Budapest, 2010. [2] G. Stoyan, Takó G., Numerikus módszerek 1., Typotex, Budapest, 2005. [3] G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, USA, 1996. [4] K. R. Rao., The Transform and Data Compression Handbook, CRC Press, New York, USA, 2001. [5] R. O. Duda, P. E. Hart, D. G. Stork, Pattern Classi cation, Second Edition, Wiley-Interscience Press, USA, 2000. [6] C. Eckart, G. Young, "The approximation of one matrix by another of lower rank",
Psychometrika, vol. 1, no. 3, pp. 211-218, Chicago, Illinois, 1936. [7] P. N. Tan, M. Steinbach, V. Kumar, "Introduction to Data Mining", Addison-
Wesley Longman Publishing Co., Inc., Boston, Massachusetts, 2005. [8] Four
face
databases
in
matlab
format,
Zhejiang
University,
China.
http://www.cad.zju.edu.cn/home/dengcai/Data/FaceData.html [9] C. Ding and J. Ye, "Two-Dimensional Singular Value Decomposition (2DSVD) for 2D Maps and Images", Proc. SIAM Int'l Conf. Data Mining (SDM'05), pp. 32-43, October 3, 2004. [10] F. Castells, P. Laguna, L. Sörnmo, A. Bollmann and J. M. Roig, "Principal Component Analysis in ECG Signal Processing", EURASIP Journal on Advances
in Signal Processing, vol. 2007, pp. 98-119, 2007.
40
