Demo 1: One-parameter scaling

In this demo we will show how to perform finite-size scaling where one free parameter is tuned to obtain the best possible data collapse. We will use the CDW transition in the square lattice Holstein model as an example, where the critical inverse temperature $\beta_c$ is to be optimized.

First, import the FiniteSizeScaling.jl package:

julia> using FiniteSizeScaling

Before we can perform finite-size scaling, we need to have data for several difference lattice sizes $L$. Suppose we have five lattice sizes $L=4, 6, 8, 10, 12$. For each lattice size, we will have an array of $X$ values (i.e. values of the inverse temperature $\beta$), an array of $Y$ values (i.e. values of the order parameter $S_{cdw}$), and an array of $E$ values which are the errors in each $Y$ value. The data should be pre-arranged into arrays like this:

X_L4 = [1.973, 2.989, 3.978, 4.513, 4.754, 4.968, 5.182, 5.476, 5.743, 5.957, 6.225, 6.492, 7.000, 7.989, 8.979]
Y_L4 = [1.250, 1.806, 2.222, 3.056, 3.333, 3.750, 4.306, 4.583, 4.861, 5.278, 5.694, 5.972, 6.250, 6.806, 7.083]
E_L4 = [0.114, 0.123, 0.132, 0.128, 0.135, 0.121, 0.130, 0.101, 0.108, 0.139, 0.133, 0.138, 0.128, 0.199, 0.137]

X_L6 = [2.000, 2.989, 4.005, 4.513, 4.754, 4.995, 5.262, 5.476, 5.743, 6.011, 6.492, 7.000, 7.455, 7.989, 9.005, 9.995, 11.010, 12.000]
Y_L6 = [1.245, 1.667, 2.778, 4.444, 5.694, 6.667, 8.056, 9.306, 10.139, 10.972, 12.222, 13.333, 13.750, 14.167, 14.444, 14.583, 14.722, 14.722]
E_L6 = [0.164, 0.153, 0.152, 0.148, 0.165, 0.151, 0.150, 0.161, 0.138, 0.169, 0.153, 0.148, 0.138, 0.149, 0.167, 0.155, 0.147, 0.160]

X_L8 = [4.487, 4.727, 4.968, 5.235, 5.503, 5.743, 6.011, 6.225, 6.492, 6.733, 6.973, 7.267, 7.481, 7.775, 8.043, 8.524, 8.979, 10.021, 10.983, 12.000]
Y_L8 = [4.861, 5.694, 7.639, 10.278, 13.472, 15.417, 17.917, 19.583, 20.694, 21.528, 22.083, 22.361, 22.778, 23.194, 23.612, 23.889, 24.028, 24.583, 24.722, 24.722]
E_L8 = [0.252, 0.289, 0.243, 0.324, 0.489, 0.524, 0.578, 0.609, 0.575, 0.589, 0.625, 0.621, 0.650, 0.623, 0.621, 0.602, 0.658, 0.665, 0.667, 0.689]

X_L10 = [4.489, 4.757, 5.000, 5.244, 5.487, 5.755, 5.971, 6.293, 6.509, 6.750, 7.018, 7.259, 7.500, 7.768, 8.035, 10.041, 12.020]
Y_L10 = [4.583, 6.111, 9.306, 15.556, 19.861, 21.806, 25.833, 27.639, 31.111, 32.361, 33.333, 33.750, 34.722, 35.972, 35.694, 37.222, 37.500]
E_L10 = [0.254, 0.289, 0.501, 1.189, 1.184, 0.790, 0.804, 1.045, 0.600, 0.584, 0.598, 0.655, 0.680, 0.640, 0.601, 0.575, 0.478]

X_L12 = [4.493, 4.734, 4.981, 5.280, 5.528, 5.637, 5.777, 6.023, 6.268, 6.539, 6.809, 7.024, 7.292]
Y_L12 = [5.432, 6.964, 11.838, 16.713, 23.120, 25.627, 31.058, 36.072, 40.669, 43.454, 45.822, 47.075, 48.050]
E_L12 = [0.321, 0.398, 0.401, 0.556, 0.531, 0.598, 0.601, 0.592, 0.625, 0.630, 0.587, 0.658, 0.688]

The next step is to create a single array of data, where each element is an array of $[X, Y, E, L]$ data for a given lattice size. The length of this new array will equal the number of different lattice sizes being used:

data_with_error = [[X_L4, Y_L4, E_L4, 4], [X_L6, Y_L6, E_L6, 6], [X_L8, Y_L8, E_L8, 8], [X_L10, Y_L10, E_L10, 10], [X_L12, Y_L12, E_L12, 12]]

If your $Y$ data does not have error bars, you can also omit the arrays of error values, for example:

data_no_error = [[X_L4, Y_L4, 4], [X_L6, Y_L6, 6], [X_L8, Y_L8, 8], [X_L10, Y_L10, 10], [X_L12, Y_L12, 12]]

At this point, you can call the function plot_data to display the raw unscaled data:

julia> plot_data(data_with_error)

As discussed in the Example page, we aim to rescale the data along new axes $X_s = (\beta - \beta_c) L$ and $Y_s = S_{cdw} L^{-7/4}$. We will take $\beta_c$ to be the one free parameter we will tune to obtain the optimal data collapse, and denote it $v_1$, i.e. our scaled axes will be $X_s = (X - v_1) L$ and $Y_s = Y L^{-7/4}$. For one-parameter scaling, the user now must define two functions which define the scaled X and Y axes:

x_scaled(X, L, v1) = (X .- v1) * L

y_scaled(Y, L, v1) = Y * (L^(-7/4))

The function defining the $X_s$ axis should take $X$, $L$, and $v_1$ as arguments, while the function defining the scaled $Y_s$ axis should take $Y$, $L$, and $v_1$ as arguments. Note that since $X$ and $Y$ are arrays, elementwise operators such as $.-$ may be necessary.

To determine the best value of $v_1$, a search is performed between an initial value $v_{1i}$ and a final value $v_{1f}$. For each candidate value of $v_1$, the data is scaled according to the user-defined formulae above, and a polynomial curve of order $p$ is fit to the data, using the method of least-squares to minimize the sum of squared residuals between the scaled data and the polynomial curve. The value of $v_1$ which yields the smallest sum of squared residuals is identified as the optimal value.

When performing the fit, one has the option to multiply each squared residual by a weight $W$. For example, to perform weighted least-squares, these weights are typically inverse variances. This means that data points with larger error bars will be deemed less important for the fit, and vice-versa. To use weights when performing the fit, define an array as follows:

fit_weights = [1 ./ (E_L4.^2), 1 ./ (E_L6.^2), 1 ./ (E_L8.^2), 1 ./ (E_L10.^2), 1 ./ (E_L12.^2)]

where each element is an array of $W$ values for a given lattice size. Here inverse variances have been chosen to perform a typical weighted least-squares. Note that the length of this array will equal the number of different lattice sizes being used. Using weights is completely optional and can be omitted if desired when calling the function fss_one_var.

The next step is to call the function fss_one_var to perform the finite-size scaling.

julia> scaled_data, residuals, min_res, best_v1 = fss_one_var(data=data_with_error, xs=x_scaled, ys=y_scaled, v1i=5.0, v1f=7.0, n1=100, p=4, weights=fit_weights);

where data is the single array of data defined previously, xs and ys are the functions previously defined for the scaled axes, v1i and v1f are the start and end points of the parameter search, and n1 is the number of values of $v_1$ in this range to check. The integer degree $p$ of the polynomial must also be specified, typically $4 \leq p \leq 8$ is sufficient.

The function fss_one_var returns four variables: an array scaled_data_array where each element is an array of $[X_s, Y_s, E_s, L]$ data for a given lattice size; an array residuals of length n1 which stores the sum of squared residuals for each value of $v_1$ checked; a scalar min_res which is the minimum value of the array residuals; and a scalar best_v1 which is the value of $v_1$ which gave the smallest overall residual. By default it will also print out the values of best_v1 and min_res:

Optimal v1 value: 6.1313131313131315 
Smallest residual: 0.12847365603386035 

At this point, you can again call the function plot_data, passing in the scaled_data_array returned by the fss_one_var function. This will produce a plot of the optimal data collapse, i.e. the scaled data with $v_1$ set to best_v1:

julia> plot_data(scaled_data)

After one-parameter finite-size scaling has been performed, we can plot the dependence of the sum of squared residuals on $v_1$. To do this use the plot_residuals function, and pass in the array of residuals returned by fss_one_var, along with the exact values of v1i, v1f, and n1 which were used:

julia> plot_residuals(residuals, v1i=5.0, v1f=7.0, n1=100)