9
$\begingroup$

I am looking for a fast (dare I say optimal?) explicit solution the 3x3 linear real problem, $\mathbf{A}\mathbf{x} = \mathbf{b}$, $\mathbf{A} \in \mathbf{R}^{3 \times 3}, \mathbf{b} \in \mathbf{R}^{3}$.

Matrix $\mathbf{A}$ is general, but close to the identity matrix with a condition number close to 1. Because $\mathbf{b}$ are actually sensor measurements with about 5 digits of precision, I do not mind losing several digits due to numerical issues.

Of course, it is not hard to come up with an explicit solution based on any number of methods, but if there is something that has been shown to be optimal in terms of FLOPS count, that would be ideal (after all, the whole problem will likely fit in the FP registers!).

(Yes, this routine is called often. I've already gotten rid of low-hanging fruit and this is next in my profiling list...)

$\endgroup$
  • $\begingroup$ Is each $A$ used only once, or are there multiple linear systems with the same matrix? This would change the costs. $\endgroup$ – Federico Poloni Oct 11 '14 at 15:06
  • $\begingroup$ In this instance, A is used only once. $\endgroup$ – Damien Oct 11 '14 at 22:54
14
$\begingroup$

You can't beat an explicit formula. You can write down the formulas for the solution $x=A^{-1}b$ on a piece of paper. Let the compiler optimize things for you. Any other method will almost inevitably have if statements or for loops (e.g., for iterative methods) that will make your code slower than any straight line code.

$\endgroup$
9
$\begingroup$

Since the matrix is so close to the identity, the following Neumann series will converge very rapidly:

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

Depending on the accuracy required it might even be good enough to truncate after 2 terms:

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

This might be slightly faster than a direct formula (as suggested in Wolfgang Bangerth's answer), though with much less accuracy.


You could get more accuracy with 3 terms: $$A^{-1} \approx I + (I - A) + (I-A)^2 = 3I - 3A + A^2$$

but if you write out the entry-by-entry formula for $(3I - 3A + A^2)b$, you are looking at a comparable amount of floating point operations as the direct 3x3 matrix inverse formula (you don't have to do a division though, which helps a little).

$\endgroup$
  • $\begingroup$ Are divisions still more expensive than the other flops? I thought it was a relic of the past. $\endgroup$ – Federico Poloni Oct 10 '14 at 10:26
  • $\begingroup$ Divisions do not pipeline well one some architectures (ARM is the contemporary example) $\endgroup$ – Damien Oct 10 '14 at 11:39
  • $\begingroup$ @FedericoPoloni With Cuda, you can see instruction throughput here, it's six times higher for multiplications/additions than for divisions. $\endgroup$ – Kirill Oct 10 '14 at 13:13
  • $\begingroup$ @Damien and Kirill I see, thanks for the pointers. $\endgroup$ – Federico Poloni Oct 10 '14 at 13:25
5
$\begingroup$

FLOPS count based on the suggestions above:

  • LU, no pivoting:

    • Mul = 11, Div/Recip = 6, Add/Sub = 11, Total = 28; or
    • Mul = 16, Div/Recip = 3, Add/Sub = 11, Total = 30
  • Gaussian Elimination with back-substitution, no pivoting:

    • Mul = 11, Div/Recip = 6, Add/Sub = 11, Total = 28; or
    • Mul = 16, Div/Recip = 3, Add/Sub = 11, Total = 30
  • Cramer's rule via cofactor expansion

    • Mul = 24, Div = 3, Add/Sub = 15, Total = 42; or
    • Mul = 27, Div = 1, Add/Sub = 15, Total = 43
  • Explicit Inverse then multiply:

    • Mul = 30, Div = 3, Add/Sub = 17, Total = 50; or
    • Mul = 33, Div = 1, Add/Sub = 17, Total = 51

MATLAB proof-of-concepts:

Cramer's Rule via Cofactor Expansion:

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 (no pivoting) and back-substitution:

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

Explicit Inverse then Multiply:

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

Gaussian Elimination:

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

Note: Please feel free to add your own methods and counts to this post.

$\endgroup$
  • $\begingroup$ Did you computed the times it took to solve with the two methods? $\endgroup$ – nicoguaro Oct 10 '14 at 22:35
  • $\begingroup$ No. The code above will not execute quickly at all. The point of it was to get an explicit FLOPS count and supply the code for review in case I missed something, $\endgroup$ – Damien Oct 10 '14 at 22:39
  • $\begingroup$ In LU, 5 divisions can be converted to 5 MULs at the expense of 2 extra reciprocal operations (i.e. 1/U11 and 1/U22). That will be arch-specific as to whether there is a gain to be made there. $\endgroup$ – Damien Oct 10 '14 at 22:41
  • 2
    $\begingroup$ Assuming I haven't mis-counted, approximating $A^{-1}b$ by $2b - Ab$ is going to require 12 multiplications, 9 addition/subtractions, and no divisions. Approximating $A^{-1}b$ by $3(b - Ab) + A^2b$ is going to require 21 multiplications and 18 addition/subtractions. Calculating $A^{-1}b$ via this explicit formula looks to be 33 multiplications, 17 additions/subtractions, and 1 division. Like I said, my numbers might be off, so you might want to double check. $\endgroup$ – Geoff Oxberry Oct 10 '14 at 23:19
  • $\begingroup$ @GeoffOxberry, I will look into it and report. $\endgroup$ – Damien Oct 10 '14 at 23:55
4
$\begingroup$

Probably Cramer's Rule. If you can avoid pivoting, maybe LU factorization; it's a 3x3 matrix, so unrolling the loops manually would be easy. Anything else will probably involve branching, and I doubt that a Krylov subspace method would converge often enough in 1 or 2 iterates for it to be worth it.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.