# Generate real random matrix with some pure imaginary eigenvalues

I want to generate a random 4x4 matrix with real entries, one or some of whose eigenvalues is the pure imaginary number I or some other imaginary number.

NOTICE The real part of one or some eigenvalues of the random real matrix is zero.

I can generate a random matrix with integer entries.

Clear["Global*"];
While[mA = RandomInteger[{-5, 5}, {4, 4}];
! MemberQ[Eigenvalues[mA], I]];

MatrixForm[mA]
Eigensystem[mA]


$$\left(\begin{array}{cccc}2 & -2 & 4 & 0 \\ -2 & -3 & -2 & -1 \\ -4 & 1 & -3 & 2 \\ 3 & 2 & 4 & 1\end{array}\right)$$

{{-3, I, -I, 0}, {{4, -12, -11, 14}, {4, -2, -3 + I, 4}, {4, -2, -3 - I, 4}, {17, -9, -13, 19}}}

But this code doesn't generate real random matrices, and I don't know if there are any errors in the code. Or is there a better way?

Clear["Global*"];
While[mA = RandomReal[{-5, 5}, {4, 4}];
! MemberQ[Eigenvalues[mA], N[I]]];

MatrixForm[mA]
Eigensystem[mA]


Some reference questions in this site.

Generate "nice" random matrix

Generate random matrix with specific eigenvalues

Update

Now I think there is an answer to this question.

Thanks for @Daniel Huber's answer. The following code is expanded according to Daniel Huber's answer. Thanks to @Nasser's answers and comments and @Yarchik's comments, the following code is to generate a matrix from pure imaginary eigenvalues according to their comments, instead of generating random matrices and then detecting whether the eigenvalues are pure imaginary.

Take the 10x10 real matrix as an example, where there are four pure imaginary eigenvalues and the other six are real eigenvalues. Since the eigenvalue of anti-symmetric matrix is pure imaginary or zero, a random anti-symmetric matrix mantisimA is generated first (there are some related problems in this station, i.e. Random real antisymmetric matrix ).

Clear["Global*"];
RandomMatrix[rank_, m_] :=
Sum[TensorProduct @@ RandomReal[{-1, 1}, {2, m}], {i, rank}];
MatrixForm[mA = RandomMatrix[4, 4]];
MatrixForm[mantisimA = Transpose[mA - Transpose[mA]]]


$$\left(\begin{array}{cccc}0 . & 0.881925 & 0.223423 & 0.329019 \\ -0.881925 & 0 . & -0.666968 & -0.746005 \\ -0.223423 & 0.666968 & 0 . & 0.136676 \\ -0.329019 & 0.746005 & -0.136676 & 0 .\end{array}\right)$$

AntisymmetricMatrixQ[mantisimA]


True

Eigenvalues[mantisimA] // Chop


{0. + 1.39773 I, 0. - 1.39773 I, 0. + 0.0484833 I, 0. - 0.0484833 I}

Then generate a random diagonal matrix mB:

MatrixForm[mB = DiagonalMatrix[RandomReal[{-1, 1}, 6]]]


$$\left(\begin{array}{cccccc}0.413869 & 0 . & 0 . & 0 . & 0 . & 0 . \\ 0 . & 0.962549 & 0 . & 0 . & 0 . & 0 . \\ 0 . & 0 . & -0.408453 & 0 . & 0 . & 0 . \\ 0 . & 0 . & 0 . & 0.278016 & 0 . & 0 . \\ 0 . & 0 . & 0 . & 0 . & 0.237626 & 0 . \\ 0 . & 0 . & 0 . & 0 . & 0 . & -0.904459\end{array}\right)$$

A block diagonal matrix mAB is generated from mantisimA and mB:

mAB = {{mantisimA, 0}, {0, mB}};

MatrixForm[mAB = ArrayFlatten[mAB]]


