Axel Ruhe
4 februari 2005		 Tillämpade numeriska metoder II, våren 2005, 2D1250 NADA bild

LABORATION 2

Utdelas 4 februari, in 25 februari.

Minsta kvadratproblem

1. Anpassning av polynom till dataserie

Polynom:

Formulera det minsta kvadratproblem man får, när man skall anpassa ett polynom till en serie mätdata y. Ta t ex m=101 ekvidistanta punkter t_i i intervallet [0<=t<=1] och försök med polynom av gradtal upp till n-1 där n=15. Vi får en m*n matris A.

Ta en titt på dess singulära värdesuppdelning (SVD)! Utred följande:
  1. Kolonnerna i A är ju de första potenserna. Plotta ut deras värden som funktioner av t. Se hur lika de blir, mycket riktigt är kolonnerna i A nära på linjärt beroende. 
  2. Jämför med de första kolonnerna i U. Dessa är ortogonala vektorer och liknar ortogonala funktioner, jämför gärna med trigonometriska funktioner.
  3. De singulära värdena blir strax små, man behöver ta till semilogy för att plotta ut de minsta. 
Samma fråga för QR faktoriseringen av A:
  1. Nu blir kolonnerna i Q ortogonala polynom, jämför med Legendrepolynomen om Du hört talas om dessa.
  2. Diagonalelementen i R avtar, men ej monotont. Är de större eller mindre till absolutbeloppet än motsvarande singulära värden?

Dataserier:

Använd följande tre dataserier:
a. y=exp(t), detta är en hel funktion som kan approximeras med polynom,  bättre desto högre gradtalet är
b. y=sqrt(t), denna är singulär i nollan och låter sig inte så gärna approximeras med ett polynom
c. y=exp(t).*(1+e*randn(m,1)), en slumpmässig störning av en hel funktion.

Jämför trunkerad QR och trunkerad SVD:

Det är inte lämpligt att ta med alla singulära värden eller alla kolonner i Q när man tar fram en lösning. Tar man för många, ökar normen på lösningen, utan att man förbättrar anpassningen något nämnvärt. Det kan studeras genom att man ser på komponenterna i det transformerade högerledet, U'*y resp Q'*y. När de blir större än motsvarande singulära värden eller diagonalelement i R, är det ingen ide att följa med längre.

När Du bestämt en rang r, räkna fram x med utnyttjande av de r första singulära värdena resp r första kolonnerna i Q. QR lösningen blir ett polynom av gradtalet r-1 medan SVD lösningen blir ett polynom av fullt gradtal som är en linjärkombination av de första kolonnerna i U. Ibland ger QR bäst anpassning, ibland SVD. Hitta gärna ett fall där QR är bäst och ett där SVD är bäst.

Ta några olika värden av störningen e i serien c, och se för hur stora e man kan känna igen exponentialfunktionen. Plotta data y och approximerande kurvor A*x. Plotta residualen y-A*x och se om den uppför sig som ett polynom dvs byter tecken lagom många gånger (vad är lagom här?) eller består av slumpvisa störningar.

Regularisering, L-kurvan:

En tredje metod att få fram en lösning med mindre norm är regularisering enl Tichonov. Då löser man minimeringsproblemet

min ||Ax-b||2+µ2||x||2
   x
Det görs lätt genom att man löser ett minsta kvadratproblem med matrisen [A ; mu*eye(n)] och högerledet [b ; zeros(n,1)], för att använda Matlabs notation. Ju större värde på µ desto större vikt lägger man vid att få liten norm ||x|| på bekostnad av en större residual ||r||.

Plotta en L kurva med dubbellogaritmisk skala loglog! Visa lösningsnormen ||x|| som funktion av residualnormen  ||r||, där r står för residualen r=A*x-b. Jämför punkterna som svarar mot olika ranger r, från SVD, gradtal k från QR och värden av regulariseringsparametern µ. Om man varierar µ så får man en kurva som förhoppningsvis ser ut som ett L, medan punkterna för de olika värdena på r och k ligger ovanför. Välj µ jämnt fördelade över ett intervall mellan det minsta och största singulära värdet på A.

