Reasonable "Random" matrices to test numerical algorithms

Question

Hello,

in numerical analysis, it is common to compare the behavior of different algorithms, and of different implementation of algorithms. This occurs not only on the theoretical level, but also on the concrete level of implementation - and not to forget, it serves the purpose of demonstration.

A prominent problem is the solution of linear systems, both general as well as various subcases.

To test, benchmark and profile numerical implementations, you run your work on several instances of the problem. However, it is difficult question to obtain a good set of these instances. You want to inspect pathological cases (diff. degrees of ill-conditionedness) as well as "real-life" examples (whatever this may mean). Ideally, you have an algorithm which puts out matrices with certain properties in a "reasonable" probability measure. A good notion of "reasonable" might be accessible, as most such LSE problems from physics or simulations have much more structure as is actually demanded by the algorithms in theory.

In so far, I wonder whether there are works in numerical analysis how to, given $n \in \mathbb N$ randomly produce

• a sequence of ${n \times n}$-matrices
• optionally constraint to be symmetric, positive definite, well-conditioned
• which is reasonable in whatever sense

This is probably an interesting topic within the theory of numerical algorithms.

Thanks, Martin

Answer 1

A good set of benchmark matrices depends on the problem being solved (sparse solvers, eigenvalue problems, special structures, et cetera). Often, you'll want to include some real-life example applications of the algorithm that you are testing, or matrices that lead to particularly difficult problems to solve. Just choosing matrices at random won't cut that. Hence good sets of benchmarks mostly contain matrices that are carefully chosen, and are publication-worthy in themselves; some get lots of citations.

Among many of them, I shall mention here:

• the all-purpose Matlab's gallery. In particular tweaking the parameters in randsvd could address quite well the requirements in your original question.

• for sparse matrices, Matrix Market and the University of Florida's sparse matrix collection.

• for control problems and Hamiltonian/symplectic matrices, carex and darex by Peter Benner.
• for image processing, the USC image database, including the long-debated Lenna (originally a cropped version of a Playboy centerfold, now deemed by many not politically correct to use).

Answer 2

One way to generate random matrices while constraining it, is to generate its LU decomposition first. That way you can restrict it to be symmetric ($L=U^T$) and gives you the control over its spectrum.

In [S.M. Rump. A Class of Arbitrarily Ill-conditioned Floating-Point Matrices. SIAM J. Matrix Anal. Appl. (SIMAX), 12(4):645-653, 1991] there is a related method for generating very ill-conditioned matrices.

Answer 3

One way to get useful information on a numerical algorithm is to choose by hand a well-conditioned or even more an ill-conditioned example and perturb it randomly. In fact it is useful to perturb it randomly such that the perturbations form a sequence approaching zero; you can then look at the sequence of solutions that you get and see what you can say about error propagation and the like.

Answer 4

This is an impossible question to answer (so maybe this should be a comment...), since (for example) the sparsity patterns of the matrices encountered are completely different in statistical applications and in finite elements -- each problem class leads to a completely different distribution on the space of matrices. As a result, all the papers I have seen in numerical analysis are of the sort: we tried heuristic X and it worked well for problem Y. This gives us reasonable confidence that it will work for problem Z, such that $|Z-Y| < \epsilon.$ People are very careful not to make general statements.

DISCLAIMER: I am not a professional numerical analyst, nor do I play one on TV.