Computation of Large-Dimension Jordan Normal Transform via Popular Platforms

We address the practical issue that popular computation platforms like Matlab©\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\copyright }$$\end{document} are unable to perform the Jordan normal (canonical) form J and the associated transform matrix P on a high-dimension matrix. Given the knowledge of the eigenvalues and eigenvectors of n-square matrix A obtained on these platforms, we present an efficient algorithm of Jordan transform with detailed computation steps. The efficiency is demonstrated by comparing the simulation results of 11 examples with varying orders of n to that computed by the symbolic-based and capacity-limited routine [P,J]=jordan(A)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[P, J] = {\textsf {jordan}}(A)$$\end{document} in Matlab©\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\copyright }$$\end{document}. Applications in stochastic dynamics are addressed.


Introduction
The similarity between square matrices is an important equivalence relationship in linear algebra.In complex field, not all square matrices can be diagonalized.But every square matrix is similar to a Jordan normal form, which is also called Jordan matrix, through a similarity transform known as the Jordan transform.Jordan matrix is an "almost diagonal" upper triangular matrix (or lower triangular matrix in some convention).Its diagonal elements are the eigenvalues of the matrix, and the form includes diagonal matrix as a special case.As long as the Jordan canonical form of a matrix is known, the linear-algebraic nature of the matrix can be seen at a glance.Jordan matrix plays an important role not only in algebraic theory, but also in physics and in engineering control [1,2] etc.
For a set of linear-differential equations of constant coefficients, the Jordan canonical form can be employed in analyzing the properties of the solutions [3].In further 1 3 Journal of Nonlinear Mathematical Physics (2023) 30:834-842 development, Frobenius first proposed one hundred years ago that a square matrix can be decomposed into the product of two symmetric matrices A n,n = S 1 S 2 with one of the two being nonsingular [4].The proposal was rediscovered by Taussky and Zassenhaus in 1959 [5].Such a significant matrix factorization can be proved within the framework of Jordan transform [6][7][8].More recently, Kwon, Ao and Thouless proved in stochastic system near fixed point that the strict Lyapunov function can be obtained with the help of Jordan transform [9].The same can be extended to the work of Chen et al. since 2020 [10,11].
Chen's work explores the rich phenomena of nonlinear quasi limit cycles near fixed points and make connections to similar properties found in the field of artificial deep neural networks [12,13].The computational analysis requires the solution of continuous Lyapunov equation.It is highly desirable that we find another method that can be fully controlled and does not rely on the command of lyap in e.g.Matlab © when performing numerical simulations.Such a task may be accomplished by having an efficient algorithm of Jordan transform to improve the computational efficiency.The ability to overcome the limitation is of great significance in the relevant field of stochastic dynamical analysis.
There have been many studies on getting the Jordan canonical form.For example, Chi and Ye [14] studied the Jordan form by the singular value decomposition (SVD).But when the difference between the individual matrix elements becomes large, the work fails to obtain satisfying accuracy.For another example, Yu and Cheng [15] provided an algorithm for computing the Jordan canonical form using elementary transforms.The pressing problem is that on popular platforms one can not obtain the Jordan normal matrix and/or the transform matrix efficiently.In the case of Matlab © , the Jordan matrix J and the transform matrix P is given by [P, J] = (A) .But the computation is actually carried out by symbolic manipulations hence is limited to very small dimensions.For an highly anisotropic n-square matrix A with n > 12 , the command could simply fail to produce a result.Even on a block-like matrix, it can barely handle a size of n > 50 .Such a capability is useless for practical applica- tions.In this paper, we present a new algorithm to compute the Jordan transform.The success is demonstrated on e.g. the Matlab © platform.

