Matrixmultiplikation mit "Einstiegsspunkt" (Householder-Transformation)

[Monophagus]

Hallo Forum,

in einem meiner Studienfächer habe ich für diese Woche eine Programmieraufgabe aufbekommen, nämlich eine Speicher-effiziente QR-Zerlegung mittels Householder-Transformation. Mein Code für eine formale QR-Zerlegung funktioniert bereits perfekt, nun würde ich das ganze gerne umschreiben, um Speicherplatz zu sparen. Im Moment scheitert es nur noch an einer Sache, nämlich einer speziellen Matrixmultiplikation wie folgt:

Schema:

1 0 0 0 0 0           * 0 0 0 
0 1 0 0 0 0           * * 0 0
0 0 * * * *     *     * * * *
0 0 * * * *           * * * *
0 0 * * * *           * * * *
0 0 * * * *           * * * *
         H      *      A

Wobei hier H für für die Householder-Matrix im dritten Schritt (k=2) steht und A meine zu zerlegende Matrix. Nun möchte ich eine Multiplikation durchführen, welche nach gewohnter Matrixmultiplikationsregel nur folgende Stellen berechnet und direkt in A speichert (gekennzeichnet mit einem a):

* 0 0 0 
     * * 0 0 
=    * * a a
     * * a a         
     * * a a     
     * * a a   
        A

Meine Idee dazu war, k als einen "Einstiegspunkt" für die Multiplikation zu wählen. Mein Codeschnipsel für diese Multiplikation sieht folgendermaßen aus, liefert aber falsche Ergebnisse:

double** HH_matrix_mult(double** H, double** A, int x, int m, int n){
	//x als Einstiegspunkt der Multiplikation
		int i, j, k;
		double z;
		for(j=x;j<n;j++){
			for(i=x;i<m;i++){
				z=0;
				for(k=0;k<m;k++) z=z+H[i][k]*A[k][j];
				A[i][j]=z;
			}
		}
	return A;
}

Liegt in dem Code irgendwo ein Fehler? Habe ich allgemein einen Denkfehler gemacht? Ich bin das ganze Schema schon dreimal durchgegangen und finde keinen Fehler - sieht jemand von euch problematische Stellen?

Das ganze ist in folgendem Code eingebettet:

(Q wird noch zur Kontrolle berechnet, fliegt zum Schluss raus. Der Hintergrund der obigen Multiplikation ist, dass ich die vorderen Spalten von A unbeeinflusst lassen will, um dort die Vektoren speichern zu können, die zur Erzeugung der Householder-Matrizen dienen. Die dürfen natürlich nicht von einer Multiplikation durcheinander gebracht werden.)

#include <stdio.h>
#include <math.h>
#include "matrix.h"
#include "QR_dec.h"

double norm_Euklid(double** V, int k, int m){
	double z=0; int i;
	for (i=k; i<m; i++) z = z + V[i][k]*V[i][k];
	return sqrt(z);
}

double skalarprodukt(double* v, int m){
	double z=0; int i;
	for(i=0; i<m; i++)  z=z+v[i]*v[i];
	return z;
}

double dyade(double* v, int i, int j){
	return v[i]*v[j];
}

int signum(double x){
	if(x>=0)
		return  1;
	else
		return -1;
}

double** HH_matrix_mult(double** H, double** A, int x, int m, int n){
	//x als Einstiegspunkt der Multiplikation
		int i, j, k;
		double z;
		for(j=x;j<n;j++){
			for(i=x;i<m;i++){
				z=0;
				for(k=0;k<m;k++) z=z+H[i][k]*A[k][j];
				A[i][j]=z;
			}
		}
	return A;
}

double** HH_matrix(double* v, int m){
	double** H=matrix_init(m,m);
	int i, j;
	for(j=0;j<m;j++){
		for(i=0;i<m;i++){
				H[i][j]=(-1)*(2/skalarprodukt(v,m))*(dyade(v,i,j));
				if(i==j) H[i][j]=H[i][j]+1;
		}
	}
	return H;
}