$$\left(\begin{array}{ccccc}0 . & 0.881925 & 0.223423 & 0.329019 & 0 \\ -0.881925 & 0 . & -0.666968 & -0.746005 & 0 \\ -0.223423 & 0.666968 & 0 . & 0.136676 & 0 \\ -0.329019 & 0.746005 & -0.136676 & 0 . & 0 \\ 0 & 0 & 0 & 0 & 0.413869 \\ 0 & 0 & 0 & 0 & 0 . \\ 0 & 0 & 0 & 0 & 0 . \\ 0 & 0 & 0 & 0 & 0 . \\ 0 & 0 & 0 & 0 & 0 . \\ 0 & 0 & 0 & 0 & 0 .\end{array}\right.\left.\begin{array}{ccccc}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 . & 0 . & 0 . & 0 . & 0 . \\ 0.962549 & 0 . & 0 . & 0 . & 0 . \\ 0 . & 0.408453 & 0 . & 0 . & 0 . \\ 0 . & 0 . & 0.278016 & 0 . & 0 . \\ 0 . & 0 . & 0 . & 0.237626 & 0 . \\ 0 . & 0 . & 0 . & 0 . & -0.904459\end{array}\right)$$

For mAB, we use a random matrix mC to do the similarity transformation to get the final matrix mat we want. Because the similarity transformation does not change the eigenvalues of the matrix, so mat looks like a random matrix. It can be seen that mat has 4 random pure imaginary number eigenvalues and 6 random real number eigenvalues, which meet our requirements:

MatrixForm[mC = RandomMatrix[10, 10]];
MatrixForm[mat = Inverse[mC] . mAB . mC]
Eigenvalues[mat] // Chop


$$\left(\begin{array}{cccccccccc}16.3385 & -41.5472 & -4.33312 & 20.4892 & 4.12538 & -6.1673 & 7.77213 & -38.2829 & 38.9472 & 23.9866 \\ 5.84092 & -14.6843 & -2.5163 & 5.77717 & 1.61527 & -1.99138 & 2.08355 & -11.9304 & 13.7347 & 9.85066 \\ 10.2421 & -27.6697 & -5.07184 & 10.6219 & 2.53251 & -7.94874 & 3.00113 & -20.4135 & 24.6155 & 20.2253 \\ -10.5361 & 24.7431 & 5.9475 & -5.92907 & -2.08304 & 7.68168 & -1.11551 & 14.7715 & -22.0036 & -16.9828 \\ 2.75025 & -7.24597 & -1.2144 & 3.09571 & 0.81912 & -2.37392 & 1.08187 & -5.75566 & 6.51882 & 4.31431 \\ -4.4096 & 11.1205 & 2.44641 & -3.8146 & -1.23057 & 2.28264 & -1.6066 & 8.30484 & -10.4105 & -6.96863 \\ 23.2956 & -59.9963 & -8.20137 & 27.3573 & 5.98617 & -10.5259 & 10.7755 & -52.3885 & 56.2111 & 35.4369 \\ -0.406186 & -1.81916 & 2.08004 & 4.85612 & 0.0409096 & -0.143971 & 1.48435 & -5.43918 & 1.31029 & 1.88427 \\ -4.71134 & 10.3438 & 0.595214 & -5.04138 & -0.942769 & 1.01936 & -2.01276 & 9.93673 & -9.89259 & -4.04313 \\ 7.18626 & -18.7654 & -2.7216 & 8.60548 & 2.05074 & -2.56859 & 3.59416 & -16.7025 & 17.8306 & 11.3804\end{array}\right)$$

{0. + 1.39773 I, 0. - 1.39773 I, 0.962549, -0.904459, 0.413869, -0.408453, 0.278016, 0.237626, 0. + 0.0484833 I, 0. - 0.0484833 I}

• You can start with some eigenvalues that you want and construct then the companion matrix mathematica.stackexchange.com/a/143181/9469 Commented Feb 11, 2023 at 22:13
• You might be better off as others are suggesting, which is start with the desired eigenvalues type and construct the matrix which would have generated these, instead of the other way around, which is keep trying different matrices until you get one which will generate the types of eigenvalues you want. Commented Feb 12, 2023 at 10:29
• Thank you. My question has been updated. Commented Feb 13, 2023 at 13:31

It is well known, that a matrix times one of its eigenvector transforms the eigenvector into a multiple of itself. Now, for a 3D rotation matrix, a vector along its axis is a eigenvector to the eigenvalue 1. However, every other real vector is NOT transformed into a real multiple of itself Therefore any other eigenvector must belong to a complex eigenvalue.

This reasoning can be carried over to any dimension: n. If n is odd, we will have one eigenvalue of: 1 and the rest are conjugate complex eigenvalues. If n is even there are only conjugate complex eigenvalues.

Note further that a generalized rotation is a matrix the leaves the length (norm) of a vector unchanged and has a determinant of 1 (no mirroring). In MMA we can get such matrices using "Orthogonalize" and check if "Det" is one.

