Programmieraufgaben zu Blatt 6: QR-Zerlegung
Die Programmieraufgabe ist in mehrere kleine Teile zerlegt. Bevor es losgeht alles von Blatt 6 (diesen Link klicken) herunterladen und in den Numerik Ordner packen. Das Makefile hat sich nicht verändert sondern ist nur zur Sicherheit dabei...
In jeder Aufgabe ist erklärt:
- Was zu tun ist.
- Welche Datei ihr editieren müsst.
- Das ganze ist so gedacht: ihr lest eine Aufgabe aufmerksam durch und legt dann los. In jeder Aufgabe sollte dabei das erklärt sein, was ihr an C++ Kenntnissen zum lösen braucht.
- Wenn etwas unkalr ist: Email schreiben oder ins Forum posten. Ich werde diese Seite dann anpassen
Und los geht's:
Aufgabe 1
Das Hauptprogramm ist bereits fix und fertig:
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
DEVector alpha(3);
GEMatrix A(4,3);
A = 2, 1, 3,
4, 4, 9,
6, 9, 16,
8, 16, 24;
VectorView v = A(_,1);
alpha(1) = norm2(v);
cout << alpha << endl;
}
Was soll dieses Programm tun?
- Es soll die euklidische Norm der ersten Spalte von A berechnen und in alpha(1) speichern.
- Zeilen, Spalten sowie rechteckige Bereiche von Matrizen kann man mit dem so genannten Range-operator "ausschneiden". Der Range-operator ist hier der Unterstrich in A(_,1).
- Mit dem Unterstrich kann man Indexbereiche angeben. A(_,1) bezieht sich auf alle Elemente in
A mit
- beliebigem Zeilenindex
- mit Spaltenindex =1
- Noch etwas ist bei diesem Programm zu beachten:
Damit legt man fest, dass v sich auf die erste Spalte bezieht. Anschliessend kann v wie ein ganz normaler Vektor benutzt werden. Wird v allerdings verändert, dann wird damit auch die erste Spalte von A verändert! Deshalb die Bezeichnung "View". Dies drück aus, dass v eben keine Kopie der ersten Spalte ist, sondern sich direkt auf diese Spalte bezieht.
VectorView v = A(_,1);
Damit dieses Programm funktioniert muss aber noch die Funktion norm2 in der Datei blatt6.h implementiert werden:
template <typename VX>
double
norm2(const DenseVector<VX> &x)
{
// Hier Programm-Code einfuegen.
// statt 0 soll die euklidische Norm von x
// zurueckgegeben werden:
return 0;
}
Das sollte nicht weiter schwer sein, denn bekanntlich gilt 
. Im
schlimmsten Fall programmiert man das mit einer for-Schleife. Allerdings kann man
das Skalarprodukt in FLENS auch einfach mit x*x erhalten. Wenn man daraus noch die
Wurzel zieht, dann ist man schon fertig. Eine einzige Programmzeile sollte also reichen.
Dabei ist folgendes zu beachten:
- FLENS kennt keine Unterscheidung zwischen Zeilen- und Spaltenvektoren.
- Sind x und y zwei Vektoren, dann ist
immer das Skalarprodukt.
double alpha = x*y; - Möchte man statt dem Skalarprodukt aber das sogenannte Tensorprodukt
berechnen (wobei x ein Spaltenvektor ist), dann geht das nicht mit dem *-operator. Es sollte aber aus der Übung und den Tutorien klar sein, dass man das Tensorprodukt für die Householder-Transformation nicht berechnen muss!
Aufgabe 2
Auch folgendes Programm ist bereits fertig. In blatt6.h muss dafür nur noch die Funktion householder programmiert werden.
Ziel der Aufgabe ist, dass
berechnet wird. Der Vektor v muss also nicht normiert sein. Das Ergebnis der Aufgabe sollte aus den Übungen bekannt sein.
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
DEVector x(3);
x = 2, 1, 2;
DEVector v(3);
v = 2, -2, 2;
householder(v, x);
cout << "Hx = " << x << endl;
}
Damit alles funktioniert muss in blatt6.h die Funktion householder programmiert werden:
template <typename VV, typename VX>
void
householder(const DenseVector<VV> &v, DenseVector<VX> &x)
{
// Hier soll x mit H*x ueberschrieben werden,
}
Dabei muss man (noch darf man) weder eine Einheitsmatrix aufstellen noch weitere Vektoren anlegen. Auch dies kann mit einer einzigen Zeile programmiert werden!
Aufgabe 3
Aufgabe 3 ist eigentlich gar keine Aufgabe sondern nur ein Test! Sogar ein doppelter Test
- Für die Funktion householder
- und für euch!
Erklärt was diese Programm tut und erklärt die Ausgabe des Programm:
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
GEMatrix A(4,3);
A = 2, 1, 3,
4, 4, 9,
6, 9, 16,
8, 16, 24;
DEVector v(4);
v = A(_,1);
v(1) = v(1) + norm2(v);
int n = A.numCols();
for (int k=1; k<=n; ++k) {
VectorView x = A(_,k);
householder(v, x);
}
cout << "A = " << A << endl;
}
Hinweis: In diesem Programm steht
DEVector v(4); v = A(_,1);
Hier ist also v keine View. Deshalb enthät v eine Kopie der ersten Spalte. Und deshalb wird A auch nicht verändert, wenn v verändert wird.
Aufgabe 4
In dieser Aufgabe muss nur das Programm aufgabe4.cc veraendert werden.
Sei 
dann wählt man im ersten Schritt der QR-Zerlegung den
Vektor 
für die Householder Transformation so:
- Falls
setze Vektor 
- sonst
.
den ersten Einheitsvektor.
Das soll nun als Programm realisiert werden:
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
DEVector alpha(3);
GEMatrix A(4,3);
A = 2, 1, 3,
4, 4, 9,
6, 9, 16,
8, 16, 24;
VectorView v = A(_,1);
// Hier v veraendern. Und zwar so:
// Sei alpha die euklidische Norm der ersten Spalte von A
// falls A(1,1) > 0, dann
// setze v := v + alpha*e1
// sonst
// setze v := v - alpha*e1
int n = A.numCols();
for (int k=2; k<=n; ++k) {
VectorView x = A(_,k);
householder(v, x);
}
cout << "alpha = " << alpha << endl;
cout << "A = " << A << endl;
}
Aufgabe 5
Was in Aufgabe 4 noch im Hauptprogramm getan wurde soll nun in die Funktion qrFactorStep (die wie der Name vage andeute einen Schritt der QR-Faktorisierung durchführt) gepackt werden:
template <typename MA>
double
qrFactorStep(GeMatrix<MA> &A)
{
typedef typename GeMatrix<MA>::VectorView VectorView;
int n = A.numCols();
VectorView v = A(_,1);
// Hier Programm-Code einfuegen.
// Satt 0 soll die Funktion alpha zurueckgeben:
return 0;
}
Damit sollte das Programm die gleiche Ausgabe liefern wie in Aufgabe 4:
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
DEVector alpha(3);
GEMatrix A(4,3);
A = 2, 1, 3,
4, 4, 9,
6, 9, 16,
8, 16, 24;
alpha(1) = qrFactorStep(A);
cout << "A = " << A << endl;
cout << "alpha = " << alpha << endl;
}
Aufgabe 6
Noch ein paar Programm Zeilen und die QR-Zerlegung ist fertig implementiert.
In blatt6.h muss dafür die Funktion qr fertig programmiert werden:
template <typename MA, typename VX>
void
qr(GeMatrix<MA> &A, DenseVector<VX> &alpha)
{
typedef typename GeMatrix<MA>::View MatrixView;
int m = A.numRows();
int n = A.numCols();
// Hier Programm-Code einfuegen
}
Dazu sollten circa 4 Zeilen Programm-Code ausreichend sein! Die Funktion soll A mit der QR-Zerlegung überschreiben und die Diagonale von R in alpha speichern.
Wie man das so leicht hinbekommt? Ganz einfach! Die Funktion qr ruft selber die Funktion qrFactorStep auf. Dabei übergibt sie verschieden Teilmatrizen von A:
- Im ersten Schritt die gesamte Matrix A
- Im zweiten Schritt A ohne die erste Zeile und ohne die erste Splate.
- Im dritten Schritt A ohne die ersten zwei Zeile und ohne die ersten zwei Spalten.
- ...
Einen Matrixausschnitt kann man so festlegen:
MatrixView B = A(_(2,m),_(2,n));
Dann bezieht sich B auf alle Elemente von A, bei denen
- der Zeilenindex zwischen 2 und m liegt,
- der Spaltenindex zwischen 2 und n liegt.
#include <blatt6.h>
typedef GeMatrix<FullStorage<double, ColMajor> > GEMatrix;
typedef DenseVector<Array<double> > DEVector;
typedef GEMatrix::VectorView VectorView;
int
main()
{
DEVector alpha(3);
GEMatrix A(4,3);
A = 2, 1, 3,
4, 4, 9,
6, 9, 16,
8, 16, 24;
qr(A, alpha);
cout << "A = " << A << endl;
cout << "alpha = " << alpha << endl;
}