Investigation of the Product of Random Matrices and Related Evolution Models

Abstract

In this paper, we study the phase structure of the product of D * D order matrices. In each round, we randomly choose a matrix from a finite set of d matrices and multiply it with the product from the previous round. Initially, we derived a functional equation for the case of matrices with real eigenvalues and correlated choice of matrices, which led to the identification of several phases. Subsequently, we explored the case of uncorrelated choice of matrices and derived a simpler functional equation, again identifying multiple phases. In our investigation, we observed a phase with a smooth distribution in steady-state and phases with singularities. For the general case of D-dimensional matrices, we derived a formula for the phase transition point. Additionally, we solved a related evolution model. Moreover, we examined the relaxation dynamics of the considered models. In both the smooth phase and the phase with singularities, the relaxation is exponential. The superiority of relaxation in the smooth phase depends on the specific case.

1. Introduction

The investigation of the product of Random Matrices is a significant area in modern statistical physics concerning disordered systems and data science [1,2,3,4,5,6,7,8]. This field of mathematics is well-established and has connections to quantum condensed matter, Hidden Markov models (an essential tool in data science) [9,10,11,12,13], genetic network population dynamics [14], evolution in dynamically changing environments [15,16,17,18,19,20,21], and quantum resetting [22]. Our recent study focused on the phase structure of the product of correlated random matrices with D = 2 [23]. We discovered a transition similar to the localized-delocalized transition previously observed in evolution models [14,24]. The choice of the matrix in the current round can either be correlated, depending on the matrix chosen in the previous round of iteration, or uncorrelated.

During the investigation of the product of random matrices, a key objective is to calculate the maximal Lyapunov index. This can be achieved using the formalism described in [12,13]: by introducing a vector and multiplying it at each round by a randomly chosen matrix type for that round. This process generates an n-th vector at the n-th round. The norm of the n-th vector grows exponentially with n, and the coefficient of this exponential growth corresponds to the maximal Lyapunov index.

In the case of correlated matrix type choice, a functional equation has recently been derived to describe the probability distribution for the vector mentioned earlier. A similar equation has also been derived for the Hidden Markov model. Further investigation of this equation revealed the identification of a phase transition in the steady-state distribution for the rotating reference vector. This transition is directly related to the phase transition found in [14]. In [14], transitions occur between different types of environments in infinitesimal periods of time, resulting in a system of ordinary differential equations (ODEs) instead of a functional equation.

In our recent investigation, we focused on the product of correlated d = 2 matrices of order D = 2, with correlation in the type choice. In Section 2, we identify the phase structure for the case of uncorrelated matrices. Additionally, we derive formulas to determine the phase transition points in the case of high-order matrices. In Section 3, we delve into the related evolution on a fluctuating fitness landscape, including its relaxation dynamics.

2. The Product of Matrices

2.1. The Uncorrelated Matrices

We have two matrices, $M_1\equiv A$ and $M_2\equiv B$. We have a probability ${\pi }_l$ for the choice of the matrix type x. We choose some initial vector $\mathbf{X}$, then multiply it by the randomly chosen matrix, obtaining the updated version of the vector $\mathbf{X}$. We define the distribution $P(n,\mathbf{X})$ at the n-th round, and $\int dXP(n,X)=1$. We have a vector $\mathbf{X}$ with a distribution $P(n,\mathbf{X})$ at the n-th round, then at the next step we obtain

$$\begin{array}{c}\hfill p(n+1,\mathbf{X})=\sum _l{\pi }_{i}P(n,M_l^{-1}\mathbf{X})\left|\mathbf{det}{\left({\mathbf{M}}_{\mathbf{l}}\right)}^{-\mathbf{1}}\right|\end{array}$$

(1)

For the case of $2\times 2$ matrices, we can parameterize the two-dimensional vector via angular parameter $\alpha$.

$$\begin{array}{c}\hfill (x_1,x_2)=d(cos\left(\alpha \right),sin\left(\alpha \right))\end{array}$$

where d is a norm of the vector.

