# Deep Linear Discriminant Analysis with Variation for Polycystic Ovary
Syndrome Classification

Raunak Joshi
Department of IT
University of Mumbai
Mumbai, 400032, India
<EMAIL_ADDRESS>
&Abhishek Gupta
Department of EXTC
University of Mumbai
Mumbai, 400032, India
<EMAIL_ADDRESS>
&Himanshu Soni
Department of EXTC
SJCEM
Palghar, 401404, India
<EMAIL_ADDRESS>
&Ronald Laban
Department of EXTC
SJCEM
Palghar, 401404, India
<EMAIL_ADDRESS>

###### Abstract

The polycystic ovary syndrome diagnosis is a problem that can be leveraged
using prognostication based learning procedures. Many implementations of PCOS
can be seen with Machine Learning but the algorithms have certain limitations
in utilizing the processing power graphical processing units. The simple
machine learning algorithms can be improved with advanced frameworks using
Deep Learning. The Linear Discriminant Analysis is a linear dimensionality
reduction algorithm for classification that can be boosted in terms of
performance using deep learning with Deep LDA, a transformed version of the
traditional LDA. In this result oriented paper we present the Deep LDA
implementation with a variation for prognostication of PCOS.

_K_ eywords Deep Learning $\cdot$ Deep LDA $\cdot$ Linear Discriminant
Analysis

##  1 Introduction

The use of medical research data for various statistical tasks has been done
from a prolonged period of time. The data after having consistent number of
records can be utilized for deriving inference using prognostication methods.
The methods that fall under the area of inferential statistics [1] which are
extended with applied areas of statistics can be used, one of which can be
used is Machine Learning [2]. Task of prognostication based on data can be
done by learning patterns from the data using machine learning. The dataset
that we have used in this paper is pertaining to polycystic ovary syndrome [3]
diagnosis, which falls under the classification [4] category. Classification
has 2 sub-divisions, viz. Binary and Multi-Class where the data used in this
paper falls under binary classification [5] precisely. The methods using
Machine Learning have already been performed on polycystic ovary syndrome
abbreviated as PCOS using logistic regression [6], bagging ensemble methods
[7], discriminant analysis [8], stacked generalization [9], boosting ensemble
methods [10] and deep neural networks [11]. The use of deep learning is
evident and we want to focus on the more variations that can be brought into
the current state-of-the-art system. Deep Learning [12, 13] provides more
depth of learning as compared to machine learning and is used when the amount
of parameters in dimensions are high. The PCOS dataset has over 41 dimensions
which are enough for instating the use of deep learning. The variation we
wanted to perform was related to some machine learning algorithm that can be
leveraged with power of deep learning. The implementation of any machine
learning algorithm using a library like scikit-learn [14, 15] meets
limitations in terms of utilization with GPU processing power. For the same
reason, using a mature framework like Tensorflow [16] can definitely bring
change to the working. This is where we decided to work with parametric
learning method which works with simple and definite procedures. The
parametric learning method we focused on using was discriminant analysis [17,
18] which has variations in it where we focused on Linear Discriminant
Analysis [19]. This actually accounted for an idea that training the linear
discriminant analysis with deep learning style will yield us Deep Linear
Discriminant Analysis [20, 21] which has been already been discovered and we
decided to proceed with out implementation using it. The variations that we
brought in the network will be explained in further sections of this paper.

##  2 Methodology

This section of the paper gives detailed insights about the implementation and
approach we have taken to solve the problem. The model considering the
implementation with respect to Deep LDA [20, 21] revolves around the idea of
the convolutional neural networks [22, 23, 24, 25]. The modification can be
done to work with numerical values. This has been implemented by various
developers and names of the developers are given in the acknowledgement
section of the paper. The Deep LDA is basically an implementation of latent
representations in linearly separable method. The Deep LDA is an extensive
implementation of the traditional Linear Discriminant Analysis which was
intended for dimensionality reduction based classification methods. The
implementation consists of 2 phases, first phase consists of linear
discriminator as the deep neural network and second phase consists of support
vector machine for detailed classification.

###  2.1 First Phase

Figure 1: First Phase with LDA Implementation

The figure 1 gives depiction of first phase of the implementation. The input
layer takes 41 dimension of features from the data. This is passed on to one
dense layer that has 1024 hidden neurons. The L2 regularization [26] is
applied as a kernel regularization for the layer. The activation function is
used sigmoid [27]. The rectified linear unit abbreviated as ReLU [28] is the
activation function commonly used for deep learning methods but using sigmoid
ensures linear based system such as developed for linear discriminant
analysis. The parameters learned from the first hidden layer are 43,008.
Similar type of hidden layers are repeated twice, where second hidden layer
learns 1,049,600 parameters and third hidden layer also learns same amount of
parameters. The output layer consists of 1 hidden neuron with sigmoid
activation function and learns 1025 parameters. The network learns total of
2,143,233 parameters where all the parameters are trainable. The loss function
used is binary cross-entropy that differs from the original implementation of
deep LDA paper. The loss optimizer used is Adam [29] optimizer with learning
rate of $1*10^{-5}$ that roughly denotes 0.00001. The implementation is done
using Keras [30] over Tensorflow [16] back-end trained for 100 epochs with 64
as batch size.

###  2.2 Second Phase

Figure 2: Second Phase with SVM Implementation

The figure 2 depicts the second phase implementation. This is done using
Support Vector Machine [31] implementation with neural network inclination.
The input layer is connected to hidden layer with 100 hidden neurons with ReLU
[28] activation function. This layer learns 200 parameters. This layer is
connected to dropout [32] layer with 50% threshold and does not learn any
parameters. The output layer has 1 hidden neuron and has sigmoid activation
function for binary classification and learns 101 parameters. The total
parameters learned are 301 and the network uses binary cross-entropy. The loss
optimizer used is Adam [29] optimizer with $1*10^{-5}$ that approximately
denotes 0.00001 learning rate. The network is trained with 100 epochs and 64
as batch size.

###  2.3 Complete Network

Figure 3: Complete Network

The complete network is accumulation of first and second phase where the
output of the first phase is the input for second phase. The output of the
first phase is 1-dimensional array from 41-dimensional output. This is given
as an input to the second phase of the network and final prediction which is
1-dimensional is achieved. The both phases are trained independently and the
output is retained from first phase and given as second phase. The results of
the network will be given in succeeding section of the paper.

##  3 Results

###  3.1 Accuracy and Loss for First Phase

$0$$20$$40$$60$$80$$0.6$$0.7$$0.8$$0.9$EpochsAccuraciesTraining AccuracyValidation Accuracy | $0$$20$$40$$60$$80$$0.2$$0.4$$0.6$$0.8$$1$EpochsLossesTraining LossValidation Loss
---|---
Figure 4: Accuracy and Loss for the First Phase

The training and validation accuracy graph can be seen from figure 4 and the
infractions between the training and validation accuracy seem a bit wider but
the values for training are 98.35% and validation 90.909% respectively.
Considering the loss, the validation loss has some infraction that can be seen
with different metrics the training loss is 6.79% whereas the validation loss
is 38.05% respectively.

###  3.2 Accuracy and Loss for Second Phase

The second phase contains support vector machine and not necessarily the graph
depiction is right measure but we have included the graph. The inference can
be drawn as there are no significant changes in the accuracy or loss. All the
metrics are learnt from the trained parameters of the first phase and it does
not make much sense to make mappings out of it. The training accuracy for the
support vector machine phase is generated as 98.354% and validation accuracy
is obtained as 90.909% which is similar to the first phase training and
validation accuracy. The training loss is 6.79% and validation loss is 38.052%
which is again similar to the first phase. The better inference can be
generated from different metrics intended for classification and not just the
graph depictions of the training and validation accuracy as well as loss.

$0$$20$$40$$60$$80$$0.75$$0.8$$0.85$$0.9$$0.95$EpochsAccuraciesTraining AccuracyValidation Accuracy | $0$$20$$40$$60$$80$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$EpochsLossesTraining LossValidation Loss
---|---
Figure 5: Accuracy and Loss for the Second Phase

###  3.3 Precision

The precision [33] is dependent on total number of samples that are predicted
to be positive among all the set of samples. This is a popular metric for
prognostication algorithms and requires basic elements of a confusion matrix
[34], viz. true positives, true negatives, false positives and false
negatives. The precision score obtained in 88.88% which is very close to 1 as
expected.

###  3.4 Recall

The recall [33] is a metrics which states all the positive elements from the
every single predicted element. The recall also utilizes every single element
from the confusion matrix same as precision. The recall generated for the
model is 80% which is again a very good score and gives inference that model
has performed adequately. The recall is not only metric that gives the final
inference and more precise metric can be obtained.

###  3.5 F-Score

This is a subtle metric that gives the overall flow of how efficiently does
the model perform. The building blocks of F-score [33, 35] are precision and
recall. The F-Score we got for the model is 84.21%, which is adequately good
and proves that model performs better on average.

##  4 Conclusion

The idea of entire paper revolves around an experimentation that can be
performed for polycystic ovary syndrome diagnosis problem. We proved a point
that simple machine learning algorithms can be leveraged using deep learning
for efficient performance based inclination. The Deep Linear Discriminant
Analysis idea was proven in this paper. We also introduced some of our
personal variations in implementation and they turned out to be effective from
the results section of the paper. The paper can definitely sum up many lost
ideas into practical implementation and enforce many young researchers for
better development perspective.

## Acknowledgements

We would genuinely like to thank Mr. Prasoon Kottarathil for making the
polycystic ovary syndrome dataset available through Kaggle platform. We would
also like to deeply thank for the contributions of VahidooX -
_https://github.com/VahidooX/DeepLDA_ and Thomas Chaton -
_https://github.com/tchaton/DeepLDA_ for providing us the basic
implementations of Deep LDA based from the original paper.

