MultiSolve 3x3

Vectorized solve of multiple 3x3 linear systems

Description

xx = multisolve3x3(AA,b)
    
Highly optimized solution of multiple 3x3 linear systems.
    
The m 3-by-3-coefficient matrices can be given in 2 forms:
(1) The argument matrix AA is 3m x 3 (i.e., the k-th linear system occupies rows 3k-2:3k).
(2) The argument AA is an 3x3xm array. (the k-th linear system is AA(:,:,k) ).
    
Similarly, the result vectors b can be given in 2 forms (independent of AA):
(1) a 3m x 1 vector, where the k-th result is in rows 3k-2:3k,
(2) a 3xm matrix, whose k-th column contains the k-th result.

Either way, the solution xx is given in size identical to that of b.

Example:

m = 10 ; % system num
AA = rand(3,3,m);
bb = rand(3,m);

xx = multisolve3x3(AA,bb);

%test results:
for jj=1:m
    max(abs( AA(:,:,jj) * xx(:,jj) - bb(:,jj) ))
end

• MATLAB release: MATLAB 7.2 (R2006a)

[John D'Errico]

Let me preface these comments with the statement that I like the idea of multisolve3x3. It has useful help. In fact, I've often enough used a similar technique myself to solve mutiple sets of small linear systems. Multisolve3x3 should also be fast.

What I worry about is the stability of multisolve3x3. Because it uses an adjoint matrix approach, dividing by 3x3 determinants to solve the blocks, it may not be stable when one or more of the systems are nearly singular. For example...

% Choose a random matrix, with moderately large
% condition number.
A = rand(3,3)^8;
cond(A)
ans =
   2.9374e+07

b0 = ones(3,1);
b = b0+10*eps*randn(3,1);

% Note that backslash amplifies the noise in b, by
% roughly the condition number of A
(A\b0) - (A\b)
ans =
   6.9122e-11
  -2.4738e-10
   1.3824e-10

% this is also true of multisolve
multisolve3x3(A,b0) - multisolve3x3(A,b)
ans =
   6.5484e-11
  -2.1828e-10
   1.1642e-10

% But there appears to be a problem with multisolve3x3
% when compared to backslash.
multisolve3x3(A,b0)-(A\b0)
ans =
    -0.058747
      0.19833
     -0.10987

% Which solution was more correct?
x0 = A\b0;
xm = multisolve3x3(A,b0);

A*x0 - b0
ans =
  -1.3601e-10
    1.047e-10
  -1.1154e-10

A*xm - b0
ans =
   5.3125e-06
   3.6284e-06
     3.87e-06

Note that multisolve3x3 was 4 orders of magnitude
worse in its residual error.

The point of my example is that while multisolve3x3 is likely to work acceptably for well conditioned systems, you should worry what happens on poorly conditioned problems. I wanted to rate this as a 3.5. I'd have preferred to see the author use a sparse block diagonal backslash to solve the problem instead of the adjoint. My rating then would be a 5.