Therefore, you can take any rotation matrix to get a matrix with complex eigenvalues. E.g. for n=4:

SeedRandom[1];
n = 4;
Until[ (Det[mat] // Chop) == 1.,
mat = RandomReal[{-1, 1}, {n, n}];
mat = Orthogonalize[mat];
];
Det[mat] // Chop
Eigenvalues[mat]

(* 1 *)
(* {-0.984172 + 0.177216 I, -0.984172 - 0.177216 I,
0.975573 + 0.219674 I, 0.975573 - 0.219674 I} *)


To get pure imaginary eigenvalues we may start from a diagonal matrix. Consider the 2x2 matrix with eigenvalues I and -I:

m1 = {{I, 0}, {0, -I}};


We may make a base change to a real matrix by:

m2 = {{I, 1}, {I, -1}};
m3 = Inverse[m2] . m1 . m2


(* {{0, 1}, {-1, 0}} *)

To get a 4x4 matrix we create a block matrix:

m41 = {{I, 0, 0, 0}, {0, -I, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 2}};{{0, 1}, {-1, 0}}


And transform it with another block matrix:

m42 = {{I, 1, 0, 0}, {I, -1, 0, 0}, {0, 0, 1, 2}, {0, 0, 2, 1}};
m43 = Inverse[m2] . m1 . m2

(* {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 7/3, 2/3}, {0, 0, -(2/3), 2/3}} *)


The eigenvalues are now:

Eigenvalues[m43]

(* {2, I, -I, 1} *)


You may further scramble the matrix m43 by a similar transform with any orthogonal matrix, because this does not change the eigenvalues:

SeedRandom[1];
m44 = Orthogonalize[RandomReal[{-1, 1}, {4, 4}]];
m45 = Inverse[m44] . m43 . m44 ;
m45 // Chop // MatrixForm


And the eigenvalues are still the same:

(* {2., 0. + 1. I, 0. - 1. I, 1.} *)

• Hi, maybe my question is not accurate. The real part of one or some eigenvalues of the random real number matrix I want to get is zero, not just complex numbers. Commented Feb 11, 2023 at 21:01
• @lotus2019 I added pure imaginary eigenvalues to my answer. Commented Feb 12, 2023 at 9:35
• Thank you. My question has been updated. Commented Feb 13, 2023 at 13:30

pure imaginary number I or some other imaginary number.

You can use FreeQ[Eigenvalues[mA], _Complex] to check for complex eigenvalues.

Clear["Global*"];
keepTrying = True;
n = 0;
While[keepTrying,
mA = RandomReal[{-5, 5}, {4, 4}];
If[Not[FreeQ[Eigenvalues[mA], _Complex]],
Print["Found one ", mA];
keepTrying = False
];
n++;
If[n > 100, keepTrying = False] (*guard*)
]


I noticed in all the trials, first one comes with complex eigenvalue. So I do not think you need a loop more than one time. But the above does a loop with a max of 100 trials.

Update

The real part of one or some eigenvalues of the random real number matrix I want to get is zero, not just complex numbers

If you mean at least one eigevalue being pure imaginary?

it could not find one matrix trying up to 100,000 iterations when using random real numbers of 4 by 4. But with random integers, it can find some.

Clear["Global*"];
keepTrying = True;
n = 0;
While[keepTrying,
(*mA=RandomReal[{-5,5},{4,4}];*)
mA = RandomInteger[{-5, 5}, {4, 4}];
If[Length[Cases[Eigenvalues[mA], x_ /; Re[x] == 0 && Im[x] != 0]] > 0,
Print["Found one at trial ", n, MatrixForm@mA, " with eigenvalues", Eigenvalues[mA]];
keepTrying = False
];
n++;
If[n > 60000, keepTrying = False] (*guard*)
]


You see the above has pure imaginary eigenvalues (3rd and 4th, and they will come as complex conjugates of each others). But when using Real numbers it could not find one. (just uncomment the code above to use RandomReal instead of RandomInteger

• Hi, maybe my question is not accurate. The real part of one or some eigenvalues of the random real number matrix I want to get is zero, not just complex numbers. Commented Feb 11, 2023 at 21:02
• @lotus2019 fyi, updated Commented Feb 12, 2023 at 3:24
• Thank you. I tried to set the WorkingPrecision -> 1 in RandomReal`, but it didn't work either. Commented Feb 12, 2023 at 4:17