Computational Logics
In complex field, an n-square matrix A has t distinct eigenvalues { 1 , 2 , … , t } , corresponding to s ≥ t linearly independent eigenvectors.Let the algebraic multi- plicity of eigenvalue i ( i ∈ [1, t] ) be m i in the Jordan form, with the geometric mul- tiplicity i that corresponds to the order of the eigenspace for the ith eigenvalue.In other words, eigenvalue i has i linearly independent eigenvector(s).We have ∑ t i=1 m i = n , m i ≥ i , and Namely, there exists an n th -order invertible matrix M such that M −1 AM becomes a diagonal matrix in which the diagonal elements are the eigenvalues of matrix A. But when s < n , matrix A has multiple eigenvalue(s) i , whose m i >  i .In this case, the n-square matrix A is not diagonalizable.There is an nth-order invertible matrix P satisfying The matrix J is the Jordan normal form, which is similar to the matrix A. In the above, the square matrix J ij ( i ) ( j = 1, 2, … , i ) is a Jordan block, which can be written as: where J ij ( i ) is a k j th-order square matrix, and In order to obtain the Jordan normal form of n-square matrix A, we first need to find the transform matrix P = {P 1 , P 2 , … , P t } , where P i = {P i1 , P i2 , … , P i i } , and P ij is a column vector array.In a single Jordan block J ij ( i ) , the corresponding eigen- vectors where we have dropped the subscript j on k for simplicity.Introducing a matrix N i = A − i I , Eq. ( 4) can be rewritten as: where we have dropped the subscript i on N i for the same reason.The relationship in Eq. ( 4) or ( 5) is called Jordan chain.Where u 1 is the eigenvector and {u 2 , … , u k j } are generalized eigenvectors of matrix A, which are respectively in the null-space of matrix N [i.e.kernel(N) ] and the generalized eigenspace [i.e.kernel(N m i ) ]. Accord- ing to Eq. ( 5), N k u k = 0 (and 1 < k ≤ m i ), so that N m i u k = 0 .Therefore, the general- ized eigenvector u k can be obtained from a basis of N m i with zero eigenvalue and then the eigenvector u 1 can be deduced inversely.Eventually the reversible matrix P can be found and the Jordan transform of square matrix A is realized.

Implementation of Jordan Transform
In accordance with the computational logics presented in Sect.2, a flow chat of the algorithm proposed for Jordan transform is shown in Fig. 1.For a given n-square matrix A, the specific steps are as follows.
(1) Compute the eigenvalues Λ = { i } and eigenvectors V = {v i } of A; (1) Journal of Nonlinear Mathematical Physics (2023) 30:834-842 (2) Record the eigenvalues and the corresponding algebraic multiplicity { i , m i } (i = 1, 2, … , t); (3) Judgment: If all the algebraic multiplicities equal 1, that is ∀ i, m i ≡ 1 , the matrix A is diagonalizable.In this case, P = V and J = Λ is a diagonal matrix; (4) Judgment: Else for each m i > 1 , compute the geometric multiplicity of the eigenvalue.If the algebraic multiplicities equal to their geometric multiplicities for all of them, that is ∀ m i > 1 , m i ≡ i , A is also diagonalizable.We too have P = V and the diagonal J = Λ; (5) Failing on the conditions set in steps ( 3) and (4), A is not diagonalizable.Single out those eigenvalues with m i = i , each of their "Jordan blocks" reduces to an Fig. 1 The workflow of the algorithm for Jordan transform m i th-order diagonal square matrix with the element i .In this case, P i = {v i } (multi- ple columns); (6) When m i >  i , there are 2 or more Jordan blocks corresponding to i .Let N = A − i I , N p = N m i , and record the corresponding null-space S = kernel(N) , G = kernel(N p ); (7) The column vector arrays S and G recorded in step ( 6) are respectively the space composed of eigenvectors and that of the generalized eigenvectors.From the two we can obtain the column vector array R that is not in S: (9) When R c ≥ 2 , it is necessary to check the validity of the column vector array p k .The leftist column in p k must be in the space of S so that Np k (∶, 1) = 0 (in Mat- lab © language, and p k (∶, 1) refers to the leftist column or the first column in the vector array p k ).Otherwise it is removed.If the dimension of the array of two first columns dim([p k (∶, 1), p k � (∶, 1)]) = 1 , the two vector groups p k and p k ′ are equiva- lent, corresponding to the same Jordan block.In this case, the longer array is kept: (10) The (generalized) eigenvector group P ik corresponding to each Jordan block can be so obtained through step (9).But it is necessary to further check if some eigenvectors in S are left out.If so add them to the array of p i : p i = {p i ∈ S, p i ∉ {P ik }}; (11) After steps ( 5)-( 10), the generalized eigenvector group for the ith eigenvalue is (12) Output the transform matrix P and the Jordan normal form J of n-square matrix A as

