Matrisinvers 8x8

[Korken]

Godagens!

Jag arbetar på det nya filtret för KFly och har stött på ett problem, men jag tänkte ta lite bakgrund först.
Mitt filter som jag utvecklar är DCM baserat, men man tar ut raderna i DCM-matrisen mha två stycken kaskad Kalman-filters.
Detta ger att mina state-vektorer blir 8 element (3 för DCM raden, 3 för vikelhastigheter och 2 för bias) och på så vis för att räkna ut inversen av kovarians-matrisen så måste jag ta en matris-invers på en generell 8x8 matris.

Så mitt problem är, hur gör jag enklast detta (ta inversen)?
Ska man använda matrisinvers-koder eller hårdkoda inversen (inversen kan lätt beräknas med tex Maple så man får formler för varje element)? Finns det några andra lösningar?
Eller attackerar jag kanske detta problem med fel vinkel?

Jag tror nästan mest på hårdkoda den, men då så kan man inte snabbt och enkelt ändra i filtret.
Däremot så kommer matris-inverskoderna bli väldigt ineffektiva då de inte tar hänsyn till alla nollor som finns i matrisen.

Detta görs på en STM32F405 med FPU och jag gör allt med 32-bitars flyttal.

[PF013]

En STM32F405? Då kan du väl nästan lika gärna porta Matlab till den?!

[Andax]

Med tanke på att matrisen innehåller uppåt kanske 2/3 nollor tror jag en hårdkodad variant är bäst.

[Korken]

Jag tog och slog igenom alla ekvation (utan bias) för att få en uppfattning på matrisen.
Något säger mig att jag måste ta den långa vägen för att ta inversen.

Det är S matrisen som ska inverteras. Tar jag inverse i Maple blir uttrycket ca 4st A4-sidor långt.
Jag ska testa och se om det går att förenkla.

eqn.png

[Andax]

Tja, matrisen du ville invertera hade visst inte så många nollor som jag föreställt mig. Däremot är den symmetrisk, vilket innebär att inversen också är det, så det reducera beräkningarna en hel del.

Så när du sätter testar maple kan du ju skapa en symbolisk matris med elementen (obersvera att jag hoppade över bokstaven o för att inte förvilla med 0 noll, och därför blev sista variabeln p istället)

Kod: Markera allt

[ a b c 0 d e
  b f g h 0 i
  c g j k l 0
  0 h k m 0 0 
  d 0 l 0 n 0
  e i 0 0 0 p ]

Du har ju fortfarande 12 st nollor i matrisen av 36 element vilket också förenklar en del.

[Korken]

Mycket bra observation! Jag hade inte sätt att den var symmetrisk...
Det blev ganska så mycket bättre, jag tog och multiplicerade med determinanten för att få ett enklare uttryck.
Det är överskådligt nu iaf. Jag behöver få matrisen på en mer överskådlig form för att kunna hitta förenklingar i den (tex se symmetrin).

eqn2.png

Determinanten:

det.png

[Korken]

Nu har jag testat att iterera systemet ett par varv och märkte en väldigt tråkig sak.
S matrisen blir en generell matris efter andra iterationen, så jag kommer köra med en matrisinverskod.

Nåväl, det kunde ha blivit snyggare iaf.

[thebolt]

Jag tror att man generellt kan säga (beror lite på processor-arkitektur osv) att 3x3 eller 4x4 så är det värt att ha en explicit form på lösningen av inversen enl. Cramers regel, vid större matriser får du färre totalt antal operationer med LU-faktorisering eller gausisk eliminering. Du får en liten förskjutning mot att Cramers regel är bra om du har mycket 0:or, men är matrisen väldigt sparse finns det andra metoder också.

[Korken]

Precis som du säger så kommer jag gör på någon av de visen. Jag lutar mer mod LU-faktorisering, det tycks vara enklast att implementera.
Har någon exempel på algoritmer som använder LU-faktorisering?

Edit: Hittade en massa. Google is your friend.

[thebolt]

Menar du algoritmer som använder LU så är det ju i princip alltid om du vill invertera, eller menar du "implementationer av LU-faktorisering"? Numerical Recipies in C (finns gratis online) har en version, den är inte snygg och inte den absolut snabbaste, men duger som grund att rensa upp

[Korken]

Om du råkar sitta på kod som gör LU-faktorisering och sedan inversen så får du gärna dela med dig.
Jag testade att implementera en LU-faktoriseringskod och den fungerar riktigt bra. Exekverar på ca 500 klockcykler för en 8x8 matris. (140 multiplikationer och 28 divisioner)

Edit: Har någon en kod som tar inversen av trianguella matriser? Jag hittar ingen sådan och har ingen direkt lust att skriva en från grunden.
Edit2: Kom precis på ett sätt att lösa det på rekursivt. Ska se om det inte fungerar.

[Korken]

Nu börjar jag få något som kan fungera. Jag gjorde en snabb rekursion för hand och teorin tycks fungera.
Det lite fina med den som jag upptäckte är att man bara behöver så många divisioner som det finns element på diagonalen.
Det är mycket bra då det tar 7ggr längre tid att utföra divisioner än multiplikationer.

Ska ta och koda den nu och se vad som händer.

Kladdtavlan:

20120212_142311.jpg