[Math] Reasonable “Random” matrices to test numerical algorithms

algorithmsapplied-mathematicsna.numerical-analysissoft-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

Best Answer

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

Related Solutions

[Math] Why does randomness work in numerical algorithms

If one had infinite computing capacity, then you would be right that randomness would not be useful (unless an adversary also had access to randomness - see below), and it would always be better to maintain full control of one's parameters. ("God does not need to play dice".) However, our computing capacity is limited, and randomness offers a way to stretch that capacity further (though it is an important open problem to quantify exactly how much randomness actually stretches our capacity; see the discussion on P=BPP in other comments).

Let's give a simple example. Suppose one is playing the following game with a computer. This computer uses some (deterministic) algorithm to generate 1000 "bad" numbers from 1 to 1 million. Meanwhile, you pick your own number from 1 to 1 million. If you pick one of the bad numbers, you lose; otherwise, you win. You get to see the computer's algorithm in advance.

With infinite computing capacity, you have a strategy which is 100% guaranteed to succeed: you run the computer's algorithm, find all the bad numbers, and pick a number not chosen by the computer.

But there is a much lazier strategy that requires no computing power: just pick a number randomly, and one has a 99.9% chance of winning.

[Note also that if the computer used a randomised algorithm instead of a deterministic one, then there is no longer any 100% winning strategy for you, even with infinite computing power, though the random strategy works just fine. This is another indication of the power of randomness.]

Many randomised algorithms work, roughly speaking, by reducing the given problem to some version of the above game. For instance, Monte Carlo integration using a set of points to sample the function might give satisfactory levels of accuracy for all but a tiny fraction of the possible choices of sample points. With infinite computing capacity, one could locate all the bad configurations of sample points and then find a configuration which is guaranteed to give a good level of accuracy, but this is a lot of work - more work, in fact, than just doing classical numerical integration. Instead, one takes the lazy way out and picks the sample points randomly. (Well, in practice, one picks the points pseudorandomly instead of randomly, because true randomness is very hard to capture by a computer; this is a subtle issue - again related to P=BPP - that I don't want to get into here.)

[Math] Pseudo-random number generation algorithms

Don't miss this wonderful post by Marsaglia. He's not a fan of the Mersenne Twister and offers some strong PRNGs with exceptionally small code footprints. One of his examples is:

static unsigned long 
x=123456789,y=362436069,z=521288629,w=88675123,v=886756453; 
      /* replace defaults with five random seed values in calling program */ 
unsigned long xorshift(void) 
{unsigned long t; 
 t=(x^(x>>7)); x=y; y=z; z=w; w=v; 
 v=(v^(v<<6))^(t^(t<<13)); return (y+y+1)*v;}

Best Answer

Related Solutions

[Math] Why does randomness work in numerical algorithms

[Math] Pseudo-random number generation algorithms

Related Question