Szybkie, wyraźne rozwiązanie dla

9

Szukam szybkiego (ośmielę się powiedzieć optymalnego?) Jawnego liniowego problemu 3x3, , . $\mathbf{A}\mathbf{x} = \mathbf{b}$ $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$

Macierz jest ogólna, ale zbliżona do macierzy tożsamości z numerem warunku bliskim 1. Ponieważ to tak naprawdę pomiary czujnika z około 5 cyfrową precyzją, nie mam nic przeciwko utracie kilku cyfr z powodu liczb problemy. $\mathbf{A}$ $\mathbf{b}$

Oczywiście nie jest trudno znaleźć jednoznaczne rozwiązanie oparte na dowolnej liczbie metod, ale jeśli istnieje coś, co okazało się optymalne pod względem liczby FLOPS, byłoby to idealne (w końcu cały problem prawdopodobnie zmieści się w rejestrach FP!).

(Tak, ta rutyna jest często nazywana . Już pozbyłem się nisko wiszących owoców, a to jest następne na mojej liście profilowania ...)

linear-solver reference-request

— Damien
źródło

Czy każde używane tylko raz, czy też istnieje wiele układów liniowych o tej samej matrycy? To zmieniłoby koszty.

A

$A$

— Federico Poloni

W tym przypadku A jest używane tylko raz.

— Damien

14

Nie możesz pobić wyraźnej formuły. Możesz zapisać formuły rozwiązania na kartce papieru. Pozwól, aby kompilator zoptymalizował dla Ciebie. Każda inna metoda będzie prawie nieuchronnie zawierać instrukcje lub pętle (np. W przypadku metod iteracyjnych), które spowodują, że Twój kod będzie wolniejszy niż jakikolwiek kod linii prostej. $x=A^{-1}b$ iffor

— Wolfgang Bangerth
źródło

9

Ponieważ matryca jest tak blisko tożsamości, następujące serie Neumanna zbiegną się bardzo szybko:

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

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

W zależności od wymaganej dokładności może być nawet wystarczająco dobra, aby obciąć po 2 warunkach:

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

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

Może to być nieco szybsze niż bezpośrednia formuła (jak sugeruje odpowiedź Wolfganga Bangertha), choć z dużo mniejszą dokładnością.

Możesz uzyskać większą dokładność dzięki 3 terminom:

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

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

ale jeśli wypiszesz formułę wpis po wpisie dla $(3I - 3A + A^2)b$ , patrzysz na porównywalną liczbę operacji zmiennoprzecinkowych jak bezpośrednia odwrotna formuła macierzy 3x3 (nie musisz jednak dokonywać podziału, co trochę pomaga).

— Nick Alger
źródło

Czy dywizje są nadal droższe niż inne klapy? Myślałem, że to relikt przeszłości.

— Federico Poloni,

Dywizje nie radzą sobie dobrze z niektórymi architekturami (ARM jest współczesnym przykładem)

— Damien

@FedericoPoloni Dzięki Cuda możesz zobaczyć przepływność instrukcji tutaj , jest sześciokrotnie wyższa dla mnożenia / dodawania niż dla podziałów.

— Kirill,

@Damien i Kirill widzę, dzięki za wskazówki.

— Federico Poloni,

5

Liczba FLOPS na podstawie powyższych sugestii:

LU, bez obrotu:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; lub
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Eliminacja Gaussa z substytucją wsteczną, bez obrotu:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; lub
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Reguła Cramera poprzez rozszerzenie kofaktora
- Mul = 24, Div = 3, Add / Sub = 15, Total = 42; lub
- Mul = 27, Div = 1, Add / Sub = 15, Total = 43
Jawne odwrotne, a następnie pomnóż:
- Mul = 30, Div = 3, Add / Sub = 17, Total = 50; lub
- Mul = 33, Div = 1, Add / Sub = 17, Total = 51

Sprawdzanie poprawności koncepcji MATLAB:

Reguła Cramera za pomocą rozszerzenia kofaktora :

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 (bez obrotu) i podstawienie wsteczne:

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

Jawne odwrotne, a następnie pomnożone:

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

Eliminacja Gaussa:

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

Uwaga: dodaj własne metody i liczenie do tego postu.

— Damien
źródło

Czy obliczyłeś czas potrzebny do rozwiązania za pomocą tych dwóch metod?

— nicoguaro

Nie. Powyższy kod w ogóle nie zostanie wykonany szybko. Chodziło o to, aby uzyskać wyraźną liczbę FLOPS i podać kod do recenzji na wypadek, gdyby coś przeoczyłem,

— Damien,

W LU 5 działów można przekształcić w 5 MUL kosztem 2 dodatkowych operacji wzajemnych (tj. 1 / U11 i 1 / U22). Będzie to zależne od łuku, czy można tam osiągnąć zysk.

— Damien

2

Zakładając, że nie pomyliłem się, przybliżam

A^{- 1} b

$A^{-1}b$ przez

2 b - A b

$2b - Ab$ będzie wymagało 12 mnożenia, 9 dodawania / odejmowania i żadnych podziałów. Przybliżenie

A^{- 1} b

$A^{-1}b$ przez

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

$3(b - Ab) + A^2b$ będzie wymagało 21 mnożenia i 18 dodawania / odejmowania. Obliczenie

A^{- 1} b

$A^{-1}b$ za pomocą tej wyraźnej formuły wygląda 33 mnożenia, 17 dodawania / odejmowania i 1 podział. Jak powiedziałem, moje numery mogą być wyłączone, więc możesz chcieć dokładnie sprawdzić.

— Geoff Oxberry

@GeoffOxberry, przyjrzę się temu i złożę raport.

— Damien

4

Prawdopodobnie Reguła Cramera. Jeśli możesz uniknąć obrotu, być może faktoryzacja LU; to macierz 3x3, więc ręczne rozwijanie pętli byłoby łatwe. Wszystko inne prawdopodobnie będzie wymagało rozgałęzienia i wątpię, aby metoda podprzestrzeni Kryłowa zbiegała się wystarczająco często w 1 lub 2 iteracjach, aby była tego warta.

— Geoff Oxberry
źródło