Trefethen & Bau & MATLAB & Julia: Lectures 20, 21, 23: Solving square systems
Three in one this time: Lecture 20, which is on Gaussian elimination / LU factorization, Lecture 21, on row pivoting, and Lecture 23, on Cholesky factorization. I mainly skipped over Lecture 22, about the curious case of the stability of pivoted LU, but the main example is dropped into the end of my coverage of pivoting.
The Julia surprises are, not surprisingly, coming less frequently. In Lecture 20 I had some fun with rational representations. I like using MATLAB’s format rat when presenting Gaussian elimination, as it allows me to recall the way the process looks when learned by hand. It’s a fun trick, but of course the underlying values are still all double precision, and the rational approximations to them are found ex post facto. By contrast, Julia offers true rational numbers, constructed and shown using the // operation.
Compare the MATLAB
format rat
I = eye(4);
L21 = I + (-5/17)*I(:,2)*I(:,1)';
L31 = I + (-9/17)*I(:,3)*I(:,1)';
L41 = I + (-4/17)*I(:,4)*I(:,1)';to the Julia
I = eye(Rational,4);
L21 = copy(I); L21[2,1] = -5//17;
L31 = copy(I); L31[3,1] = -9//17;
L41 = copy(I); L41[4,1] = -4//17;The MATLAB code requires only the format call, because it’s only the display of results that is affected. The Julia code is doing something deeper and needs more changes as a result.
Julia could use something like a format command. I almost always find MATLAB’s terminal output more readable, or at least easier to manipulate into a good form. Here’s one example using the rational output. First, MATLAB:
17 2 3 13 93
0 194/17 155/17 71/17 963/17
0 101/17 92/17 87/17 574/17
0 230/17 243/17 -18/17 1158/17And the Julia:
4x5 Array{Rational{T<:Integer},2}:
17//1 2//1 3//1 13//1 93//2
0//1 194//17 155//17 71//17 963//34
0//1 101//17 92//17 87//17 287//17
0//1 230//17 243//17 -18//17 579//17I almost never need that header line that Julia gives. The numbers are already showing themselves to be Rational, and the shape of the array is self-evident. (Though I now see that MATLAB 2016b is adding such headers to non-float output.) The zero structure also jumps out more clearly in the MATLAB case, though it’s profligate with whitespace.
Another comparison, of MATLAB (using the default format):
1 0 0 0 0 1
0 1 0 0 0 2
0 0 1 0 0 4
0 0 0 1 0 8
0 0 0 0 1 16
0 0 0 0 0 32versus Julia:
6×6 Array{Float64,2}:
1.0 0.0 0.0 0.0 0.0 1.0
0.0 1.0 0.0 0.0 0.0 2.0
0.0 0.0 1.0 0.0 0.0 4.0
0.0 0.0 0.0 1.0 0.0 8.0
0.0 0.0 0.0 0.0 1.0 16.0
0.0 0.0 0.0 0.0 0.0 32.0There’s nothing wrong per se about Julia’s. But which version would you write down, or expect to see in print? One last case, of a matrix that is supposed to be triangular but for a little roundoff. First, Julia:
4×4 Array{Float64,2}:
17.0 2.0 3.0 13.0
0.0 13.5294 14.2941 -1.05882
0.0 0.0 -2.93913 5.06957
5.55112e-16 0.0 -4.44089e-16 4.09024And MATLAB, using the default format:
17.0000 2.0000 3.0000 13.0000
0 13.5294 14.2941 -1.0588
0 0 -2.9391 5.0696
0.0000 0 -0.0000 4.0902Julia has chosen to align on the decimal point. It’s also suppressing trailing zeros, except for the first, giving an odd and false impression of values that have a precise number of significant digits. MATLAB’s choice of right alignment is visually superior, and only exact zero gets a special display. True, you might want that exponential notation for the tiny values; you can get it by changing the format.
>> format short e
>> U
U =
1.7000e+01 2.0000e+00 3.0000e+00 1.3000e+01
0 1.3529e+01 1.4294e+01 -1.0588e+00
0 0 -2.9391e+00 5.0696e+00
5.5511e-16 0 -4.4409e-16 4.0902e+00
>> format short g
>> U
U =
17 2 3 13
0 13.529 14.294 -1.0588
0 0 -2.9391 5.0696
5.5511e-16 0 -4.4409e-16 4.0902It’s nice to have options.