Titta på L kurvan och se om det finns någon bra krök. Det är dataserien c, som mest liknar ett rimligt antagande, om en kurva med slumpvisa fel. Det gäller att anpassa till kurvan och inte till störningarna! Serien a är speciell, det går bara bättre ju högre gradtal vi tar. Serien b är lite besvärligare.

2. SVD i bildbehandling: Att känna igen ansikten

Människans förmåga att känna igen ett ansikte är förunderlig. Ett litet fotografi på ett körkort skall räcka för att avgöra om det är rätt person som sitter vid ratten. Och det är inte mycket som skiljer olika ansikten åt! Finns det något medelansikte som liknar alla och ingen?

Vi skall undersöka vad som händer när vi använder singulär värdesuppdelning (SVD) för att skilja allmänna och särskilda egenskaper hos ansikten.
Vi  har lånat ett par uppsättningar inskannade bilder, lagrade som Matlab matriser. Varje ansikte består av h*b pixlar.
Vi använder bara gråskala, så vi får en m*n matris X, där var och en av n bilder bildar en kolonn av längden m=h*b. När vi skall plotta ansiktet så förvandlar vi kolonnen X(:,k) till en h*b matris FaceMatrix som plottas ut med kommandona

colormap gray
FaceMatrix=reshape(X(:,k),h,b);
imagesc(FaceMatrix)
axis image
  1. Läs in en datafil!
  2. De bilder vi får in består av heltal som anger gråskalan. Normalisera bilderna så att medelvärdet blir noll och kolonnerna får samma euklideska norm, ett.
  3. Gör SVD på den sålunda normaliserade matrisen:
    [U S V]=svd(X,0)
    Observera nollan! Den gör att man bara räknar fram de n första kolonnerna av U, annars får man en kvadratisk matris U som tar massor av plats.
  4. Plotta de singulära värdena, diag(S). Se hur snabbt de avtar! Hur många komponenter behöver man ha för att representera 70, 80 och 90% av den totala variansen hos datamaterialet? De singulära värdena är av ganska olika storlek, pröva logaritmisk skala semilogy och se om det blir bättre!
  5. Den första kolonnen i U blir ett slags medelansikte, de följande ger mera speciella egenskaper. Man brukar kalla dessa "eigenfaces". Titta på några av dessa! Hur länge liknar de ansikten över huvud taget? 
  6. Varje individuellt ansikte kan nu representeras som en linjär kombination av egenansiktena,
    FaceVector=U(:,1:r)*S(1:r,1:r)*V(k,1:r)'    
    ger en representation av individ k med de r första egenansiktena.Ta någon bild och se efter hur många termer r som behövs för att individen skall kännas igen! 
  7. Om det räcker med r termer så kan man göra ett porträttgalleri genom att spara de r kolonnerna hos U och n katalogvektorer, kolonnerna i matrisen
    Catalog=S(1:r,1:r)*V(:,1:r)';
    Undersök hur detta fungerar för identifikation av individer. Flera av bilderna i datamatrialet är på samma person. Se efter om dessas katalogvektorer ligger närmare varann än vektorer för bilder av olika individer.

Datafiler

mitfaces.mat innehåller två uppsättningar om n=34 ansikten, ganska små bilder h=51 och b=43, ger m=51*43=2193. Dessa bilder är mycket grovkorniga men noggrannt förbehandlade, man har sett till att alla blir lika placerade genom att ögonen kommer på samma plats på alla bilderna. Det finns två matriser smileFaces och neutralFaces och vi kan testa om man hittar rätt leende individ genom att jämföra dess bild med alla neutral ansiktena och se vilket som ligger närmast.

ucdfaces.mat innehåller en uppsättning om n=48 ansikten, lite större bilder, h=120 och b=128, m=15360. Dessa bilder har inte förbehandlats så mycket, bl a finns det med bakgrund här och var. I vissa fall är det flera bilder på en individ.

Använd bägge materialen och leta gärna reda på ett eget porträttgalleri!