Matrizenmultiplikation/ schneller Array-Zugriff

Hey,

ich bin grad dabei Matrizen auf verschiedene Arten zu multiplizieren und jeweils die Zeiten zu messen.
Für sehr kleine Matrizen ist der serielle Algorithmus natürlich noch schneller als der parallele aber ab ~100x100 Matrizen ändert sich dann.

Habe jetzt 4 Algorithmen:

Seriell
Parallel
Seriell transponierte Matrix
Parallel transponierte Matrix.

Um eine Matrix multiplizieren zu können muss ja die Anzahl der Spalten von Matrix A gleich der Anzahl der Zeilen von Matrix B sein.
Mit dem seriellen und dem parallelen Algorithmus kann ich auch problemlos eine (m x n)-Matrix mit einer (n x o)-Matrix multiplizieren.

Habe allerdings festgestellt, dass man mit transponierten Matrizen wesentlich schneller rechnen kann. Selbst der seriell-transponierte-Algorithmus ist für 1024x1024 Matrizen schneller als der “normal-parallele”. Der transponierte-parallele nochmal wesentlich schneller.
Im Moment kann ich leider nur quadratische Matrizen transponiert multiplizieren, weil sich ja beim transponieren die Anzahl der Zeilen und Spalten ändert.

Gibt es also eine Möglichkeit eine (m x n)-Matrix mit einer (n x o)-Matrix transponiert zu multiplizieren? Oder geht das per Definition nur mit quadratischen Matrizen?
Außerdem würde es mich mal interessieren wie dieser Performanceunterschied zu Stande kommt. Im Prinzip transponiere ich die eine Matrix nur, und führe dann einen fast-identischen Algorithmus aus, wie bei einer normalen Multiplikation? Erstaunlich für mich, dass der Algorithmus so viel schneller ist.

Hoffe mir kann da jemand helfen

Grüße

richtest du dich nach irgendeiner Quelle/ Literatur/ begründeter Erfahrung wenn du ‘transponierte-parallel’ überhaupt ausprobierst?

mein traurig-eingerostetes Mathematik-Wissen von Wiki aufgefrischt ist eine transponierte Matrix lediglich gespiegelt,
das bringt einerseits die Erkenntnis zu deiner Frage, dass
B^T * A^T = (A*B)^T ist
(schöner unter den Standardformeln in Matrix (Mathematik) – Wikipedia bei “Die transponierte Matrix”)
, also multipliziert werden kann, am Ende nochmal transponieren

siehe auch
MP-Forum: Matrizenmultiplikation / transponierte Matrix (Matroids Matheplanet)

aber bringt auch die logische Verwunderung, wie das auf die Geschwindigkeit/ auf den Arbeitsaufwand Auswirkungen haben soll?
es ist die gleiche Anzahl Arbeitsschritte, in einfachen 2x3-Matrixen gut nachzuzählen, wenn nicht eh offensichtlich

dass du bei Millionen-Rechnungen anderes beobachtest kann man ja erstmal hinnehmen, aber nochmal: wieso hast du es überhaupt ausprobiert,
aufgrund welcher Annahme, und das obwohl du anscheinend nichtmal B^T * A^T = (A*B)^T kennst?
bei Quadraten liegt der Versuch etwas näher, aber trotzdem… (edit: dort muss man aber doch genauso auf diese Formel achten?)


Performanceunterschied also, hmm, alles sieht danach aus dass es das nicht geben kann,
(wobei ich immer noch (wie du vielleicht auch) vermute/ befürchte, dass auf einmal irgendein mathematischer Geniestreich herbeigezaubert wird)
fange beim Nachprüfen wiederum mit Zählen an: wieviele einzelne Additionen/ Multiplikationen?

erhälst du
“normale Rechnung: 3548759 Rechenschritte, Zeit 1.5 sec”
“Transpo-supi Rech: 3548759 Rechenschritte, Zeit 1.1 sec”
?

achte darauf keine Fehler zu begehen wie jede Rechnung nur einmal durchzuführen, eine davon bei Programmstart schneller oder langsamer als die andere,

Ein paar mehr Infos wären nicht schlecht.

