Schnelle explizite Lösung für

9

Ich suche nach einer schnellen (wage ich zu sagen, optimal?) Expliziten Lösung für das lineare 3x3-Problem, , . $\mathbf{A}\mathbf{x} = \mathbf{b}$ $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$

Matrix ist allgemein, aber nahe an der Identitätsmatrix mit einer Bedingungsnummer nahe 1. Da tatsächlich Sensormessungen mit einer Genauigkeit von etwa 5 Stellen sind, macht es mir nichts aus, mehrere Zahlen aufgrund von Zahlen zu verlieren Probleme. $\mathbf{A}$ $\mathbf{b}$

Natürlich ist es nicht schwer, eine explizite Lösung zu finden, die auf einer beliebigen Anzahl von Methoden basiert, aber wenn sich etwas als optimal in Bezug auf die FLOPS-Anzahl erwiesen hat, wäre dies ideal (schließlich das gesamte Problem) wird wahrscheinlich in die FP-Register passen!).

(Ja, diese Routine wird oft aufgerufen . Ich habe bereits tief hängende Früchte entfernt und dies steht als nächstes in meiner Profilliste ...)

linear-solver reference-request Damien
quelle

Wird jedes nur einmal verwendet oder gibt es mehrere lineare Systeme mit derselben Matrix? Dies würde die Kosten ändern.

A

$A$

Federico Poloni

In diesem Fall wird A nur einmal verwendet.

Damien

14

Sie können eine explizite Formel nicht schlagen. Sie können die Formeln für die Lösung auf ein Blatt Papier schreiben . Lassen Sie den Compiler die Dinge für Sie optimieren. Jede andere Methode enthält fast zwangsläufig Anweisungen oder Schleifen (z. B. für iterative Methoden), die Ihren Code langsamer machen als jeder gerade Code. $x=A^{-1}b$ iffor

Wolfgang Bangerth
quelle

9

Da die Matrix der Identität so nahe kommt, konvergieren die folgenden Neumann-Reihen sehr schnell:

{EIN}^{- - 1} = \sum_{k = 0}^{\infty} (ich - - EIN)^{k}

$A^{-1} = \sum_{k=0}^\infty (I-A)^k$

Abhängig von der erforderlichen Genauigkeit kann es sogar gut genug sein, um nach 2 Begriffen abzuschneiden:

{EIN}^{- - 1} \approx ich + (ich - - EIN) = 2 ich - - EIN .

$A^{-1} \approx I + (I - A) = 2I - A.$

Dies ist möglicherweise etwas schneller als eine direkte Formel (wie in Wolfgang Bangerths Antwort vorgeschlagen), allerdings mit viel geringerer Genauigkeit.

Sie könnten mehr Genauigkeit mit 3 Begriffen erhalten:

{EIN}^{- - 1} \approx ich + (ich - - EIN) + (ich - - EIN)^{2} = 3 ich - - 3 EIN + {EIN}^{2}

$A^{-1} \approx I + (I - A) + (I-A)^2 = 3I - 3A + A^2$

Wenn Sie jedoch die Eintrag-für-Eintrag-Formel für schreiben , sehen Sie eine vergleichbare Anzahl von Gleitkommaoperationen wie die direkte inverse 3x3-Matrixformel (müssen Sie nicht mache aber eine Teilung, was ein wenig hilft). $(3I - 3A + A^2)b$

Nick Alger
quelle

Sind Divisionen immer noch teurer als die anderen Flops? Ich dachte, es sei ein Relikt der Vergangenheit.

Federico Poloni

Divisionen leiten einige Architekturen nicht gut ab (ARM ist das zeitgenössische Beispiel)

Damien

@FedericoPoloni Mit Cuda können Sie hier den Befehlsdurchsatz sehen , der für Multiplikationen / Additionen sechsmal höher ist als für Divisionen.

Kirill

@ Damien und Kirill Ich verstehe, danke für die Hinweise.

Federico Poloni

5

FLOPS zählen basierend auf den obigen Vorschlägen:

LU, kein Schwenken:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; oder
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Gaußsche Elimination mit Rücksubstitution, kein Schwenken:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; oder
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Cramer-Regel über Cofaktor-Erweiterung
- Mul = 24, Div = 3, Add / Sub = 15, Total = 42; oder
- Mul = 27, Div = 1, Add / Sub = 15, Total = 43
Explicit Inverse dann multiplizieren:
- Mul = 30, Div = 3, Add / Sub = 17, Total = 50; oder
- Mul = 33, Div = 1, Add / Sub = 17, Total = 51

MATLAB Proof-of-Concepts:

Cramer-Regel über Cofaktor-Erweiterung :

function k = CramersRule(A, m)
%
% FLOPS:
%
% Multiplications:        24
% Subtractions/Additions: 15
% Divisions:               3
%
% Total:                  42

a = A(1,1);
b = A(1,2);
c = A(1,3);

d = A(2,1);
e = A(2,2);
f = A(2,3);

g = A(3,1);
h = A(3,2);
i = A(3,3);

x = m(1);
y = m(2);
z = m(3);

ei = e*i;
fh = f*h;

di = d*i;
fg = f*g;

dh = d*h;
eg = e*g;

ei_m_fh = ei - fh;
di_m_fg = di - fg;
dh_m_eg = dh - eg;

yi = y*i;
fz = f*z;

yh = y*h;
ez = e*z;

yi_m_fz = yi - fz;
yh_m_ez = yh - ez;

dz = d*z;
yg = y*g;

dz_m_yg = dz - yg;
ez_m_yh = ez - yh;


det_a = a*ei_m_fh - b*di_m_fg + c*dh_m_eg;
det_1 = x*ei_m_fh - b*yi_m_fz + c*yh_m_ez;
det_2 = a*yi_m_fz - x*di_m_fg + c*dz_m_yg;
det_3 = a*ez_m_yh - b*dz_m_yg + x*dh_m_eg;


p = det_1 / det_a;
q = det_2 / det_a;
r = det_3 / det_a;

k = [p;q;r];

LU (kein Schwenken) und Rücksubstitution:

function [x, y, L, U] = LUSolve(A, b)
% Total FLOPS count:     (w/ Mods)
%
% Multiplications:  11    16
% Divisions/Recip:   6     3
% Add/Subtractions: 11    11
% Total =           28    30
%

A11 = A(1,1);
A12 = A(1,2);
A13 = A(1,3);

A21 = A(2,1);
A22 = A(2,2);
A23 = A(2,3);

A31 = A(3,1);
A32 = A(3,2);
A33 = A(3,3);

b1 = b(1);
b2 = b(2);
b3 = b(3);

L11 = 1;
L22 = 1;
L33 = 1;

U11 = A11;
U12 = A12;
U13 = A13;

L21 = A21 / U11;
L31 = A31 / U11;

U22 = (A22 - L21*U12);
L32 = (A32 - L31*U12) / U22;

U23 = (A23 - L21*U13);

U33 = (A33 - L31*U13 - L32*U23);

y1 = b1;
y2 = b2 - L21*y1;
y3 = b3 - L31*y1 - L32*y2;

x3 = (y3                  ) / U33;
x2 = (y2 -          U23*x3) / U22;
x1 = (y1 - U12*x2 - U13*x3) / U11;

L = [ ...
    L11,   0,   0;
    L21, L22,   0;
    L31, L32, L33];

U = [ ...
    U11, U12, U13;
      0, U22, U23;
      0,   0, U33];

x = [x1;x2;x3];
y = [y1;y2;y3];

Explizite Umkehrung und dann Multiplikation:

function x = ExplicitInverseMultiply(A, m)
%
% FLOPS count:                  Alternative
%
% Multiplications:        30            33
% Divisions:               3             1
% Additions/Subtractions: 17            17
% Total:                  50            51


a = A(1,1);
b = A(1,2);
c = A(1,3);

d = A(2,1);
e = A(2,2);
f = A(2,3);

g = A(3,1);
h = A(3,2);
i = A(3,3);

ae = a*e;
af = a*f;
ah = a*h;
ai = a*i;

bd = b*d;
bf = b*f;
bg = b*g;
bi = b*i;

