# What’s New in MATLAB R2019a and R2019b?

I didn’t have time earlier this year to write about the first MATLAB release of 2019, so in this post I will discuss R2019a and R2019b together. As usual in this series, I concentrate on the features most relevant to my work.

## Modified Akima Cubic Hermite Interpolation (R2019b)

A new function makima computes a piecewise cubic interpolant with continuous first-order derivatives. It produces fewer undulations than the existing spline function, with less flattening than pchip. This form of interpolant has been an option in interp1, interp2, interp3, interpn, and griddedInterpolant since R2017b and is now also a standalone function. For more details see Cleve Moler’s blog post about makima, from which the code used to generate the image above is taken.

## Linear Assignment Problem and Equilibration (R2019a)

The matchpairs function solves the linear assignment problem, which requires each row of a matrix to be assigned to a column in such a way that the total cost of the assignments (given by the “sum of assigned elements” of the matrix) is minimized (or maximized). This problem can also be described as finding a minimum-weight matching in a weighted bipartite graph. The problem arises in sparse matrix computations.

The equilibrate function take as input a matrix A (dense or sparse) and returns a permutation matrix P and nonsingular diagonal matrices R and C such that R*P*A*C has diagonal entries of magnitude 1 and off-diagonal entries of magnitude at most 1. The outputs P, R, and C are sparse. This matrix equilibration can improve conditioning and can be a useful preprocessing step both in computing incomplete LU preconditioners and in iterative solvers. See, for example, the recent paper Max-Balanced Hungarian Scalings by Hook, Pestana, Tisseur, and Hogg.

The scaling produced by equilibrate is known as Hungarian scaling, and its computation involves solving a linear assignment problem.

## Function Argument Validation (R2019b)

A powerful new feature allows a function to check that its input arguments have the desired class, size, and other aspects. Previously, this could be done by calling the function validateattributes with each argument in turn or by using the inputParser object (see Loren Shure’s comparison of the different approaches).

The new feature provides a more structured way of achieving the same effect. It also allows default values for arguments to be specified. Function argument validation contains many aspects and has a long help page. We just illustrate it with a simple example that should be self-explanatory.

function fav_example(a,b,f,t)
arguments
a (1,1) double               % Double scalar.
b {mustBeNumeric,mustBeReal} % Numeric type and real.
f double {mustBeMember(f,[-1,0,1])} = 0 % -1, 0 (default), or 1.
t string = "-"               % String, with default value.
end
% Body of function.
end


The functions mustBeNumeric and mustBeMember are just two of a long list of validation functions. Here is an example call:

>> a = 0; b = 2; f = 11; fav_example(a,b,f)
Error using fav_example
Invalid input argument at position 3. Value must be a member of ...
this set:
-1
0
1


One caveat is that a specification within the arguments block

v (1,:) double


allows not just 1-by-n vectors (for any n) but also n-by-1 vectors, because of the convention that in many situations MATLAB does not distinguish between row and column vectors.

For more extensive examples, see New with R2019b – Function Argument Validation and Spider Plots and More Argument Validation by Jiro Doke at MathWorks Blogs.

## Dot Indexing (R2019b)

One can now dot index the result of a function call. For example, consider

>> optimset('Display','iter').Display
ans =
'iter'


The output of the dot call can even be indexed:

>> optimset('OutputFcn',{@outfun1,@outfun2}).OutputFcn{2}
ans =
function_handle with value:
@outfun2


Of course, these are contrived examples, but in real codes such indexing can be useful. This syntax was not allowed in R2019a and earlier.

## Hexadecimal and Binary Values (R2019b)

One can now specify hexadecimal and binary values using the prefixes 0x or 0X for hexadecimal (with upper or lower case for A to F) and 0b or 0B for binary. By default the results have integer type.

>> 0xA
ans =
uint8
10
>> 0b101
ans =
uint8
5


## Changes to Function Precedence Order (R2019b)

The Release Notes say “Starting in R2019b, MATLAB changes the rules for name resolution, impacting the precedence order of variables, nested functions, local functions, and external functions. The new rules simplify and standardize name resolution”. A number of examples are given where the meaning of code is different in R2019b, and they could cause an error to be generated for codes that ran without error in earlier versions. As far as I can see, though, these changes are not likely to affect well-written MATLAB code.

