Solution explicite rapide pour

9

Je recherche une solution explicite rapide (oserais-je dire optimale?) Au problème réel linéaire 3x3, , . $\mathbf{A}\mathbf{x} = \mathbf{b}$ $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$

La matrice est générale, mais proche de la matrice d'identité avec un numéro de condition proche de 1. Parce que sont en fait des mesures de capteur avec environ 5 chiffres de précision, cela ne me dérange pas de perdre plusieurs chiffres en raison de problèmes numériques. $\mathbf{A}$ $\mathbf{b}$

Bien sûr, il n'est pas difficile de trouver une solution explicite basée sur un certain nombre de méthodes, mais s'il y a quelque chose qui s'est révélé optimal en termes de nombre de FLOPS, ce serait idéal (après tout, tout le problème rentrera probablement dans les registres FP!).

(Oui, cette routine est souvent appelée . Je me suis déjà débarrassé des fruits bas et c'est la prochaine dans ma liste de profilage ...)

linear-solver reference-request Damien
la source

Chaque

utilisé une seule fois ou existe-t-il plusieurs systèmes linéaires avec la même matrice? Cela changerait les coûts.

A

$A$

Federico Poloni

Dans ce cas, A n'est utilisé qu'une seule fois.

Damien

14

Vous ne pouvez pas battre une formule explicite. Vous pouvez écrire les formules de la solution sur une feuille de papier. Laissez le compilateur optimiser les choses pour vous. Toute autre méthode aura presque inévitablement des instructions ou des boucles (par exemple, pour les méthodes itératives) qui rendront votre code plus lent que tout code en ligne droite. $x=A^{-1}b$ iffor

Wolfgang Bangerth
la source

9

La matrice étant si proche de l'identité, les séries Neumann suivantes convergeront très rapidement:

{UNE}^{- 1} = \sum_{k = 0}^{\infty} (je - UNE)^{k}

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

Selon la précision requise, il pourrait même être suffisant de tronquer après 2 termes:

{UNE}^{- 1} \approx je + (je - UNE) = 2 je - UNE .

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

Cela pourrait être légèrement plus rapide qu'une formule directe (comme suggéré dans la réponse de Wolfgang Bangerth), mais avec beaucoup moins de précision.

Vous pouvez obtenir plus de précision avec 3 termes:

{UNE}^{- 1} \approx je + (je - UNE) + (je - UNE)^{2} = 3 je - 3 UNE + {UNE}^{2}

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

mais si vous écrivez la formule entrée par entrée pour , vous regardez une quantité comparable d'opérations en virgule flottante comme formule inverse directe de matrice 3x3 (vous n'avez pas à faire une division cependant, ce qui aide un peu). $(3I - 3A + A^2)b$

Nick Alger
la source

Les divisions sont-elles encore plus chères que les autres flops? Je pensais que c'était une relique du passé.

Federico Poloni du

Les divisions ne canalisent pas bien certaines architectures (ARM est l'exemple contemporain)

Damien

@FedericoPoloni Avec Cuda, vous pouvez voir le débit des instructions ici , il est six fois plus élevé pour les multiplications / ajouts que pour les divisions.

Kirill

@Damien et Kirill je vois, merci pour les pointeurs.

Federico Poloni du

5

Les FLOPS comptent sur la base des suggestions ci-dessus:

LU, pas de pivotement:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; ou
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Élimination gaussienne avec substitution arrière, sans pivot:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; ou
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
La règle de Cramer via l'expansion du cofacteur
- Mul = 24, Div = 3, Add / Sub = 15, Total = 42; ou
- Mul = 27, Div = 1, Add / Sub = 15, Total = 43
Inverse explicite puis multipliez:
- Mul = 30, Div = 3, Add / Sub = 17, Total = 50; ou
- Mul = 33, Div = 1, Add / Sub = 17, Total = 51

Preuve de concepts MATLAB:

Règle de Cramer via l'extension Cofactor :

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 (pas de pivotement) et substitution arrière:

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];

Inverse explicite puis multipliez:

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];

Élimination gaussienne:

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];

Remarque: n'hésitez pas à ajouter vos propres méthodes et décomptes à ce message.

Damien
la source

Avez-vous calculé le temps nécessaire pour résoudre les deux méthodes?

nicoguaro

Non. Le code ci-dessus ne s'exécutera pas rapidement du tout. Le but était d'obtenir un décompte explicite des FLOPS et de fournir le code pour examen au cas où j'aurais raté quelque chose,

Damien

En LU, 5 divisions peuvent être converties en 5 MUL au détriment de 2 opérations réciproques supplémentaires (soit 1 / U11 et 1 / U22). Ce sera spécifique à l'arc pour savoir s'il y a un gain à y faire.

Damien

2

En supposant que je n'ai pas mal compté, approximant

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$ semble être 33 multiplications, 17 additions / soustractions et 1 division. Comme je l'ai dit, mes numéros sont peut-être faux, alors vous voudrez peut-être revérifier.

Geoff Oxberry

@GeoffOxberry, je vais l'examiner et faire un rapport.

Damien

4

Probablement la règle de Cramer. Si vous pouvez éviter le pivotement, peut-être la factorisation LU; c'est une matrice 3x3, donc dérouler les boucles manuellement serait facile. Tout le reste impliquera probablement une ramification, et je doute qu'une méthode de sous-espace de Krylov converge assez souvent en 1 ou 2 itérations pour que cela en vaille la peine.

Geoff Oxberry
la source

Solution explicite rapide pour

Réponses: