Problem mit Matrix Multiplikation: Rechenzeit vs. Format (Genauigkeit)

[teo (Gast)]

Hallo an alle!

Ich versuche gerade im Rahmen meiner Diplomarbeit einen bestimmten 
Algorithmus auf dem dsPIC33 zu implementieren. Es funktioniert auch 
soweit ganz gut, bloß hat sich das Problem einer langen Rechenzeit 
ergeben, was ich auch erwartet hatte. Im Prinzip besteht der Algorithmus 
aus mehreren Matrixmultiplikationen, das sind natürlich die Dinger die 
Rechenzeit beanspruchen.
Nun bietet der dsPIC  die tolle dsp.h Library an, welche saumäßig 
schnell rechnen kann, aber eben nur im fractional-Format.

Um die MatrixMultiply() funktion nutzen zu können, müsste ich 
theoretisch alle Matrizen skalieren damit die Werte im Bereich von 
-1...1 liegen.

Problem: Meine Matrizen sind dermaßen krumm, dass eine Skalierung auf 
den größten Wert die anderen Werte ganz stark verkleinern würde und im 
Q.15 Format nicht mehr darstellbar wären.
Ich möchte euch ein Beispiel zeigen:


 const float A[3][3] = {{1,  1E-3,   0},
                {-1194172.438500119,  -5.911870074,1194172.438500119}, 
{0,                 0,                1}
              };  // A

 float P[3][3] = {{1,  0,   0},
                  {0,  1,   0},
                  {0,  0,   1};




Eine Zeile in dem Algorithmus sieht so aus:    Px = A*P*A' + Q;

Q hat sehr sehr kleine Werte im Bereich 10^-5....10^-15 oder kleiner. 
Aber auf keine Fall zu Vernachlässigen!


Ich würde gerne schnell rechnen, kann es aber nicht weil Genauigkeit 
verloren geht. Ich habe schon überlegt, ob ich z.B. A*P*A' im fractional 
Format rechne, das Ergebnis mit Fract2Float() umwandle und +Q im float 
addiere.

Eine Skalierung von A unnd P würde dennoch das Ergebnis stark verändern.

Eine weitere Überlegung war den PIC32 zu nehmen, aber der wiederum hat 
noch keine Matrix Multiplikation in der dsp library. Abgesehen davon 
liegt das Q.32 format auch nur zwischen -1 und 1, die Genauigkeit ist 
natürlich viel höher!

Ich hoffe ich habe mein Problem einigermaßen gut erklärt, falls nicht, 
fragt bitte!

Kennt sich jemand mit diesem Problem aus? Habt ihr Ideen?
Vielen vielen Dank für eure Hilfe!

Grüße,

teo

[Detlef _. (detlef_a)]

Du bist nicht der erste, der sich darüber Gedanken macht.

"kalman filter dynamic range" liefert z.B. den:

http://focus.ti.com.cnXX/cnXX/lit/an/spra249/spra249.pdf

(Die vier XX rausmachen, sonst schlägt der Spamfilter zu :-/ )

Schnelleren Prozessor nehmen.
Hast Du die Matrixsymmetrie ausgenutzt ?
Die Verstärkungen des Filters lassen sich offline vorberechnen.

Cheers
Detlef

[teo (Gast)]

Hi Detlef!

HErzlichen Dank für die Antwort. Du bist scheinbar ein Fuchs, ich hatte 
von Kalman Filter nichts erwähnt, aber du hast es scheinbar erkannt ;-) 
In der Tat geht es darum.

Schnelleren PRozessor nehmen: Ja, will ich, mach ich. Den PIC32, der ist 
schonmal doppelt so schnell. ABER: ich würde dennoch auf eine Abtastrate 
des AD-Wandlers von 2 kHz kommen (jetzt 1 kHz) und das ist deutlich zu 
wenig. Ich werde mir den Link mal ansehen, es scheint um das gleichte 
Thema zu gehen. Vielen Dank.

Was genau meinst du mit MAtrix-Symetrie ausnutzen? Zum Beispiel die 
MAtrix A... was könnte ich tun?

Bis denn!

teo

[Detlef _. (detlef_a)]

Hi,

nein, A ist ja beliebig, aber P und P* sind symmetrisch, K*R*K' macht 
auch ne symmetrische Matrix, da kannst Du dann fast Faktor 2 sparen.

Cheers
Detlef

[teo (Gast)]

Stimmt! Das muss ich ausnutzen. Das könnte schon einiges bringen. Bis 
jetzt habe ich halt ganz "normale" Matrixmultiplikationen benutzt, die 
ich selbst erstellt habe:



void mat_mult33x33(float M1[3][3],
            float M2[3][3],
     float M3[3][3]) {

 int i, j, k;

    for (i=0; i<3; i++) {
    for (j=0; j<3; j++) {
      M3[i][j]=0;
      for (k=0; k<3; k++) {   // Multiplikation
    M3[i][j] += M1[i][k] * M2[k][j];
      }
  }
    }

}


Um es jetzt auf Speed zu optimieren könnte ich für jede Multiplikation 
eine angepasste (Symmetrieausnutzende) Funktion schreiben.

Bringt es was die Kunktionen gar nicht aufzurufen, sondern direkt die 
Schleifen abzuarbeiten?


Vielen Dank für den Tip!

Gruß,

teo

[Detlef _. (detlef_a)]

>>sondern direkt die Schleifen abzuarbeiten?

ja, keine 'one serves it all' Funktion benutzen.

- Bei symmetrischen Matrizen die redundante Hälfte weder berechnen noch 
speichern.

- Vorsicht: Compiler/Plattformabhängig geht man günstiger die Spalten 
oder Zeilen der Matrix lang. Am liebsten gleich auf 2D arrays verzichten 
und die Elemente linear hintereinander speichern und bearbeiten (kann 
der DSPic prefetch ?)

- Systeme dritter Ordnung müssen doch mit >2ks/s gehen, denke ich, ohne 
das nachgerechnet zu haben.

Cheers
Detlef