Faltung (Convolution) eines riesigen Arrays optimieren

[martinibook]

Hallo,

ich habe ein großes 3D Array, das in jedem Punkt mit einem kleineren Array multipliziert werden soll. Im kleineren Arrays sind überall Nullen und eine Kugel aus Einsen in der Mitte. Letztlich entspricht das einem Volumenintegral in einem Segment des großen Arrays.

Momentan habe ich das in C++ ganz naiv gelöst, mit 6 for schleifen:
Pseudocode:

for x:
  for y:
    for z:
      summe = 0
        for i (von -radius bis +radius):
          for j (von -radius bis +radius):
            for k (von -radius bis +radius):
              if abstand i, j, k zum Mittelpunkt < radius:
                summe += 1
      resultat[x][y][z] = summe;

Das Problem ist nur leider, dass das dann 2150*2150*460 * 137*137*137 = 5.467.596.451.550.000 Iterationen sind. Und ein Datensatz besteht aus 5 solcher Segmente. Mein Programm kam nach einigen Minuten zu dem Ergebnis, dass es ungefähr 6 Jahre dauert, bis es damit fertig ist. Das kann es ja einfach nicht sein.

Bisher wurde das Problem dadurch bewältigt, dass man die Daten einfach um Faktor 5 verkleinert, das ist dann letztlich 5**6 weniger Iterationen, =15625, somit wäre man schon in so zwei Stunden fertig. Das Programm ist momentan in IDL implementiert, einer Arraybasierten Sprache wie Matlab. Und dort hat der ganze Algorithmus, der aus drei Teilen dieser Art besteht, nur eine Stunde gebraucht.

Hat jemand vielleicht eine Idee, wie man das etwas optimieren könnte?

Das Ziel der ganzen Sache ist es im großen Array Kugeln zu finden, als Ergebnis der Faltung erhält man Maxima an den Punkten, wo die Kugelmaske am besten gepasst hat.


Grüße,

martinibook

 

[Marduk]

Kurze Zwischenfrage damit ich das Problem ganz verstehe: Sollte statt "summe += 1" nicht "summe += 3dArray[x][y][z]" heissen?
Zudem: Enthält das 3dArray beliebige Zahlen oder nur '0' und '1'?

 

[anegmawad]

Du solltest die i,j,k-Schleife nach außen setzen, da der Vergleich in der Schleife nur davon abhängt. Die dann inneren x,y,z-Schleifen machst du dann ins if.

 

[Dalgliesh]

So, wie Du es beschreibst, ist das keine Faltung.

Irgendwas stimmt nicht, denn für jedes (x,y,z) bekommst Du den gleichen Wert.

Wenn es wirklich eine Faltung ist, kannst Du durch eine Fast-Fourier-Transformation die Laufzeit von quadratisch nach linear reduzieren.

 

[martinibook]

Marduk: Ja, absolut. Ich habe den Fehler auch irgendwann im Code entdeckt.

Den Code habe ich ein kleines bisschen optimiert, in dem in der inneren (i, j, k) Schleife nicht der komplette Würfel durchgelaufen wird, sondern nur die wirkliche Kugel. Sofern ich mich in der Mathematik nicht verrannt habe, spart das ungefähr 45% der Iterationen ein und man kommt auch um das Berechnen des Abstandes in jeder Iteration herum.

for (int z = 0; z < sz; z++) {
	for (int y = 0; y < sy; y++) {
		for (int x = 0; x < sx; x++) {
			pos = x + sx * y + sx * sy * z;
			unsigned int sum = 0;

			int limit_k = MIN(sz, z + radius);
			for (int k = MAX(0, z - radius); k < limit_k; k++) {

				int limit_j = MIN(sy, y + sqrt(pow(radius, 2) - pow(k, 2)));
				for (int j = MAX(0, y - radius); j < limit_j; j++) {

					int limit_i = MIN(sx, x + sqrt(pow(radius, 2) - pow(k, 2) - pow(j, 2)));
					for (int i = MAX(0, x - radius); i < limit_i; i++) {
						unsigned int pos2 = i + sx * j + sx * sy * k;
						sum += threshold[pos2];
					}
				}
			}
			// apply a threshold, this does not save the data if it is too low
			if (sum >= thresholdValue) {
				conv[pos] = sum;
			}
		}
	}
}

Die Daten im Array threshold[] sind entweder 0 oder 1, da auf sie ein Grenzwert angewandt wurde.

@Dalgliesh: Wie ist quadratisch zu verstehen? Meinst du, dass ich letztlich O(N⁶) auf O(N³) kommen könnte, wenn man N als die Kantenlänge des großen und kleinen Würfels nimmt?

FFT müsste ich mir mal anschauen. Ich könnte mir vorstellen, dass der Algorithmus in IDL entsprechend so funktioniert und eben schneller ist.

 

[dmssmd]

Falls du nicht über FFT gehen willst könntest du noch “SIMD Within A Register” (SWAR) ausprobieren. Durch das geschickte ausnutzten der MMX registers könntest du ein Speed up von ca. 3-4 bekommen.
Ist zwar nicht so schön zu programmieren aber rentiert sich.

Hier mal ein Header file in der die Befehle sind: Anhang anzeigen mmx.h.zip

Du könntest auch deine interne for loop aufdröseln. Sodass du jeweils 4 schritte ausführst, dann nutzen die meisten Compiler die MMX register, wenn du die Flag -O setzt.

Deine inner loop würde dann so aussehen. Man muss dann nur aufpassen, dass die Grenzen passen.

unsigned int pos2;
for (int i = MAX(0, x - radius); i < limit_i; i+=4) {
	tpos2 = i + sx * j + sx * sy * k;
	sum += threshold[pos2];
        tpos2 = i +1+ sx * j + sx * sy * k;
	sum += threshold[pos2];
        tpos2 = i+2 + sx * j + sx * sy * k;
	sum += threshold[pos2];
        tpos2 = i+3 + sx * j + sx * sy * k;
	sum += threshold[pos2];
}

Du solltest auch vielleicht darauf achten alle Operationen wie sqrt und pow möglichst wenig oft in tiefen loops zu machen, sqrt ist sehr teuer. Und pow(x,2) kann man einfach durch x*x ersetzten. Auch solltest du variablen wie z.B. pos2 nach Möglichkeit nicht in der Loop deklarieren, in der du sie benötigst sondern einen Ebene höher.

Die Compiler sind meistens leider sehr dumm.

 

[martinibook]

Für die Limits brauche ich das sqrt() wahrscheinlich schon. Dass ich pow durch x*x ersetzen sollte erscheint logisch, da pow(double, double) ja viel allgemeiner ist und ja das ja mit e^ln() macht, und ln() ist sicher noch teurer als sqrt().

Ich hatte angenommen, dass die Variablen nicht immer neu initialisiert werden, wenn sie im Loop sind, da der Compiler das vielleicht irgendwie merkt. Gut, dann schreibe ich die alle nach ganz außen. Wahrscheinlich wäre ich dann ein klein wenig schneller, aber ein Jahr braucht das Programm vielleicht immer noch.

Danke für die Tipps bis hier her!

 

[Marduk]

Du solltest die innersten drei Loops ganz weg machen. Dann wird es einiges schneller ;-) Das ist möglich, wenn man das Resultat der vorangegangen Berechnung wieder benutzt. Du subtrahierst einfach die Werte vom alten Resultat die du nicht mehr willst und addierst die neuen Werte.
Ausserdem solltest du jegliche 'if' im innern vermeiden. Das ist so ziemlich das teuerste das du machen kannst. Im Moment hast du ja ausser bei threshold keine aber wenn du die Loops weg machst hast du da Probleme weil du die Randbedingungen speziell behandeln musst. Um diese 'if' statments loszuwerden alloziert du das Input-Array einfach um den Faktor 2*Radius grösser und füllst es am Anfang und am Ende (jeweils 1*Radius) einfach mit nullen. Die die Loops starten dann beim Index 'Radius'. So hast du keine speziellen Randbedinungen, sondern du addierst wenn du halt am Rand bist einfach nullen.

Wichtig ist auch, dass du deine Datentypen möglichst klein wählst. Wenn du z.B. weisst, dass das Array nur Werte zwischen 0 und 256 haben kann, dann sollte das Array vom Typ unsigned char sein.

Also für einen einfachen 2D Fall stelle ich mir das ungefähr so vor(Randbedingungen ignoriert (ich mache z.B. x-2 und wenn x dann 0 ist...)):
Radius=3

unsigned int oldSum=[initialisieren mit einer vorangegangen 'vollen' Berechnung]
[array wie beschrieben mit zusätzlichen Speicher initialisieren]
unsigned int newSum;
for (int x = 0; x < sz; x++) {
  for (int y = 0; y < sy; y++) {
    //bei einem radius von 3 muss man genau 5 Werte subtrahieren/addieren
    //fünf Werte hat man schon von der alten Summe
    newSum=oldSum-array[x-1][y+2]-array[x-2][y+1]-array[x-3][y]-array[x-2][y-1]-array[x-1][y-2]+array[x][y+2]+array[x+1][y+1]+array[x][y+2]+array[x+1][y-2]+array[x][y-2];
    [irgendwas damit machen/in ein result array schreiben ect.]
    oldSum=newSum;
  }
}

Am besten nimmst du dir ein Stück Papier und zeichnest es auf. Dann merkst du was ich meine.
Wenn ich mich nicht irre musst du für dein inneres Array von der Grösse 137 und 3 Dimensionen, 410 (3*137) Werte addieren/subtrahieren. Bin da aber nicht sicher. Keine Wurzeln/Potenzen rechnen nur ziemlich viel/hässlicher Code.

Das macht dan 2.126.350.000 Iterationen in denen du 410 einfache Operationen machst => 871.803.500.000. Bei einem 1Ghz Rechner und optimalen Cacheverhältnissen würde das noch ca. 15 Minuten brauchen (was natürlich niemals so sein wird...)

 

[martinibook]

Klar, letztlich muss man ja nur die Kugelhülle abgehen und entsprechend Daten addieren und subtrahieren, dadurch spart man sich im inneren sehr viel Arbeit. Und die Oberfläche skaliert quadratisch und die kubisch wie das Volumen. Das wäre schon mal ein Gewinn. Nur wird das dann interessant, wenn man in eine neue Zeile mit der Kugel geht, aber selbst wenn man es dann noch mal neu berechnet, lohnt sich das ja immer noch.

 

[Dalgliesh]

Ja, von N^6 nach N^3 durch FFT.

Beim Kugelhüllen-Ansatz mußt Du übrigens regelmäßig die ganze Kugel neu berechnen, sonst bekommst Du haufenweise Rundungsfehler, da für den a + b + c - a nicht genau b + c liefert, sondern b + c + Rundungsfehler, und nach 100 Schritten hast Du den Rundungsfehler * 100.

 

[apple][c]

Gebt es endlich zu: Ihr habt von dem, was ihr da schreibt genausoviel Ahnung wie ich.

Nämlich keine.

Wenn ihr das hier wirklich drauf habt, dann Chapeau!

 

[martinibook]

Ich habe mich heute da in die DFT und FFT eingelesen ... schon harter Stoff. Allerdings glaube ich, dass ich es mit der fftw Bibliothek vielleicht irgendwie hinbekommen kann. Die Daten habe ich schon mal entsprechend verkleinert, damit das ganze in wenigen Sekunden fertig sein kann.

Aber mit FFT werde ich wohl nicht unter N³ * log(N) kommen, aber das ist ja deutlich besser als N⁶ mit der naiven Variante oder so N⁵, die man mit diesen Schalen vielleicht hinbekommen könnte.

 

[Dalgliesh]

Aber mit FFT werde ich wohl nicht unter N³ * log(N) kommen

Natürlich nicht. Du mußt jeden Punkt mindestens einmal anschauen, dann hast Du schon N^3. Und log N kannst Du vergessen.