How we made a better backslash in MATLAB R2024a

Posted by Mike Croucher, April 2, 2024

35 views (last 30 days) | 0 Likes | 4 comments

Few things say 'MATLAB' better than the backslash. So good, we put it on a T-shirt.

As most people will be aware, the backslash operator solves a system of linear equations. That is, it finds the x that solves A*x = b. For example

A = [2 1 1;

    -1 1 -1;

    1 2 3];

b = [2;3;-10];

x = A\b

So simple and yet there is a lot of numerical linear algebra behind backslash. You see, exactly how one writes an algorithm to solve a linear system of equations depends very much on the structure and type of the matrix A. Is A Hermitian? If so, is its diagonal all positive or all negative and so on? Then you have to ask is A single, double or complex? Is the matrix full or sparse? The answers to these questions determine which algorithm is best to use. Get it right and you get the correct answer very quickly, making maximal use of both your computer hardware and decades of numerical linear algebra research. 

In the older MATLAB R2023b, the logic that backslash uses for full (i.e. not sparse) matrices is as follows

In the latest version of MATLAB, R2024a, we rearranged this flow-chart a little and added a test for tridiagonal matrices. This isn't a particularly expensive test since we are already some of the way there by doing the upper Hessenberg check.

As such, the backslash flow chart for full matrices is as follows (These flowcharts are from the documentation for mldivide by the way)

The practical upshot is that backslash is now significantly faster when the input matrix is dense and tridiagonal. In general, you are better off using sparse matrices when dealing with tridiagonal matrices (see later in this post) but sometimes you happen to have a dense one (e.g. from the output of ldl) and we wanted to ensure that backslash would deal with it efficiently. 

Tridiagonal matrices are zero everywhere apart from the main diagonal and the two diagonals either side of it (the so-called superdiagonal and subdiagonal).  Here's an example generated using MATLAB's gallery function

full(gallery("tridiag",8))

ans = 8×8
  -1     0     0     0     0     0     0
    -1     2    -1     0     0     0     0     0
  -1     2    -1     0     0     0     0
   0    -1     2    -1     0     0     0
   0     0    -1     2    -1     0     0
   0     0     0    -1     2    -1     0
   0     0     0     0    -1     2    -1
   0     0     0     0     0    -1     2

Tridiagonal matrices are particularly important because they appear in many areas of science and engineering including quantum mechanics, partial differential equations and splines to name a few.

To see how much faster we've made things, let's construct a problem involving a random tridiagonal matrix

n = 1e4;

mainDiag =  rand(n, 1);

upperDiag = rand(n-1, 1);

lowerDiag =  rand(n-1, 1);

triDiag = diag(mainDiag) + diag(upperDiag, 1) + diag(lowerDiag, -1);

Just to make sure that I haven't made a mistake while generating triDiag, I turn to the isbanded function. The lower and upper bandwidths of a tridiagonal matrix are both 1 so I can test as follows

isbanded(triDiag,1,1)  % If this returns true, D is tridiagonal

ans = logical
   1

OK, so now I'm convinced that D is tridiagonal and that I can use it to time this update to backslash

b = rand(n,1); % Construct the RHS of the equation 

f = @() triDiag\b;

runtime = timeit(f)

runtime = 0.0716

 On the same machine, R2023b did this in 0.5637 seconds so we have

R2024a: 0.0716 seconds
R2023b: 0.5637 seconds

Almost 7.9x speed-up! This is a little more than the 6.5x speed-up discussed in the release notes. The author of that section of the release notes used a Windows 11, AMD EPYC 74F3 24-Core Processor @ 3.19 GHz test system whereas my system is rather more modest.

cpuinfo

ans = struct with fields:
        CPUName: '11th Gen Intel(R) Core(TM) i7-11700 @ 2.50GHz'
          Clock: '2501 MHz'
          Cache: 4194304
    TotalMemory: 3.4010e+10
        NumCPUs: 1
     TotalCores: 8
         OSType: 'Windows'
      OSVersion: 'Microsoft Windows 11 Pro'
       Hostname: 'DESKTOP-N3C7C6Q'

Clearly, how much faster the new version of backslash is depends on your system but whatever machine you use, it will definitely be faster in R2024a.

Going sparse for maximum speed

In linear algebra, the path to speed often relies on taking advantage of matrix structure and there is one aspect of the structure of tridiagonal matrices that is important -- the fact they are very sparse. We haven't sped this up in 24a but it's worth pointing out that if you can use sparse matrices for these problems then you really should. It's so much faster, even taking into account the speed-ups shown above. The memory benefits are huge too! 

Let's convert triDiag to a sparse matrix to see 

spTriDiag = sparse(triDiag);

f = @() spTriDiag\b;

runtime = timeit(f)

runtime = 7.2862e-04

That's almost 100x faster than the full-matrix version in R2024a and almost 800x faster than the full-matrix version in R2023b.

More on backslash

Backslash - A historical post where Cleve tells us how backslash came to represent the solution of linear systems in MATLAB
What is A\A? - The answer: A\A is always I, except when it isn't.