How To Run An Analysis
This example shows how to run VecchiaMLE for a given VecchiaMLEInput. There are some necessary parameters the user needs to input, which are the following:
using VecchiaMLE
using SparseArrays # To show matrix at the end
# Things for model
n = 20
k = 10
Number_of_Samples = 100100These are the three major parameters: n (dimension), k (conditioning for Vecchia), and the NumberofSamples. they should be self-explanatory, but the documentation for the VecchiaMLEInput should clear things up. Next, we should generate the samples. At the moment, the code has only been verified to run for generated samples by the program, but there should be no difficulty in inputting your own samples. Just make sure the length of each sample should be a perfect square. This is a strange restriction, and will be removed when we allow user input locations. Nevertheless, let's move on.
The sample generation requires a covariance matrix. Here, we will generate a matern like covairance matrix which VecchiaMLE has a predefined function for. Although not exposed directly to the user, we allow the generation of samples in the following manner:
params = [5.0, 0.2, 2.25, 0.25]
MatCov = VecchiaMLE.generate_MatCov(n, params)
samples = VecchiaMLE.generate_Samples(MatCov, n, Number_of_Samples)100×400 Matrix{Float64}:
4.42603 2.17307 0.498556 -1.13273 … -0.0894268 -1.96642
4.5491 5.61946 6.83208 6.56174 -3.01765 -2.392
1.77223 0.161599 -2.02022 -3.21684 -1.8481 1.19667
-2.72575 -1.93194 0.171668 3.23915 3.87761 4.52429
-1.13836 -1.03013 -1.8493 -2.1122 -12.8495 -11.4842
8.97171 8.74703 9.72335 10.4759 … 1.47777 3.78649
-3.25559 -1.32086 -0.242984 -0.189268 -7.28878 -4.65416
4.05125 2.5488 0.557156 -0.425997 -5.96266 -7.07601
2.10653 -0.352877 -2.93825 -4.89028 0.517696 1.68152
3.22881 4.04061 3.88867 2.02811 -5.5524 -4.44156
⋮ ⋱
3.3266 3.33989 3.73334 2.96618 2.47901 2.88396
-4.09188 -1.79362 -0.212817 -0.293876 -2.31641 -1.81008
4.32659 2.41783 0.449138 -1.9389 0.0472859 -0.21798
2.76224 1.93997 0.959987 1.98891 -5.7332 -5.44557
2.00094 2.98639 2.89258 1.88634 … 5.56521 3.41733
4.00389 4.36053 2.91871 -0.438409 3.83293 0.636207
4.87534 4.42061 4.14321 2.95393 -4.3462 -4.94394
5.5964 7.63256 8.62329 6.57248 -4.79785 -6.14058
-4.53222 -0.870149 2.41098 5.65378 0.755557 4.05335Next, we create and fill in the VecchiaMLEInput struct. This is done below.
input = VecchiaMLE.VecchiaMLEInput(n, k, samples, Number_of_Samples, 5, 1)VecchiaMLEInput{Matrix{Float64}}(20, 10, [4.426032299453231 2.1730703121048442 … -0.0894268174816274 -1.9664176646874885; 4.549100003194665 5.619461884193633 … -3.0176500921474787 -2.3920046741307788; … ; 5.59639617697812 7.632555092782814 … -4.7978510366926574 -6.140577504916721; -4.532215782292826 -0.8701487424501329 … 0.7555574001333263 4.053349294068503], 100, MadNLP.ERROR, VecchiaMLE.cpu, [[0.0, 0.0], [0.10526315789473684, 0.0], [0.21052631578947367, 0.0], [0.3157894736842105, 0.0], [0.42105263157894735, 0.0], [0.5263157894736842, 0.0], [0.631578947368421, 0.0], [0.7368421052631579, 0.0], [0.8421052631578947, 0.0], [0.9473684210526315, 0.0] … [1.0526315789473684, 2.0], [1.1578947368421053, 2.0], [1.263157894736842, 2.0], [1.368421052631579, 2.0], [1.4736842105263157, 2.0], [1.5789473684210527, 2.0], [1.6842105263157894, 2.0], [1.7894736842105263, 2.0], [1.894736842105263, 2.0], [2.0, 2.0]])Here, we set (though not necessary) the optimizer (MadNLP) print level to 5 (ERROR), which keep it silent. Furthermore, we set the compute mode to 1 (cpu). The other option would be to set it to 2 (gpu), but since not all machines have a gpu connected, we opt for the sure-fire approach.
All that's left is to run the analysis. This is done in one line:
d, o = VecchiaMLE_Run(input)(VecchiaMLE.Diagnostics(2.63786402, 0.36384893599999996, 14.508672361, -1575.7179026138365, 4.072575098211863e-12, 3.473312452207548e-9, 5, VecchiaMLE.cpu), [1.7809436947855384 0.0 … 0.0 0.0; -2.2594566938128766 2.2110108254054204 … 0.0 0.0; … ; 0.0 0.0 … 0.5893620501942315 0.0; 0.0 0.0 … -0.5605425757910939 0.20794347253597226])We can see the structure of the output of the program has an approximately banded structure. This is due to the generation of the sparsity pattern depends on the euclidean distance between any given point and previously considered points in the grid.
sparse(o)400×400 SparseArrays.SparseMatrixCSC{Float64, Int64} with 3955 stored entries:
⎡⢷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠻⣯⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠻⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠁⠈⠙⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠁⠈⠙⣮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠁⠈⠻⣮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠁⠈⠻⣮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠁⠈⡻⢮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⡙⢦⡻⣦⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⡙⢆⠻⣦⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢀⠻⢆⠻⣦⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢀⠳⢆⠻⣦⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢀⡓⢦⡛⣦⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡙⢦⡹⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡹⣬⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⡻⣮⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⡱⣬⡻⣮⡳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢈⡛⢮⡛⣮⡳⣄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡙⢮⡻⣮⡳⣄⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠹⣮⠻⣮⠳⣄⎦The function VecchiaMLE_Run also returns some diagnostics that would be difficult otherwise to retrieve, and the result of the analysis. This is the cholesky factor to the approximate precision matrix. You can check the KL Divergence of this approximation, with a function inside VecchiaMLE, though this is not recommended for larger dimensions (n >= 30). This obviously requires the true covariance matrix, which should be generated here if not before.
println("Error: ", VecchiaMLE.KLDivergence(MatCov, o))Error: 54.606012298076706