## References

  * [1] Gill Marshall and Leon Jonker.  An introduction to inferential statistics: A review and practical guide.  Radiography, 17, 2011.
  * [2] Iqbal H Sarker.  Machine learning: Algorithms, real-world applications and research directions.  SN Computer Science, 2(3):1–21, 2021.
  * [3] Uche Anadu Ndefo, Angie Eaton, and Monica Robinson Green.  Polycystic ovary syndrome: a review of treatment options with a focus on pharmacological approaches.  Pharmacy and therapeutics, 38(6):336, 2013.
  * [4] R. M. Cormack.  A review of classification.  Journal of the Royal Statistical Society. Series A (General), 134(3):321–367, 1971.
  * [5] Roshan Kumari and Saurabh Kr. Srivastava.  Machine learning: A review on binary classification.  International Journal of Computer Applications, 160(7):11–15, Feb 2017.
  * [6] Preeti Chauhan, Pooja Patil, Neha Rane, Pooja Raundale, and Harshil Kanakia.  Comparative analysis of machine learning algorithms for prediction of pcos.  In 2021 International Conference on Communication information and Computing Technology (ICCICT), pages 1–7, 2021.
  * [7] Nandan Kanvinde, Abhishek Gupta, and Raunak Joshi.  Binary classification for high dimensional data using supervised non-parametric ensemble method.  arXiv preprint arXiv:2202.07779, 2022.
  * [8] Abhishek Gupta, Himanshu Soni, Raunak Joshi, and Ronald Melwin Laban.  Discriminant analysis in contrasting dimensions for polycystic ovary syndrome prognostication.  arXiv preprint arXiv:2201.03029, 2022.
  * [9] Sruthi Nair, Abhishek Gupta, Raunak Joshi, and Vidya Chitre.  Combining varied learners for binary classification using stacked generalization.  arXiv preprint arXiv:2202.08910, 2022.
  * [10] Abhishek M Gupta, Sannidhi S Shetty, Raunak M Joshi, and Ronald Melwin Laban.  Succinct differentiation of disparate boosting ensemble learning methods for prognostication of polycystic ovary syndrome diagnosis.  In 2021 International Conference on Advances in Computing, Communication, and Control (ICAC3), pages 1–5. IEEE, 2021.
  * [11] Abhishek Gupta, Sruthi Nair, Raunak Joshi, and Vidya Chitre.  Residual-concatenate neural network with deep regularization layers for binary classification.  arXiv preprint arXiv:2205.12775, 2022.
  * [12] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [13] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  2016\.
  * [14] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [15] Lars Buitinck, Gilles Louppe, Mathieu Blondel, Fabian Pedregosa, Andreas Mueller, Olivier Grisel, Vlad Niculae, Peter Prettenhofer, Alexandre Gramfort, Jaques Grobler, Robert Layton, Jake VanderPlas, Arnaud Joly, Brian Holt, and Gaël Varoquaux.  API design for machine learning software: experiences from the scikit-learn project.  In ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108–122, 2013.
  * [16] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [17] Leslie J. King.  Discriminant analysis: A review of recent theoretical contributions and applications.  Economic Geography, 46:367–378, 1970.
  * [18] Somesh Das Gupta.  Discriminant analysis.  In Stephen E. Fienberg and David V. Hinkley, editors, R.A. Fisher: An Appreciation, pages 161–170, New York, NY, 1980. Springer New York.
  * [19] Alaa Tharwat, Tarek Gaber, Abdelhameed Ibrahim, and Aboul Ella Hassanien.  Linear discriminant analysis: A detailed tutorial.  AI Commun., 30:169–190, 2017.
  * [20] Matthias Dorfer, Rainer Kelz, and Gerhard Widmer.  Deep linear discriminant analysis.  arXiv preprint arXiv:1511.04707, 2015.
  * [21] Qing Tian, Tal Arbel, and James J. Clark.  Deep lda-pruned nets for efficient facial gender classification.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 512–521, 2017.
  * [22] Yann LeCun, Yoshua Bengio, et al.  Convolutional networks for images, speech, and time series.
  * [23] Laith Alzubaidi, Jinglan Zhang, Amjad J Humaidi, Ayad Al-Dujaili, Ye Duan, Omran Al-Shamma, José Santamaría, Mohammed A Fadhel, Muthana Al-Amidie, and Laith Farhan.  Review of deep learning: Concepts, cnn architectures, challenges, applications, future directions.  Journal of big Data, 8(1):1–74, 2021.
  * [24] Abhishek Gupta, Raunak Joshi, and Ronald Laban.  Detection of tool based edited images from error level analysis and convolutional neural network.  arXiv preprint arXiv:2204.09075, 2022.
  * [25] Raunak M Joshi and Deven Shah.  Refactoring faces under bounding box using instance segmentation algorithms in deep learning for replacement of editing tools.  In Intelligent Computing and Networking, pages 236–247. Springer, 2022.
  * [26] Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh.  L2 regularization for learning kernels.  arXiv preprint arXiv:1205.2653, 2012.
  * [27] Ali A. Minai and Ronald D. Williams.  Original contribution: On the derivatives of the sigmoid.  Neural Netw., 6(6):845–853, jun 1993.
  * [28] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  arXiv preprint arXiv:1803.08375, 2018.
  * [29] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [30] François Chollet et al.  Keras.  https://github.com/fchollet/keras, 2015.
  * [31] M.A. Hearst, S.T. Dumais, E. Osuna, J. Platt, and B. Scholkopf.  Support vector machines.  IEEE Intelligent Systems and their Applications, 13(4):18–28, 1998\.
  * [32] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(56):1929–1958, 2014.
  * [33] David MW Powers.  Evaluation: from precision, recall and f-measure to roc, informedness, markedness and correlation.  arXiv preprint arXiv:2010.16061, 2020.
  * [34] Kai Ming Ting.  Confusion matrix.  Encyclopedia of machine learning and data mining, 260, 2017.
  * [35] Cyril Goutte and Eric Gaussier.  A probabilistic interpretation of precision, recall and f-score, with implication for evaluation.  In European conference on information retrieval, pages 345–359. Springer, 2005.


# Finite free point processes

Adam W. Marcus
École Polytechnique Fédérale de Lausanne Research supported by NSF CAREER
Grant DMS-1552520.

###### Abstract

We use techniques from finite free probability to analyze matrix processes
related to eigenvalues, singular values, and generalized singular values of
random matrices. The models we use are quite basic and the analysis consists
entirely of expected characteristic polynomials. A number of our results match
known results in random matrix theory, however our main result (regarding
generalized singular values) seems to be more general than any of the standard
random matrix processes (Hermite/Laguerre/Jacobi) in the field. To test this,
we perform a series of simulations of this new process that, on the one hand,
confirms that this process can exhibit behavior not seen in the standard
random matrix processes, but on the other hand provides evidence that the true
behavior is captured quite well by our techniques. This, coupled with the fact
that we are able to compute the same statistics for this new model that we are
for the standard models, suggests that further investigation could be both
interesting and fruitful.

##  1 Introduction

For integers $m,n$, let $\mathcal{M}_{m,n}$ denote the collection of real
$m\times n$ matrices and let $X_{m\times n}\in\mathcal{M}_{m,n}$ have entries
consisting of independent, identically distributed
symmetric111$\mathbb{P}\left[{Y}\right]=\mathbb{P}\left[{-Y}\right]$. random
variables with variance $1$ (and all other moments bounded). We consider the
following questions:

1\. Eigenvalues:

If $A\in\mathcal{M}_{n,n}$, how do the eigenvalues of

$A+\sqrt{\theta}X_{n\times n}$

evolve for $\theta\geq 0$?

2\. Hermitian Eigenvalues:

If $B\in\mathcal{M}_{n,n}$, how do the eigenvalues of

$(B+\sqrt{\theta}X_{n\times n})+(B+\sqrt{\theta}X_{n\times n})^{\intercal}$

evolve for $\theta\geq 0$?

3\. Singular Values:

If $C\in\mathcal{M}_{m,n}$, how do the (squares of the) singular values of

$C+\sqrt{\theta}X_{m,n}$

evolve for $\theta\geq 0$?

4\. Generalized Singular Values:

If $D_{1}\in\mathcal{M}_{n_{1},k}$ and $D_{2}\in\mathcal{M}_{n_{2},k}$, how do
the (squares of the) generalized singular values of the pair

$\left\\{D_{1}+\sqrt{\theta}X_{n_{1}\times
k},D_{2}+\sqrt{\theta}X_{n_{2}\times k}\right\\}$

evolve for $\theta\geq 0$?

In this paper we investigate topics 2., 3., 4. using techniques from finite
free probability222Topic 1. has the added complication of having a
decomposition whose matrices may or may not be unitary, which is something
that must be addressed by alternative means (see [21] for developments in this
direction).. Our aim is to show that these simple models and the information
we can get about them connect deeply to the standard models in random matrix
theory.

We will first show that 2. is related to a problem that has already been
written about by Tao [22], so (in fact) the real focus in this paper will be
on 3. and 4. and the majority of the work will lie in developing a tool that
is robust enough to handle both cases at once (Section 4). Our main tool will
be finite free probability — or, in other words, expected characteristic
polynomials. Despite only containing a fraction of the information about a
distribution, we will find that expected characteristic polynomials capture
essentially all of the “non-random” elements of random matrix processes. In
particular, the point process derived from 2. will reproduce the Hermite (or
Gaussian) process (Section 2) and the point process derived from 3. will
reproduce the Laguerre (or Wishart) process (Section 3).

While those familiar with random matrix processes might naturally expect the
point process derived from 4. to reproduce the Jacobi process, we find that is
not the case (though the two are related — see Section 5.2). Most notably, the
dynamics of 4. depend very heavily on the initial matrix, a feature not
present in the Hermite, Laguerre, and Jacobi processes. As far as we can tell,
the dynamics displayed by 4. are (in this sense) more complex than any of the
processes typically studied in random matrix theory, and so may suggest a
fruitful direction of further research.

In order to get a better understanding of the point process derived from 4.,
we provide the results of a simulation that produces two vastly different
behaviors from the same starting points and then compare this to the
information provided by the expected characteristic polynomials. This, along
with some other experimental results, is presented in Section 6. Finally, we
end with some concluding comments and acknowledgments in Section 7.

###  1.1 Conventions

As mentioned in the introduction, we will let $\mathcal{M}_{m,n}$ denote the
collection of $m\times n$ real matrices and will always use $X_{m\times
n}\in\mathcal{M}_{m,n}$ to denote a matrix with entries consisting of
independent, identically distributed symmetric random variables with variance
$1$ (and all other moments bounded)333It is likely this condition can be
weakened substantially, but that is not something we will be concerned with in
this paper.. The choice of reals is arbitrary — one could substitute the
complex (and conjugate transpose for transpose) and get the same results
(except for Section 6).

For an operator
$L:\mathbb{R}[x_{1},\dots,x_{n}]\to\mathbb{R}[x_{1},\dots,x_{n}]$, we will
write

$L\left\\{{p}\right\\}$

to denote $L$ applied to the polynomial $p$.

We will use $\mathcal{SPM}_{n}\subseteq\mathcal{M}_{n,n}$ to denote the set of
signed permutation matrices of size $n$. We will say that a random matrix
$Y\in\mathcal{M}_{m\times n}$ is $\mathcal{SPM}$–bi-invariant if

$\mathbb{P}\left[{Y}\right]=\mathbb{P}\left[{Q_{m}Y}\right]=\mathbb{P}\left[{YQ_{n}}\right]$

for $Q_{m}\in\mathcal{SPM}_{m}$ and all $Q_{n}\in\mathcal{SPM}_{n}$ and we
will say that a random matrix $Z\in\mathcal{M}_{n\times n}$ is
$\mathcal{SPM}$–conjugation-invariant if

$\mathbb{P}\left[{Y}\right]=\mathbb{P}\left[{Q_{n}YQ_{n}^{T}}\right]$

for all $Q_{n}\in\mathcal{SPM}_{n}$. It should be clear that (by construction)
$X_{m,n}$ is $\mathcal{SPM}$–bi-invariant for all $n,m$. It is then easy to
check that $X_{n,n}+X_{n,n}^{\intercal}$ is
$\mathcal{SPM}$-conjugation–symmetric.

Lastly, we note that the convention in random matrix theory is to consider the
squares of the (generalized) singular values, as this can then be viewed as an
“eigenvalue” problem and we will do the same (in order to compare results
effectively).

###  1.2 Some Random Matrix Theory

Here we list the relevant random matrix processes along with the known results
concerning the evolution of their derived point processes. In each case, the
behavior of the point process differs depending on whether the Brownian motion
involved is real valued ($\beta=1$) or complex valued ($\beta=2$)
valued444Technically, quaternion valued ($\beta=4$) is possible, but that
presents complications we wish to avoid here..

####  1.2.1 Hermite Process

The real (complex) matrix Hermite process was introduced and studied by Dyson
[7]. It can be constructed as a Wiener process $(H_{t})_{t\in 0,\infty}$ with
$H_{0}$ a fixed symmetric (Hermitian) matrix and with each increment
$H_{t}-H_{s}$ equal to $\sqrt{t-s}$ times a Gaussian orthogonal (unitary)
matrix. Assuming $H_{0}$ has distinct eigenvalues, one can show that the
eigenvalues of $H_{t}$ form a point process $\lambda=\lambda(t)$ which is the
unique strong solution to the stochastic differential equation (SDE)

$d\lambda_{i}=\,\mathrm{d}W_{i}+\frac{\beta}{2}\left[\sum_{j\neq
i}\frac{1}{\lambda_{i}-\lambda_{j}}\right]\,\mathrm{d}t$ (1)

where the $W_{i}$ are (1-dimensional) independent Brownian motions.

####  1.2.2 Laguerre Process

The real (complex) matrix Laguerre process — also known as the Wishart process
— was introduced and studied in [3]. It can be constructed as a Wiener process
$(L_{t})_{t\in 0,\infty}$ with $L_{0}$ a fixed $m\times n$ real (complex)
matrix with $m\leq n$ and with each increment $L_{t}-L_{s}$ equal to an
$m\times n$ matrix with entries independent real (complex) Gaussian random
variables with variance $t-s$. Assuming $L_{0}$ has distinct singular values,
the $m$ eigenvalues of $L_{t}L_{t}^{\intercal}$ (so the squares of the
singular values of $L_{t}$) then form a point process
$\lambda=\lambda(t)\in\mathbb{R}^{m}$ which is the unique strong solution to
the SDE

$d\lambda_{i}=2\sqrt{\lambda_{i}}\,\mathrm{d}W_{i}+\beta\left[n+\sum_{j\neq
i}\frac{\lambda_{i}+\lambda_{j}}{\lambda_{i}-\lambda_{j}}\right]\,\mathrm{d}t$
(2)

where the $W_{i}$ are (1-dimensional) independent Brownian motions.

####  1.2.3 Jacobi Process

The real (complex) matrix Jacobi process was introduced and studied in [6].
The method of construction given in [6] is to first construct a Brownian
motion $Q_{n}$ on the group of $n\times n$ orthogonal (unitary) using a heat
kernel as in [1]. Then for fixed parameters $m,p,q\leq n$ with $q=n-p$, one
forms the $m\times m$ matrix

$J_{m}:=P_{m}Q_{n}P_{p}^{\intercal}P_{p}Q_{n}^{\intercal}P_{m}^{\intercal}$

where $P_{k}$ denotes the $k\times n$ projection matrix

$P_{k}=\begin{bmatrix}I_{k}&0_{k\times(n-k)}\end{bmatrix}.$

It was shown in [6] that when $J_{m}$ is real (complex) and $\max\\{p,q\\}\geq
m-1+2/\beta$ for $\beta=1$ ($\beta=2$), the eigenvalue process
$\lambda=\lambda(t)$ is the unique strong solution of the stochastic
differential equation

$d\lambda_{i}=2\sqrt{\lambda_{i}(1-\lambda_{i})}\,\mathrm{d}W_{i}+\beta\left[p(1-\lambda_{i})-q\lambda_{i}+\sum_{j\neq
i}\frac{\lambda_{i}(1-\lambda_{j})+\lambda_{j}(1-\lambda_{i})}{\lambda_{i}-\lambda_{j}}\right]\,\mathrm{d}t$
(3)

where the $W_{i}$ are (1-dimensional) independent Brownian motions.

######  Remark 1.1.

It should be noted that the construction of the Jacobi process is quite
different from the previous two. We suspect that this is because the natural
way to generalize the constructions of the Hermite and Laguerre processes to
generalized singular values does not reduce to a simple SDE (this is discussed
in Section 5.1).

###  1.3 Some remarks

As we claimed earlier, the goal is to show (and attempt to exploit) a
seemingly strange connection between the point processes discussed in the
previous section and ones derived from expected characteristic polynomials. A
common theme in random matrix theory is the fact that certain distributions
depend on the number field over which the entries are drawn. This dependence
on the number field is typically denoted $\beta$ and the general goal is to
express solutions as a function of $\beta$ in such a way that $\beta=1,2,4$
corresponds to the real, complex, and quaternion number systems. Once one has
an object (say, for example, a distribution) expressed in terms of $\beta$ in
this way, one can consider the collection of objects one would get by plugging
in any $\beta\geq 0$ as a parameter. While these may not necessarily
correspond to random matrix, understanding the dependence of the objects on
the parameter $\beta$ (apart from being of interest on its own) can often give
insights into the original random matrix problems corresponding to the
specific cases $\beta=1,2,4$. The idea that random matrices represent three
distinct points on a continuum of probability models was called the “Threefold
Way” by Dyson in his original paper. The author shares a belief first put
forth by Edelman555At least as far as the author knows. that there is at least
one other “Way” — that is, a fourth distinct point on the continuum of
probability models that can be modeled using random matrices. This distinct
“point” is the limit $\beta\to\infty$, with the corresponding model being the
algebra of expected characteristic polynomials known as “finite free
probability”.

It is worth noting that a common interpretation of the $\beta$ parameter is as
an inverse temperature (stemming from Dyson’s original paper [7] where he
noted that a number of the distributions he derived would be similar to what
one would get from an abstract “log-gas”). In that respect, one would expect
that a model for $\beta\to\infty$ would consist of the “non-random”
contributions to a given random process. While we will not attempt to define
“non-random” rigorously (in general), it is clear what this would mean in the
specific context of the processes listed in Section 1.2. Each of these
processes is a drift–diffusion process, with the parts that integrate against
$\,\mathrm{d}t$ being referred to as the “drift” and the parts that integrate
against Brownian motions being referred to as the “diffusion.” Hence in this
respect, we should expect a $\beta\to\infty$ model to exhibit similar drift
behavior to the $\beta=1,2$ models, but without any diffusion term.

##  2 Hermitian eigenvalues

We start by discussing the case of Hermitian eigenvalues. Our main tool will
be a theorem that first appeared (in spirit) in [19] but was proved in the
generality of $\mathcal{SPM}$–conjugation-invariance in [18].

######  Theorem 2.1.

Let $A,B\in\mathcal{M}_{n,n}$ be independent random matrices with $B$ being
$\mathcal{SPM}$–conjugation-invariant. Furthermore, let $P,Q$ be power series
for which the operators $P(\partial)$ and $Q(\partial)$ satisfy

$\mathbb{E}_{A}\left\\{{\mathrm{det}\left[{xI+A}\right]}\right\\}=P(\partial)\left\\{{x^{n}}\right\\}\quad\text{and}\quad\mathbb{E}_{B}\left\\{{\mathrm{det}\left[{xI+B}\right]}\right\\}=Q(\partial)\left\\{{x^{n}}\right\\}.$

Then

$\mathbb{E}_{A,B}\left\\{{\mathrm{det}\left[{xI+A+B}\right]}\right\\}=(PQ)(\partial)\left\\{{x^{n}}\right\\}$

where $PQ$ denotes the multiplication of the two power series.

As mentioned in Section 1.1, $X_{n,n}+X_{n,n}^{\intercal}$ is
$\mathcal{SPM}$-conjugation–invariant, and so Theorem 2.1 gives us (in theory)
a way to compute

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-(A+\sqrt{\theta}X_{n,n})-(A+\sqrt{\theta}X_{n,n})^{\intercal}}\right]}\right\\}$

for any $A\in\mathcal{M}_{n,n}$ simply by knowing

$\mathbb{E}_{A}\left\\{{\mathrm{det}\left[{xI-(A+A^{\intercal})}\right]}\right\\}\quad\text{and}\quad\mathbb{E}\left\\{{\mathrm{det}\left[{xI-\sqrt{\theta}(X_{n,n}+X_{n,n}^{\intercal}}\right])}\right\\}$

To do so, however, we will need to find a power series $Q$ which satisfies

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-\theta(X_{n,n})+X_{n,n}^{\intercal}}\right])}\right\\}=Q(\partial)\left\\{{x^{n}}\right\\}$

This can be derived from a result of Edelman which showed that the expected
characteristic polynomial of the Gaussian Unitary Ensemble is a Hermite
polynomial [8]. We give an alternative derivation here:

######  Lemma 2.2.

For all integers $n>0$, we have

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-\sqrt{\theta}(X_{n\times
n}+X_{n\times
n}^{\intercal})}\right]}\right\\}=e^{-\theta\partial^{2}}\left\\{{x^{n}}\right\\}$

where we are writing $e^{-\theta\partial^{2}}$ as shorthand666This also has
the benefit of revealing the underlying semigroup. for the operator

$e^{-\theta\partial^{2}}=\sum_{i=0}^{\infty}\frac{(-\theta)^{i}}{i!}\partial^{2i}.$

###### Proof.

Let $Y=\frac{X_{n\times n}+X_{n\times n}^{\intercal}}{\sqrt{2}}$. This is a
symmetric matrix with random variable entries, each of which has expectation
$0$ and variance $1$. Furthermore, all of the entries of $Y$ are independent
except for the imposed symmetry $Y_{i,j}=Y_{j,i}$. We now consider the
expansion of the polynomial

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI+Y}\right]}\right\\}$ (4)

using the formula

$\mathbb{E}\left\\{{\mathrm{det}\left[{A}\right]}\right\\}=\sum_{\sigma\in
S_{n}}(-1)^{\mathrm{sgn}(\sigma)}\mathbb{E}\left\\{{A_{1,\sigma(a)}\dots
A_{n,\sigma(n)}}\right\\}.$

Of particular note is the fact that the diagonal entries of $Y$ do not
contribute in any way to (4) — the expansion is affine in terms of each
diagonal entry and so when the expectation is taken, these terms will all
become $0$. Hence we can replace the diagonal entries of $Y$ with
(deterministic) $0$ entries and not change (4) — we call the remaining matrix
$Z$.

The utility of this transformation is that $Z$ can now be viewed as the
weighted adjacency matrix of the complete graph with independent random
variables on the edges, each of which has mean $0$ and variance $1$. It is now
a somewhat well-known result of Godsil and Gutmann that the expected
characteristic polynomial of such an adjacency matrix is the matching
polynomial of the graph [12]. It is also well known (see [11]) that the
matching polynomial of the complete graph on $n$ vertices is the Hermite
polynomial, which has the formula [24]

$H_{n}(x)=e^{-\frac{\partial^{2}}{2}}\left\\{{x^{n}}\right\\}.$

From this, it is easy to check that

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-
tZ}\right]}\right\\}=t^{n}H_{n}(x/t)=e^{-t^{2}\frac{\partial^{2}}{2}}\left\\{{x^{n}}\right\\}.$

The theorem follows by setting $t=\sqrt{2\theta}$. ∎

A combination of Theorem 2.1 and Lemma 2.2 imply that for a given matrix
$B\in\mathcal{M}_{n,n}$ with

$p_{B}(x)=\mathrm{det}\left[{xI-B-B^{\intercal}}\right]$

we have the formula

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-(B+X_{n,n})-(B+X_{n,n})^{\intercal}}\right]}\right\\}=e^{-\theta\partial^{2}}\left\\{{p_{B}(x)}\right\\}.$
(5)

At this point we could more or less refer to the reader to [22], where an
exploration of this process was conducted for other reasons (and contains a
number of interesting observations that are beyond the scope of this paper).
However it will be instructive to outline relevant parts of the analysis
briefly since we take a slightly different approach which will be somewhat
instructive for the results in Section 5. To analyze (5), we start by finding
a differential equation that it satisfies. Let $p(x)$ be an arbitrary
polynomial and let

$q(x,t)=e^{-\theta\partial^{2}}\left\\{{p(x)}\right\\}.$