If we have a vector with angle $\alpha$ after multiplying it by matrices $M_l$, then the matrix A will move it to the point $f^{-1}\left(\alpha \right)$, where $f^{(-1)}$ is the inverse function of f. It means that matrix A acts as a transformation that takes the vector from the angle $\alpha$ to its inverse according to the function f.

$$\begin{array}{c}\hfill \alpha ->f_l^{-1}\left(\alpha \right)\end{array}$$

(2)

where

$$\begin{array}{c}\hfill f_1^{-1}\left(\alpha \right)=arccos(\frac{(A_{11}cos\left(\alpha \right)+A_{21}sin\left(\alpha \right))}{\sqrt{{(A_{11}cos\left(\alpha \right)+A_{21}sin\left(\alpha \right))}^2+{(A_{12}cos\left(\alpha \right)+A_{22}sin\left(\alpha \right))}^2}})\\ \hfill f_2^{-1}\left(\alpha \right)=arccos(\frac{(B_{11}cos\left(\alpha \right)+B_{21}sin\left(\alpha \right))}{\sqrt{{(B_{11}cos\left(\alpha \right)+B_{21}sin\left(\alpha \right))}^2+{(B_{12}cos\left(\alpha \right)+B_{22}sin\left(\alpha \right))}^2}})\end{array}$$

(3)

then the functional Equation (1) transforms to

$$\begin{array}{c}\hfill P(n+1,\alpha )=\sum _l\left|f_l^{\prime }\left(\alpha \right)\right|P\left(f_l\left(\alpha \right)\right){\pi }_l\end{array}$$

(4)

Equation (4) is derived for general assumptions.

For the numerical investigation (see the Appendix B for the details of the numerical algorithm) it is convenient to look at the probability distribution with fixed types $P_l(n,\alpha )$,

$$\begin{array}{c}\hfill P_l(n+1,\alpha )=f_l^{\prime }\left(\alpha \right)|P(n,f_l\left(\alpha \right)){\pi }_l,\\ \hfill P(n,\alpha )=P_1(n,\alpha )+P_2(n,\alpha )\end{array}$$

(5)

For the real eigenvalues case, we assume a steady state solution $P\left(\alpha \right)=P_1\left(\alpha \right)+P_2\left(\alpha \right)$ for Equation (4)

$$\begin{array}{c}\hfill P\left(\alpha \right)=\sum _l\left|f_l^{\prime }\left(\alpha \right)\right|P\left(f_l\left(\alpha \right)\right){\pi }_l\end{array}$$

(6)

First, we need to identify the support $[{\alpha }_1,{\alpha }_2]$ where the solution of Equation (4) is defined. We define ${\alpha }_1,{\alpha }_2$ as the angles of the eigenvectors that are associated with the maximal eigenvalues:

$$\begin{array}{c}\hfill {\alpha }_1=f_1\left({\alpha }_1\right),\\ \hfill {\alpha }_2=f_2\left({\alpha }_2\right)\end{array}$$

(7)

We assume that the angles $[{\alpha }_1,{\alpha }_2]$ are the special points for the functions functions $P_1,P_2$.

For the steady state case of Equation(4), we expand the function $f_l\left(\alpha \right)$ near the point ${\alpha }_l$

$$\begin{array}{c}\hfill f_1\left(\alpha \right)\approx {\alpha }_1+k_1(\alpha -{\alpha }_1),f_2\left(\alpha \right)\approx {\alpha }_2+k_2(\alpha -{\alpha }_2)\end{array}$$

(8)

We look an Ansatz

