{"id":2192,"date":"2024-04-02T06:54:20","date_gmt":"2024-04-02T10:54:20","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=2192"},"modified":"2024-04-02T07:07:23","modified_gmt":"2024-04-02T11:07:23","slug":"how-we-made-a-better-backslash-in-matlab-r2024a","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2024\/04\/02\/how-we-made-a-better-backslash-in-matlab-r2024a\/","title":{"rendered":"How we made a better backslash in MATLAB R2024a"},"content":{"rendered":"
Few things say 'MATLAB' better than the backslash. So good, we put it on a T-shirt.<\/span><\/div>
\"\"<\/img><\/div>
As most people will be aware, the backslash operator solves a system of linear equations. That is, it finds the <\/span>x<\/span> that solves <\/span>A*x = b<\/span>. For example<\/span><\/div>
A = [2 1 1;<\/span><\/span><\/div><\/div>
-1 1 -1;<\/span><\/span><\/div><\/div>
1 2 3];<\/span><\/span><\/div><\/div>
 <\/div><\/div>
b = [2;3;-10];<\/span><\/span><\/div><\/div>
 <\/div><\/div>
x = A\\b<\/span><\/span><\/div>
x = <\/span>3\u00d71<\/span><\/div><\/div>
3\r\n 1\r\n -5\r\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
So simple and yet there is a <\/span>lot<\/span> of numerical linear algebra behind backslash. You see, exactly <\/span>how<\/span> one writes an algorithm to solve a linear system of equations depends very much on the structure and type of the matrix <\/span>A<\/span>. Is A Hermitian? If so, is its diagonal all positive or all negative and so on? Then you have to ask is A <\/span>single<\/span>, <\/span>double<\/span> or <\/span>complex? <\/span>Is the matrix <\/span>full<\/span><\/a> or <\/span>sparse<\/span><\/a>? 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. <\/span><\/div>
In the older MATLAB R2023b, the logic that backslash uses for full (i.e. not sparse) matrices is as follows<\/span><\/div>
\"\"<\/img><\/div>
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.<\/span><\/div>
As such, the backslash flow chart for full matrices is as follows (These flowcharts are from the documentation for <\/span>mldivide<\/span><\/a> by the way<\/span>)<\/span><\/div>
\"\"<\/img><\/div>
The practical upshot is that backslash is now significantly faster when the input matrix is dense and <\/span>tridiagonal<\/span><\/a>. In general, you are better off using sparse matrices when dealing with tridiagonal matrices<\/span> (see later in this post) but sometimes you happen to have a dense one (e.g. from the output of <\/span>ldl)<\/span> and we wanted to ensure that backslash would deal with it efficiently. <\/span><\/div>
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 <\/span>gallery<\/span><\/a> function<\/span><\/div>
full(gallery(<\/span>\"tridiag\"<\/span>,8))<\/span><\/span><\/div>
ans = <\/span>8\u00d78<\/span><\/div><\/div>
2 -1 0 0 0 0 0 0\r\n -1 2 -1 0 0 0 0 0\r\n 0 -1 2 -1 0 0 0 0\r\n 0 0 -1 2 -1 0 0 0\r\n 0 0 0 -1 2 -1 0 0\r\n 0 0 0 0 -1 2 -1 0\r\n 0 0 0 0 0 -1 2 -1\r\n 0 0 0 0 0 0 -1 2\r\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
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.<\/span><\/div>
To see how much faster we've made things, let's construct a problem involving a random tridiagonal matrix<\/span><\/div>
n = 1e4;<\/span><\/span><\/div><\/div>
mainDiag = rand(n, 1);<\/span><\/span><\/div><\/div>
upperDiag = rand(n-1, 1);<\/span><\/span><\/div><\/div>
lowerDiag = rand(n-1, 1);<\/span><\/span><\/div><\/div>
triDiag = diag(mainDiag) + diag(upperDiag, 1) + diag(lowerDiag, -1);<\/span><\/span><\/div><\/div><\/div>
Just to make sure that I haven't made a mistake while generating <\/span>triDiag<\/span>, I turn to the <\/span>isbanded<\/span><\/a> function. The lower and upper bandwidths of a tridiagonal matrix are both 1 so I can test as follows<\/span><\/div>
isbanded(triDiag,1,1) <\/span>% If this returns true, D is tridiagonal<\/span><\/span><\/div>
ans = logical<\/span><\/span><\/div>
1\r\n<\/div><\/div><\/div><\/div><\/div><\/div>
OK, so now I'm convinced that D is tridiagonal and that I can use it to time this update to backslash<\/span><\/div>
b = rand(n,1); <\/span>% Construct the RHS of the equation <\/span><\/span><\/div><\/div>
 <\/div><\/div>
f = @() triDiag\\b;<\/span><\/span><\/div><\/div>
runtime = timeit(f)<\/span><\/span><\/div>
runtime = 0.0716<\/div><\/div><\/div><\/div>
On the same machine, R2023b did this in 0.5637 seconds so we have<\/span><\/div>