Now define the root path $r(\theta)$ so that $q(r(\theta),\theta)=0$. Taking a
derivative in $\theta$, we have

$0=\frac{\partial}{\partial\theta}q(r(\theta),\theta)=q_{1}(r(\theta),\theta)r^{\prime}(\theta)+q_{2}(r(\theta),\theta)$

so (setting $\theta=0$)

$0=p^{\prime}(r)r^{\prime}-p^{\prime\prime}(r)$

or

$r^{\prime}(0)=\frac{p^{\prime\prime}(r)}{p^{\prime}(r)}.$

Now if we set

$p(x)=\prod_{i=0}^{d}(x-\lambda_{i})$

then at $x=\lambda_{i}$ we have

$\frac{p^{\prime\prime}(\lambda_{i})}{p^{\prime}(\lambda_{i})}=\sum_{j\neq
i}\frac{2}{\lambda_{i}-\lambda_{j}}.$ (6)

Hence evolution of the roots of (5) satisfy the differential equations

$\frac{\partial\lambda_{i}}{\partial\theta}=\sum_{j\neq
i}\frac{2}{\lambda_{i}-\lambda_{j}}.$

At this point in his presentation, Tao [22] notes “Curiously, this system
resembles that of Dyson Brownian motion (except with the Brownian motion part
removed, and time reversed777His comment about time reversal stems from the
fact that he was using the operator $e^{t\partial^{2}}$ instead of
$e^{-t\partial^{2}}$.).” Specifically, he is noting that this evolution
matches (up to a constant) the drift term in the Hermite process (1).

The “curiosity” comes from the fact that one would not expect this to be the
case — in some sense, this correspondence can be interpreted as asserting that
a random variable $Y$ (in our case, the characteristic polynomial of a certain
random matrix) satisfies

$f(\mathbb{E}\left\\{{Y}\right\\})=\mathbb{E}\left\\{{f(Y)}\right\\}$

for some nonlinear function $f$ (in our case, the “$k^{th}$ root” function
applied to a polynomial). This is of course not true in general, but (as we
shall see in Section 3) seems to occur more generally in this context.

##  3 Singular Values

For the sake of observing the parallels with the previous section, we will
quickly discuss the case of singular values. Our main tool will be a theorem
that first appeared (in spirit) in [15] but was proved in the generality of
$\mathcal{SPM}$–bi-invariance in [17].

######  Theorem 3.1.

Let $A$ and $B$ be $m\times n$ independent random matrices with $B$ being
$\mathcal{SPM}$–bi–invariant. Furthermore, let $P,Q$ be power series for which
the operators $P(\partial_{x}\partial_{y})$ and $Q(\partial_{x}\partial_{y})$
satisfy

$y^{n-m}\mathrm{det}\left[{xyI-
AA^{T}}\right]=P(\partial_{x}\partial_{y})\left\\{{x^{m}y^{n}}\right\\}\quad\text{and}\quad
y^{n-m}\mathrm{det}\left[{xyI-
BB^{T}}\right]=Q(\partial_{x}\partial_{y})\left\\{{x^{m}y^{n}}\right\\}.$

Then

$\mathbb{E}_{Q,R}\left\\{{y^{n-m}\mathrm{det}\left[{xyI-(A+QBR)(A+QBR)^{T}}\right]}\right\\}=(PQ)(\partial_{x}\partial_{y})\left\\{{x^{m}y^{n}}\right\\}$
(7)

where $PQ$ denotes the multiplication of the two power series.

As mentioned in Section 1.1, $X_{n,m}$ is $\mathcal{SPM}$-bi–invariant, and so
Theorem 3.1 gives us (in theory) a way to compute

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-(A+\sqrt{\theta}X_{m,n})(A+\sqrt{\theta}X_{m,n})^{\intercal}}\right]}\right\\}$

for any $A\in\mathcal{M}_{m,n}$ simply by knowing

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-
AA^{\intercal}}\right]}\right\\}\quad\text{and}\quad\mathbb{E}\left\\{{\mathrm{det}\left[{xI-\theta
X_{m,n}X_{m,n}^{\intercal}}\right]}\right\\}.$

This can again be derived from a result of Edelman which showed that the
expected characteristic polynomial of the Wishart Ensemble is a Laguerre
polynomial [8]. In Section 4.2 we will give an alternative derivation of the
following fact (see Lemma 4.2):

######  Lemma 3.2.

For all integers $m,n>0$, we have

$\mathbb{E}\left\\{{y^{n-m}\mathrm{det}\left[{xyI-\theta
X_{m,n}X_{m,n}^{\intercal}}\right]}\right\\}=e^{-\theta\partial_{x}\partial_{y}}\left\\{{x^{m}y^{n}}\right\\}$

where we are writing $e^{-\theta\partial_{x}\partial_{y}}$ as shorthand for
the operator

$e^{-\theta\partial_{x}\partial_{y}}=\sum_{i=0}^{\infty}\frac{(-\theta)^{i}}{i!}\partial_{x}^{i}\partial_{y}^{i}.$

A combination of Theorem 3.1 and Lemma 3.2 imply that for a given matrix
$B\in\mathcal{M}_{m,n}$ with

$p_{BB^{\intercal}}(x)=\mathrm{det}\left[{xI-BB^{\intercal}}\right]$

we have the formula

$y^{n-m}\mathbb{E}\left\\{{\mathrm{det}\left[{xyI-(B+X_{m,n})(B+X_{m,n})^{\intercal}}\right]}\right\\}=e^{-\theta\partial^{2}}\left\\{{y^{n-m}p_{BB^{\intercal}}(xy)}\right\\}.$
(8)

We now proceed similarly to Section 2. For the sake of simplicity, we will
assume $n\geq m$ and then find a differential equation that (8) satisfies. Let
$p(x,y)$ be an arbitrary polynomial and let

$q(x,y,t)=e^{-\theta\partial_{x}\partial_{y}}\left\\{{p(x,y)}\right\\}.$

Now define the root path $r(\theta)$ so that $q(r(\theta),1,\theta)=0$. Taking
a derivative in $\theta$, we have

$0=\frac{\partial}{\partial\theta}q(r(\theta),1,\theta)=q_{1}(r(\theta),1,\theta)r^{\prime}(\theta)+q_{3}(r(\theta),1,\theta)$

so (setting $\theta=0$)

$0=p_{1}(r,1)r^{\prime}-p_{12}(r,1)$

or

$r^{\prime}(0)=\frac{p_{12}(r,1)}{p_{1}(r,1)}.$ (9)

We can now use the fact that $x$ and $y$ play similar roles in (8), so that
for $p(x,y)=y^{n-m}p_{BB^{\intercal}}(xy)$ we have

$p_{1}=y^{n-m+1}p^{\prime}_{BB^{\intercal}}(xy)\quad\text{and}\quad
p_{2}=(n-m)y^{n-m-1}p_{BB^{\intercal}}(xy)+xy^{n-m}p^{\prime}_{BB^{\intercal}}(xy).$

Hence we can write

$p_{2}=(n-m)\frac{p}{y}+\frac{x}{y}p_{1}\quad\text{and}\quad
p_{12}=(n-m+1)\frac{p_{1}}{y}+\frac{x}{y}p_{11}$

Now if we set $y=1$, we see that (9) becomes

$r^{\prime}(0)=(n-m+1)+r\frac{p_{11}(r,1)}{p_{1}(r,1)}.$

and so for $x=\lambda_{i}$, (6) implies the evolution of the roots of
$p_{BB^{\intercal}}$ satisfy the differential equations

$\frac{\partial\lambda_{i}}{\partial\theta}=(n-m+1)+\lambda_{i}\sum_{j\neq
i}\frac{2}{\lambda_{i}-\lambda_{j}}=n+\sum_{j\neq
i}\frac{\lambda_{i}+\lambda_{j}}{\lambda_{i}-\lambda_{j}}$

So, in particular, we see that the evolution again exhibits behavior similar
to the drift term in the Brownian motion system (2). Furthermore, one can
consider this correspondence to be “curious” for the same reasons as the
Hermite case.

##  4 A tool for generalized singular values

We start by briefly reviewing the generalized singular value decomposition and
and its associated “characteristic polynomial” as motivation.

###  4.1 The generalized singular value decomposition

The generalized singular value decomposition (GSVD) is often described as a
decomposition on a pair of matrices, however we prefer the approach of [23] of
using a single block matrix. For simplicity, we will restrict to the case when
this block matrix has full column rank, however we will find a way to
eliminate this requirement shortly.

For a matrix $M\in\mathcal{M}_{n_{1}+n_{2},k}$ with rank $k$ and block form

$M=\begin{bmatrix}M_{1}\\\ M_{2}\end{bmatrix}\begin{array}[]{c}n_{1}\\\
n_{2}\end{array}$

the GSVD provides a simultaneous decomposition of $M_{1}$ and $M_{2}$ into

$M_{1}=U_{1}CH\quad\text{and}\quad M_{2}=U_{2}SH$

where

  * •

$U_{1}\in\mathcal{M}_{n_{1},n_{1}}$ satisfies
$U_{1}^{\intercal}U_{1}=U_{1}U_{1}^{\intercal}=I_{n_{1}}$

  * •

$U_{2}\in\mathcal{M}_{n_{2},n_{2}}$ satisfies
$U_{2}^{\intercal}U_{2}=U_{2}U_{2}^{\intercal}=I_{n_{2}}$

  * •

$C\in\mathcal{M}_{n_{1},k}$ and $S\in\mathcal{M}_{n_{2},k}$ are pseudo-
diagonal with $C^{\intercal}C+S^{\intercal}S=I_{k}$, and

  * •

$H$ is an $k\times k$ invertible matrix.

The diagonal entries of $C$ and $S$ satisfy $c_{i}^{2}+s_{i}^{2}=1$, and
because of this, these matrices are often referred to as cosine and sine
matrices.

The standard singular value decomposition can then be recovered by letting
$M_{2}=I$ and noting that this implies that $M_{1}$ has the decomposition

$M_{1}=U_{1}CS^{-1}U_{2}.$

When $M$ has rank $k$, there is an easy way to calculate the generalized
singular values without needing to form the entire decomposition. Letting
$W_{1}=M_{1}^{\intercal}M_{1}$ and $W_{2}=M_{2}^{\intercal}M_{2}$, it is easy
to check that

$W=(W_{1}+W_{2})^{-1/2}W_{1}(W_{1}+W_{2})^{-1/2}$ (10)

is a positive semidefinite Hermitian matrix that has the same
eigenvalues888What one actually calls the “generalized singular values”
($c_{i}$?, $c_{i}/s_{i}$?) is subject to interpretation, but all
interpretations are functions of the eigenvalues of $C^{\intercal}C$, which
will always be our subject of interest in this paper. as $C^{\intercal}C$.
Letting

$\mathrm{det}\left[{xI-(W_{1}+W_{2})^{-1/2}W_{1}(W_{1}+W_{2})^{-1/2}}\right]=\mathrm{det}\left[{W_{1}+W_{2}}\right]^{-1}\mathrm{det}\left[{x(W_{1}+W_{2})-W_{1}}\right].$
(11)

and so we can use the characteristic polynomial

$\mathrm{det}\left[{x(W_{1}+W_{2})-W_{1}}\right]=\mathrm{det}\left[{(x-1)W_{1}+xW_{2}}\right]$
(12)

which has the added benefit that it does not require $W_{1}+W_{2}$ to be
invertible.

###  4.2 A generalized characteristic polynomial

For matrices $A\in\mathcal{M}_{n_{1},k}$ and $B\in\mathcal{M}_{n_{2},k}$ we
will consider the polynomial

$p_{A,B}(x,y,z)=\mathrm{det}\left[{xI+yA^{\intercal}A+zB^{\intercal}B}\right].$
(13)

Note that $p_{A,B}(x,0,-1)$ reduces to the polynomial $p_{BB^{\intercal}}(x)$
from Section 3 and $p_{A,B}(0,z-1,z)$ reduces to (12). Hence understanding the
evolution of the polynomial

$p_{A+\sqrt{\theta}X_{n_{1},k},B+\sqrt{\theta}X_{n_{2},k}}$ (14)

will allow us to derive the information we want regarding both decompositions.

The main property of $X_{m,n}$ that we will exploit in our analysis is the
fact that it is $\mathcal{SPM}$–bi-invariant. In particular, it was shown in
[17] that for any $\mathcal{SPM}$–bi-invariant matrices $C,D$, the
coefficients of

$\mathbb{E}\left\\{{p_{A+C,B+D}(x,y,z)}\right\\}$

are each multilinear functions of the coefficients of

$\mathbb{E}\left\\{{p_{A,B}(x,y,z)}\right\\}\quad\text{and}\quad\mathbb{E}\left\\{{p_{C,D}(x,y,z)}\right\\}.$

To state the correspondence, it will be beneficial (for the moment) to instead
work with a slight transformation of the polynomials $p_{A,B}$:

$q_{A,B}(x,y,z)=y^{n_{1}}z^{n_{2}}p(x,1/y,1/z).$

The results of [17] then imply the following:

######  Theorem 4.1.

Let $A,C\in\mathcal{M}_{n_{1},k}$ and $B,D\in\mathcal{M}_{n_{2},k}$ be random
matrices with $C$ and $D$ $\mathcal{SPM}$–bi-invariate. Now let $F$ and $G$ be
bivariate polynomials for which

$\mathbb{E}\left\\{{q_{A,B}(x,y,z)}\right\\}=F(\partial_{x}\partial_{y},\partial_{x}\partial_{z})\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}\quad\text{and}\quad\mathbb{E}\left\\{{q_{C,D}(x,y,z)}\right\\}=G(\partial_{x}\partial_{y},\partial_{x}\partial_{z})\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}.$

Then

$\mathbb{E}\left\\{{q_{A+C,B+D}(x,y,z)}\right\\}=F(\partial_{x}\partial_{y},\partial_{x}\partial_{z})G(\partial_{x}\partial_{y},\partial_{x}\partial_{z})\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}.$

Given Theorem 4.1, an obvious next step to computing the polynomial in (14)
would be to find the value of the polynomial
$\mathbb{E}\left\\{{q_{aX_{n_{1},k},bX_{n_{2},k}}(x,y,z)}\right\\}$ for
general $a,b\in\mathbb{R}$. This can be done combinatorially (as in Theorem
2.2), but we will find it easier (and more instructive) to use a central limit
theorem argument.

######  Lemma 4.2.

Let $\\{A_{i}\\}_{i=1}^{\infty}\subseteq\mathcal{M}_{n_{1},k}$ and
$\\{B_{i}\\}_{i=1}^{\infty}\subseteq\mathcal{M}_{n_{2},k}$ be sequences of
independent random $\mathcal{SPM}$–bi-invariant matrices for which

$\mathbb{E}\left\\{{\mathrm{Tr}\left[{A_{i}A_{i}^{\intercal}}\right]}\right\\}=\theta_{1}n_{1}k\quad\text{and}\quad\mathbb{E}\left\\{{\mathrm{Tr}\left[{B_{i}B_{i}^{\intercal}}\right]}\right\\}=\theta_{2}n_{2}k$

for all $i$. Define the random matrices

$C_{m}=\frac{\sum_{i=1}^{n}A_{i}}{\sqrt{m}}\quad\text{and}\quad
D_{m}=\frac{\sum_{i=1}^{n}B_{i}}{\sqrt{m}}.$

Then

$\lim_{m\to\infty}\mathbb{E}\left\\{{q_{C_{m},D_{m}}(x,y,z)}\right\\}=e^{\theta_{1}\partial
x\partial
y+\theta_{2}\partial_{x}\partial_{z}}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$
(15)

###### Proof.

Due to the boundedness of the underlying random variables, we can write

$\displaystyle\mathbb{E}\left\\{{q_{A_{i}/\sqrt{m},B_{i}/\sqrt{m}}(x,y,z)}\right\\}$
$\displaystyle=x^{k}y^{n_{1}}z^{n_{2}}+\frac{\theta_{1}n_{1}k}{m}x^{k-1}y^{n_{1}-1}z^{n_{2}}+\frac{\theta_{2}n_{2}k}{m}x^{k-1}y^{n_{1}}z^{n_{2}-1}+O\left(\frac{1}{m^{2}}\right)$
$\displaystyle=1+\frac{\theta_{1}}{m}\partial_{x}\partial_{y}+\frac{\theta_{2}}{m}\partial_{x}\partial_{z}+O\left(\frac{1}{m^{2}}\right)\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

as $m\to\infty$. So by Theorem 4.1, we have

$\mathbb{E}\left\\{{q_{C_{m},D_{m}}(x,y,z)}\right\\}=\left(1+\frac{\theta_{1}}{m}\partial_{x}\partial_{y}+\frac{\theta_{2}}{m}\partial_{x}\partial_{z}+O\left(\frac{1}{{m}^{2}}\right)\right)^{m}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

which converges to the claimed polynomial as $m\to\infty$. ∎

We now show that (15) is (in fact) the same polynomial as
$\mathbb{E}\left\\{{q_{aX_{n_{1},k},bX_{n_{2},k}}(x,y,z)}\right\\}$.

######  Corollary 4.3.

$\mathbb{E}\left\\{{q_{aX_{n_{1},k},bX_{n_{2},k}}(x,y,z)}\right\\}=e^{a^{2}\partial
x\partial
y+b^{2}\partial_{x}\partial_{z}}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

###### Proof.

Let $G_{1}\in\mathcal{M}_{n_{1},k}$ and $G_{2}\in\mathcal{M}_{n_{2},k}$ be
matrices of independent standard Gaussians and let $a$ and $b$ be real
numbers. We will apply Lemma 4.2 in the case where the $A_{i}$ are independent
copies of $aG_{1}$ and the $B_{i}$ are independent copies of $bG_{2}$.

On the one hand, it is easy to see using the properties of standard Gaussians
that

$C_{m}=\frac{\sum_{i=1}^{n}A_{i}}{\sqrt{m}}=aG_{1}\quad\text{and}\quad
D_{m}=\frac{\sum_{i=1}^{n}B_{i}}{\sqrt{m}}=bG_{2}$

for all $m$. Hence we have

$\lim_{m\to\infty}\mathbb{E}\left\\{{q_{C_{m},D_{m}}(x,y,z)}\right\\}=\mathbb{E}\left\\{{q_{aG_{1},bG_{2}}(x,y,z)}\right\\}.$

On the other hand, it is easy to calculate that

$\mathbb{E}\left\\{{\mathrm{Tr}\left[{a^{2}G_{1}^{\intercal}G_{1}}\right]}\right\\}=a^{2}n_{1}k\quad\text{and}\quad\mathbb{E}\left\\{{\mathrm{Tr}\left[{b^{2}G_{2}^{\intercal}G_{2}}\right]}\right\\}=b^{2}n_{2}k$

so Lemma 4.2 implies that

$\lim_{m\to\infty}\mathbb{E}\left\\{{q_{C_{m},D_{m}}(x,y,z)}\right\\}=e^{\theta_{1}\partial
x\partial
y+\theta_{2}\partial_{x}\partial_{z}}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

Equating the two gives us

$\mathbb{E}\left\\{{q_{aG_{1},bG_{2}}(x,y,z)}\right\\}=e^{a^{2}\partial
x\partial
y+b^{2}\partial_{x}\partial_{z}}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

The remainder of the proof follows from the fact that $q_{A,B}$ is at most
quadratic in the entries of the matrices $A$ and $B$. So when those entries
are random variables, the only contributions to the polynomial come from the
first two moments. Hence having a matrix of independent standard Gaussians
will give the same result as any set of independent mean $0$, variance $1$
random variables, and so we have

$\mathbb{E}\left\\{{q_{aX_{n_{1},k},bX_{n_{2},k}}(x,y,z)}\right\\}=\mathbb{E}\left\\{{q_{aG_{1},bG_{2}}(x,y,z)}\right\\}=e^{a^{2}\partial
x\partial
y+b^{2}\partial_{x}\partial_{z}}\left\\{{x^{k}y^{n_{1}}z^{n_{2}}}\right\\}$

as claimed. ∎

Theorem 4.1 combined with Corollary 4.3 gives us the following formula for any
matrices $A\in\mathcal{M}_{n_{1},k}$ and $B\in\mathcal{M}_{n_{2},k}$:

$\mathbb{E}\left\\{{q_{A+aX_{n_{1},k},B+bX_{n_{2},k}}(x,y,z}\right\\}=e^{a^{2}\partial
x\partial y+b^{2}\partial_{x}\partial_{z}}\left\\{{q_{A,B}(x,y,z)}\right\\}.$
(16)

We end this subsection by showing how one can translate this back to the
polynomials we are truly interested in. For the sake of notation, for a
polynomial $f=f(y)$, let us define the operator $U^{y}_{k}$ as

$U^{y}_{k}\left\\{{p(y)}\right\\}=y^{k}p\left(\frac{1}{y}\right).$ (17)

######  Lemma 4.4.

Let $i,j$ be integers with $j\leq s$. Then

$U_{y}^{s}\circ e^{\theta\partial_{x}\partial_{y}}\circ
U_{y}^{s}\left\\{{x^{i}y^{j}}\right\\}=(1+\theta
y\partial_{x})^{s-j}\left\\{{x^{i}y^{j}}\right\\}$

###### Proof.

We compute:

$\displaystyle U_{y}^{s}\circ e^{\theta\partial_{x}\partial_{y}}\circ
U_{y}^{s}\left\\{{x^{i}y^{j}}\right\\}$ $\displaystyle=U_{y}^{s}\circ
e^{\theta\partial_{x}\partial_{y}}\left\\{{x^{i}y^{s-j}}\right\\}$
$\displaystyle=U_{y}^{s}\left\\{{\sum_{k}\frac{\theta^{k}}{k!}{\partial_{x}^{k}\partial_{y}^{k}}\left\\{{x^{i}y^{s-j}}\right\\}}\right\\}$
$\displaystyle=U_{y}^{s}\left\\{{\sum_{k}\theta^{k}\partial_{x}^{k}\binom{s-j}{k}x^{i}y^{s-j-k}}\right\\}$
$\displaystyle=(1+\theta y\partial_{x})^{s-j}\left\\{{x^{i}y^{j}}\right\\}.$

∎

###  4.3 Returning to the decompositions

We now return the respective decompositions to give a proof of Lemma 3.2 and
the corresponding result for generalized singular values. Consider formula
(16) in the case when $A$ and $B$ are $0$ matrices and $a=b=1$. Then
$q_{A,B}(x,y,z)=x^{k}$ and so Lemma 4.4 implies

$\mathbb{E}\left\\{{p_{X_{n_{1},k},X_{n_{2},k}}(x,y,z)}\right\\}=(1+y\partial_{x})^{n_{1}}(1+z\partial_{x})^{n_{2}}\left\\{{x^{k}}\right\\}.$

As noted at the beginning of Section 4.2, we can get the singular value
polynomial by setting $y=-1$ and $z=0$:

$\mathbb{E}\left\\{{\mathrm{det}\left[{xI-X^{\intercal}_{n_{1},k}X_{n_{1},k}}\right]}\right\\}=(1-\partial_{x})^{n_{1}}\left\\{{x^{k}}\right\\}=(-1)^{k}k!L^{(n_{1}-k)}_{k}(x)$

where $L_{n}^{(\alpha)}(x)$ is a Laguerre polynomial [24]. Setting $x=0$ and
$y=z-1$, on the other hand, gives us the generalized singular value
polynomial:

$\mathbb{E}\left\\{{\mathrm{det}\left[{(z-1)X^{\intercal}_{n_{1},k}X_{n_{1},k}+zX^{\intercal}_{n_{2},k}X_{n_{2},k}}\right]}\right\\}=(1+(z-1)\partial_{x})^{n_{1}}(1+z\partial_{x})^{n_{2}}\left\\{{x^{k}}\right\\}=P^{(n_{2}-k,n_{1}-k)}_{k}(2x-1)$

where $P_{n}^{(\alpha,\beta)}(x)$ is a Jacobi polynomial [24]. For those
familiar with random matrix theory, these two results are perhaps not
surprising as the random matrix ensembles one would get by using i.i.d.
standard normals random variables in the random matrices $X_{k,k},X_{n_{1},k}$
and $X_{n_{2},k}$ are often called the Hermite, Laguerre, and Jacobi ensembles
(for this reason) [10].

##  5 Brownian Motion

Using the results in Section 4.2, we have reduced the problem of interest to
understanding the effects of a certain differential operator on a generalized
characteristic polynomial. The operator of interest is $Q^{\theta}_{s,t}$
defined by

$Q^{\theta}_{s,t}[x^{i}y^{j}z^{k}]\mapsto(1+\theta
y\partial_{x})^{s-j}(1+\theta
z\partial_{x})^{t-k}\left\\{{x^{i}y^{j}z^{k}}\right\\}$ (18)

and then extended to general polynomials linearity.

To understand the evolution of these polynomials, we first derive a
differential equation that this operator satisfies. For a fixed polynomial
$p(x,y,z)$, we define the polynomial

$\tilde{p}(x,y,z,\theta)=Q^{\theta}_{s,t}\left\\{{p(x,y,z)}\right\\}$

and then consider variables
$\tilde{x}=\tilde{x}(\theta),\tilde{y}=\tilde{y}(\theta),\tilde{z}=\tilde{z}(\theta)$
for which

$\tilde{p}(\tilde{x}(\theta),\tilde{y}(\theta),\tilde{z}(\theta),\theta)=0.$

In particular, we can take a derivative in $\theta$ to get

$\tilde{p}_{1}\tilde{x}^{\prime}+\tilde{p}_{2}\tilde{y}^{\prime}+\tilde{p}_{3}\tilde{z}^{\prime}+\tilde{p}_{4}=0.$

######  Lemma 5.1.

$\tilde{p}_{4}(x,y,z,0)=(n_{1}y+n_{2}z)\tilde{p}_{1}(x,y,z,0)-y^{2}\tilde{p}_{12}(x,y,z,0)-z^{2}\tilde{p}_{13}(x,y,z,0)$

###### Proof.

Write

$p(x,y,z)=\sum_{i,j,k}c_{i,j,k}x^{i}y^{j}z^{k}$

where the degree of $y$ is at most $s$ and the degree of $z$ is at most $t$.
Then

$\displaystyle Q^{\theta}_{s,t}\left\\{{p(x,y,z)}\right\\}$
$\displaystyle=\sum_{i,j,k}c_{i,j,k}(1+\theta
y\partial_{x})^{n_{1}-j}(1+\theta
z\partial_{x})^{n_{2}-k}\left\\{{x^{i}y^{j}z^{k}}\right\\}$
$\displaystyle=\sum_{i,j,k}c_{i,j,k}\left(1+\theta((n_{1}-j)y\partial_{x}+(n_{2}-k)z\partial_{x})+O(\theta^{2})\right)\left\\{{x^{i}y^{j}z^{k}}\right\\}$
$\displaystyle=\sum_{i,j,k}c_{i,j,k}\left(x^{i}y^{j}z^{k}+i(n_{1}-j)\theta
x^{i-1}y^{j+1}z^{k}+i(n_{2}-k)\theta x^{i-1}y^{j}z^{k+1}\right)+O(\theta)^{2}$

and so

$\frac{\partial}{\partial\theta}Q^{\theta}_{n_{1},n_{2}}\left\\{{p(x,y,z)}\right\\}\bigg{|}_{\theta=0}=\sum_{i,j,k}ic_{i,j,k}\left((n_{1}-j)x^{i-1}y^{i+1}z^{k}+(n_{2}-k)x^{i-1}y^{j}z^{k+1}\right).$

Now we rewrite

$i(n_{1}-j)x^{i-1}y^{j+1}z^{k}=in_{1}x^{i-1}y^{j+1}z^{k}-ijx^{i-1}y^{j+1}z^{k}=n_{1}y\partial_{x}\left\\{{x^{i}y^{j}z^{k}}\right\\}-y^{2}\partial
x\partial y\left\\{{x^{i}y^{j}z^{k}}\right\\}$

and similarly

$i(n_{2}-k)x^{i-1}y^{j}z^{k+1}=in_{2}x^{i-1}y^{j}z^{k+1}-ikx^{i-1}y^{j}z^{k+1}=n_{2}z\partial_{x}\left\\{{x^{i}y^{j}z^{k}}\right\\}-z^{2}\partial
x\partial z\left\\{{x^{i}y^{j}z^{k}}\right\\}.$

Hence

$\tilde{p}_{4}(x,y,z,0)=\frac{\partial}{\partial\theta}Q^{\theta}_{n_{1},n_{2}}\left\\{{p(x,y,z)}\right\\}\bigg{|}_{\theta=0}=(n_{1}y+n_{2}z)p_{1}(x,y,z)-y^{2}p_{12}(x,y,z)-z^{2}p_{13}(x,y,z)$

and the result follows from the fact that $\tilde{p}(x,y,z,0)=p(x,y,z)$. ∎

It is worth mentioning that the equation we have derived:

$\tilde{p}_{1}\tilde{x}^{\prime}+\tilde{p}_{2}\tilde{y}^{\prime}+\tilde{p}_{3}\tilde{z}^{\prime}=-(n_{1}y+n_{2}z)\tilde{p}_{1}+y^{2}\tilde{p}_{12}+z^{2}\tilde{p}_{13}$

still has a bit of freedom to be manipulated due to the homogeneity if $p$.
How we do that manipulation will depend on how we intend to use the
polynomial. For example, if our goal is going to be to plug in
$\tilde{y}=\tilde{w}-1$ and $\tilde{z}=\tilde{w}$ at some point, then the
terms

$y^{2}p_{12}\quad\text{and}\quad z^{2}p_{13}$

will not reduce to functions of $\tilde{w}$ easily. We can fix this as
follows:

######  Lemma 5.2.

Let $p(x,y,z)$ be a $k$-homogeneous polynomial. Then

$y^{2}p_{12}+z^{2}p_{13}=(k-1)(y+z)p_{1}-x(y+z)p_{11}-yz(p_{12}+p_{13}).$

###### Proof.

Since $p$ is $k$ homogeneous, $p_{1}$ is $k-1$ homogeneous, and so

$xf_{11}+yf_{12}+zf_{13}=(k-1)f_{1}$

Solving for $f_{13}$ (and separately for $f_{12}$) then gives

$z^{2}f_{13}=-xzf_{11}-yzf_{12}+z(d-1)f_{1}\quad\text{and}\quad
y^{2}f_{12}=-xyf_{11}-yzf_{13}+y(d-1)f_{1}.$

and the lemma follows by adding them together. ∎

Hence we have derived the following equation for $p$ at the point
$(\tilde{x},\tilde{y},\tilde{z})$ for $\theta=0$:

$p_{1}\tilde{x}^{\prime}+\tilde{p}_{2}\tilde{y}^{\prime}+\tilde{p}_{3}\tilde{z}^{\prime}=-(n_{1}-k+1)\tilde{y}p_{1}-(n_{2}-k+1)\tilde{z}p_{1}-\tilde{x}(\tilde{y}+\tilde{z})p_{11}-\tilde{y}\tilde{z}(p_{12}+p_{13}).$
(19)

###  5.1 Jacobi Process

We now make the substitution described in (12), letting

$h(x,u)=p(x,u-1,u)$

and setting $\tilde{y}=\tilde{w}-1$ and $\tilde{z}=\tilde{w}$ in (19) we get

$\displaystyle\tilde{x}^{\prime}h_{1}(\tilde{x},\tilde{w})+\tilde{w}^{\prime}h_{2}(\tilde{x},\tilde{w})=$
$\displaystyle-(n_{1}-k+1)(\tilde{w}-1)h_{1}(\tilde{x},\tilde{w})-(n_{2}-k+1)\tilde{w}h_{1}(\tilde{x},\tilde{w})$
$\displaystyle-\tilde{x}(2\tilde{w}-1)h_{11}(\tilde{x},\tilde{w})-\tilde{w}(\tilde{w}-1)h_{12}(\tilde{x},\tilde{w})$

which for $\tilde{x}=0$ simplifies to

$\tilde{w}^{\prime}\frac{h_{2}}{h_{1}}=-(n_{1}-k+1)(\tilde{w}-1)-(n_{2}-k+1)\tilde{w}-\tilde{w}(\tilde{w}-1)\frac{h_{12}}{h_{1}}.$
(20)

We would like to compare this to the evolution of the random matrix version of
the Jacobi process given in (3) with parameters $p=n_{1}$ and $q=n_{2}$. At
first glance, these are not the same, but we claim that this is because (20)
is actually significantly more general than (3).

Before addressing this, however, we wish to point out that not only are (20)
and (3) different, they are fundamentally different. In particular, (3) is an
eigenvalue evolution that only depends on the original eigenvalues. That is,
if we let

$A_{\theta}=A+\sqrt{\theta}X_{n_{1},k}\quad\text{and}\quad
B_{\theta}=B+\sqrt{\theta}X_{n_{2},k}$

then the evolution of (3) would imply that one can compute the eigenvalues of

$(A_{\theta}A_{\theta}^{T}+B_{\theta}B_{\theta}^{T})^{-1/2}A_{\theta}A_{\theta}^{T}(A_{\theta}A_{\theta}^{T}+B_{\theta}B_{\theta}^{T})^{-1/2}$
(21)

as a function of the eigenvalues of

$(A_{0}A_{0}^{T}+B_{0}B_{0}^{T})^{-1/2}A_{0}A_{0}^{T}(A_{0}A_{0}+B_{0}B_{0})^{-1/2}.$

This is not the case for (20), and for good reason. As we will see in Section
6, the eigenvalues of (21) will depend on (among other things) the
relationship between the eigenvectors of $A_{0}A_{0}^{T}$ and
$B_{0}B_{0}^{T}$, and so to some extent it is actually surprising that we can
get all of the necessary information even from the complete polynomial
$h(\tilde{x},\tilde{w})$.

###  5.2 The matrix Jacobi process, revisited

To see why the two solution differ, let us first review how the standard
matrix Jacobi process is defined. As we saw in Section 4.1, the method for
obtaining the generalized singular values of

$M=\begin{bmatrix}M_{1}\\\ M_{2}\end{bmatrix}\begin{array}[]{c}n_{1}\\\
n_{2}\end{array}$

was via the formula (10) where $W_{1}=M_{1}^{\intercal}M_{1}$ and
$W_{2}=M_{2}^{\intercal}M_{2}$. Note that $W_{1}+W_{2}=M^{\intercal}M$, and
that

$M_{1}=\begin{bmatrix}I_{n_{1}}&0\end{bmatrix}M$

and so another way to write (10) would be as

$W=(M^{\intercal}M)^{-1/2}M^{\intercal}\begin{bmatrix}I_{n_{1}}\\\
0\end{bmatrix}\begin{bmatrix}I_{n_{1}}&0\end{bmatrix}M(M^{\intercal}M)^{-1/2}M^{\intercal}.$

One can simplify this using the (normal) singular value decomposition: letting
$M=U\Sigma V^{\intercal}$, and simplifying, we get that

$W=\begin{bmatrix}I_{k}&0\\\
\end{bmatrix}U^{\intercal}\begin{bmatrix}I_{n_{1}}\\\
0\end{bmatrix}\begin{bmatrix}I_{n_{1}}&0\end{bmatrix}U\begin{bmatrix}I_{k}\\\
0\end{bmatrix}=C^{\intercal}C$

where $C$ is the upper $n_{1}\times k$ corner of the unitary matrix
$U\in\mathcal{M}_{n_{1}+n_{2},n_{1}+n_{2}}$.

When $M$ is a matrix of Gaussians, the derived distribution on $U$ is the
Haar-distribution. This shows the equivalence of two methods for creating a
Jacobi ensemble (a fact that has been used numerous times in the literature
[4, 9]):

  1. 1.

as the generalized singular values of a random Gaussian matrix

  2. 2.

as the singular values of the upper corner of a Haar-distributed unitary
matrix

The issue occurs when one tries to turn this into a process. The standard
matrix Jacobi process uses the “upper corner” model, replacing the Haar
distribution with a unitary Brownian motion (starting, possibly at a fixed
unitary matrix) [6]. And this is precisely where the two processes diverge —
for general matrices $A_{\theta}$, approximating the left eigenvectors of
$A_{\theta}$ with a unitary Brownian motion centered at the left eigenvectors
of $A_{0}$ is only that (an approximation). The actual distribution will
depend on the singular values of $A$, which (in turn) will depend on the
relationship between the two blocks.

On the other hand, if one was able to artificially force the singular values
of $A_{\theta}$ to remain the same throughout the process, then this extra
effect (while not necessarily zero) will not change, and so one can then hope
to derive an evolution equation which holds for all $A_{\theta}$. And this is
(essentially) what is being done in the “upper corner process” — by using a
unitary Brownian motion, one is forcing all of the singular values of the
$A_{\theta}$ to be $1$ (for all $\theta$)999This is often how Brownian motion
is defined — as the process on the space of unitary matrices which has
instantaneous change equal to that of the additive process [20].. And it is
not hard to show that for any $A_{\theta}$ where this is the case, dynamics in
(20) at that $\theta$ can be computed explicitly in terms of the current point
configuration:

######  Lemma 5.3.

Let $A\in\mathcal{M}_{(n_{1}+n_{2}),k}$ be such that $A^{\intercal}A=I_{k}$.
Then $h_{1}(0,u)=h_{2}(0,u)$ and

$\frac{h_{12}(0,\lambda_{i})}{h_{1}(0,\lambda_{i})}=\sum_{j\neq
i}\frac{2}{\lambda_{i}-\lambda_{j}}$

where $\\{\lambda_{i}\\}_{i=1}^{k}$ are the roots of $h(0,u)$.

###### Proof.

Letting

$A=\begin{bmatrix}A_{1}\\\ A_{2}\end{bmatrix}\begin{array}[]{c}n_{1}\\\
n_{2}\end{array}$

we have that $A^{\intercal}A=I_{k}$ is equivalent to
$A_{1}^{\intercal}A_{1}+A_{2}^{\intercal}A_{2}=I$, so plugging these into the
generalized characteristic polynomial (13), you get

$h(x,u)=p(x,u-1,u)=\mathrm{det}\left[{xI+(u-1)A_{1}^{\intercal}A_{1}^{\intercal}+u(I-A_{1}^{\intercal}A_{1})}\right]=\mathrm{det}\left[{(x+u)I-A_{1}^{\intercal}A_{1}}\right].$

Hence derivatives in $x$ and derivatives in $u$ are the same, so (in
particular) $h_{1}(0,u)=h_{2}(0,u)$. For the same reason, we have

$\frac{h_{12}(0,u)}{h_{1}(0,u)}=\frac{h_{22}(0,u)}{h_{2}(0,u)}$

where if we let $g(u)=h(0,u)$, then we can again use (6) to get

$\frac{h_{22}(0,\lambda_{i})}{h_{2}(0,\lambda_{i})}=\frac{g^{\prime\prime}(\lambda_{i})}{g^{\prime}(\lambda_{i})}=\sum_{j\neq
i}\frac{2}{\lambda_{i}-\lambda_{j}}$

as claimed. ∎

Hence whenever $A_{\theta}^{\intercal}A_{\theta}=I_{k}$, (20) reduces to

$\displaystyle\frac{\partial\lambda_{i}}{\partial\theta}$
$\displaystyle=-(n_{1}-k+1)(\lambda_{i}-1)-(n_{2}-k+1)\lambda_{i}-\sum_{j\neq
i}\frac{2\lambda_{i}(\lambda_{i}-1)}{\lambda_{i}-\lambda_{j}}$
$\displaystyle=-n_{1}(\lambda_{i}-1)-n_{2}\lambda_{i}+\sum_{j\neq
i}\frac{\lambda_{i}(1-\lambda_{j})+\lambda_{j}(1-\lambda_{i})}{\lambda_{i}-\lambda_{j}}$

matching the drift term in (3) exactly. To end the section, we find it worth
mentioning that the GSVD of a unitary matrix is a decomposition that is well-
studied in its own right (known as the CS-decomposition) [13]. This fact was
noted in [9] as well, however only in the context of the Jacobi ensemble
itself (not the process). Hence this might be the more appropriate
“decomposition” to associate directly to the matrix Jacobi process and, if
this is the case, then it suggests that a random matrix process related to the
true GSVD has yet to be considered in full generality.

##  6 Simulation

The purpose of this section is 3-fold. Firstly, we hope to convince the reader
of the assertion made in Section 5.1 that the point process described by 4. is
fundamentally different than the point processes defined by 2. and 3. in the
respect that future point configurations are not determined by the current
point configuration alone.

Secondly, we hope to explore the degree to which the $k$th root of the
expected characteristic polynomial correspond to the expected value of the
$k$th root. The use of free probability as a tool for estimating random matrix
statistics has a long history, hence there is hope that finite free
probability has similar potential. The obvious upside of finite free
probability is that it holds in fixed dimensions, and so one might hope that
in situations where a fixed dimensions ensemble has features that do not
remain in large dimensions that finite free probability could give better
results. The obvious downside is that, unlike in the case of free probability,
it is unclear how well the statistics that one is able to compute (the
expected characteristic polynomials) relate to the statistics one might want
to compute (moments or probabilities). Thus we hope to convey some indication
as to what the potential is in this respect.

Thirdly, we wish to explore the recurring pattern in this paper of the point
process defined by the polynomial convolution matching the “non-random” or “as
$\beta\to\infty$” part of some random matrix process. While it is only a
limited set, we wish to examine the $\beta=1$ cases and $\beta=2$ cases to see
whether we can find support for such a claim experimentally and also to see
whether there is some further intuition we can gain in this regard (like the
appearance of monotonicity, for example).

###  6.1 The experiment

We first describe the experiment. We set $k=4,n_{1}=5,n_{2}=10$, and use the
following $k\times k$ diagonal matrices:101010One would not be wrong to
consider the choice of matrices $H$ and $K$ to be somewhat extreme and not
representative of a typical application. One can show that even small
differences between $H$ and $K$ result in different behaviors, but “how
different” is unclear. The choice of these particular matrices was made with
the hopes of making various features visually recognizable.

$D_{1}=\begin{bmatrix}2&&&\\\ &2&&\\\ &&5&\\\ &&&5\end{bmatrix}\quad
D_{2}=\begin{bmatrix}10&&&\\\ &10&&\\\ &&10&\\\ &&&10\end{bmatrix}\quad
H=\begin{bmatrix}1&&&\\\ &1&&\\\ &&1&\\\ &&&500\end{bmatrix}\quad
K=\begin{bmatrix}500&&&\\\ &1&&\\\ &&1&\\\ &&&1\end{bmatrix}$

We add $(n_{1}-k)$ rows of $0$ to $D_{1}$ to form
$F_{1}\in\mathcal{M}_{n_{1},k}$ (and similarly to form
$F_{2}\in\mathcal{M}_{n_{2},k}$ from $D_{2}$). We then pick constants
$\delta,\epsilon$ and form

$A=\begin{bmatrix}A_{1}\\\ A_{2}\end{bmatrix}=\delta\begin{bmatrix}F_{1}H\\\
F_{2}H\end{bmatrix}\quad\text{and}\quad B=\begin{bmatrix}B_{1}\\\
B_{2}\end{bmatrix}=\epsilon\begin{bmatrix}F_{1}K\\\ F_{2}K\end{bmatrix}.$

with $\delta,\epsilon$ chosen so that
$\mathrm{Tr}\left[{A^{\intercal}A}\right]=\mathrm{Tr}\left[{B^{\intercal}B}\right]=1$.

Using (10), one can easily check that multiplying on the right by an
invertible matrix does not change the eigenvalues of $C^{\intercal}C$ and so
the squares of the generalized singular values of $A$ and $B$ will be the
same:

$\left(\frac{4}{13},\frac{4}{13},\frac{25}{34},\frac{25}{34}\right)\approx(0.308,0.308,0.735,0.735).$

For each trial, we consider two independent point processes: one starting at
$A^{0}=A$ and the other at $B^{0}=B$. At time step $i$, we generate a random
matrix $Z_{i}\in\mathcal{M}_{(n_{1}+n_{2}),k}$ with independent normal entries
($\mu=0,\sigma^{2}=10^{-8}$) to add to both $A$ and $B$111111The reason for
adding the same matrix is to ensure that any differences that will appear are
not a result of random effects. However numerous tests were done using
independent matrices and each gave similar results.:

$A^{i}\leftarrow A^{i-1}+Z_{i}\quad\text{and}\quad B^{i}\leftarrow
A^{i-1}+Z_{i}.$

We calculate the squares of the generalized singular values
$c_{A^{i}}^{2},c_{B^{i}}^{2}\in\mathbb{R}^{k}$ (both in increasing order). We
then average over $500$ trials, obtaining the (average) squared generalized
singular values
$\overline{c}_{A^{i}}^{2},\overline{c}_{B^{i}}^{2}\in\mathbb{R}^{k}$. These
will be compared to the roots of the expected generalized singular value
polynomials

$\tilde{p}_{A^{i}}(x,y,z)=Q^{(i\sigma^{2})}_{n_{1},n_{2}}\left\\{{p_{A}(x,y,z)}\right\\}\quad\text{and}\quad\tilde{p}_{B^{i}}(x,y,z)=Q^{(i\sigma^{2})}_{n_{1},n_{2}}\left\\{{p_{B}(x,y,z)}\right\\}.$

To simplify the wording slightly, rather than writing (for example) “notice
that $\tilde{p}_{A^{i}}$ (respectively $\tilde{p}_{B^{i}}$) are”, we will
write simply “notice that the $\tilde{p}_{M^{i}}$ are” (essentially anywhere
there is an $M$, it should be taken as being a statement for both $A$ and $B$,
separately).

###  6.2 Results

We first hope to convey the fact that the evolution of $A$ and $B$ — despite
starting at the same generalized singular values and using the same “random”
matrices — are fundamentally different. For this purpose we can examine the
evolution of the point process over various time scales (Figure 1). There are,
of course, similarities, it is easy to see that both $A$ and $B$ exhibit
repulsion between the paths early on in addition to converging to the roots of
a Jacobi polynomial in the long term (both of which should be expected).
However $B$ seems to converge very quickly to the asymptotic limit whereas $A$
takes a (very) long time. Of particular interest is the path of the largest
value (in red); the path of $A$ actually moves away from the eventual
asymptotic limit for some time during the middle range. This is in direct
contradiction to the other paths of $A$ and all of the paths of $B$ which seem
to move in the direction of the asymptotic limit.

First $10^{3}$ steps.

First $10^{5}$ steps.

First $10^{7}$ steps.

Figure 1: Evolution of roots of $\tilde{p}_{A^{i}}$ (+) vs $\tilde{p}_{B^{i}}$
(-) vs asymptotes (dashed lines).

For the second goal, we will focus on the part of the process where the most
action happens. Figure 2 shows the first 100 steps of the process for each
matrix process (in colors) as well as the paths taken by the roots of the
$\tilde{p}_{M^{i}}$, in black). The third plot shows the the two processes
together. Figure 2 suggests the $\tilde{p}_{M^{i}}$ are capturing the general
behavior of each $\overline{c}^{2}_{M^{i}}$ quite well. Furthermore, the
distance between the two seems to be rather consistent. This is corroborated
by Figure 3, which shows the actual distributions of $c^{2}_{A^{i}}$
normalized so that the root of $\tilde{p}_{A^{i}}$ is at the center.