cd = c*d;
ce = c*e;
cg = c*g;
ch = c*h;

dh = d*h;
di = d*i;

eg = e*g;
ei = e*i;

fg = f*g;
fh = f*h;

dh_m_eg = (dh - eg);
ei_m_fh = (ei - fh);
fg_m_di = (fg - di);

A = ei_m_fh;
B = fg_m_di;
C = dh_m_eg;
D = (ch - bi);
E = (ai - cg);
F = (bg - ah);
G = (bf - ce);
H = (cd - af);
I = (ae - bd);

det_A = a*ei_m_fh + b*fg_m_di + c*dh_m_eg;

x1 =  (A*m(1) + D*m(2) + G*m(3)) / det_A;
x2 =  (B*m(1) + E*m(2) + H*m(3)) / det_A;
x3 =  (C*m(1) + F*m(2) + I*m(3)) / det_A;

x = [x1;x2;x3];

Gaußsche Eliminierung:

function x = GaussianEliminationSolve(A, m)
%
% FLOPS Count:      Min   Alternate
%
% Multiplications:  11    16
% Divisions:         6     3
% Add/Subtractions: 11    11
% Total:            28    30
%

a = A(1,1);
b = A(1,2);
c = A(1,3);

d = A(2,1);
e = A(2,2);
f = A(2,3);

g = A(3,1);
h = A(3,2);
i = A(3,3);

b1 = m(1);
b2 = m(2);
b3 = m(3);

% Get to echelon form

op1 = d/a;

e_dash  = e  - op1*b;
f_dash  = f  - op1*c;
b2_dash = b2 - op1*b1;

op2 = g/a;

h_dash  = h  - op2*b;
i_dash  = i  - op2*c;
b3_dash = b3 - op2*b1; 

op3 = h_dash / e_dash;

i_dash2  = i_dash  - op3*f_dash;
b3_dash2 = b3_dash - op3*b2_dash;

% Back substitution

x3 = (b3_dash2                  ) / i_dash2;
x2 = (b2_dash        - f_dash*x3) / e_dash;
x1 = (b1      - b*x2 -      c*x3) / a;

x = [x1 ; x2 ; x3];

Hinweis: Sie können diesem Beitrag auch Ihre eigenen Methoden und Zählungen hinzufügen.

Damien
quelle

Haben Sie die Zeit berechnet, die zum Lösen mit den beiden Methoden benötigt wurde?

Nicoguaro

Nein. Der obige Code wird überhaupt nicht schnell ausgeführt. Der Punkt war, eine explizite FLOPS-Zählung zu erhalten und den Code zur Überprüfung bereitzustellen, falls ich etwas verpasst habe

Damien

In LU können 5 Divisionen auf Kosten von 2 zusätzlichen wechselseitigen Operationen (dh 1 / U11 und 1 / U22) in 5 MULs umgewandelt werden. Das wird bogenspezifisch sein, ob dort ein Gewinn erzielt werden kann.

Damien

2

A^{- 1} b

$A^{-1}b$

2 b - A b

$2b - Ab$

A^{- 1} b

$A^{-1}b$

3 (b - A b) + A^{2} b

$3(b - Ab) + A^2b$

A^{- 1} b

$A^{-1}b$ über diese explizite FormelEs sieht aus wie 33 Multiplikationen, 17 Additionen / Subtraktionen und 1 Division. Wie ich bereits sagte, sind meine Nummern möglicherweise nicht korrekt, daher sollten Sie dies überprüfen.

Geoff Oxberry

@GeoffOxberry, ich werde es untersuchen und berichten.

Damien

4

Wahrscheinlich Cramers Regel. Wenn Sie ein Schwenken vermeiden können, möglicherweise eine LU-Faktorisierung. Es ist eine 3x3-Matrix, daher wäre es einfach, die Schleifen manuell abzuwickeln. Alles andere wird wahrscheinlich eine Verzweigung beinhalten, und ich bezweifle, dass eine Krylov-Subraummethode in 1 oder 2 Iterationen oft genug konvergiert, damit es sich lohnt.

Geoff Oxberry
quelle

Schnelle explizite Lösung für

Antworten: