Randomness in zgetu0w()

This is about the function zgetu0w() in the Fortran code. It uses random numbers, but they're not as random as you might think. The Fortran version of zgetu0w() uses the same sets of random numbers every time the code is run.

zgetu0w() constructs a random array in order to do…whatever is that function does. (I don't know; it has no comments.) It calls the LAPACK function zlarnv() to generate that array.

zlarnv() requires as input an array of four integers (why four?) that are used to seed the random number generator. The RNG is a pseudo-RNG, that is, it will always return the same values given the same seeds. zgetu0w() takes the unusual approach of hardocding the four seeds, so it receives the same output from the RNG every time it runs.

There's evidence in commented-out code that other seeds were tried and discarded, but no information as to why some were deemed better than others.

Effects on the Fortran Code

I added the function generate_seeds() to zget0w.f as an experiment. It generates random seeds (which aren't particularly random themselves, but they're better than hardcoding) to feed to zlarnv(). That revealed something interesting: using different seeds affects the outcome.

One easy-to-observe difference is that the RMSE reported by my test harness varies depending on the seeds used. For instance, when using the "standard" hardocded seeds, the data file press_cp3 generates an RMSE of 11.73. I get the same RMSE for most sets of random seeds. But some sets of seeds (~10%) generate an RMSE of 11.54, which is a slightly better result.

From this I conclude the following –

  • The algorithm is not entirely stable, which is probably no surpise given its complexity.
  • The hardcoded seeds are not optimal for all inputs.

I would also be willing to bet the following –

  • Seeds that produce an optimal output for one set of input may not produce optimal output for a different set of input. In other words, there is no set of seeds that's optimal for all inputs.

Effects on the Python Code

While porting, I copied the array generated by zlarnv() to a text file and read that into Python to ensure that Python used the same "random" numbers as Fortran. Once my code produced the same results, I switched to using "truly" random numbers from Python's RNG.

The results are similar but different. For press_cp3, I still get RMSEs of 11.54 and 11.74, but the former is the result of about ~95% of runs.

One very disappointing result is that when using truly random numbers in the Python code, the data file laser_2010102901.xml occasionally (~3%) causes the code to generate a garbage result. Whether this is a result of my mistake in my port or a revelation of an algorithmic weakness, I don't know.

The Fortran code consistently produces a good result with laser_2010102901.xml even when using random seeds. However, in a small percentage of cases (< 1%) it produces an RMSE of 1.15 instead of the usual 0.90.

Interestingly enough, the Fortran code calls zgetu0w() an unusually large number of times for some seeds, but that apparently isn't enough to drive it to a bad outcome. It might be interesting to capture the random array(s) that cause the Fortran code to struggle, feed them to the Python code and see what happens.

Neither the Python nor the Fortran produces a garbage result with any data file other than laser_2010102901.xml.