Solução rápida e explícita para , , número de condição baixo

9

Estou procurando uma solução explícita rápida (ouso dizer ideal?) Para o problema real linear 3x3, , . Ax=bA∈R3×3,b∈R3

Matrix A é geral, mas próximo à matriz de identidade com um número de condição próximo a 1. Como b são na verdade medições de sensores com cerca de 5 dígitos de precisão, não me importo de perder vários dígitos devido a números. problemas.

Obviamente, não é difícil encontrar uma solução explícita com base em vários métodos, mas se há algo que se mostrou ideal em termos de contagem de FLOPS, isso seria ideal (afinal, todo o problema provavelmente caberá nos registros FP).

(Sim, essa rotina é chamada com frequência . Eu já me livrei das frutas baixas e isso é o próximo na minha lista de perfis ...)

linear-solver reference-request Damien
fonte
Cada usado apenas uma vez ou existem vários sistemas lineares com a mesma matriz? Isso mudaria os custos. A
Federico Poloni
Nesse caso, A é usado apenas uma vez.
Damien

Respostas:

14

Você não pode vencer uma fórmula explícita. Você pode anotar as fórmulas da solução em um pedaço de papel. Deixe o compilador otimizar as coisas para você. Qualquer outro método quase inevitavelmente terá instruções ou loops (por exemplo, para métodos iterativos) que tornarão seu código mais lento do que qualquer código de linha reta.x=A−1biffor

Wolfgang Bangerth
fonte
9

Como a matriz está tão próxima da identidade, as seguintes séries Neumann convergirão muito rapidamente:

A−1=∑k=0∞(I−A)k

Dependendo da precisão necessária, pode até ser bom o suficiente para truncar após 2 termos:

A−1≈I+(I−A)=2I−A.

Isso pode ser um pouco mais rápido que uma fórmula direta (como sugerido na resposta de Wolfgang Bangerth), embora com muito menos precisão.

Você pode obter mais precisão com três termos:

A−1≈I+(I−A)+(I−A)2=3I−3A+A2

mas se você escrever a fórmula entrada por entrada para , estará analisando uma quantidade comparável de operações de ponto flutuante como a fórmula inversa direta da matriz 3x3 (não é necessário faça uma divisão, o que ajuda um pouco).(3I−3A+A2)b

Nick Alger
fonte
As divisões ainda são mais caras que os outros fracassos? Eu pensei que era uma relíquia do passado.
Federico Poloni 10/10
Divisões não faça dutos bem um algumas arquiteturas (ARM é o exemplo contemporâneo)
Damien
@FedericoPoloni Com o Cuda, você pode ver o rendimento das instruções aqui , é seis vezes maior para multiplicações / adições do que para divisões.
Kirill
@ Damien e Kirill vejo, obrigado pelos ponteiros.
Federico Poloni 10/10
5

Os FLOPS contam com base nas sugestões acima:

  • LU, sem articulação:

    • Mul = 11, Div / Recip = 6, Adição / Sub = 11, Total = 28; ou
    • Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30

  • Eliminação Gaussiana com substituição traseira, sem articulação:

    • Mul = 11, Div / Recip = 6, Adição / Sub = 11, Total = 28; ou
    • Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30

  • Regra de Cramer via expansão de cofator

    • Mul = 24, Div = 3, Adição / Sub = 15, Total = 42; ou
    • Mul = 27, Div = 1, Adição / Sub = 15, Total = 43

  • Inverso explícito e multiplique:

    • Mul = 30, Div = 3, Add / Sub = 17, Total = 50; ou
    • Mul = 33, Div = 1, Adição / Sub = 17, Total = 51

Prova de conceitos do MATLAB:

Regra de Cramer via expansão de cofator :

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 (sem rotação) e substituição traseira:

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 e multiplicado:

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

Eliminação 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: sinta-se à vontade para adicionar seus próprios métodos e contagens a esta postagem.

Damien
fonte
Você calculou o tempo que levou para resolver com os dois métodos?
nicoguaro
Não. O código acima não será executado rapidamente. O ponto de que era para ter uma explícita ALETA contar e fornecer o código de revisão no caso eu perdi alguma coisa,
Damien
Na LU, 5 divisões podem ser convertidas em 5 MULs à custa de 2 operações recíprocas extras (ou seja, 1 / U11 e 1 / U22). Isso será específico do arco para saber se há um ganho a ser feito lá.
Damien
2
UMA-1 1b2b-UMAbUMA-1 1b3(b-UMAb)+UMA2bUMA-1 1bparece ser 33 multiplicações, 17 adições / subtrações e 1 divisão. Como eu disse, meus números podem estar desativados, por isso, verifique novamente.
Geoff Oxberry 10/10
@ GeoffOxberry, analisarei e reportarei.
Damien
4

Provavelmente a regra de Cramer. Se você pode evitar a rotação, talvez a fatoração da LU; é uma matriz 3x3, então seria fácil desenrolar os loops manualmente. Qualquer outra coisa provavelmente envolverá ramificação, e duvido que um método de subespaço de Krylov convergisse com frequência suficiente em 1 ou 2 iterações para que valesse a pena.

Geoff Oxberry
fonte