## Program File Size (R2019b)

Program files larger than “approximately 128 MB” will no longer open or run. This is an interesting change because I cannot imagine writing an M-file of that size. Good programming practice dictates splitting a code into separate functions and not including all of them in one program file. Presumably, some MATLAB users have produced such long files—perhaps automatically generated.

For a modern introduction and reference to MATLAB, see MATLAB Guide (third edition, 2017)

# What’s New in MATLAB R2018b?

The MATLAB R2018b release notes report a number of performance improvements, including faster startup and faster calls to built-in functions. I pick out here a few other highlights from the release (excluding the toolboxes) that are of interest from the numerical analysis point of view.

The new xline and yline functions add vertical or horizontal lines to a plot—something that I have needed from time to time and have previously had to produce with a line or a plot command, together with hold on. For example, xline(pi) plots a vertical line at x = pi.

The stackedplot function is mainly of interest for plotting multiple variables from a table in a single plot, but it can also plot the columns of a matrix. In this example, A is the symmetric eigenvector matrix for the second difference matrix:

A = gallery('orthog',10,1); stackedplot(A);


The resulting plot clearly shows that the number of sign changes increases with the column index.

String arrays, introduced in R2016b, are now supported throughout MATLAB. In the previous example I could have typed A = gallery("orthog",10,1).

A new option 'all' to the functions all, any, max, min, prod, sum (and a few others) makes them operate on all the dimensions of the array. For example:

>> A = pascal(3)
A =
1     1     1
1     2     3
1     3     6
>> max(A,[],'all')
ans =
6
>> [prod(A,'all'), sum(A,'all')]
ans =
108    19


The empty second argument in the max call is needed because max(x,y) is also a supported usage. This is a useful addition. I have often written norm(A(:),inf) to avoid max(max(abs(A))) (which does not work for arrays with more than two dimensions), but now I can write max(abs(A),[],'all') without incurring the overhead of forming A(:).

New functions sinpi and cospi plot the sine and cosine functions at the specified multiple of $\pi$. Thus sinpi(x) is the same as sin(pi*x) except that it does not explicitly compute pi*x and so returns exact answers for integer arguments:

>> [sin(pi)         sinpi(1)]
ans =
1.2246e-16            0


# What’s New in MATLAB R2018a?

MATLAB R2018a was released in March 2018. With each biannual release I try to give a brief overview of the changes in MATLAB (not the toolboxes) that are of most interest to me. These are not comprehensive summaries of what’s new and you should check the release notes for full details.

Complex empty arrays, such as that generated with complex([]) now have an (empty) imaginary part instead of being real.

% R2017b
>> complex([])
ans =
[]

% R2018a
>> complex([])
ans =
0×0 empty complex double matrix


This is a consequence of a change in the way MATLAB stores complex arrays. Prior to R2018a it stored the real parts together and the imaginary parts together. In R2081 it uses an interleaved format in which the real and imaginary parts of a number are stored together; this is the storage format used in the C and C++ languages. For more details see MATLAB Support for Interleaved Complex API in C MEX Functions. Unless you write MEX files or create complex empty arrays, this change should have no effect on you.

The Live Editor continues to gain improved functionality, including embedded sliders and drop-down menus.

Legends can now have multiple columns, specified with the NumColumns property for Legend objects. I must admit that I had not realized that is was already possible to arrange legends in a horizontal (1-by-n) rather than vertical orientation (n-by-1) using the Orientation property.

Tall arrays have several new capabilities, including the ability to compute a least squares solution to an overdetermined linear system Ax = b with a tall A, which is done by QR factorization.

MATLAB starts up faster and has further optimizations in the execution engine.

The GraphPlot Object has some additional options for the 'force', 'force3', and 'circle' layouts.

I noticed that the recent Windows 10 Spring Creators Update broke the Symbolic Toolbox. An update to fix this (and various other issues) is available at this page (MathWorks login required).

For a concise but wide-ranging introduction and reference to MATLAB, see MATLAB Guide (third edition, 2017)

# 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 $1$s and $-1$s:

>> 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.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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: 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:

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?

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

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.

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.

## 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

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.

## 1987

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

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

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. 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. 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. 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.