Wiki

Clone wiki

NumericalUniversality / Home

Numerical Universality

This Mathematica library contains the software used to perform all the computations in

P. Deift, G. Menon, S. Olver and T. Trogdon, arXiv:1407.3829 [math.NA], 2014 pdf

The package also has some plotting functionality and should make any extensions of the software easy to perform. The following algorithms are represented. A tutorial, by example, can be found in Tutorial.nb.

Jacobi Algorithm

The package implements a simple-minded database. The Available command allows one to explore what is in the database:

#!Mathematica
In[3]:= Available["Jacobi", "GOE"]

Out[3]= {"10-10.m", "10-8.m", "11-8.m", "30-10.m", "30-8.m", 
"50-10.m", "50-8.m", "70-10.m", "70-8.m", "90-10.m", "90-8.m"}

In[4]:= Available["Jacobi", "GOE"]
Available["Jacobi", "BE"]

Out[4]= {"10-10.m", "10-8.m", "11-8.m", "30-10.m", "30-8.m", 
"50-10.m", "50-8.m", "70-10.m", "70-8.m", "90-10.m", "90-8.m"}

Out[5]= {"10-10.m", "10-8.m", "30-10.m", "30-8.m", "50-10.m", 
"50-8.m", "70-10.m", "70-8.m", "90-10.m", "90-8.m"}
Typically, with the exception of Curie-Weiss, the filenames are of the form N-n.m where N is the dimension (i.e. matrix size) and n is the exponent in the tolerance: 10^(-n).

Use the HaltingDatabase command to access or create data:

#!Mathematica
In[8]:= HaltingDatabase["Jacobi", "GOE", 50, 10, 16000][[1 ;; 10]]

Out[8]= {4510, 4565, 4531, 4584, 4560, 4622, 4530, 4594, 4514, 4558}

In[9]:= HaltingDatabase["Jacobi", "GOE", 51, 10, 16000]

Getting 16000 more samples.

By default the code uses the ParallelTable Mathematica routine. Use ParallelFalse[] to turn this off and perform the calculations serially. Both the PlotSamplesPaper and PlotSamplesWide commands can be used to produce histograms:

PlotSamplesPaper["Jacobi", {"GOE", "BE"}, {30, 50, 70, 
  90}, 10, 16000, 4, {{-3, 3}, {0, .6}}, "Matrix size(s)"]

It takes in multiple ensembles and multiple dimensions and overlays all histograms once normalized to mean zero and variance one.

Tutorial1.jpg

QR Algorithm

Here is some sample code to access the database for the QR algorithm

In[11]:= Available["QR"]
Available["QR", "GUE"]
Available["QR", "QUE"]
Available["QR", "CoshUE"]

Out[11]= {"CoshUE", "GOE", "GUE", "QUE"}

Out[12]= {"10-10.m", "110-10.m", "130-10.m", "150-10.m", "20-10.m", 
"30-10.m", "40-10.m", "50-10.m", "70-10.m", "90-10.m"}

Out[13]= {"10-10.m", "110-10.m", "150-10.m", "20-10.m", "30-10.m", 
"50-10.m", "70-10.m", "90-10.m"}

Out[14]= {"10-10.m", "10-8.m", "110-10.m", "150-10.m", "30-10.m", 
"70-10.m"}

The spectra for the CoshUE and QUE matrices is in the SpectraDatabase folder. It is not generated as it is called. For this reason SampleHaltingTime for these ensembles should not be called on its own.

SampleHaltingTime["QR", "CoshUE", 110, 10]

Define CoshCounter

# of requested samples exceeded # in SpectraDatabase, counter undefined or no SpectraDatabase entry exists

If it is called after HaltingDatabase it may produce dependent results.

To plot:

PlotSamplesWide["QR", {"CoshUE", "GUE", "QUE"}, {150, 110}, 10, 16000,
  1]

Tutorial2.jpg

Conjugate Gradient Algorithm

Here is sample code for the conjugate gradient algorithm.

In[4]:= Available["CG"]
Available["CG", "cLOE"]
Available["CG", "cPBE"]

Out[4]= {"cLOE", "cPBE"}

Out[5]= {"100-8.m", "200-8.m", "300-8.m"}

Out[6]= {"100-8.m", "200-8.m", "300-8.m"}

PlotSamplesPaper["CG", {"cPBE", 
  "cLOE"}, {300}, 8, 16000, 1, {{-3, 3}, {0, .5}}, "Matrix size(s)"]
Tutorial3.jpg

GMRES Algorithm

Here is sample code for the GMRES algorithm

Available["GMRES"]
Available["GMRES", "BDE"]
Available["GMRES", "cSBE"]
Available["GMRES", "cSGE"]
Available["GMRES", "UDE"]

{"BDE", "cSBE", "cSGE", "UDE"}

{"100-8.m", "200-8.m", "300-8.m", "400-8.m", "500-8.m", "600-8.m"}

{"100-8.m", "200-8.m", "300-8.m", "400-8.m", "500-8.m", "550-8.m", 
"600-8.m"}

{"100-8.m", "200-8.m", "300-8.m", "400-8.m", "500-8.m", "550-8.m", 
"600-8.m"}

{"100-8.m", "150-8.m", "200-8.m", "300-8.m", "400-8.m", "500-8.m", 
"50-8.m", "600-8.m"}
PlotSamplesPaper["GMRES", {"BDE", "cSBE"}, {100, 
  200}, 8, 16000, 1, {{-3, 3}, {0, .5}}, "Matrix size(s)"]

Tutorial4.jpg

Genetic Algorithm

The syntax for the Genetic algorithm is modified slightly because there are two sub-methods.

In[12]:= Available["Genetic"]
Available[{"Genetic", "Gaussian"}]
Available[{"Genetic", "Split"}]
Available[{"Genetic", "Gaussian"}, "Bernoulli"]
Available[{"Genetic", "Gaussian"}, "Uniform"]
Available[{"Genetic", "Split"}, "Bernoulli"]
Available[{"Genetic", "Split"}, "Uniform"]

Out[12]= {"Gaussian", "Split"}

Out[13]= {"Bernoulli", "Exact", "Uniform"}

Out[14]= {"Bernoulli", "Exact", "Uniform"}

Out[15]= {"10-2.m", "20-2.m", "30-2.m", "40-2.m"}

Out[16]= {"10-2.m", "20-2.m", "30-2.m", "40-2.m"}

Out[17]= {"10-2.m", "20-2.m", "30-2.m", "40-2.m"}

Out[18]= {"10-2.m", "20-2.m", "30-2.m", "40-2.m"}
In[20]:= HaltingDatabase[{"Genetic", "Gaussian"}, "Bernoulli", 20, 2, 
  16000][[1 ;; 10]]

Out[20]= {281, 746, 220, 558, 433, 426, 402, 489, 447, 434}
PlotSamplesWide[{"Genetic", "Gaussian"}, {"Bernoulli", 
  "Uniform"}, {10, 20}, 2, 16000, 10]

Tutorial5.jpg

Curie-Weiss System

Finally, here are the commands for the Curie-Weiss decision model.

In[5]:= Available["CurieWeiss"]
Available["CurieWeiss", "Original"]
Available["CurieWeiss", "Cubic"]

Out[5]= {"Cubic", "Eight", "Original", "Remove", "Shift"}

Out[6]= {"100-0.5-1.3.m", "10-0.5-1.3.m", "110-0.5-1.3.m", \
"120-0.5-1.3.m", "130-0.5-1.3.m", "150-0.5-1.3.m", "170-0.5-1.3.m", \
"190-0.5-1.3.m", "20-0.5-1.3.m", "210-0.5-1.3.m", "230-0.5-1.3.m", \
"250-0.5-1.3.m", "30-0.5-1.3.m", "40-0.5-1.3.m", "50-0.5-1.3.m", \
"60-0.5-1.3.m", "70-0.5-1.3.m", "80-0.5-1.3.m", "90-0.5-1.3.m"}

Out[7]= {"10-0.5-1.3.m", "110-0.5-1.3.m", "130-0.5-1.3.m", \
"150-0.5-1.3.m", "170-0.5-1.3.m", "190-0.5-1.3.m", "210-0.5-1.3.m", \
"30-0.5-1.3.m", "50-0.5-1.3.m", "70-0.5-1.3.m", "90-0.5-1.3.m"}
Notice that there is the extra argument 1.3 in the filename. This is the inverse temperature in the model. It is hard-coded at this value but can be changed for further simulations. The other difference here is that n is no longer the exponent of the tolerance. It is now the threshold for decision. Once the average over all spins, in absolute value, is greater than this number a "decision" is made.

PlotSamplesPaper["CurieWeiss", {"Cubic", "Original"}, {210, 
  150}, .5, 16000, .2, {{-2, 5}, {0, .8}}]

Tutorial6.jpg

Updated