$\overline{c}^{2}_{A^{i}}$ vs $\tilde{p}_{A^{i}}$

$\overline{c}^{2}_{B^{i}}$ vs $\tilde{p}_{B^{i}}$

$\overline{c}^{2}_{A^{i}}$ (+) vs $\overline{c}_{B^{i}}$ (-)

Figure 2: Experimental values (colors) versus polynomial predictions (black)
for the first 100 steps.

blue (smallest)

yellow

green

red (largest)

Figure 3: The distributions of $c^{2}_{A^{i}}$ shifted so that the centers are
the roots of $\tilde{p}_{A^{i}}$ (for $i=20,100$).

To address the third goal we will find it more useful to examine the moments
induced by the $\overline{c}^{2}_{M^{i}}$ instead of the points themselves.
Figure 4 and Figure 5 show the relationship between the moments generated by
the $\overline{c}_{M^{i}}$ (in yellow) and the moments generated by the roots
of the $\tilde{p}_{M^{i}}$ (in blue). An additional data point (in green)
consists of the moments of the $\overline{c}^{2}_{M^{i}}$ when we change the
simulation slightly to use complex Gaussians when forming the random $Z_{i}$
(all transposes become conjugate transposes as well). Note that these are the
moments of $\overline{c}_{M^{i}}$ and not $\overline{c}^{2}_{M^{i}}$, so the
first image is (in some sense) the $1/2$ moment of $\overline{c}^{2}_{M^{i}}$.
We found this view to be more compelling as it shows the concave functions
monotone increasing in $\beta$ and the convex ones monotone decreasing (signs
of possible majorization)121212One might also notice in Figure 2 that the
errors in the largest value (red) are consistently positive, whereas the
errors in the smallest value (blue) are negative (also signs of possible
majorization).. We will discuss conjectures related to this relationship in
Section 7.