$$\begin{array}{c}\hfill {\rho }_1\left(\alpha \right)=a_1{|\alpha -{\alpha }_1|}^{-1+n_1},\\ \hfill {\rho }_2\left(\alpha \right)=a_2\left(\right|\alpha -{\alpha }_2{|}^{-1+n_2}\end{array}$$

(9)

near the points ${\alpha }_1$ and ${\alpha }_2$.

Let us define $n_1,n_2$. Looking the solutions with the scaling behavior near the corresponding critical points

$$\begin{array}{c}\hfill P_1\left(\alpha \right)\sim (n_1+1){(\alpha -{\alpha }_1)}^{n_1},\\ \hfill P_2\left(\alpha \right)\sim (n_2+1){(\alpha -{\alpha }_2)}^{n_2}\end{array}$$

(10)

When the transition rates tend to 0, then our Ansatz Equation (10) transforms to $\delta (\alpha -{\alpha }_l)$

We obtain from Equation (7)

$$\begin{array}{c}\hfill 1={\pi }_1{\left(k_1\right)}^{n_1},\\ \hfill 1={\pi }_2{\left(k_2\right)}^{n_2}\end{array}$$

(11)

Thus

$$\begin{array}{c}\hfill n_1=-\frac{log\left({\pi }_1\right)}{log\left(k_1\right)},n_2=-\frac{log\left({\pi }_2\right)}{log\left(k_2\right)}\end{array}$$

(12)

The critical cases are when $n_1=1$.

$$\begin{array}{c}\hfill {\pi }_1{k}_1=1\end{array}$$

(13)

and

$$\begin{array}{c}\hfill {\pi }_2{k}_2=1\end{array}$$

(14)

We derived Equations (13) and (14) by assuming a scaling Ansatz as given in Equation (10). Although the scaling Ansatz is not observed in all situations according to our numerical simulations, Equations (13) and (14) are still well-confirmed and hold true.

For simplicity, we use the following representation for the matrices

$$\begin{array}{c}\hfill A=|{\alpha }_1>|{\alpha }_1{>}^{*}+{\lambda }_1|\frac{\pi }{2}+{\alpha }_1>|\frac{\pi }{2}+{\alpha }_1{>}^{*},\\ \hfill B=|{\alpha }_2>|{\alpha }_2{>}^{*}+{\lambda }_2|\frac{\pi }{2}+{\alpha }_2>|\frac{\pi }{2}+{\alpha }_2{>}^{*}\end{array}$$

(15)

where $|{>}^{*}$ means a transposition of the vector. We denote also $k_1=1/{\lambda }_1,k_2=1/{\lambda }_2$

Then we have

(16)

Figure 1 and Figure 2 illustrate the phase transition. In our investigation, we choose both matrices with equal probabilities.

2.2. The Relaxation Dynamics

During our investigation, we also studied the relaxation dynamics and confirmed that there is an exponential fast relaxation to the steady state in both phases. In these phases, the square distance to the steady state decreases approximately as $exp[-2rt]$. Although we were not able to derive a general expression for the relaxation dynamics in all cases, we observed that when we choose ${\lambda }_1={\lambda }_2=\lambda$, it is reasonable to consider r as the gap between the maximal and next eigenvalues.

$$\begin{array}{c}\hfill r=1-\lambda \end{array}$$

(17)

where the maximal eigenvalue is 1 (the probability balance condition) and r is the gap between the maximal and the next to the maximal eigenvalues, $\lambda \equiv max[{\lambda }_1,{\lambda }_2]$ see Equation (15). Figure 3 gives the results of numerics for the r.

We also investigated the relaxation dynamics in case of correlated choice of the matrices: at the next round, we have the current choice of the matrix type with a probability q and an alternative choice with a probability $1-q$. Figure 4 gives the results of the numerics for the r.

2.3. The Identification of the Phase Transition in Higher Dimensions

Consider the case when there are $d>2$ order matrices. For the sake of simplicity, we take the following representation for the matrices

$$\begin{array}{c}\hfill A=|V_1>|V_1{>}^{*}{\lambda }_1+{\lambda }_{11}|V_{11}>|V_{11}|{>}^{*}+{\lambda }_{12}|V_{12}>|V_{12}{>}^{*}\\ \hfill B=|V_2>|V_2{>}^{*}{\lambda }_2+{\lambda }_{21}|V_{21}>|V_{21}{>}^{*}+{\lambda }_{22}|V_{22}>|V_{22}{>}^{*}\end{array}$$

(18)

The Matrix A has a maximal eigenvalue ${\lambda }_1$, the next eigenvalues are ${\lambda }_{11},{\lambda }_{12}$, and ${\lambda }_1>{\lambda }_{11}>{\lambda }_{12}$. The Matrix B has a maximal eigenvalue ${\lambda }_2$, the next eigenvalues are ${\lambda }_{21},{\lambda }_{22}$ and ${\lambda }_2>{\lambda }_{21}>{\lambda }_{22}$.

If we take a vector $|{\mathbf{V}}_{\mathbf{1}}>+\delta |{\mathbf{V}}_{\mathbf{12}}>$, then multiplying by A we obtain ∼$|{\mathbf{V}}_{\mathbf{1}}>+\delta \frac{{\lambda }_{12}}{{\lambda }_1}|{\mathbf{V}}_{\mathbf{12}}>$, thus the distance to the ${\mathbf{V}}_{\mathbf{1}}$ decreased ${\lambda }_1/{\lambda }_{12}$ times. We have that there is

$$\begin{array}{c}\hfill \frac{{\lambda }_1}{{\lambda }_{12}}\end{array}$$

(19)

times attraction. Thus for d = 2 case we identify $f_1^{\prime }\left(\mathbf{ff}\right)={\lambda }_1/{\lambda }_{12}$.

Exact critical point corresponds to the limit

$$\begin{array}{c}\hfill \pi |\frac{{\lambda }_1}{{\lambda }_{12}}|=1\end{array}$$

(20)

When

$$\begin{array}{c}\hfill \pi |\frac{{\lambda }_1}{{\lambda }_{12}}|>1\end{array}$$

(21)

The steady-state distribution is expected to be focused on an eigenvector related to the maximal eigenvalue. We derived Equations (20) and (21) without making any scaling assumption. Our derivation is on a qualitative level, in contrast to the strict derivation of Equations (13) and (14) after using the Ansatz (10). Currently, the state of the system is described by a d-dimensional unit vector $\mathbf{v}$.

In the Figure 5 we take $V_1\sim (1,0,0),V_{11}=(0,1,0),V_{12}=(0,0,1)$, $V_2\sim (1,1,1),V_{21}\sim (51,-1,0),V_{22}\sim (1,1,-2)$. Our numeric confirms the existence of phase transition according to Equations (21) and (22). We calculate numerically the distribution via $x=(1+\mathbf{v}.{\mathbf{V}}_{\mathbf{1}})/\mathbf{2}.$ for the case with the A matrix choice, and via $x=(1+\mathbf{v}.{\mathbf{V}}_{\mathbf{2}})/\mathbf{2}.$ for the alternative case.

In the d = 2 case, we observed non-ergodicity with several separated domains. However, in the current case, we do not observe such non-ergodic behavior with separated domains. This difference in behavior is attributed to the one-dimensional character of the functional equation in the current scenario. In summary, the change in the dimensionality of the functional equation has an impact on the presence or absence of non-ergodicity and separated domains in the system.

Although we are unable to derive a steady-state distribution for $\alpha$ similar to the one described by Equation (5), we can still measure some distribution numerically during the iteration process. By doing so, we can identify the phase transition by observing the singularities of the distribution. Even without an analytical expression for the steady-state distribution, numerical methods can provide valuable insights into the system’s behavior and help identify important transitions and patterns.

2.4. Complex Eigenvalues Case

We assume that the support of the steady state distribution can be larger than $[0,2\pi ]$. We take the following matrix

(22)

and

$$\begin{array}{c}\hfill f_3^{-1}\left(\alpha \right)=arccos(\frac{cos\left(\alpha \right)+csin\left(\alpha \right)}{\sqrt{1+c^2}})\end{array}$$

(23)

In the case where one matrix has a real eigenvalue ${\alpha }_1$, while another matrix has complex eigenvalues, we observe a transition in both distributions, as described by Equation (13), as shown in Figure 6.

The presence of complex eigenvalues in one matrix introduces additional dynamics and can lead to interesting transitions in the distributions. The interplay between real and complex eigenvalues influences the behavior of the system and can result in distinct phase transitions.

3. Investigation of Related Evolution Model

In the previous sections, we have identified various qualitative behaviors for the steady-state distributions. These behaviors arise due to the interplay between random transitions between different types of matrices and forces of selection (attraction towards the attractor state). Now, in this section, we proceed to investigate a related model involving the evolution of a population with two genotypes, where there exist random transitions between two landscapes.

The Crow–Kimura evolution model [25,26] considers a population with a genome length L and two alleles for each gene. The mutation rate is represented by µper genome. In this model, the growth rate of a genotype is equal to its fitness. We choose J as the growth rate of the “peak” genotype, which is the genotype with a single peak fitness [27]. The fitnesses of other genotypes are set to zero. In this model, the evolution of the population is driven by mutation and natural selection, where genotypes with higher fitness (higher growth rates) have a greater chance of survival and reproduction. The model allows us to investigate the dynamics of population evolution in the presence of mutation and fitness differences among genotypes.

We have the following equation [27] for the fraction of replicators with the high fitness sequence

$$\begin{array}{c}\hfill dx/dt=f\left(x\right)\\ \hfill f\left(x\right)=((J-\mu )x-J{x}^2)\end{array}$$

(24)

where x is the fraction of the peak genotype in the population.

Let us consider two fitness landscapes with $\mu =1$ and $J=k_1$ for the first case and $J=k_2$ for the second case.

Let us consider two fitness landscapes where the mutation rate is $\mu =1$. In the first case, the growth rate $J=k_1$, and in the second case, $J=k_2$. We now introduce random transitions between the landscapes with the rate $\alpha$. To describe the process stochastically, we introduce the probability distributions $p_1(x,t),p_2(x,t)$, where x represents the genotype and t represents time.

We can derive the following system of equations, similar to the one found in [14]:

$$\begin{array}{c}\hfill \frac{d{p}_1(x,t)}{dt}=-d\frac{p_1(x,t)f_1\left(x\right)}{dx}-\alpha p_1(x,t)+\alpha p_2(x,t),\\ \hfill f_1\left[x\right]=((k_1-1)x-k_1{x}^2)\\ \hfill \frac{d{p}_2(x,t)}{dt}=-d\frac{p_2(x,t)f_2\left(x\right)}{dx}-\alpha p_2(x,t)+\alpha p_1(x,t)\end{array}$$

(25)

The support where the steady state solution is nonzero is $[x_1,x_2]$, $x_1=1-1/k_1,$$x_2=1-1/k_2$.

Consider the steady-state distribution. We again assume a scaling Ansatz similar to the one in Equation (10):

$$\begin{array}{c}\hfill p_1\left(x\right)\sim (n_1+1){(x-x_1)}^{n_1},\\ \hfill p_2\left(x\right)\sim (n_2+1){(x-x_2)}^{n_2}\end{array}$$

(26)

By applying the Ansatz to Equation (25) and near the critical point, we make the replacement $f_l\left(x\right)=-(k_l-1)(x-x_l)$, where $k_l$ is the growth rate of genotype l and $x_l$ is the steady state position of genotype l. This approximation allows us to simplify the equation and analyze the behavior of the system near the critical point.

$$\begin{array}{c}\hfill n_1+1=\frac{\alpha }{k_1-1},\\ \hfill n_2+1=\frac{\alpha }{k_2-1}\end{array}$$

(27)

If we put $n_l=0$, we can derive the transition points for the system.

$$\begin{array}{c}\hfill 1=\frac{\alpha }{k_1-1},\\ \hfill 1=\frac{\alpha }{k_2-1}\end{array}$$

(28)

When $n_1>0$ or $n_2$, we observe singularities in the distribution of $p_l\left(x\right)$ near the point $x_l$. This indicates that the system is in a localized phase. On the other hand, if $n_1=n_2=0$, we do not have singularities, and the system is in a de-localized phase. These phases, characterized by the presence or absence of singularities in the distributions of $p_l\left(x\right)$, were first discovered in [14]. The presence of singularities indicates localization, while the absence of singularities signifies the de-localization of the population in the respective phases.

Let us investigate the model numerically. Maintaining probability balance is crucial in the numerical algorithm at each step of the iteration. This means that the probabilities of all possible events must sum up to 1. It ensures the correctness and accuracy of the algorithm and prevents numerical errors or discrepancies. Please refer to the Appendix B section for further details on how the probability balance is upheld throughout the numerical iteration process. For the numerical algorithm, it is principal to hold the probability balance at any step of iteration, see the Appendix B.

Instead of discretization of the equation, a related random walk model on the strip is formulated, with transitions between neighbors on each chain and between two chains. The model considers L points in the interval, and the probabilities at the points $x_l$ in the grid are defined using two sets of probabilities, $p_l^a,p_l^b$. The numerical simulations have been performed, and the results indicate that the relaxation dynamics strongly depend on the phase of the system, as shown in Figure 7.

4. Conclusions

In this current article, our investigation focused on the product of random matrices and its related evolution model. The study aimed to understand the behavior and dynamics of the system under consideration. In this article, we formally define the “localized phase” as a phase characterized by a singularity in the steady-state distribution and the “de-localized phase” as a phase with a smooth distribution. This definition provides a rigorous classification of the phases based on the behavior of the steady-state distributions. Upon analysis, we identify that the de-localized phase exhibits a wide distribution, indicating a broader range of possible outcomes for the system. In our investigation, we numerically studied the product of uncorrelated discrete random matrices and identified a phase transition between localized and de-localized phases. The derived master Equation (3) for this case was found to be simpler than the one for correlated random matrices. By considering matrices with real eigenvalues and complex eigenvalues separately, we were able to obtain analytical formulas for the transition points, including a formula (Equation (19)) for the appearance of singularities in the phase space.

For the case of 2 × 2 matrices, we observed a series of transitions, as shown in Figure 1. However, we did not find any indication that different phases in the statics have different expressions for the Lyapunov indices. Regarding the relaxation dynamics, we investigated both the correlated and uncorrelated cases. We found that there is no transition in the dynamics, but relaxation is faster in the singular phase. Despite obtaining answers to some questions related to the random matrix problem, certain issues remain unsolved. Two significant problems are determining the conditions under which relaxation in the localized phase is faster than in the de-localized phase and identifying cases of explicit non-ergodic behavior in the steady-state distribution. Furthermore, we investigated the related evolution model, which resembles the case of random matrix products when matrices have real eigenvalues. In this model, we verified that the relaxation dynamics are somehow faster in the non-singular phase. Although our research focused on the Master equation, we acknowledge that the existence of the phase transition may have practical implications for related models, such as Hidden Markov Models. In conclusion, our work contributes to a better understanding of the random matrix problem, but some challenges remain unresolved. We hope that future research will address these difficult problems and provide further insights into the practical implications of the phase transition in related models.

For Hidden Markov Models (HMMs), in the functional equations such as Equation (5), the transition probabilities between states depend on the current state of the system. In HMMs, the system’s state is hidden and can only be inferred from the observed output (emission) probabilities associated with each state. The transition probabilities in the functional equations determine the probabilities of transitioning between different hidden states based on the current observed state. The existence of phase transition here could have a serious impact. Indeed, the view that learning is especially successful near the phase transition point is supported by a vast amount of research and experience over the last decades. In various complex systems, such as statistical physics, machine learning, neural networks, and optimization algorithms, it has been observed that the system’s performance improves significantly when operating near a phase transition point [28].

The Kalman filter [29] is another interesting model related to the product of correlated random matrices. The Kalman filter is an essential tool in control theory and estimation, used for state estimation in dynamic systems. It deals with estimating the state of a system based on noisy measurements. In the linear case, the Kalman filter has a well-defined and efficient solution. However, in the nonlinear case, the correct solution remains an unsolved problem [30]. The nonlinear Kalman filter is an active area of research in control theory and estimation. The challenge of solving the nonlinear Kalman filter is due to the complexity introduced by the nonlinearity of the system dynamics and the uncertainty in the process and measurement noise. Nonlinear Kalman filters, such as the Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF), are commonly used as approximations, but they may have limitations and are not guaranteed to provide the globally optimal solution. Researchers continue to work on developing better algorithms and approaches to handle the nonlinear case effectively and efficiently. We hope that the methods developed here and in [31] could work for the nonlinear Kalman filter case as well.