# Fast explicit solution for $\mathbf{A}\mathbf{x} = \mathbf{b}$, $\mathbf{b} \in \mathbf{R}^3$, low condition number

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

• Is each $A$ used only once, or are there multiple linear systems with the same matrix? This would change the costs. Oct 11, 2014 at 15:06
• In this instance, A is used only once. Oct 11, 2014 at 22:54

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.

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

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

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


• 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. Oct 10, 2014 at 23:19