Figure 4:  The first four moments derived from $\overline{c}_{A^{i}}$ for
$Z_{i}$ real valued (yellow) and complex valued (green). The moments derived
from the roots of $\tilde{p}_{A^{i}}(0,w-1,w)$ are in blue.

Figure 5:  The first four moments derived from $\overline{c}_{B^{i}}$ for
$Z_{i}$ real valued (yellow) and complex valued (green). The moments derived
from the roots of $\tilde{p}_{B^{i}}(0,w-1,w)$ are in blue.

##  7 Conclusion

In many ways, it seems that the work in this paper suggests more interesting
problems than it solves. A number of these questions were discussed at the
beginning of Section 6, however we would like to comment a bit more on the
theme of the polynomial convolution matching the “non-random” part of some
random matrix process. One intriguing part of the polynomial convolutions is
that (unlike many other things in random matrix theory) they do not depend on
$\beta$, a fact that has been discovered in various contexts a number of times
[2, 14]. It is also quite common for $\beta$-ensembles on $n\times n$ matrices
to converge as $\beta\to\infty$ to a uniform distribution on a set of $n$
points. Any such distribution (finite, on $n$ points) is completely determined
by the values of its expected elementary symmetric polynomials, and so one can
hope to find a polynomial convolution which captures this behavior completely.
The typical way to prove this, however, would be to express the behavior of
interest as a function of $\beta$ and then take the limit. In situations where
it is computationally infeasible to find such a function, are there other
methods one could use to prove such a correspondence? If so, Section 6
suggests that there may be a measurable relationship between distributions as
$\beta$ increases, and so such a result could lead to improved estimates for
classical ensembles.

###  7.1 Random matrix theory

The most obvious question in this direction is whether a Gaussian point
process that comes from 4. can be solved. That is, given matrices
$A,G\in\mathcal{M}_{m,n}$ where $A$ is fixed and the entries of $G$ are i.i.d.
Gaussians, can we find the exact distribution of the (squares of the)
generalized singular values of $A+\theta G$? It is certainly understandable if
any previous attempt seemed overly daunting — as we have mentioned,
instantaneous behavior of the Hermite, Laguerre, and Jacobi matrix processes
depend only on the current point configuration and a small number of
parameters, whereas it should be clear from Section 6 that the point process
derived from 4. depends on a much larger number of parameters. However the
results of Section 5.1 suggest that these parameters can be captured by
natural relationships between the two matrices, in which case an explicit
formula could be possible. We would think that such a formula would certainly
be an interesting development in field.

###  7.2 Finite free probability

Finding a polynomial convolution for a matrix operation (when it exists) tends
to be a fairly reasonable task; the results of [17], for example, reduce a
number of possible convolution combinations to a straightforward calculation.
The opposite direction — given an operation on polynomials, trying to find
matrices (and matrix operations) that they correspond to — seems much harder.
Even in the case where we know the operation, it tends to be hard to prove
that there exist concrete (non-random) matrices that behave in this way. The
major tool in this respect is a result of Helton and Vinnikov concerning real
stable polynomials [16], however their result is quite quickly reaching the
end of its utility (it is not true when the homogeneous polynomials in
question have more than $3$ variables). That said, it would not be surprising
if polynomial convolutions that are based on matrix operations did not also
maintain determinantal representation, and it would be interesting to find new
ways to prove such statements.

The condition of the random matrix $X_{m,n}$ being symmetric turns out to be
far more than we need — as we have seen, the expected characteristic
polynomial (13) is (at most) quadratic in the entries of the individual
matrices. The utility of having the symmetry condition is that it makes the
resulting random matrix distribution $\mathcal{SPM}$–bi-invariant. This
suggests that a weaker condition than $\mathcal{SPM}$–bi-invariance might be
sufficient for gaining the required amount of symmetry to be able to use
Theorem 4.1. One natural candidate to replace the signed permutation matrices
is the collection of matrices in the standard representation of $S_{n+1}$ (a
set of size $(n+1)!$ instead of $2^{n}n!$). Furthermore, the validity of this
replacement would follow easily from a conjecture in [17] that Theorem 4.1
holds in a slightly more general context (we refer the interested reader to
[17] for more details).

###  7.3 Acknowledgments

The author would like to recognize the IPAM program in Quantitative Linear
Algebra, without which this paper would likely not exist. The results
presented here all came, in one way or another, from discussions that were
started at IPAM. The author would also like to thank Benno Mirabelli, who
helped to continue these discussions after the program had ended (and
eventually led to the ideas in this paper).

## References

  * [1] P. Biane. Free Brownian motion, free stochastic calculus and random matrices. Fie. Inst. Comm. 12, Amer. Math. Soc. Providence, RI, 1997
  * [2] A. Borodin, Stochastic higher spin six vertex model and Macdonald measures. J. of Math. Physics 59.2 (2018): 023301.
  * [3] M. F. Bru, Wishart Process. J. Theor. Prob. 4, (1991) 725–751. MR113213
  * [4] B. Collins, Intégrales matricielles et probabilités non-commutatives, Ph.D. thesis, Université Paris 6, 2003.
  * [5] N. Demni. Free Jacobi process. J. of Theor. Prob. 2008 Mar 1; 21(1): 118-143.
  * [6] Y. Doumerc, Matrices, aléatoires, processus stochastiques et groupes de réflexions, Ph.D. thesis, Université Toulouse III.
  * [7] F. J. Dyson, A Brownian-motion model for the eigenvalues of a random matrix, J. Math. Phys. 3 1191-1198 (1962)
  * [8] A. Edelman, Eigenvalues and Condition Numbers of Ranom Matries. Thesis.
  * [9] A. Edelman, B. D. Sutton, The Beta-Jacobi Matrix Model, the CS Decomposition, and Generalized Singular Value Problems. Found. of Comp. Math., Vol 8 (2008), pp. 259–285.
  * [10] P. J. Forrester, Log-gases and random matrices (LMS-34). Princeton University Press. 2010.
  * [11] C. D. Godsil. Algebraic combinatorics. Chapman & Hall/CRC, 1993\.
  * [12] C. D. Godsil, I. Gutman, On the matching polynomial of a graph. In L. Lovász and V. T. Sós, ed., Algebraic Methods in graph theory, Vol. I of Colloq. Math. Soc. János Bolyai, 25, pages 241-249. János Bolyai Math. Soc., 1981.
  * [13] G. H. Golub and C. F. Van Loan, Matrix computations, 4th Edition, JHU Press, 2013.
  * [14] V. Gorin, A. W. Marcus, Crystallization of random matrix orbits. Intern. Math. Res. Notices 2020.3 (2020): 883-913.
  * [15] A. Gribinski, A. W. Marcus, Existence and polynomial time construction of biregular, bipartite Ramanujan graphs of all degrees, arXiv: 2108.02534.
  * [16] W. Helton, V. Vinnikov, Linear matrix inequality representation of sets. Comm. on Pure and Appl. Math., 60(5), 2007: 654-–674.
  * [17] A. W. Marcus, A class of multivariate polynomial convolutions (and applications), arXiv: 2010.08996
  * [18] A. W. Marcus, Discrete Unitary Invariance, arXiv: 1607.06679
  * [19] A. W. Marcus, Polynomial convolutions and (finite) free probability, preprint.
  * [20] Mehta, M. L., Random Matrices. 3rd ed. Vol. 142. Amsterdam: Elsevier Academic, 2004.
  * [21] B. P. Mirabelli, Hermitian, Non-Hermitian and Multivariate Finite Free Probability, Ph.D. Thesis (2020).
  * [22] T. Tao, public communication.
https://terrytao.wordpress.com/2017/10/17/heat-flow-and-zeroes-of-polynomials/

  * [23] C. F. Van Loan, Generalizing the Singular Value Decomposition. SIAM J. Numer. Anal. 13.1 (1976): 76–83.
  * [24] G. Szegő. Orthogonal polynomials. Amer. Math. Soc.; 1939.