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

  • $\begingroup$ I think the ultimate goal of the question is valid and useful, though it requires more thought to first state more precise questions (several concerns exist, including those shared with random number generation, statistical tests of randomness, sampling from matrix valued distributions, complexity theoretic concerns, etc. etc.) $\endgroup$ – Suvrit Jan 1 '11 at 9:57
  • $\begingroup$ Another remark: it is often not a good idea to pick random matrices as examples, as the experiments cannot be reproduced and compared to other algorithms. Also, no one can tell if you restarted your experiments several times until you got a "good-looking" graph. If you still use random matrices, e.g. for completeness, it is a good idea to set the random seed to a fixed number before you start your experiments, and to use well-known and reproducible generators (such as the one in MATLAB, which, albeit commercial code, is the de-facto standard in numerical linear algebra). $\endgroup$ – Federico Poloni Jan 1 '11 at 18:13

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:


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.


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.


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.

  • $\begingroup$ I would like to inquire about which television show has a numerical analyst as a character. $\endgroup$ – Michael Lugo Feb 10 '11 at 5:44
  • $\begingroup$ Sadly, there is not one -- otherwise I would be able to fill that much-needed gap in my resume. $\endgroup$ – Igor Rivin Feb 10 '11 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.