Быстрое явное решение для

9

Я ищу быстрое (смею сказать, оптимальное?) Явное решение линейной вещественной задачи 3x3, $\mathbf{A}\mathbf{x} = \mathbf{b}$ , $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$ ,

матрица $\mathbf{A}$ является общим, но близким к единичной матрице с номером условия, близким к 1. Потому что $\mathbf{b}$ на самом деле измерения датчика с точностью около 5 цифр, я не против потерять несколько цифр из-за числовых проблем.

Конечно, нетрудно придумать явное решение, основанное на любом количестве методов, но если есть что-то, что оказалось оптимальным с точки зрения количества FLOPS, это было бы идеально (в конце концов, вся проблема скорее всего будет вписываться в регистры ФП!).

(Да, эту процедуру часто называют . Я уже избавился от низко висящих фруктов, и это следующее в моем списке профилирования ...)

linear-solver reference-request Damien
источник

Каждый

A

$A$ используется только один раз, или существует несколько линейных систем с одинаковой матрицей? Это изменило бы расходы.

Федерико Полони

В этом случае A используется только один раз.

Дэмиен

14

Вы не можете превзойти явную формулу. Вы можете записать формулы для решения $x=A^{-1}b$ на листе бумаги. Позвольте компилятору оптимизировать вещи для вас. Любой другой метод почти неизбежно будет иметь ifоператоры или forциклы (например, для итеративных методов), которые сделают ваш код медленнее, чем любой прямой код.

Вольфганг Бангерт
источник

9

Поскольку матрица настолько близка к единице, следующие ряды Неймана будут очень быстро сходиться:

A^{- 1} знак равно Σ_{К знак равно 0}^{\infty} (я - A)^{К}

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

В зависимости от требуемой точности это может быть даже достаточно хорошо для усечения после 2 сроков:

A^{- 1} \approx я + (я - A) знак равно 2 я - A,

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

Это может быть немного быстрее, чем прямая формула (как предложено в ответе Вольфганга Бангерта), но с гораздо меньшей точностью.

Вы можете получить больше точности с 3 терминами:

A^{- 1} \approx я + (я - A) + (я - A)^{2} знак равно 3 я - 3 A + A^{2}

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

но если вы напишите формулу для каждой записи $(3I - 3A + A^2)b$ вы смотрите на сопоставимое количество операций с плавающей запятой в качестве прямой формулы обратной матрицы 3х3 (хотя вам не нужно делать деление, что немного помогает).

Ник Алджер
источник

Дивизии все еще дороже, чем другие флопы? Я думал, что это пережиток прошлого.

Федерико Полони

Подразделения плохо прорабатывают некоторые архитектуры (ARM - современный пример)

Дамьен,

@FedericoPoloni С Cuda вы можете видеть пропускную способность инструкций здесь , для умножения / сложения она в шесть раз выше, чем для делений.

Кирилл

@ Дэмиен и Кирилл, я вижу, спасибо за указатели.

Федерико Полони

5

FLOPS рассчитывать на основе предложений выше:

LU, без поворота:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; или
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Гауссовское исключение с обратной заменой, без поворота:
- Mul = 11, Div / Recip = 6, Add / Sub = 11, Total = 28; или
- Mul = 16, Div / Recip = 3, Add / Sub = 11, Total = 30
Правило Крамера через расширение кофактора
- Mul = 24, Div = 3, Add / Sub = 15, Total = 42; или
- Mul = 27, Div = 1, Add / Sub = 15, Total = 43
Явное обратное, затем умножить:
- Mul = 30, Div = 3, Add / Sub = 17, Total = 50; или
- Mul = 33, Div = 1, Add / Sub = 17, Total = 51

Подтверждение концепций MATLAB:

Правило Крамера через расширение кофактора :

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 (без поворота) и обратная замена:

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

Явное обратное, затем умножение:

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

Гауссово исключение:

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

Примечание: пожалуйста, не стесняйтесь добавлять свои собственные методы и подсчет к этому сообщению.

Damien
источник

Вы рассчитали время, которое потребовалось для решения двумя методами?

Никогуаро

Нет. Код выше не будет выполняться быстро вообще. Смысл этого состоял в том, чтобы получить явный счетчик FLOPS и предоставить код для проверки на случай, если я что-то пропустил,

Дамьен

В LU 5 делений могут быть преобразованы в 5 MUL за счет 2 дополнительных взаимных операций (т.е. 1 / U11 и 1 / U22). Это будет зависеть от того, есть ли здесь какая-то прибыль.

Дэмиен

2

Предполагая, что я не ошибся, приближение

A^{- 1} b

$A^{-1}b$ по

2 b - A b

$2b - Ab$ собирается потребовать 12 умножений, 9 сложений / вычитаний и без делений. Аппроксимация

A^{- 1} b

$A^{-1}b$ по

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

$3(b - Ab) + A^2b$ собирается потребовать 21 умножения и 18 сложения / вычитания. расчета

A^{- 1} b

$A^{-1}b$ через эту явную формулу выглядит как 33 умножения, 17 сложений / вычитаний и 1 деление. Как я уже сказал, мои номера могут быть выключены, так что вы можете проверить еще раз.

Джефф Оксберри

@ GeoffOxberry, я посмотрю на это и сообщу.

Дэмиен

4

Вероятно, правило Крамера. Если вы можете избежать поворота, возможно, LU факторизация; это матрица 3х3, поэтому развернуть петли вручную будет легко. Все остальное, вероятно, будет включать ветвление, и я сомневаюсь, что метод подпространства Крылова достаточно часто сходится за 1 или 2 итерации, чтобы это стоило того.

Джефф Оксберри
источник

Быстрое явное решение для

Ответы: