The Strange Case of the Determinant of a Matrix of 1s and -1s

By Nick Higham and Alan Edelman (MIT)

In a 2005 talk the second author noted that the MATLAB det function returns an odd integer for a certain 27-by-27 matrix composed of 1s and -1s:

>> A = edelman; % Set up the matrix.
>> format long g, format compact, det(A)
ans =
           839466457497601

However, the determinant is, from its definition, a sum of an even number (27 factorial) of odd numbers, so is even. Indeed the correct determinant is 839466457497600.

At first sight, this example is rather troubling, since while MATLAB returns an integer, as expected, it is out by 1. The determinant is computed as the product of the diagonal entries of the U factor in the LU factorization with partial pivoting of A, and these entries are not all integers. Standard rounding error analysis shows that the relative error from forming that product is bounded by nu/(1-nu), with n=27, where u \approx 1.1 \times 10^{-16} is the unit roundoff, and this is comfortably larger than the actual relative error (which also includes the errors in computing U) of 6 \times 10^{-16}. Therefore the computed determinant is well within the bounds of roundoff, and if the exact result had not been an integer the incorrect last decimal digit would hardly merit discussion.

However, this matrix has more up its sleeve. Let us compute the determinant using a different implementation of Gaussian elimination with partial pivoting, namely the function gep from the Matrix Computation Toolbox:

>> [Lp,Up,Pp] = gep(A,'p'); det(Pp)*det(Up)
ans =
           839466457497600

Now we get the correct answer! To see what is happening, we can directly form the products of the diagonal elements of the U factors:

>> [L,U,P] = lu(A); 
>> d = diag(U); dp = diag(Up); 
>> rel_diff_U_diags = norm((dp - d)./d,inf)
rel_diff_U_diags =
      7.37206353875273e-16
>> [prod(d), prod(dp)]
ans =
          -839466457497601          -839466457497600
>> [prod(d(end:-1:1)), prod(dp(end:-1:1))]
ans =
          -839466457497600          -839466457497600

We see that even though the diagonals of the two U factors differ by a small multiple of the unit roundoff, the computed products differ in the last decimal digit. If the product of the diagonal elements of U is accumulated in the reverse order then the exact answer is obtained in both cases. Once again, while this behaviour might seem surprising, it is within the error bounds of a rounding error analysis.

The moral of this example is that we should not be misled by the integer nature of a result; in floating-point arithmetic it is relative error that should be judged.

Finally, we note that numerical evaluation of the determinant offers other types of interesting behaviour. Consider the Frank matrix: a matrix of integers that has determinant 1. What goes wrong here in the step from dimension 24 to 25?

>> A = gallery('frank',24); det(A)
ans =
         0.999999999999996
>> A = gallery('frank',25); det(A)
ans =
          143507521.082525

The Edelman matrix is available in the MATLAB function available in this gist, which is embedded below. A Julia notebook exploring the Edelman matrix is available here.

function A = edelman
%EDELMAN Alan Edelman's matrix for which det is computed as the wrong integer.
% A = EDELMAN is a 27-by-27 matrix of 1s and -1s for which the
% MATLAB det function returns an odd integer, though the exact
% determinant is an even integer.
A = [%
1 1 1 1 -1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 -1 -1 -1
1 -1 -1 1 -1 1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1 1
-1 1 -1 -1 1 1 1 1 1 -1 1 -1 1 1 1 1 -1 -1 1 -1 -1 1 1 -1 1 1 -1
-1 -1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1 -1 -1 -1
-1 1 1 -1 -1 -1 -1 -1 1 -1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 1 -1 1 1 1 -1
1 1 1 1 1 1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 1 1 1 -1 1 -1
-1 -1 -1 -1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 1 -1 -1 1 -1
1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 -1 -1 1 1 -1
-1 1 -1 1 1 1 -1 -1 1 -1 1 1 1 1 -1 1 1 1 -1 1 1 -1 1 1 1 -1 -1
-1 1 -1 1 1 -1 -1 -1 -1 1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1
-1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 -1 -1 1 1 1 1 -1 1 1 1 -1 -1 1
1 1 -1 -1 1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 1 -1 1 1 -1 1 1 1 -1 -1
-1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1
1 1 -1 1 1 -1 1 1 -1 -1 -1 1 1 1 1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 -1
-1 1 -1 1 -1 1 1 1 -1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1 1 -1 1 1 -1 -1
-1 -1 1 -1 1 -1 1 -1 -1 1 -1 -1 -1 -1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1
1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 1 -1 -1 1 1 -1 1 -1 1 1 -1 1 1 -1
1 1 -1 1 1 -1 1 1 1 -1 1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 1
1 -1 -1 -1 1 -1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 1 1 -1
-1 -1 -1 1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1
1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 -1 -1 1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1
-1 -1 1 1 -1 -1 1 -1 -1 -1 1 -1 1 -1 -1 1 1 -1 1 -1 1 1 -1 1 -1 1 -1
-1 1 -1 1 1 1 1 1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 1 1 1 -1 -1 -1 1
-1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1 1 1 1 -1
1 1 1 1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 1 1 1 -1 1 -1 -1 1 1 -1 1 1
1 -1 1 -1 1 -1 1 1 -1 1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 1 1 -1 1 -1 -1
1 -1 1 1 -1 -1 1 1 1 -1 1 1 -1 1 1 1 1 1 -1 1 -1 1 -1 1 1 -1 1];
view raw edelman.m hosted with ❤ by GitHub

What’s New in MATLAB R2017b?

Following my earlier posts What’s New in MATLAB R2016b? and What’s New in MATLAB R2017a? I take a look here at the R2017b release of MATLAB. As before, this is not a comprehensive treatment (for which see the Release Notes), but rather a brief description of the changes in MATLAB (not the toolboxes) that are the most notable from my point of view.

Decomposition Object

MATLAB now has a way of avoiding unnecessarily repeating the factorization of a matrix. In the past, if we wanted to solve two linear systems Ax = b and Ay = c involving the same square, general nonsingular matrix A, writing

x = A\b;
y = A\c;

was wasteful because A would be LU factorized twice. The way to avoid this unnecessary work was to factorize A explicitly then re-use the factors:

[L,U] = lu(A);
x = U\(L\b);
y = U\(L\c);

This solution lacks elegance. It is also less than ideal from the point of view of code maintenance and re-use. If we want to adapt the same code to handle symmetric positive definite A we have to change all three lines, even though the mathematical concept is the same:

R = chol(A);
x = R\(R'\b);
y = R\(R'\c);

The new decomposition function creates an object that contains the factorization of interest and allows it to be re-used. We can now write

dA = decomposition(A);
x = dA\b;
y = dA\c;

MATLAB automatically chooses the factorization based on the properties of A (as the backslash operator has always done). So for a general square matrix it LU factorizes A, while for a symmetric positive definite matrix it takes the Cholesky factorization. The decomposition object knows to use the factors within it when it encounters a backslash. So this example is functionally equivalent to the first two.

The type of decomposition can be specified as a second input argument, for example:

dA = decomposition(A,'lu');
dA = decomposition(A,'chol');
dA = decomposition(A,'ldl');
dA = decomposition(A,'qr');

These usages make the intention explicit and save a little computation, as MATLAB does not have to determine the matrix type. Currently, 'svd' is not supported as a second input augment.

This is a very welcome addition to the matrix factorization capabilities in MATLAB. Tim Davis proposed this idea in his 2013 article Algorithm 930: FACTORIZE: An Object-Oriented Linear System Solver for MATLAB. In Tim’s approach, all relevant functions such as orth, rank, condest are overloaded to use a decomposition object of the appropriate type. MATLAB currently uses the decomposition object only in backslash, forward slash, negation, and conjugate transpose.

The notation dA used in the MathWorks documentation for the factorization object associated with A grates a little with me, since dA (or \Delta A) is standard notation for a perturbation of the matrix A and of course the derivative of a function. In my usage I will probably write something more descriptive, such as decompA or factorsA.

String Conversions

Functions written for versions of MATLAB prior to 2016b might assume that certain inputs are character vectors (the old way of storing strings) and might not work with the new string arrays introduced in MATLAB R2016b. The function convertStringsToChars can be used at the start of a function to convert strings to character vectors or cell arrays of character vectors:

function y = myfun(x,y,z)
[x,y,z] = convertStringToChars(x,y,z);

The pre-2016b function should then work as expected.

More functions now accept string arguments. For example, we can now write (with double quotes as opposed to the older single quote syntax for character vectors)

>> A = gallery("moler",3)
A =
     1    -1    -1
    -1     2     0
    -1     0     3

Tall Arrays

There is more support in R2017b for tall arrays (arrays that are too large to fit in memory). They can be indexed, with sorted indices in the first dimension; functions including median, plot, and polyfit now accept tall array arguments; and the random number generator used in tall array calculations can be independently controlled.

Word Clouds

The new wordcloud function produces word clouds. The following code runs the function on a file of words from The Princeton Companion to Applied Mathematics.

words = fileread('pcam_words.txt');
words = string(words);
words = splitlines(words);
words(strlength(words)<5) = [];
words = categorical(words);
figure
wordcloud(words)

I had prepared this file in order to generate a word cloud using the Wordle website. Here is the MATLAB version followed by the Wordle version: wordcloud_pcam.jpg wordle-edited2.jpg Wordle removes certain stop words (common but unimportant words) that wordcloud leaves in. Also, whereas Wordle produces a different layout each time it is called, different layouts can be obtained from wordcloud with the syntax wordcloud(words,'LayoutNum',n) with n a nonnegative integer.

Minima and Maxima of Discrete Data

New functions islocalmin and islocalmax find local minima and maxima within discrete data. Obviously targeted at data manipulation applications, such as analysis of time series, these functions offer a number of options to define what is meant by a minimum or maximum.

The code

A = randn(10); plot(A,'LineWidth',1.5); hold on
plot(islocalmax(A).*A,'r.','MarkerSize',25); 
hold off, axis tight

finds maxima down the columns of a random matrix and produces this plot: islocalmax_plot.jpg

Whereas min and max find the smallest and largest elements of an array, the new functions mink(A,k) and maxk(A,k) find the k smallest and largest elements of A, columnwise if A is a matrix.

How Fast is Quadruple Precision Arithmetic?

170724-2024-49-7501.jpg
The 24th IEEE Symposium on Computer Arithmetic included several talks on multiprecision arithmetic.

When I am testing an algorithm running in double precision arithmetic I often want to compare the computed solution with a reference solution: a solution that is fully accurate to double precision. To obtain one I solve the same problem in quadruple precision and round the result back to double precision. If the conditioning of the problem leads me to doubt that that this reference solution will be correct to double precision, I compute it at a precision even higher than quadruple. Whether it is feasible to compute a reference solution in this way depends on the size of the problem and the speed of quadruple precision arithmetic. This is just one scenario where quadruple precision arithmetic is needed; others are discussed in my earlier post The Rise of Mixed Precision Arithmetic.

Roughly speaking, a quadruple precision number x can be represented by a pair of double precision numbers (x_1,x_2), where x = x_1 + x_2 and x_1 and x_2 represent the leading and trailing significant digits of x, respectively. Then a product x y can be written (x_1 + x_2)(y_1 + y_2) = x_1y_1 + x_1y_2 + x_2y_1 + x_2y_2, which involves four double precision multiplications and three double precision additions. We can therefore expect quadruple precision multiplication to be at least seven times slower than double precision multiplication. Almost all available implementations of quadruple precision arithmetic are in software (an exception is provided by the IBM z13 mainframe system), so we should expect quadruple precision to be slower than the theoretical prediction would suggest. Bailey and Borwein (2015) state that quadruple precision implemented in software typically runs five to ten times slower than hardware double precision arithmetic.

Various aspects need to be taken into consideration when we compare double precision and quadruple precision arithmetics. First, the relative speed of the arithmetics may depend on the type of operations being performed (scalar versus vectorizable operations), memory traffic considerations, and to what extent the implementations exploit multicore processors. Second, are we to compare the same code running in different precisions or computational kernels coded specially for each precision? I will take the second approach, as it is more relevant to my usage of quadruple precision. However, this does mean that I am comparing algorithms and their implementations as well as the underlying arithmetic.

I have done some very simple comparisons in MATLAB using the VPA arithmetic of the Symbolic Math Toolbox and the MP arithmetic of the third party Multiprecision Computing Toolbox from Advanpix.

The Multiprecision Computing Toolbox supports IEEE-compliant quadruple precision arithmetic, which is invoked by using the function mp.Digits to set the precision to 34 decimal digits (in fact, this is the default). For the VPA arithmetic in the Symbolic Math Toolbox 34 decimal digit precision is specified with the command digits(34) (the default is 32 digits). According to the documentation, VPA arithmetic uses a few guard digits, that is, it computes with a few extra digits of precision. The number of guard digits cannot be set by the user. In the Multiprecision Computing Toolbox the number of guard digits can be set with the function mp.GuardDigits, and the default is to use no guard digits.

These experiments were performed on an Intel Broadwell-E Core i7-6800K CPU @ 3.40GHz, which has 6 cores. The software details are:

MATLAB Version: 9.2.0.538062 (R2017a)
Operating System: Microsoft Windows 10 Home Version 10.0
                  (Build 15063)
Advanpix Multiprecision Computing Toolbox   Version 4.4.1.12580   
Symbolic Math Toolbox                       Version 7.2  (R2017a)

Further details on the Multiprecision Computing Toolbox can be seen by typing mp.Info, which lists details of open source libraries that are used.

In these experiments I do not check that the computed results are correct. Such checks are done in more extensive tests available in the test scripts at https://www.advanpix.com/2015/02/12/version-3-8-0-release-notes/. I note that the timings vary from run to run, but it is the order of magnitudes of the ratios that are of interest and these are correctly reflected by the reported results.

Test 1: LU Factorization

This code compares the speed of the arithmetics for LU factorization of a 250-by-250 matrix.

%QUAD_PREC_TEST1  LU factorization.
rng(1), n = 250;
A = randn(n); A_mp = mp(A,34); A_vpa = vpa(A,34);
t = clock; [L,U] = lu(A);     t_dp = etime(clock, t)
t = clock; [L,U] = lu(A_mp);  t_mp = etime(clock, t)
t = clock; [L,U] = lu(A_vpa); t_vpa = etime(clock, t)
fprintf('mp/double: %8.2e, vpa/double: %8.2e, vpa/mp: %8.2e\n',....
         t_mp/t_dp, t_vpa/t_dp, t_vpa/t_mp)

The output is

t_dp =
   1.0000e-03
t_mp =
   9.8000e-02
t_vpa =
   2.4982e+01
mp/double: 9.80e+01, vpa/double: 2.50e+04, vpa/mp: 2.55e+02

Here, the MP quadruple precision is 98 times slower than double precision, but VPA is 25,000 times slower than double precision.

Test 2: Complete Eigensystem of Nonsymmetric Matrix

%QUAD_PREC_TEST2  Complete eigensystem of nonsymmetric matrix.
rng(1), n = 125;
A = randn(n); A_mp = mp(A,34); A_vpa = vpa(A,34);
t = clock; [V,D] = eig(A);     t_dp = etime(clock, t)
t = clock; [V,D] = eig(A_mp);  t_mp = etime(clock, t)
t = clock; [V,D] = eig(A_vpa); t_vpa = etime(clock, t)
fprintf('mp/double: %8.2e, vpa/double: %8.2e, vpa/mp: %8.2e\n',....
         t_mp/t_dp, t_vpa/t_dp, t_vpa/t_mp)

The output is

t_dp =
   1.8000e-02
t_mp =
   1.3430e+00
t_vpa =
   1.0839e+02
mp/double: 7.46e+01, vpa/double: 6.02e+03, vpa/mp: 8.07e+01

Here, MP is 75 times slower than double, and VPA is closer to the speed of MP.

Test 3: Complete Eigensystem of Symmetric Matrix

%QUAD_PREC_TEST3  Complete eigensystem of symmetric matrix.
rng(1), n = 200;
A = randn(n); A = A + A'; A_mp = mp(A,34); A_vpa = vpa(A,34);
t = clock; [V,D] = eig(A);     t_dp = etime(clock, t)
t = clock; [V,D] = eig(A_mp);  t_mp = etime(clock, t)
t = clock; [V,D] = eig(A_vpa); t_vpa = etime(clock, t)
fprintf('mp/double: %8.2e, vpa/double: %8.2e, vpa/mp: %8.2e\n',....
         t_mp/t_dp, t_vpa/t_dp, t_vpa/t_mp)

The output is

t_dp =
   1.1000e-02
t_mp =
   3.5800e-01
t_vpa =
   1.2246e+02
mp/double: 3.25e+01, vpa/double: 1.11e+04, vpa/mp: 3.42e+02

Note that there are at least three different algorithms that could be used here (the QR algorithm, divide and conquer, and multiple relatively robust representations), so the three eig invocations may be using different algorithms.

Test 4: Componentwise Exponentiation

%QUAD_PREC_TEST4  Componentwise exponentiation.
rng(1), n = 1000;
A = randn(n); A_mp = mp(A,34); A_vpa = vpa(A,34);
t = clock; X = exp(A);     t_dp = etime(clock, t)
t = clock; X - exp(A_mp);  t_mp = etime(clock, t)
t = clock; X - exp(A_vpa); t_vpa = etime(clock, t)
fprintf('mp/double: %8.2e, vpa/double: %8.2e, vpa/mp: %8.2e\n',....
         t_mp/t_dp, t_vpa/t_dp, t_vpa/t_mp)

The output is

t_dp =
   7.0000e-03
t_mp =
   8.5000e-02
t_vpa =
   3.4852e+01
mp/double: 1.21e+01, vpa/double: 4.98e+03, vpa/mp: 4.10e+02

Both MP and VPA come closer to double precision on this problem of computing the scalar exponential at each matrix element.

Summary

Not too much emphasis should be put on the precise timings, which vary with the value of the dimension n. The main conclusions are that 34-digit MP arithmetic is 1 to 2 orders of magnitude slower than double precision and 34-digit VPA arithmetic is 3 to 4 orders of magnitude slower than double precision.

It is worth noting that the documentation for the Multiprecision Computing Toolbox states that the toolbox is optimized for quadruple precision.

It is also interesting to note that the authors of the GNU MPFR library (for multiple-precision floating-point computations with correct rounding) are working on optimizing the library for double precision and quadruple precision computations, as explained in a talk given by Paul Zimmermann at the 24th IEEE Symposium on Computer Arithmetic in July 2017; see the slides and the corresponding paper.

[Updated September 6, 2017 to clarify mp.Info and April 6, 2018 to change “Multiplication” to “Multiprecision” in penultimate paragraph.]

Elements of MATLAB Style

140415-1443-12-4032-cropped.jpg

Style is an important aspect of writing, and also of programming. While MATLAB is a quick and easy language in which to program, style should not be neglected. Good style aids readability, which in turn makes it easier to debug and maintain code. It also fosters confidence in the code. Here are six style tips, in which the examples given are all inspired by real-life code that I have come across.

For more on MATLAB style see the book MATLAB Guide, and in particular section 16.1, “Elements of Coding Style”.

1. Omit Unnecessary Parentheses

% Bad style.
x = (A')\b;
y = [1];
D = (diag(x));

The parentheses around A' are unnecessary. Without them, the statement is legal and unambiguous. But of course the parentheses would be necessary in, for example, x = (A*B)'\c;

y is a scalar. There is no need to enter it as a 1-by-1 matrix by putting it in square brackets.

There is no need for parentheses around the diag function.

% Good style.
x = A'\b;
y = 1;
D = diag(x);

2. Use Spaces Consistently

% Bad style.
q=sqrt(2);
a = 1;
Z =zeros(3);
for i=1:10, q=q+1/q; end

Stick to a convention about spacing around around equals signs and plus and minus signs. I think spaces make the code more readable.

% Good style.
q = sqrt(2);
a = 1;
Z = zeros(3);
for i = 1:10, q = q + 1/q; end

3. Omit Unnecessary Semicolons and Commas

% Bad style.
figure;
for i = 2:n,
    x(i) = x(i-1)^2*y(i);
end;

Semicolons suppress output, but in this example there is no need for the them since no output is returned by the figure or end functions. The comma would only be necessary if the for loop was collapsed onto one line.

% Good style.
figure
for i = 2:n
    x(i) = x(i-1)^2*y(i);
end

4. Use the Appropriate Construct

The following if statement is perfectly correct.

% Bad style.
if flag == 1
   fprintf('Failed to converge.\n')
elseif flag == 2
   fprintf('Overflow encountered.\n')
elseif flag == 3
   fprintf('Division by zero.\n')
elseif flag == 4
   fprintf('Tolerance too small.\n')
end

But given its regular structure, it is better written as a switch statement, which is shorter, brings out the parallelism in the logic, and is easier to check.

% Good style.
switch flag
   case 1, fprintf('Failed to converge.\n')
   case 2, fprintf('Overflow encountered.\n')
   case 3, fprintf('Division by zero.\n')
   case 4, fprintf('Tolerance too small.\n')
end

5. Omit Trailing Tildes

In the list of output arguments in a function call a tilde signifies that the output argument in that position is not wanted and so can be discarded (and hence need not be computed). This notation was introduced in MATLAB R2009b. Generally, there is no point in providing a tilde as the last output argument, as a well written function will check the number of requested outputs and not compute any trailing outputs that are not requested.

The singular value decomposition (SVD) of a matrix A is computed by the call [U,S,V] = svd(A). In the following example we wish to compute the left singular vectors of A but not the singular values or right singular vectors.

% Bad style.
[U,~,~] = svd(A);

% Good style
[U,~] = svd(A);

An obvious question is why the second example was not shortened to

U = svd(A);

The reason is that with one output argument the svd function has a special behavior: it returns a matrix of singular values, not the matrix U of left singular vectors.

6. Include an H1 line in Program Files

The H1 line in a MATLAB program file is a comment line that is the second line in the file and contains the program name followed by a one-line description of the program. It is good practice to provide every program with an H1 line, as MATLAB itself does. Several MATLAB functions make use of H1 lines. For example, when the help function is invoked with a directory name and the directory in question does not contain a contents.m file, a list of contents is created from the H1 lines of the program files in the directory and displayed.

Here is an example of a function with an H1-line (this is an updated version of a function from The Matrix Computation Toolbox).

function show(x)
%SHOW   Display signs of matrix elements.
%   SHOW(X) displays X in `FORMAT +' form, that is,
%   with `+', `-' and  blank representing positive, negative
%   and zero elements respectively.

f = get(0,'Format'); % Get current numerical format.
format +
disp(x)
format(f)            % Restore original numeric format.

What’s New in MATLAB R2017a?

MATLAB R2017a was released last week. Many of the changes reported in the release notes are evolutionary, building on and extending major new features introduced previously. For example, the Live Editor continues to gain expanded capabilities. In this post I pick out a few new features that caught my eye. This is very much a personal selection. For full details of what’s new, see the release notes, and see also my previous post What’s New in MATLAB R2016b if you are not familiar with R2016b.

Parula

The parula color map has been modified slightly. The difference is subtle, but as the following example illustrates, the R2017a parula is a bit more vibrant and has a bit less cyan and yellow in the blues and greens.

parula-2016b-2017a.jpg

I note, however, that the difference between the old and the new parula is smaller when the images are converted to the CMYK color space, as they must be for printing.

Heatmap

The new heatmap function plots a heatmap of tabular data. Although it is intended mainly for use with the table data type, I think heatmap will be useful for getting insight into the structure of matrices, as illustrated by the following examples.

heatmaps.jpg

String Arrays

String arrays, introduced in MATLAB R2016b, can now be formed using double quotes:

>> s = string('This is a string') % R2016b
s = 
    "This is a string"
>> t = "This is a string"         % R2017a
t = 
    "This is a string"
>> isequal(s,t)
ans =
  logical
   1
>> whos
  Name      Size            Bytes  Class     Attributes

  s         1x1               166  string              
  t         1x1               166  string

However, there is one major caveat. Many MATLAB functions that take a char as input argument have not yet been adapted to accept strings. Hence

>> A = gallery('moler',3)
A =
     1    -1    -1
    -1     2     0
    -1     0     3
>> A = gallery("moler",3)
Error using nargin
Argument must be either a character vector or a function handle.
Error in gallery (line 191)
nargs = nargin(matname);

I expect that such functions will be updated in future releases.

Missing

A new function missing creates missing values appropriate to the data type in question.

>> A = ones(2); A(2,2) = missing
A =
     1     1
     1   NaN

> d = datetime('2014-05-26')
d = 
  datetime
   26-May-2014
>> d(2) = missing
d = 
  1×2 datetime array
   26-May-2014   NaT

Until converted to the target datatype, a missing value has class missing:

>> m = missing
m = 
  missing
    
>> class(m)
ans =
    'missing'

Performance Improvements

The release notes report performance improvements under a variety of headings, including execution engine, scripts, and mathematics functions. These are very welcome, as the user automatically benefits from them. One comment that caught my eye is “The backslash command A\B is faster when operating on negative definite matrices”. I think this means that MATLAB checks whether the matrix, A, is symmetric with all negative diagonal elements and, if it is, attempts a Cholesky factorization of -A.

Tracing the Early History of MATLAB Through SIAM News

A recent blog post by Ned Gulley points out that the new mathematics gallery (“Mathematics: The Winton Gallery”) at the Science Museum, London, contains a copy of the disk and manual for MATLAB 1.3, from 1985, sitting next to a trial assembly of Charles Babbage’s analytical engine.

This set me thinking about the early history of MATLAB and The MathWorks. Items such as manuals and disks must be quite rare nowadays. What other traces are there of early MATLAB history? I have recently been looking through some back issues of SIAM News and spotted a number of adverts for MATLAB over the period 1985–1991. Let’s see what historical insight these early adverts give.

1985

sinews-1985-18-2-mathworks_ad.jpg PDF file

The first advert I can find is from the March 1985 SIAM News, and it is for PC-MATLAB, priced at $695, running on an IBM-PC or compatible computer. The advert features the now famous MathWorks logo, which represents an eigenfunction of the L-shaped membrane. It quotes a time of 10.1 seconds for a 50-by-50 real matrix multiplication on a machine with an Intel 8087 coprocessor. This floating-point coprocessor was a useful add-on to the IBM-PC, which used the Intel 8086 chip.

MATLAB benefited from being launched at a time when the IBM PC with 8087 coprocessor was just starting to become popular. The 8086-8087 combination made it possible to carry out computations on a desktop PC that had previously required a minicomputer—and they could now be done with the interactive MATLAB interface.

The advert mentions “mainframe MATLAB”, which it says is written in Fortran, runs on “larger computers”, and is in use in several hundred organizations worldwide.

PC-MATLAB had been rewritten in C, and it supported graphics, IEEE arithmetic (as implemented in the 8087), and “user-defined functions” (M-files).

Note that at this time MathWorks was located at 124 Foxwood Road, Portola Valley, California. In his article The Growth of MATLAB and The MathWorks Over Two Decades, Cleve Moler explains that this first mailing address was “a rented A-frame cabin where Jack [Little] lived in the hills above Stanford University in Portola Valley, California”.

The September 1985 issue of SIAM News contains a rather different advert that now refers to “the original mainframe version of MATLAB” and emphasizes the ease of use of MATLAB.

sinews-1985-18-5-mathworks_ad.jpg PDF file

1987

sinews-1987-20-1-mathworks_ad.jpg PDF file

The next advert is from the January 1987 SIAM News. MATLAB is now also available as Pro-MATLAB running on a Sun workstation or a VAX computer. M-files are now mentioned by that name and LINPACK benchmark figures are stated, the largest figure being 98 Kflops on the MicroVAX II.

The MathWorks has now moved to Massachusetts.

1990

sinews-1990-23-3-mathworks_ad.jpg PDF file

The advert in the May 1990 SIAM News gives a version number, MATLAB 3.5, and it announces the Signal Processing Toolbox. It boasts that MATLAB has over 400 built-in functions and supports an enlarged range of computers, which include Cray supercomputers. The benchmark for a 50-by-50 real matrix multiplication is now down to 0.71 seconds on a 20 Mhz 386-based PC: a reduction by a factor of 14 from the figure quoted in 1985.

A strange feature of the advert is that the toolbox is only explicitly mentioned in the box at top left and the meaning of “toolbox” is not stated.

The MathWorks address has changed again, to 21 Eliot St., South Natick, Massachusetts. In the article mentioned above, Cleve explains, “When the company reached about a dozen employees, we moved several miles east to take over the second floor of a lovely building in South Natick, Massachusetts”.

1991

sinews-1991-24-5-mathworks_ad.jpg PDF file

The advert in the May 1991 issue of SIAM News focuses on two new toolboxes: the Spline Toolbox and the Optimization Toolbox. The gray box defines toolboxes as “sets of routines written in the MATLAB programming language for specialized applications”.

MathWorks is now located at Cochituate Place on Prime Park Way in Natick. The street was so-named because it had been the home of the Prime Computer Corporation, which produced minicomputers from 1972 to 1992.

These adverts give some insight into the development of MATLAB in its early years. They also show how the rapid growth of MathWorks necessitated frequent relocation of the company. Indeed, in the article mentioned above Cleve Moler notes that the number of employees roughly doubled every year for the first seven years.

Preparing CMYK Figures for Book Printing

All printing is done in CMYK, the color space based on the four colors cyan, magenta, yellow, and black. Figures that we generate in MATLAB and other systems are invariably saved in RGB format (red, green, blue). When we send a paper for publication we submit RGB figures and they are transformed to CMYK somewhere in the production process, usually without us realizing that it has been done.

If you are one of those people who writes books and prefers to generate the final PDF yourself then you will need to convert any color figures to CMYK. The generation and use of CMYK files is something of a dark art. Here I report what I found out about it when producing the third edition (2017) of MATLAB Guide, co-authored with Des Higham. This is the first edition of the book to use color.

CMYK produces a different range of colors than RGB. Since it is a subtractive color space, designed for inks, it cannot produce some of the brilliant colors that RGB can, especially in the blues. In other words, some RGB colors are out of gamut in CMYK. (Less importantly, the converse is true: some CMYK colors cannot be produced in RGB.) This is something we are used to, but usually do not notice. Whenever we print a document on a laser printer we view a CMYK representation of the colors. In many cases, a figure will look very similar in print and on screen, but there are plenty of exceptions.

The following image shows an RGB image on the left, the result of converting that image to CMYK and then back to RGB in the middle, and a scan of a laser printer’s reproduction of the RGB image on the right. circle-rgb-cmyk-scan.jpg The differences between the RGB version and the other two versions may be shocking! Fortunately, when the CMYK version is viewed on a printed page in isolation from the RGB version (necessarily displayed on a monitor) it does not look so bad. [Of course, if you are reading a printed version of this post then the first two images will look essentially the same.]

After some experimentation, I settled on the following procedure, which I describe for a PDF workflow. For a PostScript workflow, PDF files need to be replaced by EPS files.

  • Print the whole document from an RGB file.
  • Find figures that look very different in print than on screen and select from those any where the difference is not acceptable. In most cases an unacceptable difference will be one where contrast between colors, or saturation of colors, has been reduced.
  • Edit the selected figures in Photoshop, or some other image manipulation package that can save in CMYK form, in order to produce a better looking CMYK file. In Photoshop, Image-Mode-CMYK Color converts an image to CMYK form, but before using this command you should set the correct CMYK working space under Edit-Color Settings (see the Color Space section below). You can edit the RGB image (making use of View-Proof Color to see how the image will look in CMYK, and View-Gamut Warning to see colors that will be out of gamut in CMYK) and then convert to CMYK, or you can convert to CMYK and then edit the file. Save the CMYK image as a PDF file.
  • Generate the PDF file of the book using the edited CMYK figures and the RGB forms of all the other figures, that is, those that did not need special treatment.
  • Load the PDF file into Adobe Acrobat Pro and issue the command
    Tools - Print Production - Convert Colors
    

    This converts all objects in the file to CMYK.

Of the well over 100 figures in MATLAB Guide, only two needed special treatment.

However, there was one problem with this procedure. A handful of images are screen dumps that show a MATLAB window, and these are necessarily low resolution, at 72dpi. When converted to CMYK in Adobe Acrobat these images degraded badly. The solution was to resample the images to 150dpi in Photoshop via Image-Image Size and save to PDF; conversion in Acrobat then worked fine. (Such resampling is probably good practice anyway.)

This is not quite the full story. Here are some more gory details, which are worth reading if you are actually producing a book.

Color Spaces

There is no unique CMYK or RGB color space: these spaces are device dependent. RGB spaces can have different white points and CMYK spaces are designed for particular combinations of paper and ink. Adobe Acrobat Pro offers “U.S. Web Coated (SWOP) v2”, “Uncoated FOGRA29 (ISO 12647-2:2004)” and many other inscrutably named CMYK profiles. So “convert to CMYK” is an ill-defined instruction unless the precise profile is specified. The conversion settings I used are shown in this screen dump of the dialog box, and were provided by the company who printed the book. Your settings might need to be different. acrobat-cymk-conversion.jpg There are a couple of things to note. First, Adobe Acrobat does not remember the settings in the dialog box, so they need to be re-entered every time. Second, since “Embed profile” is not selected with the settings shown, there is no way to tell which profile has been used once the conversion has been done.

Does it really matter what CMYK profile you use? Yes, if you want the best possible results. Consider the following figure. circle-rgb_cmtk-2convs.jpg On the left is the result of converting the original RGB image to “U.S. Web Coated (SWOP) v2” and on the right is the result of converting to “Photoshop 5 default CMYK” (and converting back to RGB in both cases, the latter conversion producing no visual difference in Photoshop). There is a significant difference between the two conversions in every color except white! In principle, both conversions will produce the same result when printed with the target ink and paper, but if you choose the profile without knowing the target you have little idea what the printed result will be.

Photography Books

An important assumption in the process described above is that the accuracy of the color reproduction is not vital. For photography books this assumption clearly does not hold and conversion to CMYK enters a new and frightening realm where images might need to be individually fine-tuned. Indeed high quality photography books typically undergo one or more proof stages using actual galley proofs from the ultimate printing device, sometimes with the author present at the printing press.

Generating CMYK Files in MATLAB

Most of the figures in MATLAB Guide are (naturally) produced in MATLAB. Figures saved with the print function are by default in RGB mode. CMYK is supported for EPS and PS files only, using the ‘-cmyk’ option. This is not very useful if you are using a PDF workflow, as I do.

Instead of using the print function I experimented with export_fig, which is available on MathWorks File Exchange and can save in CMYK format to PDF, EPS, or TIFF files. export_fig creates a PDF file via an EPS file using Ghostscript.

This sounds straightforward, but I found that many of the PDF files generated by export_fig were not fully in CMYK Mode. A PDF file can contain objects in different color spaces and I had files where a color plot was in CMYK but the colorbar was in RGB! The question arises of how one can check whether a given PDF file has any RGB objects. Unfortunately, there seems to be no simple way to do so in Adobe Acrobat Pro. What one can do is invoke

Tools - Print Production - Preflight - PDF analysis - 
List objects using ICCbased CMYK - Analyze

then open up

Overview - color spaces

(I am using Adobe Acrobat Pro XI, 11.0.18; usage could differ in other versions). This should contain “DeviceCMYK color space”, but if it contains “DeviceRGB color space” or “DeviceGray color space” then RGB or Grayscale objects are present. If the latter terms appear in a multi-page PDF file, double clicking on the corresponding icon should take you to the last such object.

You can also open up

Overview - Images

to get a list of pages and object types on those pages. Double clicking on a page icon takes you to that page.

Good Times in MATLAB: How to Typeset the Multiplication Symbol

The MATLAB output

>> A = rand(2); whos
  Name      Size            Bytes  Class     Attributes

  A         2x2                32  double

will be familiar to seasoned users. Consider this, however, from MATLAB R2016b:

>> s = string({'One','Two'})

s = 
  1×2 string array
    "One"    "Two"

At first sight, you might not spot anything unusual, other than the new string datatype. But there are two differences. First, MATLAB prints a header giving the type and size of the array. It does so for arrays of type other than double precision and char. Second, the times symbol is no longer an “x” but is now a multiplication symbol: “×”.

The new “times” certainly looks better. There are still remnants of “x”, for example in whos s for the example above, but I presume that all occurrences of “x” will be changed to the new symbol in the next release

However, there is a catch: the “×” symbol is a Unicode character, so it will not print correctly when you include the output in LaTeX (at least with the version provided in TeX Live 2016). Moreover, it may not even save correctly if your editor is not set up for Unicode characters.

Here is how we dealt with the problem in the third edition (published in January 2017) of MATLAB Guide. We put the code

\usepackage[utf8x]{inputenc}
\DeclareUnicodeCharacter{0215}{\ensuremath{\times}}

in the preamble of the master TeX file, do.tex. We also told our editor, Emacs, to use a UTF-8 coding, by putting the following code at the end of each included .tex file (we have one file per chapter):

%%% Local Variables:
%%% coding: utf-8
%%% mode: latex
%%% TeX-master: "do"
%%% End:

With this setup we can cut and paste output including “×” into our .tex files and it appears as expected in the LaTeX output.

MATLAB Guide, Third Edition (2017)

The third edition of MATLAB Guide, which I co-wrote with Des Higham, has just been published by SIAM. It is a major update of the second edition (2005) to reflect the many changes in MATLAB over the last twelve years, and is 25 percent longer. There are new sections and chapters, and almost every page has changed.

170106-1123-14-5528.jpg

170104-1243-13-5487.jpg

The new chapters are

  • Object-Oriented Programming: presents an introduction to object-oriented programming in MATLAB through two examples of classes.
  • Graphs: describes the new MATLAB classes graph and digraph for representing and manipulating undirected graphs and directed graphs.
  • Large Data Sets: describes MATLAB features for handling data sets so large that they do not fit into the memory of the computer.
  • The Parallel Computing Toolbox: describes this widely used and increasingly important toolbox.

The chapter The Symbolic Math Toolbox has been revised to reflect the change of the underlying symbolic engine from Maple (at the time of the second edition) to MuPAD.

New sections include Empty Matrices, Matrix Properties, Argument Checking and Parsing, Fine Tuning the Display of Arrays, Live Editor, Unit Tests, String Arrays, Categorical Arrays, Tables and Timetables, and Timing Code.

Two other big changes are that figures are now printed in color and there are thirteen “Asides”, highlighted in gray boxes, which contain discussions of MATLAB-related topics, such as anonymous functions, reproducibility, and color maps.

The book was launched with a reception hosted by The Mathworks and SIAM at the SIAM booth at the Joint Mathematics Meetings in Atlanta on January 6, 2017. Jim Rundquist (Senior Education Technical Evangelist) represented MathWorks, and several SIAM staff, including SIAM Publisher David Marshall, were present.

Two delicious cakes, one containing a representation of the cover of the book, were enjoyed by reception attendees. Inspired by MATLAB, the cakes were served using slice, deal, and input, and an occasional reshape or rotate, with a pool of workers consuming them asynchronously.

170106-1547-18-5574-Edit.jpg
Elizabeth Greenspan (Executive Editor, SIAM), Jim Rundquist (Senior Education Technical Evangelist, The MathWorks), Nick Higham (University of Manchester), Bruce Bailey (Book Marketing Associate, SIAM), and David Marshall (SIAM Publisher).

Hyphenation Question: Row-wise or Rowwise?

Sam Clark of T&T Productions, the copy editor for the third edition of MATLAB Guide (co-authored with Des Higham and to be published by SIAM in December 2016), recently asked whether we would like to change “row-wise” to “rowwise”.

A search of my hard disk reveals that I have always used the hyphen, probably because I don’t like consecutive w’s. Indeed, in 1999 I published a paper Row-Wise Backward Stable Elimination Methods for the Equality Constrained Least Squares Problem .

A bit more searching found recent SIAM papers containing “rowwise”, so it is clearly acceptable usage to omit the hyphen..

My dictionaries and usage guides don’t provide any guidance as far as I can tell. Here is what some more online searching revealed.

  • The Oxford English Dictionary does not contain either form (in the entry for “row” or elsewhere), but the entry for “column” contains “column-wise” but not “columnwise”.
  • The Google Ngram Viewer shows a great prevalence of the hyphenated form, which was about three time as common as the unhyphenated form in the year 2000.
  • A search for “row-wise” and “rowwise” at google.co.uk finds about 724,000 and 248,00 hits, respectively.
  • A Google Scholar search for “row-wise” and “rowwise” finds 31,600 and 18,900 results, respectively. For each spelling, there are plenty of papers with that form in the title. The top hit for “rowwise” is a 1993 paper The rowwise correlation between two proximity matrices and the partial rowwise correlation, which manages to include the word twice for good measure!

Since the book is about MATLAB, it also seemed appropriate to check how the MATLAB documentation hyphenates the term. I could only find the hyphenated form:

doc flipdim:
When the value of dim is 1, the array is flipped row-wise down

But for columnwise I found that MATLAB R2016b is inconsistent, as the following extracts illustrate, the first being from the documentation for the Symbolic Math Toolbox version of the function.

doc reshape:
The elements are taken column-wise from A ...
Reshape a matrix row-wise by transposing the result.

doc rmmissing:
1 for row-wise (default) | 2 for column-wise

doc flipdim:
When dim is 2, the array is flipped columnwise left to right.

doc unwrap:
If P is a matrix, unwrap operates columnwise.

So what is our conclusion? We’re sticking to “row-wise” because we think it is easier to parse, especially for those whose first language is not English.