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...)

Improve this question
$\endgroup$
2
  • $\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, 2014 at 15:06
  • $\begingroup$ In this instance, A is used only once. $\endgroup$
    – Damien
    Oct 11, 2014 at 22:54

4 Answers 4

Reset to default
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.

Improve this answer
$\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).

Improve this answer
$\endgroup$
4
  • $\begingroup$ Are divisions still more expensive than the other flops? I thought it was a relic of the past. $\endgroup$
    – Federico Poloni
    Oct 10, 2014 at 10:26
  • $\begingroup$ Divisions do not pipeline well one some architectures (ARM is the contemporary example) $\endgroup$
    – Damien
    Oct 10, 2014 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, 2014 at 13:13
  • $\begingroup$ @Damien and Kirill I see, thanks for the pointers. $\endgroup$
    – Federico Poloni
    Oct 10, 2014 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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Did you computed the times it took to solve with the two methods? $\endgroup$
    – nicoguaro ♦
    Oct 10, 2014 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, 2014 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, 2014 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, 2014 at 23:19
  • $\begingroup$ @GeoffOxberry, I will look into it and report. $\endgroup$
    – Damien
    Oct 10, 2014 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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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