Aber allgemein: Erstmal “geht” nur AxB*BxC. Ob eine Matrix transponiert werden soll, entscheidet sich ja nicht willkürlich, sondern abhängig davon, was man modelliert. Und wenn eine Matrix nicht quadratisch ist, kann man sie nicht “transponiert multiplizieren, damit es schneller wird”.

Bei BLAS ist es im allgemeinen so, dass Matrizen nicht transponiert werden. Stattdessen gibt es bei den jeweiligen Funktionen einen Parameter, der angibt, ob eine der Matrizen als transponiert angenommen werden soll. Stark vereinfacht gesagt: Man übergibt zwar weiterhin die UN-transponierte Matrix, aber durch den Parameter werden intern andere Schleifen verwendet, so dass das Ergebnis dem entspricht, was rausgekommen wäre, wenn man eine normale Multiplikation mit der transponierten Matrix gemacht hätte.

Das hat einige Vorteile (ist aber ggf. aufwändiger zum implementieren). Ein Vorteil ergibt sich ggf. daraus, dass man nicht jede Menge Speicher umherschieben muss, “nur” um eine Matrix zu transponieren.

Das hängt auch mit der Frage zusammen, warum denn das eine oder andere schneller ist:

Ich (lehne mich in Anbetracht der Tatsache, dass du nichts über die Implementierung gesagt hast, mal weit aus dem Fenster, und) gehe davon aus, dass die Performanceunterschiede praktisch ausschließlich auf caching zurückzuführen sind. Du kannst ja zur Illustration mal diesen Test starten:

public class CachingTest
{
    public static void main(String args[])
    {
        for (int n=1000; n<=5000; n+=1000)
        {
            int input[] = new int[n*n];
            for (int i=0; i<input.length; i++)
            {
                input** = i;
            }

            long before = 0;
            long after = 0;
            int runs = 5;
            for (int i=0; i<runs; i++)
            {
                before = System.nanoTime();
                long result0 = runRowMajor(input, n);
                after = System.nanoTime();
                System.out.println("size "+n+" run "+i+" row    major "+
                result0+" time "+(after-before)/1e6+" ms");

                before = System.nanoTime();
                long result1 = runColumnMajor(input, n);
                after = System.nanoTime();
                System.out.println("size "+n+" run "+i+" column major "+
                result1+" time "+(after-before)/1e6+" ms");
            }
        }
    }


    private static long runRowMajor(int input[], int n)
    {
        long result = 0;
        for (int c=0; c<n; c++)
        {
            for (int r=0; r<n; r++)
            {
                result += input[r+c*n];
            }
        }
        return result;
    }

    private static long runColumnMajor(int input[], int n)
    {
        long result = 0;
        for (int r=0; r<n; r++)
        {
            for (int c=0; c<n; c++)
            {
                result += input[r+c*n];
            }
        }
        return result;
    }


}

Er durchläuft einen 1D Array “als 2D-Array” und summiert die Elemente auf. Ziemlich trivial. Einmal ist der Durchlauf row-major und einmal column-major. Das ist ein Pseudo-Microbenchmark, und man sollte die Ergebnisse nicht überbewerten, aber… bei mir sind sie rund um sowas wie


...
size 5000 run 4 row    major 312499987500000 time 14.58925 ms
size 5000 run 4 column major 312499987500000 time 237.098033 ms

Man sieht, dass der Unterschied locker Faktor 15 sein kann - NUR abhängig von der Durchlaufrichtung! Entscheidend ist also nicht so sehr, ob die Matrix transponiert ist oder nicht, sondern “wie ihr Inhalt für die Multiplikation verwendet wird”.

(Related: http://forum.byte-welt.net/sonstiges/spielwiese/10283-javaone-shirt-contest-2.html#post69723 - und, weil der Foreninterne Link bedauerlicherweise vermutlich nicht so dauerhaft ist, wie er sein sollte, noch der direkte: https://github.com/javagl/HazelcastMatMul )

Du musst berücksichtigen, dass eine Matrixmultiplikation nicht kommutativ ist.
Es gilt [tex]\left( A \cdot B \right)^T = B^T \cdot A^T[/tex]. Du musst also die Argumente vertauschen, wenn du transponiert multiplizierst.

Edit: ok, wenn ich mir @Marco13 s Beitrag ansehe, dann geht mein Beitrag wohl am Thema vorbei…

interessanter Test, mit ‚Cache‘ der meist gleich bleibenden Multiplikation von c*n in runRowMajor() gibts bei mir nochmal 10% auf dessen Zeit

bei Matrixenmultiplikation mit zwei Matrixen, entweder beide normal oder beide transponiert eingerechnet, sollten aber derartige Vorteile allgemein nicht möglich sein, oder?
heben sich ja wenn überhaupt gegenseitig auf

edit: man kann natürlich ohne konkret an Transponierung zu denken beide ‚Datenmassen‘ optimal vorbereitet ablegen,
ob sich dieser Aufwand gegenüber Berechnung lohnt?

sag das nicht, das wäre dann ja auch eine Aussage über mein Posting :wink:

Hey, danke erstmal für eure ausführlichen Antworten.
Hätte den Sourcecode direkt mitposten sollen.

Also ich hab die Algorithmen allesamt mit kleineren 2x2, 3x3 getestet, und händisch ausgerechnet ob das Ergebnis stimmt, deswegen denke ich auch, dass es für größere geht. Die Algorithmen erscheinen mir auch sehr logisch.
Nun zum Code:

[SPOILER]```public Matrix multiply(Matrix bval) {
double [][] res = new double [val.length][bval.val[0].length];
double z,x;
if (val[0].length == bval.val.length){
for (int i=0; i< val.length; i++){
for (int j=0; j< bval.val[0].length; j++){
x=0;
for (int k=0; k<bval.val.length;k++){
z = val**[k] * bval.val[k][j];
x += z;
}
res**[j]=x;
}
}

		return new Matrix(res);
	}
	throw new IllegalArgumentException();

}```[/SPOILER]

Das ist der Code für eine ganz normale serielle Multiplikation.
Komplexität sollte nach meinem Verständnis bei O(n³) liegen.

Der parallele Code ist dementsprechend sehr ähnlich:

[SPOILER]public ParallelMatrix multiply(Matrix B) { double[][] res = new double[val.length][B.val[0].length]; ParallelMatrix C = new ParallelMatrix(res); Thread[] threads = new Thread[val.length]; for (int i = 0; i < val.length; i++) { threads** = new MultThread(C,this,B,i); threads**.start(); } for (int i =0; i < val.length; i++ ) { try { threads**.join(); } catch (InterruptedException e) { e.printStackTrace(); } } return C; }[/SPOILER]

Die Thread run-Methode führt dann praktisch nur die inneren beiden for-Schleifen aus, wie sie auch im seriellen Code sind.

Soo. Und jetzt das für mich immer noch Erstaunliche:
Die Multiplikation mit einer transponierten Matrix.
Zuerst transponiere ich sie:
[SPOILER]public Matrix transpose() { double[][] res = new double[val[0].length][val.length]; for (int i = 0; i < val.length; i++) { for (int j = 0; j < val[0].length; j++) { res**[j] = val[j]**; } } return new Matrix(res); }[/SPOILER]

Dann multipliziere ich wieder. Praktisch genau nach dem gleichen Schema, wie schon im seriellen Code.
[SPOILER]```public Matrix multiply(Matrix B) {
double [][] res = new double [val.length][B.val[0].length];
double z,x;
Matrix Bt = B.transpose();
double[][] Bval = Bt.val;
if (val[0].length == Bval.length){
for (int i=0; i< val.length; i++){
for (int j=0; j< Bval[0].length; j++){
x=0;
for (int k=0; k<Bval.length;k++){
z = val**[k] * Bval[j][k];
x += z;
}
res**[j]=x;
}
}

		return new Matrix(res);
	}
	return new Matrix(res);
    }```[/SPOILER]

Komplexität sollte auch hier bei O(n³) liegen. Im Prinzip wurden nur die Indizes leicht angepasst, nachdem transponiert wurde.
Parallel sind das ganze dann wieder so aus wie schon zuvor. Nur eben leicht angepasst.

Lasse ich nun eine Zeitmessung laufen komme ich auf folgende Zeiten:

Seriell: 11064
Parallel: 3092
Transponiert Parallel: 237
Transponiert Seriell: 978

Hier mal der Source zur Zeitmessung. Da sollte alles mit rechten Dingen zugehen:

[SPOILER]```public void timing() {

	Matrix aSeriell = new Matrix(a);
	Matrix bSeriell = new Matrix(b);
	ParallelMatrix aParallel = new ParallelMatrix(a);
	ParallelMatrix bParallel = new ParallelMatrix(b);
	ParallelTransMatrix aTrans = new ParallelTransMatrix(a);
	ParallelTransMatrix bTrans = new ParallelTransMatrix(b);
	TransMatrix aT = new TransMatrix(a);
	TransMatrix bT = new TransMatrix(b);
	
	
	long ts = System.currentTimeMillis();
	aSeriell.multiply(bSeriell);
	ts = System.currentTimeMillis() - ts;
	System.out.println("Seriell: " + ts);
	
	long tp = System.currentTimeMillis();
	aParallel.multiply(bParallel);
	tp = System.currentTimeMillis() - tp;
	System.out.println("Parallel: " + tp);

	long ttp = System.currentTimeMillis();
	aTrans.multiply(bTrans);
	ttp = System.currentTimeMillis() - ttp;
	System.out.println("Transponiert Parallel: " + ttp);

	long tts = System.currentTimeMillis();
	aT.multiply(bT);
	tts = System.currentTimeMillis() - tts;
	System.out.println("Transponiert Seriell: " + tts);

}```[/SPOILER]

Die Messungen sind für zwei 1024x1024 Matrizen. Ich finde man sieht auch ganz schön, dass der Verhältnis gleich ist. Der normal-parallele ist ~3x so schnell wie der normal-serielle.
Der parallel-transponierte wiederum ~3x so schnell wie der transponiert-serielle.

Naja, man sieht ja deutlich, dass der transponiert-paralle “alles in den Schatten stellt”. Deswegen auch meine eigentliche Frage. Es geht ja im Prinzip nur mit quadratischen Matrizen, da sonst das Verhältnis der Zeilen bzw Spalten zueinander nicht mehr passt. Da ich mathematisches auch nicht wirklich bewandert bin, hab ich mir die Frage gestellt ob es nicht eine Möglichkeit gibt, diese gute Performance auch für nicht-quadratische Matrizen zu nutzen.

Grüße

#Den Test werde ich sofort mal laufen lassen.

nur eine Matrix transponiert,
gut, da kommt dann natürlich genau die Array-Durchlauf-Frage heraus, die erstaunlicherweise anscheinend so viel ausmacht, muss ich mir merken,

aber A * B^T hilft dir nicht zur Berechnung von A * B, auch im quadratischen Falle, das sieht man ja an kleinen Beispielen,
im anderen Fällen in der Tat gar nicht erst berechnenbar, was aber auch nicht schlechter ist

die Hoffnung mag ja bestehen aus dem völlig anderen Zahlen von A * B^T noch auf A * B zu kommen, ich fürchte aber vergebens

womöglich bzw. offensichtlich lohnt sich aber wie gesagt, ohne Transponierung die Werte in neuen Arrays besser abzulegen für linearen Durchlauf,
dann auch Geschwindigkeitsgewinn + richtiges Ergebnis + auch für nicht quadratische,
einfach nur dafür sorgen dass immer Arrays der Reihe nach durchlaufen

Das richtige Ergebnis kommt schon raus:
Hier kleines Beispiel:

MatrixScript script = new MatrixScript();
		double a[][] = script.setArray(4, 4);
		double b[][] = script.setArray(4, 4);
		TransMatrix am = new TransMatrix(a);
		TransMatrix bm = new TransMatrix(b);
		am.print();
		bm.print();

		Matrix cm = am.multiply(bm);
		cm.print();

Bringt mich zur folgender Ausgabe:

Und das sollte die richtige Matrix sein.

Das kann schon allein deswegen nicht stimmen, weil ja jede Matrix die transponierte der transponierten ist - d.h. wenn du eine „Statistik“ über alle Matrizen machen würdest, dann würde das Transponieren überhaupt keinen Unterschied machen…

  1. Wie sind deine Matrizen überhaupt gefüllt - mit Zufallszahlen?

  2. Hast du die Ergebnisse verglichen, d.h. rechnen alle deine Verfahren korrekt?

  3. Warum die verschiedenen Konstruktoren? eine Matrix sollte nichts davon wissen, wie sie multipliziert wird?

  4. du musst rücktransponieren (wenn du AB ausrechnen willst, dann musst du ja (B^t)(A^t)^t berechnen?)

*** Edit ***

Glaubst du das rechnet einer von uns nach?

Multipliziere doch mal zwei Matrizen A und B mit allen vier Verfahren - erhältst du immer das gleiche Ergebnis?

ok, bei der Transpo-Multiplikation wird ja gar nicht nach Vorschrift multipliziert (etwa Zeilen mit Spalten je nach Interpretation), sondern Zeilen mit Zeilen,
also alles richtig gemacht:

die Werte so umsortiert dass für Berechnung ideal die Arrays der Reihenfolge nach durchlaufen werden,
von mir ein Sternchen

und immer noch die Frage, wie du auf die Idee kamst, falls nicht bereits Wissen darum :wink:


edit: noch die Fragen nach nicht-quadratischen, hmm, mal sehen ob noch vor Mittagspause was zu sagen

Ja genau, die werden mit Zufallszahlen gefüllt. Im Beispiel von oben waren die gleich, weil die direkt hintereinander gefüllt wurden und der Systemtime als Seed wohl gleich war.

Hab mal ein Sleep eingefügt und hatte zwei unterschiedliche Arrays, wo aber immer noch das richtige Ergebnis rauskam.

Geprüft hab ich das hier:
[video]http://www.jetzt-rechnen.de/Mathematik/Matrixmultiplikation.html#[/video]

Der Konstruktor ist grundätzlich gleich. Die Klassen leiten alle von Matrix ab und rufen dessen Konstruktor mir super() auf.
Mehrere Klassen habe ich, da ich verschiedene Threads benutze. Hab die Threads, die ich erzeuge für jede Matrix als innere Klasse implementiert.

*** Edit ***

So hab mich mal registriert, jetzt geht das editieren auch :stuck_out_tongue:
Ich mache das zusammen mit einem Partner für die Uni. Hatten eine Art vereinfachten Pseudo-Code, den wir versuchen umzusetzen.

Was mich dann noch zu meiner ursprünglichen Frage leitet, ob es überhaupt möglich ist, den transponierten Algorithmus für nicht-quadratische Matrizen anzuwenden.
Der Speedup gegenüber den anderen Algorithmen scheint nämlich enorm zu sein. Dann noch einen Threadpool nutzen und man kriegt vllt noch ein bisschen mehr raus :o

Mich würde es schon wundern, ob a**[zugriff] vs. a[zugriff]** wirklich so viel bringen können im Vergleich zu den Additionen/Multiplikationen?

Wenn die Matrizen nicht quadratisch sind, kannst du die ja mit 0en auffüllen zu quadratischen Matrizen.

*** Edit ***

Versuch?

leg doch die Matrizen einfach in ein laaaaaaaaanges 1-dim Array double[n^2]

und berechne die nötigen Zugriffsindizes beim Multiplizieren von hand - sollte noch schneller sein

[quote=Tarrew]Die Klassen leiten alle von Matrix ab und rufen dessen Konstruktor mir super() auf.
Mehrere Klassen habe ich, da ich verschiedene Threads benutze. Hab die Threads, die ich erzeuge für jede Matrix als innere Klasse implementiert.[/quote]

Das scheint mir ein recht seltsames Vorgehen zu sein. Darüber würde ich nochmal nachdenken, ob die Threads wirklich innerhalb von Matrix-Klassen verwaltet werden sollen. Ich würde das immer strikt voneinander trennen und den Threadkram von außen “drüberstülpen”.

also

            double [][] res = new double [val.length][B.val[0].length];
            double z,x;
            Matrix Bt = B.transpose();
            double[][] Bval = Bt.val;
            if (val[0].length == Bval.length){
                for (int i=0; i< val.length; i++){
                    for (int j=0; j< Bval[0].length; j++){
                        x=0;
                        for (int k=0; k<Bval.length;k++){
                        z = val**[k] * Bval[j][k];
                        x += z;
                        }
                        res**[j]=x;
                    }
                }
               
                return new Matrix(res);
            }
            return new Matrix(res);
            }

hier schlicht die if-Bedingung und eines der sowieso überflüssig doppelten returns weglassen?
wichtig ist nur, dass A.val[0].length und B.val.length bzw. dann ‘transponiert’ Bval[0].length gleich sind, denn die beiden werden gleichmäßig mit der k-Schleife durchlaufen,

die anderen für die i und j-Schleife dürfen unterschiedlich sein, gibt dann nicht-quadratisches res-double[][]

sinnvoll ist vielleicht noch ein Test, ob nicht B.val.length > A.val[0].length ist, was in der k-Schleife nicht auffallen würde,
dass die beiden gleich sind ist die einzige Bedingung (neben an sich Matrix-Eigenschaften)

oder bei welchem Gegenbeispiel läuft es nicht?
ich poste gleich noch meine angefangene double[][]-pur-Variante

edit: die transpo-Methode war bisher noch falsch, statt j,i auf i,j zuzuweisen muss es andersrum geschehen, fällt bei quadratisch natürlich nicht auf

hier rein mit double[]:

    public static void main(String[] args)  {
        double[][] a =       {
                {1, 2, 3, 4},
                {5, 6, 7, 8}};
        double[][] b =   {
                {1, 2},
                {1, 2},
                {1, 2},
                {1, 2}};
        double[][] res = multiply(a, b);
        for (int i = 0; i < res.length; i++)  {
            System.out.println(Arrays.toString(res**));
        }
    }

    public static double[][] multiply(double[][] a, double[][] b) {
        if (a[0].length != b.length) throw new RuntimeException("sorry");
        double[][] bt = transpose(b);

        double[][] res = new double[a.length][bt.length];
        double z, x;
        for (int i = 0; i < a.length; i++)   {
            for (int j = 0; j < bt.length; j++)      {
                x = 0;
                for (int k = 0; k < a[0].length; k++)       {
                    z = a**[k] * bt[j][k];
                    x += z;
                }
                res**[j] = x;
            }
        }
        return res;
    }

    public static double[][] transpose(double[][] m)  {
        double[][] res = new double[m[0].length][m.length];
        for (int i = 0; i < m.length; i++)    {
            for (int j = 0; j < m[0].length; j++)  {
                res[j]** = m**[j];
            }
        }
        return res;
    }
}```
Ergebnis:

[10.0, 20.0]
[26.0, 52.0]

Teilen der Diskussion bin ich jetzt auf die Schnelle nicht gefolgt (und der Code… den hätte ich wohl mal ausprobiert, wenn es EIN Codeblock gewesen wäre, mit allen Klassen und einer main(), den man direkt in eine IDE Copy+Pasten könnte - einschließlich der ja nicht ganz unwichtigen Klasse “Matrix” (auch wenn die nur ein “Behälter” für einen double[][]-Array zu sein scheint))

Deswegen noch ein bißchen high-level-Gelaber:
@Bleiglanz : “Mich würde es schon wundern, ob a**[zugriff] vs. a[zugriff]** wirklich so viel bringen können im Vergleich zu den Additionen/Multiplikationen?”

In der Tat (siehe den “CachingTest” von oben). Einen Speedup von 15 durch die “richtige” Zugriffsreihenfolge kann einfach sein. Einen Speepdup von 15 durch Parallelisierung (selbst auf einem 16-Kerner) muss man erstmal hinkriegen. Wir sind schon lange gegen die Speicher-Durchsatz-Wand gefahren.

Diesen Effekt des Cachings hatte ich schonmal zur Diskussion gestellt, aber den Link grabe ich jetzt aus politischen Gründen mal nicht aus. Es ist ja so: Man hat zwei Matrizen, die man multiplizieren will. Bei der einen passiert der Zugriff “row-major”, und bei der anderen “column-major”. D.h. egal, ob die Matrizen row-major oder column-major gespeichert sind, eins von beidem ist immer kagge und bremst das ganze aus. Ich hatte deswegen auch schon “absurde” Optionen in Betracht gezogen - nämlich ob es effizienter sein könnte, die “falsche” Matrix in eine neue Matrix zu kopieren, die für die Multiplikation das richtige Format hat. (Aber getestet habe ich das AFAIK dann nicht…). (Vielleicht ist das ähnlich zu dem, was Tarrew vielleicht indirekt (!) mit dem “tranponiert multiplizieren” meinte…!?)


Zum Thema “Parallele Matrizenmultiplikation” an sich: Es gibt ja verschiedene Optionen, worüber man parallelisieren kann. Man könnte einfach bei den 3 verschachtelten Schleifen die äußerste parallel abarbeiten lassen (so ist es im Moment gemacht, wenn ich das richtig gesehen habe). Für “große” Matrizen und verteiltes Rechnen gibt es aber ein Verfahren, das wohl “DAS Standardverfahren” ist. Das heißt “SUMMA”, hatte ich für die oben verlinkte Parallele Matrix Multiplikation verwendet, und ist in diesem Paper beschrieben: http://www.netlib.org/lapack/lawnspdf/lawn96.pdf (Darin wird auch die Frage nach den transponierten Varianten behandelt). Die Kernidee ist, dass man alle Einträge der Ergebnismatrix parallel berechnet. Das ist dann jeweils ein dot-Product aus Zeilen/Spalten von A/B, die aufsummiert werden. Welche Auswirkungen DAS bei shared memory auf den Cache hat (im Vergleich zum anderen Parallelisierungsansatz) müßte man sich aber genauer ansehen.

Nö, wenn ich das richtig verstandenhabe dann macht er das so, dass er eine der beiden transponiert und dann „falsch multipliziert“ - d.h. er geht von

 z = val**[k] * bval.val[k][j];//k läuft
                    x += z;

zu

                    z = val**[k] * Bval[j][k]; //k läuft
                    x += z;

Aufgrund des 2d-Zugriffs ist das wohl schneller, aber dass es so schnell ist, das ist schon merkwürdig.

man mag es kaum glauben, die Schleifen müssen sein, die Multiplikationen und Additionen alle exakt gleich, die Arrayzugriffe sind da,
sind nicht wegzuoptimieren (ok, in der Ausführung vielleicht doch, aber im Code real vorhanden),
nur dass zwischen Arrays gewechselt wird soll 60% sämtlichen Zeitverbrauches der ganzen Aktion ausmachen?

aber doch sieht es so aus, hier ein double-Testprogramm mit Berechnung der Matrix-Multiplikation,
wenn der einfachere Benchmark (dort Faktor 15, hier nur 3) nicht ausreicht:

{
    static double[][] a = new double[800][805];
    static double[][] b = new double[a[0].length][900];

    public static void main(String[] args) {
        fill(a);
        fill(b);

        check(true);
        check(false);
        check(true);
        check(false);
        check(true);
        check(false);
    }

    static void check(boolean fast)  {
        long time = System.currentTimeMillis();
        double[][] res = fast ? multiply(a, b) : multiplySlow(a, b);
        time = System.currentTimeMillis() - time;
        long sum = 0;
        for (int i = 0; i < res.length; i++)
            for (int j = 0; j < res**.length; j++)
                sum += (i + 1) * (j + 10) * res**[j];
        System.out.println(sum + " - " + time);
    }

    static Random r = new Random(42);

    static void fill(double[][] m) {
        for (int i = 0; i < m.length; i++)
            for (int j = 0; j < m[0].length; j++)
                m**[j] = r.nextInt(15);
    }

    public static double[][] multiply(double[][] a, double[][] b) {
        if (a[0].length != b.length) throw new RuntimeException("sorry");
        double[][] bt = transpose(b);

        double[][] res = new double[a.length][bt.length];
        double z, x;
        for (int i = 0; i < a.length; i++)  {
            for (int j = 0; j < bt.length; j++) {
                x = 0;
                for (int k = 0; k < a[0].length; k++) {
                    z = a**[k] * bt[j][k];
                    x += z;
                }
                res**[j] = x;
            }
        }
        return res;
    }

    public static double[][] multiplySlow(double[][] a, double[][] b)   {
        if (a[0].length != b.length) throw new RuntimeException("sorry");

        double[][] res = new double[a.length][b[0].length];
        double z, x;
        for (int i = 0; i < a.length; i++) {
            for (int j = 0; j < b[0].length; j++)  {
                x = 0;
                for (int k = 0; k < a[0].length; k++)    {
                    z = a**[k] * b[k][j];
                    x += z;
                }
                res**[j] = x;
            }
        }
        return res;
    }

    public static double[][] transpose(double[][] m) {
        double[][] res = new double[m[0].length][m.length];
        for (int i = 0; i < m.length; i++){
            for (int j = 0; j < m[0].length; j++) {
                res[j]** = m**[j];
            }
        }
        return res;
    }
}

Ausgabe bei mir, schnelle und langsames Matrixprodukt 3x im Wechsel, jeweils Kontrollsumme + Zeit:


5239993975677955 - 2292
5239993975677955 - 5627
5239993975677955 - 1749
5239993975677955 - 5473
5239993975677955 - 1744
5239993975677955 - 5471

[QUOTE=SlaterB]man mag es kaum glauben, die Schleifen müssen sein, die Multiplikationen und Additionen alle exakt gleich, die Arrayzugriffe sind da,
sind nicht wegzuoptimieren (ok, in der Ausführung vielleicht doch, aber im Code real vorhanden),
[/quote]

Im einfachsten Fall stimmt das, man hat O(n^3). Mit dem Coppersmith–Winograd algorithm - Wikipedia, the free encyclopedia hätte man O(n^2.37…). (Theoretisch…)

aber doch sieht es so aus, hier ein double-Testprogramm mit Berechnung der Matrix-Multiplikation,
wenn der einfachere Benchmark (dort Faktor 15, hier nur 3) nicht ausreicht:

Der erste Test macht ja wirklich „nichts“ außer Daten lesen, vermutlich nicht zuletzt deswegen ist der Unterschied dort deutlich größer. Aber selbst Faktor 3 ist ja nicht unerheblich (und deckt sich mit den Messungen von Tarrew, auch wenn ich den Code in beiden Fällen nicht 100% nachvollzogen habe)

na, dort eher Faktor 10+, ich hoffe mal wegen seriell unnötig langsamer, warum weiß ich aber nicht und will Code dann auch nicht ausführen

"Seriell: 11064
vs Transponiert Seriell: 978

Parallel: 3092
vs Transponiert Parallel: 237"

wobei die Arbeit letztlich immer nur in 2-3 Zeilen steckt, hmm…,
vielleicht rechnerspezifisch :wink:

*** Edit ***

So ich melde mich mal wieder.
Hab den Code inzwischen angepasst mit dem Beispiel von SlaterB. Funktioniert tadellos. :slight_smile: Danke dafür.

Verfolge das gerade ziemlich interessiert und diesen wahnsinnigen Speedup zwischen transponiert und normal find ich ziemlich spannend.

Wenn ich normal-sequentiell und normal-parallel vergleiche, komme ich auf einen maximalen Speedup von ~4. Denke mal das hängt dann von der Anzahl der Cores ab? Auf meinem Xeon hab ich halt einen maximalen Speedup von ungefähr 4-5. Auf meinem Laptop dementsprechend von ~2.
Beim Vergleich von transponiert-sequentiell und transponiert-parallel siehts ähnlich auch. Erscheint mir eigentlich auch logisch.

Für zwei 2048x2048 Matrizen komme ich zB auf:

Groesse der Matrix 2048

Seriell: 145974
Parallel: 35808
Speedup 4.0765750670241285
Transponiert Parallel: 1756
Transponiert Seriell: 7734
Speedup 4.404328018223235

Der Unterschied zwischen transponiert seriell und normal-seriell liegt hier ja schon bei Faktor 20. Ziemlich beeindruckend ;D

Btw: Wusste nicht in welchem Forum ich eine Frage über Java am besten stelle, dachte am naheliegensten wäre vllt java-forum.org. Hab dann zufällig diesen Artikel gelesen und bin hier gelandet. Muss sagen ihr seid echt super :stuck_out_tongue: Großes Lob. Bin so viel Aktivität und Hilfsbreitschaft aus anderen Foren eher nicht gewohnt.