Solución explícita rápida para

9

Estoy buscando una solución explícita rápida (¿me atrevo a decir que es óptima?) El problema real lineal de 3x3, , . $\mathbf{A}\mathbf{x} = \mathbf{b}$ $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$

La matriz es general, pero está cerca de la matriz de identidad con un número de condición cercano a 1. Debido a que son en realidad mediciones de sensores con aproximadamente 5 dígitos de precisión, no me importa perder varios dígitos debido a los valores numéricos cuestiones. $\mathbf{A}$ $\mathbf{b}$

Por supuesto, no es difícil encontrar una solución explícita basada en cualquier número de métodos, pero si hay algo que se ha demostrado que es óptimo en términos de conteo de FLOPS, sería ideal (después de todo, todo el problema probablemente encajará en los registros FP!

(Sí, esta rutina se llama con frecuencia . Ya me he librado de la fruta baja y esta es la siguiente en mi lista de perfiles ...)

linear-solver reference-request Damien
fuente

¿Cada

usa solo una vez, o hay múltiples sistemas lineales con la misma matriz? Esto cambiaría los costos.

A

$A$

Federico Poloni

En este caso, A se usa solo una vez.

Damien

14

No se puede superar una fórmula explícita. Puede escribir las fórmulas para la solución en una hoja de papel. Deje que el compilador optimice las cosas por usted. Cualquier otro método tendrá casi inevitablemente sentencias o bucles (por ejemplo, para métodos iterativos) que harán que su código sea más lento que cualquier código de línea recta. $x=A^{-1}b$ iffor

Wolfgang Bangerth
fuente

9

Dado que la matriz está tan cerca de la identidad, las siguientes series de Neumann convergerán muy rápidamente:

{UNA}^{- 1} = \sum_{k = 0 0}^{\infty} (yo - UNA)^{k}

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

Dependiendo de la precisión requerida, incluso podría ser lo suficientemente bueno para truncar después de 2 términos:

{UNA}^{- 1} \approx yo + (yo - UNA) = 2 yo - UNA .

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

Esto podría ser un poco más rápido que una fórmula directa (como se sugiere en la respuesta de Wolfgang Bangerth), aunque con mucha menos precisión.

Puede obtener más precisión con 3 términos:

{UNA}^{- 1} \approx yo + (yo - UNA) + (yo - UNA)^{2} = 3 yo - 3 UNA + {UNA}^{2}

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

pero si escribe la fórmula de entrada por entrada para , está viendo una cantidad comparable de operaciones de coma flotante como la fórmula inversa directa de matriz 3x3 (no tiene que hacer una división, lo que ayuda un poco). $(3I - 3A + A^2)b$

Nick Alger
fuente

¿Las divisiones son aún más caras que los otros fracasos? Pensé que era una reliquia del pasado.

Federico Poloni

Las divisiones no se canalizan bien en algunas arquitecturas (ARM es el ejemplo contemporáneo)

Damien

@FedericoPoloni Con Cuda, puede ver el rendimiento de la instrucción aquí , es seis veces mayor para las multiplicaciones / adiciones que para las divisiones.

Kirill

@Damien y Kirill Ya veo, gracias por los consejos.

Federico Poloni

5

Los FLOPS cuentan según las sugerencias anteriores:

LU, sin pivote:
- Mul = 11, Div / Recip = 6, Agregar / Sub = 11, Total = 28; o
- Mul = 16, Div / Recip = 3, Agregar / Sub = 11, Total = 30
Eliminación gaussiana con sustitución de espalda, sin pivote:
- Mul = 11, Div / Recip = 6, Agregar / Sub = 11, Total = 28; o
- Mul = 16, Div / Recip = 3, Agregar / Sub = 11, Total = 30
La regla de Cramer a través de la expansión del cofactor
- Mul = 24, Div = 3, Agregar / Sub = 15, Total = 42; o
- Mul = 27, Div = 1, Agregar / Sub = 15, Total = 43
Inverso explícito y luego multiplica:
- Mul = 30, Div = 3, Agregar / Sub = 17, Total = 50; o
- Mul = 33, Div = 1, Agregar / Sub = 17, Total = 51

Prueba de conceptos de MATLAB:

Regla de Cramer a través de la expansión de 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 (sin pivote) y sustitución posterior:

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

Inverso explícito y luego multiplicar:

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

Eliminación gaussiana:

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

Nota: No dude en agregar sus propios métodos y recuentos a esta publicación.

Damien
fuente

¿Calculaste los tiempos que tomó resolver con los dos métodos?

nicoguaro

No. El código anterior no se ejecutará rápidamente en absoluto. El objetivo era obtener un recuento explícito de FLOPS y proporcionar el código para su revisión en caso de que me perdiera algo,

Damien

En LU, 5 divisiones se pueden convertir en 5 MUL a expensas de 2 operaciones recíprocas adicionales (es decir, 1 / U11 y 1 / U22). Eso será específico del arco en cuanto a si hay una ganancia que hacer allí.

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$ través de esta fórmula explícitaparece ser 33 multiplicaciones, 17 sumas / restas y 1 división. Como dije, mis números pueden estar apagados, por lo que es posible que desee verificar dos veces.

Geoff Oxberry

@ GeoffOxberry, lo investigaré e informaré.

Damien

4

Probablemente la regla de Cramer. Si puede evitar pivotar, quizás factorización LU; Es una matriz de 3x3, por lo que desenrollar los bucles manualmente sería fácil. Cualquier otra cosa probablemente implique ramificación, y dudo que un método de subespacio de Krylov converja con la frecuencia suficiente en 1 o 2 iteraciones para que valga la pena.

Geoff Oxberry
fuente

Solución explícita rápida para

Respuestas: