path: root/report/report.tex

\documentclass{article}

\usepackage[margin=1in]{geometry}
\usepackage{parskip}
\usepackage{float}

\usepackage{xcolor}
\usepackage{graphicx}
\graphicspath{{./figures/} {./figures/tracy-widom-approx/}, {./figures/results/} {./figures/wall-times/}}
\usepackage{amsmath}
\usepackage{booktabs}
\usepackage{longtable}
\usepackage{caption}
\captionsetup{margin=0.5cm} % Sets a 1cm margin on both sides
\captionsetup{labelfont=bf}

\usepackage{hyperref}
\hypersetup{colorlinks=true,
            linkcolor=black,
            urlcolor=blue,
            citecolor=black,
            pdftitle={Reconstructing the Tracy-Widom Distribution},
            pdfpagemode=FullScreen}

\usepackage{listings}
\usepackage{algorithm}
\usepackage{algpseudocode}

\title{Reconstructing the Tracy-Widom Distribution\vspace{-0.75em}}
\author{Connor Moore, 100826701. \today{}}
\date{}

\begin{document}
\maketitle

\section{Introduction}
Random numbers form the basis of various studies in mathematics and the sciences. One strong example of this is \emph{Random Matrix Theory}, or the study of the construction and properties of random matrices. One of the first applications of random matrix theory was proposed by Eugene Wigner (1902-1995) in his application to model the nuclei and spectrum of various heavy atom \cite{anderson2009An_Int, wigner1955Charac}. It has since found different applications including work in computational mechanics, control theory, systems engineering, and others.

Random matrix theory borrows some distributions from statistics to construct matrices. This includes the Gaussian (normal) distribution, which is used to make various matrices depending on which number system is used. Wigner's work included Orthogonal (real symmetric), Unitary (Complex Hermitian), and Symplectic (Quaternion Hermitian) matrices and their eigenvalues. 

This work considers Gaussian orthogonal ensembles (GOEs) specifically and attempts to reconstruct the Tracy-Widom distribution directly by computing eigenvalues. An overview of relevant theory, the necessary methods, and implementation details is provided before results are presented. The distribution is examined as a function of matrix dimension and sample size, and statistical testing is performed to support a Tracy-Widom being observed as opposed to a normal distribution.

\section{Relevant Theory}

\subsection{Random Ensembles}
The GOE-based random matrix has some properties that go beyond simply being normally distributed. Namely, the matrix is square, symmetric, and the standard deviation $\sigma$ is unity except for elements on the diagonal where $\sigma=2$. The mean $\mu$ is considered as zero for all elements. Consider the normal distribution with notation $N(\mu, \sigma$). Then, the matrix is constructed using:
\begin{equation*}
    \mathbf{A}(i,j) \sim
    \begin{cases}
    N(0, 1) & \forall (i \neq j) \\
    N(0, 2) & \forall (i = j)
    \end{cases}
\end{equation*}
The diagonal standard deviation of 2 is what makes this ensemble orthogonal. Note that the representation above does not explicitly state the symmetry, but $\mathbf{A}(i,j) = \mathbf{A}(j,i)\, \forall(i,j)$.

\subsection{The Tracy-Widom Distributions}
Various distributions can be reconstructed from the eigenvalues of a GOE. For example, the case of the global eigenvalues yields a semicircle distribution \cite{anderson2009An_Int}. Consider the $N$ eigenvalues of a vector $|\lambda_1| \ge |\lambda_2| \ge |\lambda_3\| \ge ... \ge |\lambda_N|$. In this case, $\lambda_1$ is considered the \emph{dominant} eigenvalue so long as $|\lambda_1| > |\lambda_2|$. Taking the distribution of the dominant eigenvalues $\lambda_1$ yields an entirely different result that is known as the \emph{Tracy-Widom} distribution. This distribution can be parametrized in terms of which division algebra is used in its definition. Wigner's work considered real, complex, and Quaternion cases, shown below in Figure \ref{fig:tracy-widoms}:

\begin{figure}[H]
    \centering
    \input{figures/tracy-widom-approx/tracy_widom.tex}
    \caption{Various Tracy-Widom distributions for $\beta=1,2,4$. The different values of $\beta$ depend on which division algebra is used for the Gaussian ensemble (real/complex/quaternion). This work considers only the TW\textsubscript{1} distribution defined for real numbers, shown in red as the shallowest peak compared with the complex and quaternion versions.}
    \label{fig:tracy-widoms}
\end{figure}

Compared with a normal distribution, it is also difficult to distinguish between the two visually. Figure \ref{fig:tracy-widom-vs-normal} shows a comparison between a Tracy-Widom and a normal distribution:

\begin{figure}[H]
    \centering
    \input{figures/tracy-widom-approx/tracy_widom_compare.tex}
    \caption{The TW\textsubscript{1} distribution compared with a normal distribution of arbitrary mean and standard deviation. Note that while the peaks are centered, the Tracy-Widom has a `fatter' tail on the $+x$ side and a slimmer one on the $-x$ side. Both distributions are quite similar visually, so this asymmetry will play an important role in quantifying their differences.}
    \label{fig:tracy-widom-vs-normal}
\end{figure}
Note that the Tracy-Widoms plotted above have no closed-form expression and must be approximated, typically by using a scaled Gamma function \cite{chiani2012Distri}, or more cleverly by downloading some else's finished approximation off of the \textsc{Matlab} file exchange \cite{marco2013}. Because the distributions are for visual purposes only, a smart approach to validating the generation of a Tracy-Widom should include some type of quantifiable, statistical assurance not based on the simple value of the distribution itself.

\subsection{Statistical Differentiation}
There are various ways to go about discerning a Tracy-Widom distribution from a normal one. Perhaps the most elegant is to exploit some basic differences in properties between the two. As mentioned earlier, the Tracy-Widom has a fatter side and a slimmer side, which makes it asymmetric. This is not the case with the normal distribution which is identical on either side of the mean value $\mu$. It turns out that the symmetry of a distribution is readily quantifiable using existing statistical tools, and as such that will be the backbone of the statistical analysis. First consider the \emph{expected value} of a probability density function (PDF) $p(x)$, defined as \cite{siegrist2021}:
\begin{equation}
    E[x] = \int_s x p(x)\, dx \approx \left(\sum_i x_i\,p(i)\right) \Delta x
    \label{eqn:expected}
\end{equation}
where the approximation applies for a discrete PDF with equi-distance bins along $x$. The expected value has \emph{moments} that give quantifiable insights into the behaviour and shape of a PDF. The general equation for the $n$th central moment of $E(x)$ is given as:
\begin{equation}
    E\left[(x-a_n)^n\right]
\end{equation}
where $a_n$ is some constant. The first moment (taken as a so-called \emph{raw moment} with $a=0$) simply reduces to \ref{eqn:expected}. The second moment is the variance of the function and is taken with $a=\mu$:
\begin{equation}
    \sigma^2 = E\left[(x - E[x])^2\right] = E\left[ (x - \mu)^2\right]
\end{equation}
The next two moments are of particular value for this project as they quantify the symmetry of the distribution. The third moment is the \emph{skewness} of the distribution, which deviates from the central moments above and is normalized by the standard deviation $\sigma$:
\begin{equation}
    \mathrm{Skewness} = E\left[\frac{(x-\mu)^3}{\sigma}\right]
\end{equation}
Lastly, the 4th moment is the kurtosis of the distribution. Kurtosis comes from the Greek word for \emph{bulging}, and as a parameter it quantifies the `sharpness' of the PDF. A larger value of kurtosis implies that the distribution has fatter tails \cite{siegrist2021}. Note that kurtosis is always a positive value and is equal to 3 for a standard normal distribution. Therefore, the  \emph{excess} kurtosis is considered by subtracting 3. The relation is given as:
\begin{equation}
    \mathrm{Excess\ Kurtosis} = E\left[ \frac{(x - \mu)^4}{\sigma} \right] - 3
\end{equation}
These tests, while perhaps simple, will do a fine job in the strict context of differentiating a Tracy-Widom distribution from a normal distribution. In the case of the normal distribution, the skewness and excess kurtosis are both zero. For comparison, the approximate Tracy-Widom has a skewness of $\approx 0.2935$ and an excess kurtosis of $\approx 0.1292$.This was calculated using an Octave script\footnote{\href{https://git.hhmoore.ca/mcsc-6030g/p3-tracy-widom/tree/report/figures/tracy-widom-approx/plot_tracy_widom.m}{Available online} as part of the Git repository for the project.} that sampled 50,000 bins in the range $[-9,9]$ for both distributions. Performing calculations with (1,3,4,5) reconstructed the mean and standard deviation for the normal distribution sampled from its the closed-form equations down to an error of $\sim10^{-5}$.

\section{Numerical Methods}
The natural extension of defining the prerequisite statistical tests to confirm a Tracy-Widom is to actually sample a Tracy-Widom. The issue with this, however, is that it can be quite computationally demanding. The Tracy-Widom is quite similar to the normal except for near the very tails, meaning that most samples (occurring near the peaks) will show agreement between the two. Additionally, it is theorized that the dimension of the matrix $N$ has to be on the order of thousands to get a very Tracy-Widom distribution. The combination of these two points motivates finding a convenient way to calculate the maximum eigenvalue of a large matrix repeated such that thousands of calculations can be performed in a meaningfully short amount of time.
\subsection{The Power Method}
One relatively cheap way to find the maximum eigenvalue is through the so-called power method. The power method is an iterative approach to finding the maximum eigenvalue of a matrix that exploits the spectral gap ($\lambda_1 - \lambda_2$). Consider an initial guess in the form of a vector of length $N$ given as $\mathbf{x_0}$. Then, perform the multiplication:
\begin{equation}
    \mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k
\end{equation}
After multiplying, then normalize the vector such that $\left| \mathbf{x}_{k+1} \right|$ is unity:
\begin{equation}
    \mathbf{x}_{k+1} = \frac{\mathbf{x}_{k+1}}{\left| \mathbf{x}_{k+1} \right|}
\end{equation}
Or, more compactly written as one step:
\begin{equation}
    \mathbf{x}_{k+1} = \frac{\mathbf{A}\mathbf{x}_{k}}{\left|\mathbf{A}\mathbf{x}_{k}\right|}
\end{equation}
Repeating this will converge to the dominant eigenvalue $\lambda_1$ so long as $\mathbf{x}_0$ is not orthogonal to the dominant eigenvector \cite{novak2022}. The convergence of the method is $\mathcal{O}\left( \left| \lambda_1 / \lambda_2 \right|^k \right)$ and is largely dependent on the spectral gap. The main advantage of this method is that computationally it consists of only a matrix multiplication ($\mathcal{O}(n^2)$) and normalization ($\mathcal{O}(n)$) which is rather inexpensive compared with other methods per iteration \cite{watkins2002}. Worth noting, however, is that in cases where $\mathbf{A}$ has a small spectral gap it is not guaranteed to converge in an acceptable number of iterations. This poses an issue as it would bias the calculations for reconstructing the Tracy-Widom by not including random ensembles specifically with lower spectral gaps. To ratify this, more expensive methods exist which guarantee a solution that are used as a fallback.
\subsection{Schur Decomposition}
One method that guarantees a solution $\forall \mathbf{A}$ is Schur Decomposition. Named after the late Russian mathematician Issai Schur (1875-1941), Schur Decomposition involves decomposing the matrix $\mathbf{A}$ into a combination of some unitary matrix $\mathbf{U}$ and upper-triangular matrix $\mathbf{T}$ \cite{watkins2002}:
\begin{equation}
    \mathbf{A} = \mathbf{U} \mathbf{T} \mathbf{U}^\mathrm{T}
\end{equation}
In this case, the diagonal entries in $\mathbf{T}$ are then the eigenvalues of the original matrix $\mathbf{A}$. The proof in Schur's theorem is non-constructive and does not provide a solution without knowing the eigenvectors of $\mathbf{A}$, but it can be still be used in combination with the so-called QR algorithm developed independently by John Francis (1934-present) and Vera Kublanovskaya (1920-2012) starting from 1959 \cite{golub2009The_QR}. Essentially, the QR algorithm is the computational scheme that provides the matrices $\mathbf{U}$ and $\mathbf{T}$ from the input $\mathbf{A}$. It does so with a complexity of $\mathcal{O}(n^3)$, which for larger values of $N$ deviates significantly from the $\mathcal{O}(n^2)$ of the power method. For this reason it is used only when the power method fails.
\section{Implementation and MPI Parallelization}
The last remaining hurdle to computing the eigenvalues is the magnitude of memory required to store such large matrices. Because the requirements will be in the gigabytes range for most calculations, shared-memory parallelism can have issues supporting large-scale calculations needed to accurately define the Tracy-Widom. Because of this, the solution is designed around MPI parallelism with distributed and independent memory between ranks. There are fewer opportunities for race conditions and threads `overstepping' each other in MPI, however care must be taken that the now independent workers can communicate with each other properly.

\subsection{Worker-Manager Structure}
The program used in this project employed a worker-manager structure to drive the eigenvalue calculations. The first worker is assigned to be a manager and handle the seeding and results collection for the other workers. Everyone else is focused solely on calculating eigenvalues and reporting their results to the manager when finished. This setup allows for the independent calculation of as many eigenvalues as MPI ranks are supported at once, which can approach the order of hundreds or even thousands if compute resources are available on a respectable cluster.

\subsection{The Manager Process}
The manager process is relatively straightforward and focuses more on the `book keeping' side of the eigenvalue calculations. To demonstrate this, a pseudocode explaining the behaviour of the manager MPI rank is presented below in Algorithm \ref{alg:manager}:

\begin{algorithm}[H]
\caption{MPI Manager Process}
\begin{algorithmic}[1] % [1] enables line numbers
    \State N,L $\gets$ Matrix size and number of data points
    \State Open logfile and write $N$, $L$ for reference
    \While{Workers not finished}
        \For{$i=1:N_{workers}$}
            \If{Worker $i$ is ready and needed}
                \State seed $\gets$ \texttt{/dev/urandom}
                \State \texttt{MPI Broadcast}: send the worker seed
            \EndIf
        \EndFor
        \If{All workers are finished}
            \State Break loop and finish managing
        \EndIf
        \State \texttt{MPI Receive}: wait for a worker to finish calculating...
        \State Identify which worker is finished calculating and their result
        \If{Eigenvalue was calculated successfully}
            \State Tally to the number of successful eigenvalue calculations
            \State Write the worker number and eigenvalue to logfile
        \Else
            \State Tally to the number of failed eigenvalue calculations
            \State Print a warning message
        \EndIf
        \If{Worker is still needed}
            \State Set worker flag to `Ready'
        \Else
            \State Tell worker to conclude
        \EndIf
    \EndWhile
\end{algorithmic}
\label{alg:manager}
\end{algorithm}

Note that the actions labelled with involving MPI communication are blocking, i.e., the worker/manager will wait to receive a broadcast before continuing at all.

\subsection{The Calculation (Worker) Process}
The worker process involves calculating the eigenvalues and is more complicated than the manager process. Various sub-processes are included in different algorithms, namely for building the matrix, computing the eigenvalue, and applying the power method. The overview of the whole process is outlined in Algorithm \ref{alg:worker}:

\begin{algorithm}[H]
\caption{MPI Worker Process}
\begin{algorithmic}[1] % [1] enables line numbers
    \State size $\gets$ \texttt{MPI Receive}: matrix size
    \State Initialize matrix space in memory
    \While{Worker not Finished}
        \State seed $\gets$ \texttt{MPI Receive}: new matrix seed
        \State Call \texttt{open\_prng\_seed}(seed)
        \State matrix $\gets$ \texttt{build\_random\_matrix}(size)
        \State eig $\gets$ \texttt{calculate\_eigenvalue}(matrix)
        \State \texttt{MPI Send}: eig and relevant information
    \EndWhile
\end{algorithmic}
\label{alg:worker}
\end{algorithm}
The algorithm for assembling the random matrix given a size is provided below in Algorithm \ref{alg:random-matrix}:
\begin{algorithm}[H]
\caption{MPI Worker: Random Matrix Building}
\begin{algorithmic}[1] % [1] enables line numbers
\Procedure{build\_random\_matrix}{size}
\For{$i = 1:$ size}
    \State \textbf{A}(i,i) = N(0,2)
    \For{$j = 1:i-1$}
        \State \textbf{A}(i,j) = N(0,1)
        \State \textbf{A}(j,i) = \textbf{A}(i,j)
    \EndFor
\EndFor
\State \Return \textbf{A}
\EndProcedure
\end{algorithmic}
\label{alg:random-matrix}
\end{algorithm}

And likewise for the eigenvalue calculation in Algorithm \ref{alg:eigenvalue-calculation}:

\begin{algorithm}[H]
\caption{MPI Worker: Eigenvalue Calculation}
\begin{algorithmic}[1] % [1] enables line numbers
\Procedure{calculate\_eigenvalue}{matrix}
    \State Initialize convergence tolerance $\varepsilon$ and maximum iterations for power iteration
    \State Call \texttt{power\_iteration}(matrix, shift=0)
    \If{power iteration converged}
        \If{$\lambda < 0$}
            \State Call \texttt{power\_iteration}(matrix, shift=$\lambda$)
        \EndIf
    \Else
        \State call \texttt{MKL\_schur\_decomposition}(matrix)
        \If{decomposition successful}
            \State eig $\gets$ Largest eigenvalue from decomposition
            \State Set `converged' flag to `true'
        \Else
            \State Print warning message \Comment{Not touching the converged flag leaves it `false'}
        \EndIf
    \EndIf
\EndProcedure
\end{algorithmic}
\label{alg:eigenvalue-calculation}
\end{algorithm}

So far the structure has involved a minimal amount of math. This is only because it has been hidden behind abstraction and by using the MKL library to call pre-made functions for computing eigenvalues. The power iteration method was implemented by hand and involves more math as shown in Algorithm \ref{alg:power-iteration}:

\begin{algorithm}[H]
\caption{MPI Worker: Power Iteration}
\begin{algorithmic}[1] % [1] enables line numbers
\Procedure{power\_iteration}{matrix, shift=0}
\State Initialize a random vector $\mathbf{x}_o$ = N(0,1)
\State $\mathbf{x}_0 \gets \frac{\mathbf{x}_0}{\left|\mathbf{x}_0\right|}$ \Comment{Normalize}
\State mem $\gets$ $-100$

\For{$i = 1:$ max\_iterations}
    \State $\mathbf{y} \gets \mathbf{A}\mathbf{x}_{i-1}$
    \State $\mathbf{y} \gets \mathbf{y} + $shift$\times\mathbf{x}_{i-1}$
    \State eig $\gets$ \texttt{dot\_product}($\mathbf{x}_{i}$,$\mathbf{x}_{i-1}$) - shift
    \State error $\gets\ \frac{\left| \mathrm{mem} - \mathrm{eig}\right|}{\left| \mathrm{mem} \right|}$
    \If{err $< \varepsilon$}
        \State Set `converged' flag to `true'
        \State \Return{eig}
    \Else
        \State $\mathbf{x}_i \gets \frac{\mathbf{x}_i}{\left| \mathbf{x}_i \right|}$
        \State mem $\gets$ eig
    \EndIf
\EndFor

\EndProcedure
\end{algorithmic}
\label{alg:power-iteration}
\end{algorithm}

These algorithms compose the entirety of the program for reconstructing the Tracy-Widom distributions discussed earlier. It is worth noting that the expected complexity is a function of the required calculation complexity $n$, the number of MPI ranks used $m$, and the sample size $l$. Assuming a total FLOPs count of $\sim\mathcal{O}(n^3)$ in a worst-case scenario for Schur decomposition, repeated $l$ times, in reality each rank will only have to contribute 1-$m$th of the work. Note that while the QR algorithm is not done in parallel, repeating in $l$ times is embarrassingly parallel. So the total time complexity should be roughly:
\begin{equation}
    \sim\mathcal{O}\left(\frac{n^3\, l}{m}\right)
\end{equation}
which will be tested in the next section.

\section{Results}
Before the eigenvalues can be directly compared to a Tracy-Widom they must be scaled to match the existing distribution. It is possible to perform the scaling empirically, but a closed-form solution is available that makes things easier \cite{konig2005}:
\begin{equation}
    \lambda^\prime = \left(\lambda - 2\sqrt{N} \right)\times N^{1/6}
\end{equation}
This specific transform is derived from the Wigner semicircular distribution mentioned earlier in section 2. Although it describes the global distribution, taking the extrema yields a bound for the maximum eigenvalue.

\subsection{Tracy-Widom Convergence}
Following some software issues related to the OneAPI suite, the author was unable to produce enough samples in time to complete this report. However, a moment is appropriate to thank all of the teammates for being kind enough to share their data for analysis in this report. The following numbers are from Joe's runs but the analysis was done individually. Samples here are presented with respect to their matrix dimension $N$, number of samples $L$, skewness $S$ and excess kurtosis $K$. Taking, for example, different sample sizes for $N=1750,500$ yields a half-dozen examples of a clear trend emerging where the histograms tend towards a normal-\emph{ish} shape, shown in Figure \ref{fig:randoms}:

\begin{figure}[H]
    \centering
    \input{./figures/results/size_pts.tex}
    \caption{Progression of the Tracy-Widom convergence as a function of matrix size ($N$) and sample size ($L$). Also included are the skewness ($S$) and excess kurtosis ($K$) for each sample. An increase in both sample size and matrix size yields a distribution comparable with but not necessarily closer to the Tracy-Widom reference. Asymmetry is noted especially for the cases with $N=1750$, but with an incorrect sign that implies a fatter tail in the $-x$ direction.}
    \label{fig:randoms}
\end{figure}

Despite the visuals, however, the metrics tell a bleaker story. Recall that the Tracy-Widom is authentically constructed with a fatter tail on the right and thinner on the left, which seems to be the opposite of what is presented above. This is also reflected in the skewness assuming a negative value for the curves associated with $N=1750$, implying that the distribution is more favourable to the left of the mean value. As the sample size increases, the shape of the distribution becomes more refined and tends to that of the expected profile, which seems intuitive. It is possible to consider the maximum number of samples available for each matrix size $N$ to confirm whether or not this can yield better results, shown in Figure \ref{fig:big-ones}:


\begin{figure}[H]
    \centering
    \input{./figures/results/maxcomp.tex}
    \caption{Reconstructed distributions from the maximum number of samples $L_{max}$ for each size $N$. Table \ref{tbl:results} provides more information on the behaviour of skewness and excess kurtosis for the high-sample runs. Qualitatively, under filling seems prevalent at larger $N$ and a strong $-x$ tail is visible.}
    \label{fig:big-ones}
\end{figure}
The parameters calculated for the distributions are presented in Table \ref{tbl:results} for convenience:

\begin{table}[H]
    \centering
    \caption{Parameters of the Larger-$L$ Distributions for Various $N$}
    \begin{tabular}{cccc}
        \toprule
        N & L & Skewness & Excess Kurtosis \\
        \midrule
        500 & 50,104 & $-8.993\times10^{-3}$ & $8.588\times10^{-1}$ \\
        750 & 55,397 & $-5.745\times10^{-2}$ & $7.666\times10^{-1}$ \\
        1,000 & 55,536 & $-9.220\times10^{-2}$ & $6.924\times10^{-1}$ \\ 
        1,250 & 55,726 & $-1.185\times10^{-1}$ & $6.616\times10^{-1}$ \\
        1,500 & 55,948 & $-1.188\times10^{-1}$ & $4.567\times10^{-1}$ \\
        1,750 & 56,097 & $-1.471\times10^{-1}$ & $4.227\times10^{-1}$ \\
        \bottomrule
    \end{tabular}
    \label{tbl:results}
\end{table}

Despite the significant increase in sample size ($\sim10\times$ at a minimum), the skewness and excess kurtosis do not converge to their expected values. Rather, the skewness diverges from the expected value of 0.2935, and the excess kurtosis has a minimum error of at least 200\% relative to the expected 0.1292 described earlier. These trends are shown below in Figure \ref{fig:errors}:

\begin{figure}[H]
    \centering
    \input{./figures/results/parameters.tex}
    \caption{Behaviour of the two parameters skewness $S$ and excess kurtosis $K$ as function of matrix size $N$ for a maximal number of samples $L$. Observed responses are shown in solid lines, with the expected values being horizontal dotted lines and the percent-error represented with dashed lines. Neither of the trends observed in the skewness or expected kurtosis are consistent with theory, and their behaviours cannot confirm the existence of a Tracy-Widom distribution from the calculations.}
    \label{fig:errors}
\end{figure}

The tendency of the skewness to diverge from zero as the matrix size increases does point to a non-Gaussian distribution being observed, but it cannot be concluded that it is a Tracy-Widom.

\subsection{Wall Times}

Lastly, the wall time of running the simulations was investigated. These results were more aligned with the expected values of complexity. Two different perspectives are considered: the variation of sample size with a constant matrix size, and the vice-versa case. The results for the constant sample sizes is discussed first as it should correspond to the complexity mentioned in Section 4. Figure \ref{fig:N_walltimes} shows a complexity tending to $\mathcal{O}(n^3)$ (with slight variations over the matrix size domain) for all sample sizes:


\begin{figure}[H]
    \centering
    \input{./figures/wall-times/N_walltimes.tex}
    \caption{Wall times for different sample sizes $L$ as a function of the matrix size $N$. All runs were conducted using 128 MPI ranks for this figure. Wall times show a roughly $\mathcal{O}(n^3)$ scaling consistent between all sample sizes tested, although there is a magnitude difference between the fastest ($L=128$) and slowest ($L=4096$) runs.}
    \label{fig:N_walltimes}
\end{figure}

The vice-versa case shows a different story, with a complexity closer to $\mathcal{O}(n)$. This makes intuitive sense, as MPI ranks can calculate different matrices of a common sample in parallel, explaining the rather weak coupling. Compare this to increasing the matrix size itself, which is solved in serial by one rank at time, making the scaling appear much more sensitive to changes in the independent variable. The trend is visualized below in Figure \ref{fig:walltimes}:

\begin{figure}[H]
    \centering
    \input{./figures/wall-times/walltimes.tex}
    \caption{Wall times for different matrix sizes $N$ as a function of the number of sampled eigenvalues $L$. All runs were conducted using 128 MPI ranks for this figure. Wall times show a roughly $\mathcal{O}(n)$ scaling consistent between all matrix sizes tested, although there was a $\sim2$ magnitude difference between the fastest ($N=1$) and slowest ($N=1750$) runs. Small deviations exist for low sample sizes.}
    \label{fig:walltimes}
\end{figure}

\section{Conclusions}
This work set out to reconstruct a Tracy-Widom distribution via the direct sampling of calculated eigenvalues of Gaussian orthogonal ensembles. A code leveraging MPI for distributed-memory parallelism and the MKL BLAS library for high-performance numerical methods was used to sample various numbers of ensembles, sizes of matrices, and various other factors to generate data and draw conclusions. While the worker-manager structure worked well for the batch computing of eigenvalues, the statistical results proved inconclusive. It can be said with confidence that the resulting distribution was non-Gaussian but it was not necessarily a Tracy-Widom distribution. This is because of the statistical testing, consisting of the skewness and excess kurtosis, were not in support of a Tracy-Widom distribution. The discrepancies showed at least one trend divergent to the expected value and therefore it is unlikely that further sampling of eigenvalues would mitigate this issue. Rather, it is likely that there is an inconsistency related to transforming the eigenvalues to fit the standard Tracy-Widom. The wall time scaling, however, was consistent with expectations and was invested with respect to both matrix and ensemble size.

\bibliographystyle{ieeetr}
\bibliography{refs.bib}

\end{document}