Exploiting the MATLAB Language in Implementing a Matrix Collection
“Concise code” may not be commonly used or defined, but as in languages, conciseness—the quality of communicating information in as few words as is feasible—can be a good guide in programming. Conciseness of code, which could mean short and self-explanatory lines or blocks of code, is easier to achieve in higher level languages such as Python or MATLAB, and more difficult when writing, for example, in RISC assembly languages with short instructions each doing very little.
For developing the recently released matrix collection, Anymatrix [1,2] (with Nick Higham), we used MATLAB and exploited its string handling routines to do complicated string manipulations in a few lines of code. This post describes some details of the implementation and further information can be found in the user’s guide [3].
Groups of Matrices
Anymatrix enforces the organization of matrices into groups, which means a new matrix has to be placed either into one of the existing groups or a new group. For example, the built-in groups are contest, core, gallery, hadamard, matlab, nessie, and regtools. Each of these come with around 20 matrices, which can be fixed-sized matrices or matrix generators.
Matrices are identified by
Around the time we started working on this collection, MATLAB R2020b introduced pattern, to search and match strings conveniently. An example usage in Anymatrix is as follows.
% Matrix ID pattern
matrix_ID_pat = ...
asManyOfPattern(alphanumericsPattern | ...
characterListPattern("_") | ...
characterListPattern("-")) + ...
'/' + ...
asManyOfPattern(alphanumericsPattern | ...
characterListPattern("_") | ...
characterListPattern("-"));
This code creates a pattern for matrix IDs, which matches sequences of letter, digit, underscore or hyphen characters, followed by a single forward slash, followed again by the sequence. This sort of pattern can then be used, for example, to match strings (using matches) and extract substrings (using extract):
>> S = '...core/beta...';
>> extract(S, matrix_ID_pat)
ans =
1×1 cell array
{'core/beta'}
The pattern functionality helped in writing concise code. We used it in various places of Anymatrix, for example in detecting matrix IDs in the commands and extracting substrings in the Boolean search expressions to separate properties from brackets and operators and/or/not.
Properties of Matrices
Matrix properties can be specified in two ways, either as an assignment properties = {...} inside the matrix .m files, or in the am_properties.m files in group folders. With the latter option we may specify properties of multiple matrices in one assignment. Anymatrix contains a list of recognized properties, such as ill conditioned or banded. When Anymatrix scans the file structure, it gives a warning if a matrix contains an unrecognized property.
Properties can be searched by Boolean expressions with operators and/or/not and brackets for specifying precedence. For example:
>> anymatrix('p', 'integer and ill conditioned')
ans =
3×1 cell array
{'core/wilson' }
{'gallery/dramadah'}
{'matlab/pascal' }
The following function that transforms the supplied Boolean expression into an expression that can be evaluated in MATLAB, demonstrates a lot of MATLAB’s features that we made use of.
function IDs = search_by_properties(expression)
IDs = {};
% Replace 'and', 'or', and 'not' by corresponding MATLAB symbols.
expression = replace(expression, ' and ', ' & ');
expression = replace(expression, ' or ', ' | ');
expression = replace(expression, ' not ', ' ~');
expression = replace(expression, '(not ', '(~');
if startsWith(expression, 'not')
expression = expression(4:length(expression));
expression(1) = '~';
end
% Assume properties are made up letters, can include a hyphen
% or a white space character, and there is no case sensitivity.
pat = (lettersPattern + whitespacePattern + lettersPattern) ...
| (lettersPattern + characterListPattern('-') ...
+ lettersPattern) ...
| lettersPattern;
% Extract properties from the logical expression and replace
% them by function calls to test for membership.
ex_props = extract(expression, pat);
ind = 1;
new_expression = '';
for p = ex_props.'
mod_prop = strcat('ismember(''', ...
strrep(lower(p{1}), '-', ' '), ...
''', strrep(lower(matrix_properties{1}), ''-'', '' ''))');
trunc_exp = expression(ind:end);
% Find where the property is in the expression.
prop_index = strfind(trunc_exp, p{1});
% Take anything before the property and append the modified
% version of the property.
new_expression = strcat(new_expression, ...
trunc_exp(1:prop_index(1)-1), ...
mod_prop);
% Move the index after the property that was replaced.
ind = ind + prop_index(1) + length(p{1}) - 1;
end
new_expression = strcat(new_expression, expression(ind:end));
% Find matrices whose properties satisfy the specified logical
% expression.
k = 0;
for matrix_properties = properties.'
k = k + 1;
% Test if the expression is true for this matrix
% and add it's ID.
if eval(new_expression)
IDs = [IDs; matrix_IDs{k}];
end
end
end
On lines 4 to 11 we replace the operators and/or/not by the corresponding MATLAB ones &/|/~. We then extract properties from the expression, making use of patterns, modify them to test for membership of the properties and place them back into the expression (lines 14 to 38). Finally, on lines 42 to 50 we use eval to run the Boolean expression as MATLAB code for every matrix, and return a list of suitable matrices.
In the example command anymatrix('p', 'integer and ill conditioned'), the expression 'integer and ill conditioned' is transformed internally to 'ismember('integer', strrep(lower(matrix_properties{1}), '-', ' ')) & ismember('ill conditioned', strrep(lower(matrix_properties{1}), '-', ' '))' which can then be evaluated using eval. Notice that Anymatrix replaces hyphens by white spaces to allow two ways to specify properties, and converts properties to lower case to avoid case sensitivity. This way, ill conditioned and Ill-conditioned, for example, are equivalent properties in Anymatrix.
Testing in Anymatrix
There are two ways to implement testing in Anymatrix. One way is to add a function in file test_run.m in the group folder which can be invoked by an appropriate Anymatrix command. Another way is to test matrices for their properties. This is provided in the directory testing/ in the root folder of Anymatrix.
It is worth noting that not all supported properties can be tested, so only a subset in the built-in collection are. Each property that is feasible to test has a corresponding test_
function test_core_fourier(testcase)
A = anymatrix('core/fourier',3);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',5);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',8);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',10);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',15);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',24);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',25);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',30);
anymatrix_check_props(A, 'core/fourier', testcase);
A = anymatrix('core/fourier',31);
anymatrix_check_props(A, 'core/fourier', testcase);
end
Above, the function anymatrix_check_props looks for files test_
for prop = P.'
test_func_name = strcat('test_', ...
strrep(strrep(lower(prop{1}), ...
'-', '_'), ' ', '_'));
if isfile(strcat(root_path, '/private/', test_func_name, '.m'))
handle = str2func(test_func_name);
verifyTrue(testcase, handle(M), ...
strcat("Matrix ", matrix_ID, " is not ", prop{1}, "."));
end
end
Using this, Anymatrix adds a level of almost automatic testing invoked by adding new matrices or properties to matrices. A good example of this is the MATLAB built-in gallery group, which we made available in Anymatrix. By appending each matrix with properties we can check that MATLAB’s gallery matrices have the expected properties.
Summary
Anymatrix is both a collection of matrices and a tool to organize them. It is worth noting that much of the infrastructure is not specific to organizing matrices and so it can be reused to organize into groups, search, and access through IDs any kind of objects appended with properties.
Anymatrix 1.0 was released in October 2021, but the development continues. We welcome contributions to the collection either in the form of remote groups that can be downloaded into Anymatrix given a git url (an example is matrices-mp-cosm group by X. Liu), or suggestions to include matrices or groups in the built-in collection. We also welcome reports of bugs or recommendations for improvements to the clarity, correctness and conciseness of the source code.
References
[1] N. J. Higham and M. Mikaitis. Anymatrix: An Extensible MATLAB Matrix Collection. MIMS EPrint 2021.16, Manchester Institute for Mathematical Sciences, The University of Manchester, UK. Oct. 2021 (to appear in Numer. Algorithms).
[2] N. J. Higham. Anymatrix: An Extensible MATLAB Matrix Collection. Nov. 2021.
[3] N. J. Higham and M. Mikaitis. Anymatrix: An Extensible MATLAB Matrix Collection. Users’ Guide. MIMS EPrint 2021.15, Manchester Institute for Mathematical Sciences, The University of Manchester, UK. Oct. 2021.
and 
, the solution of this ODE lies on the unit circle in the 
-plane [p. 51, 5]. When solving it numerically our is aim to obtain the solution that lies as close as possible to the unit circle. A similar example is approximating orbits of asteroids in Physics [Sec. II A, 8], although it is more complicated since it additionally includes velocity and acceleration.
is a step size. Now, solving this with 
in bfloat16 arithmetic [6] with round to nearest (RN) and stochastic rounding in arithmetic operations, we obtain the solution in Figure 1 (dotted line—exact unit circle, dashed—RN, solid—SR), which spirals out as previously shown [p. 51, 5]. To improve this we could use a different solver, such as Trapezium or Leapfrog [p.51, 5], but for the purposes of testing arithmetics rather than solvers we stick with Euler and try different step sizes.
-axis represents 
and the 
-axis represents 
. Dotted line—exact unit circle, dashed—RN, and solid—SR.
and the plot is shown on the left hand side of Figure 2. Visually, this timestep size for the bfloat16 arithmetic with RN looks ideal since it quite closely approximates the unit circle and the solution does not spiral out as before. SR performs quite well here as well but noticeably worse than RN.
and the solution with that is plotted on the right hand side of Figure 2. In this case the solution with RN has been affected by rounding errors—the approximation of the unit circle looks visually as an octagon rather than a circle!
and 
, we expect that both 
the term 
is too small initially to change 
.
(left) and 
.
, when there is no higher precision available in hardware (therefore 
stage of Gaussian elimination are generated from the formula 
, where 
is a multiplier. If 
is a dense matrix of normalized floating-point numbers with norm of order 1, it is extremely unlikely that any of the 
will become subnormal. However, for sparse matrices we can identify a mechanism whereby fill-in cascades down a column and small multipliers combine multiplicatively. This mechanism is described, with an illustration for better understanding, in the EPrint, as well as strategies to automatically flush subnormals to zero, to avoid the performance penalties.
when single precision is used in place of double precision.
for a dense 
matrix with random entries. The figure on the right displays the performance in GigaFLOPS. Unlike the standard algorithm, the new one does not incur a significant loss of accuracy by storing the matrix in half precision. In doing so, it greatly reduces data movement which leads to significant speedup: the new algorithm can be up to twice faster than the standard one.
of a matrix, with the function dependence restricted to the evaluation of 
function, which implements the Schur-Parlett algorithm, is not able to compute. Our new algorithm greatly expands the class of readily computable matrix functions, and its MATLAB code is available here.
of order 
to binary64 accuracy by complementing a low-precision direct solver with a high-precision refinement strategy. The code generates 
and 
in binary64, computes the LU factorization without pivoting of 
using forward and back substitution, and finally solves the linear system with GMRES in binary64 arithmetic, using 
, it should also be easy to generate at scale on a distributed memory environment. Furthermore, we require that 
-norm condition number just above 4 for large 
for which LU factorization without pivoting is stable. The parameters 
and 
that yield a specified 
floating-point operations.
. More recently (
years ago) the FMA (fused multiply-add) operation also became prevalent in general purpose devices such as CPUs and GPUs. Software builds on top of these operations, for example compilers use a series of these and other hardware instructions to implement more complex functions and algorithms such as the exponential function, matrix multiply, or the algorithms for pseudo random number generation. These arithmetic operations have been standardized by the IEEE 754 floating-point standard since 1985 and most current systems are compliant with it.
in mixed-precision of fp16 and fp32, as shown in Figure 1. The most recent revision of tensor cores in the A100 updated both the precisions available and the dimensions of input/output matrices. See the NVIDIA V100 and NVIDIA A100 whitepapers for more details.
refers to the matrix multiply-accumulate operation where 
and 
are 
matrices, and 
have size 
and 
, respectively.
and a value from 
.
.
.
and 
minus the distance of 
that are random variables satisfying
,
.
denotes the expectation. A key consequence of these results is that we can replace the worst case error bound proportional to 
, ubiquitous in backward error analyses, by a more informative probabilistic bound proportional to 
. What is a rule of thumb for round to nearest becomes a rule for stochastic rounding: it is proved that our rounding errors satisfy the above properties.
of vectors sampled uniformly from ![[0,1]](https://s0.wp.com/latex.php?latex=%5B0%2C1%5D&bg=ffffff&fg=111111&s=0&c=20201002)
in both round to nearest and stochastic rounding in IEEE half precision arithmetic. Shown is the backward error for each value of 
and 
. 
bound. This is due to the phenomenon we call stagnation, and low precisions are particularly susceptible to it. As we compute recursively 
, eventually the partial sum 
can grow so large that the update 
is less than half of the spacing of floating-point numbers around 
, that is, the expected value of the computed result under stochastic rounding is the exact result.
Numerical Linear Algebra Group