Results
The effectiveness of the new routine ([P, J] = (A)) proposed in this paper may be verified by comparing it with [P 1 , J 1 ] = (A) , the Jordan trans- form routine in Matlab © .Due to the ordering of eigenvalues, the sequence of the resulting Jordan blocks may be different.The check the difference, it is sufficient (6) Journal of Nonlinear Mathematical Physics (2023) 30:834-842 , where a ij and b ij are the matrix elements of J 1 and J respectively.Take the 6 low-order square matrices A i ( i = 1, 2, … , 6 ) as examples 1 to exam- ple 6, Then we combine them to get Example 7 and Example 8 as ) and A 8 = diag(A 7 , A 7 ) .Finally, make three ran- dom rows A � 9,10,11 = round(rand(n) * 10) with order n = 5, 10, 1000 , and let ) be examples 9~11 (We will not show A ′ 11 in this paper as n = 1000 is too large).On a Matlab © platform, the computation times for the two routines jordan and newjordan and the errors between them are shown in Table 1 (not recorded if the times exceeded 10 min).They are averaged over five calculations for each sample.

Discussions
We conclude this work with exemplary discussions relevant to real applications.Consider a set of stochastic differential equation ( 8) where x is a multidimensional vector and is a random noise.According to the structural analysis of stochastic dynamics near a fixed point published in PNAS [9], the equations can be decomposed into In the above, the matrix D is symmetric and semi-positive definite, the matrix Q is antisymmetric, and U is symmetric.On a stable fixed point and setting D to the diffusion matrix by (t) , U becomes a semi-positive matrix for the thermodynamic potential, a crucial functional for various physical analyses.
Based on the symmetry and antisymmetry of the matrices, one can obtain This is a type of continuous Lyapunov equation.For real symmetric F, Eq. ( 12) can be rewritten as Here , are the eigenvalues of F. Q are D are the matrix elements of Q and D under the representation of F.
But for general cases it is necessary to consider the Jordan transform of the F, where the matrix L −1 is the transform matrix computed in this paper and Λ is the corresponding Jordan matrix.Similarly, in the representation of matrix F, Q = LQL T , D = LDL T , the upper triangular matrix elements (  <  ) of matrix Q can be expressed as: Complete details can be found in [9].To summarize, this paper discusses how to effectively calculate the Jordan canonical form and the transform matrix of nth-order square matrix A on popular platforms like Matlab © .The eigenvalues, eigenvectors and generalized eigenvectors of an n-square A provides new insights into the Jordan transform.On a very general consideration, the new algorithm proposed realizes Jordan transform without the limitation by the order of the matrix.It quickly gives the Jordan normal form J and the corresponding transform matrix P. When dealing with practical problems, such as the high-dimensional complex, nonlinear systems presented in [10,11,16,17], the method may play a crucial role and improve the computational efficiency.are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.
where R c is the number of columns in R and the columns in array R are mutually orthogonal; (8) Construct an array of columns for each column in the column vector array R: Let p k = [r k ] , compute Nr k and add the result to p k , i.e. p k = [Nr k , p k ] (in Mat- lab © language).Repeat the procedure with r k = Nr k until Nr k = 0 (maximum m i iterations);

Table 1
The results of commands of jordan.m and newjordan.mon Matlab ©