double** QR_dec(double** A, int m, int n){
// QR-Zerlegung nach HH
// A mxn, m größergleich n, sinnvoll bei rg(A)=n
		int i, k;
		double a;
		double*  s=vector_init(m);
		double** H;
		double** Q=matrix_init(m,m);
//Initialisieren von Q_0:
		for(k=0;k<m;k++){
			for(i=0;i<m;i++){
				if(i==k)
					Q[i][k]=1;
				else
					Q[i][k]=0;
			}
		}
		for(k=0; k<n; k++){
				printf("\nQR_Dec Iterationsschritt: %i\n\n",k);
			for(i=0;i<k;i++)  A[i][k]=0;
			for(i=0;i<m;i++)  s[i]=A[i][k];
			a=norm_Euklid(A,k,m)*signum(A[k][k]);
			s[k]=s[k]+a;
			for(i=k+1;i<m;i++) s[i]=s[i]/s[k]; s[k]=1;
			H=HH_matrix(s,m);
			s[k]=a;
				printf("Q_%i:\n",k);
			Q=matrix_mult(Q,H,m,m,m);
				printf("Q_%i*A_%i:\n",k,k);
			A=matrix_mult(H,A,m,n,m);
			//A=HH_matrix_mult(H,A,k,m,n);
			for(i=0;i<m;i++) A[i][k]=s[i];
		}
	return A;
}

Vielen Dank für einen Hinweis oder etwas Multiplikations-Theorie!

[Wutz]

Kann es sein, dass du nach zuweisen des Ergebnisses in einen A-Wert diesen nochmal bei folgenden Berechnungen verwendest (und dann natürlich nicht mehr den Originalwert liest)?

[Monophagus]

Kann es sein, dass du nach zuweisen des Ergebnisses in einen A-Wert diesen nochmal bei folgenden Berechnungen verwendest (und dann natürlich nicht mehr den Originalwert liest)?

Whoa, Tatsache! Ich dachte, mit der Deklaration von z hätte ich bereits entgegengewirkt; aber stimmt, für die nachfolgende Werte wird ja dennoch ein verändertes Element in A miteingerechnet. Vielen Dank für den Hinweis, scheint wirklich der Fehler zu sein. Scheinbar komme ich um eine Teilallokierung nicht vorbei.

[Monophagus]

Sorry, Doppelpost

[*john 0]

Monophagus schrieb:

Meine Idee dazu war, k als einen "Einstiegspunkt" für die Multiplikation zu wählen. Mein Codeschnipsel für diese Multiplikation sieht folgendermaßen aus, liefert aber falsche Ergebnisse:

double** HH_matrix_mult(double** H, double** A, int x, int m, int n){
	//x als Einstiegspunkt der Multiplikation
		int i, j, k;
		double z;
		for(j=x;j<n;j++){
			for(i=x;i<m;i++){
				z=0;
				for(k=0;k<m;k++) z=z+H[i][k]*A[k][j];
				A[i][j]=z;
			}
		}
	return A;
}

Liegt in dem Code irgendwo ein Fehler? Habe ich allgemein einen Denkfehler gemacht? Ich bin das ganze Schema schon dreimal durchgegangen und finde keinen Fehler - sieht jemand von euch problematische Stellen?

Warum werden "double**" verwendet? Das ergibt keinen Sinn, auch wenn es Dir in den Funktionen das Leben vereinfacht. Als Hinweis schau Dir mal die BLAS bzw cBLAS an. Einfache "double*" sind besser (das ist flexibler und kompatibel mit fix deklarierten Arrays) und sparen schon einiges an Speicherplatz ein. Dann kann man R natürlich im packed format speichern, da es sich um eine obere Dreiecksmatrix handelt, wenn das denn erlaubt ist.

Dann stellt sich die Frage, warum die Funktion eine "double**" zurückliefert. A ist ein Parameter und wird überschrieben, dann sollte man das auch so dokumentieren.