Eigenvalues of pair of two large matrices

Hi All,

In my FEM code I have two square matrices (Mass and Stiffness matrices). Both are symmetric and more populated near the diagonal. Since they are quite large, the EIG command consumes a lot of time.

Is there any command that can shorten the time and compute first 10 or 30 frequencies only.

Best regards, Rehan Jamshed

Answer 1

[Matt J]

See EIGS

[Rehan Jamshed]

Thanks but the link to eigs does not give answer specifically for a pair of matrices (as asked in my original question)

I can only see one line ...

[V,D,flag] = eigs(A) also returns a convergence flag. If flag is 0 then all the eigenvalues converged; otherwise not all converged.

But with no example. Its too brief. Requesting help again.

Best regards Rehan

[Matt J]

Thanks but the link to eigs does not give answer specifically for a pair of matrices (as asked in my original question)

Does "pair of matrices" mean you are trying to solve a generalized eigenvalue problem? The syntax

eigs(A,B,k)

will find the k-largest generalized eigenvalues of sparse matrices A and B. If that's not what you're after, please clarify the problem.

[Rehan Jamshed]

My requirement is eigenvalues E and eigenvectors V out of the pair A and B for k initial few eigenvalues.

[John D'Errico]

Then read Matt's response. What about eigs is not satisfactory? Read the help.

[Rehan Jamshed]

I clarify more:

I have mass matrix, M (1000x1000) and Stiffness matrix K (1000x1000) Then I issue the command :

 [A,B]=eig(K,M)   % This results fine but with huge time

 [A,B]=eigs(K,M,10) % to get first smallest 10 frequencies but results in the following error...

 Generalized matrix B must be the same
size as A and either a symmetric
positive (semi-)definite matrix or its
Cholesky factor.

I dont understand what should be the syntax

Best regards

[Matt J]

The syntax is fine. The error message means that K and M are not both 1000x1000 as you've been assuming. Check what they really are.

[Rehan Jamshed]

The same 1000x1000 K and M matrices produce perfectly correct eigenvalues if simple eig command used. The above mentioned error appears only with eigs command.

[John D'Errico]

But how well is M conditioned? Is M nearly singular? My bet is that it is so. What is the condition number of M?

cond(M)

I can only guess that M is close enough to the edge that eigs fails, while the eig algorithm can still survive. A guess since you have not provided the matrices K and M as attachments. You can do that and save much guessing on our parts.

[Rehan Jamshed]

Here are the two matrices K and M ... ...And we need eigenvalues and eigenvector matrices

[Matt J]

On my machine eigs() handles the data with no difficulty,

>> tic; [A,B]=eigs(K_struct,M_struct,10); toc
 Elapsed time is 4.824523 seconds.

It goes even faster if you actually tell MATLAB that K and M are sparse matrices,

    >> K=sparse(K_struct); M=sparse(M_struct);

    >> tic; [A,B,flag]=eigs(K,M,10); toc, flag
    Elapsed time is 0.868296 seconds.

    flag =

         0

[Rehan Jamshed]

Well on my machine I am getting the same error...

??? Error using ==> eigs>checkInputs at 805
Generalized matrix B must be the same size as A and
either a symmetric positive (semi-)definite matrix
or its Cholesky factor.

What is wrong in my case?

P.S. I appreciate everybody here for their continued support for sparing time

[Rehan Jamshed]

Has it to do with the version ? Mine is 7.6 (2008) and the machine is 64 bit core i5

[Matt J]

Only way to tell is to run on a more recent MATLAB version...

I would first try the following, to ensure that EIGS sees a perfectly symmetric M matrix,

   [A,B, flag]=eigs(K,(M+M.')/2,10);

Answer 2

[Matt J]

Your K and M matrices, in addition to having a very pretty snow-flake pattern, both appear to be Kronecker products of a 5x5 matrix with a 243x243 matrix.

>> spy(M)

If that is indeed the case, you can make the computation even more efficient by using the EIG method in my KronProd class ( Download ). The class doesn't have an eigs() method, so you are forced to find all the eigenvectors and eigenvalues. Nevertheless, the whole thing is 10 times faster on my machine than [A,B]=eigs(sparse(K),sparse(M),10).

    Kx=rand(5); Ky=rand(243);
    Mx=rand(5); My=rand(243);

     K=KronProd({Ky,Kx});
     M=KronProd({My,Mx});

     tic;
       [A,B]=eig(K,M);
     toc
     %Elapsed time is 0.097657 seconds.

[Rehan Jamshed]

The KronProd is indeed impressive and thanks for sharing it.

I used it for random Kx, Ky, Mx and My and it worked fine. Now, as a final result, I need Natural Frequencies (Eigenvalues) and Mode Shapes (Eigenvectors) for K_struct and M_struct (Both are 375x375 this time). Here, they are Kronecker products of a 5x5 matrix with a 75x75 matrix.

The natural freq(eigenvalues) need be sorted in ascending order and first few (e.g. first 20) are important.

In the following is how I get Freq_Hz without KronProd (Vector Freq_Hz attached)

tic; [A,B]=eig(K_struct,M_struct); toc

Freq1=diag(B);

Freq1=abs(Freq1);

Freq_Hz=((sqrt(Freq1)))/(2*pi);
Freq_Hz=sort(Freq_Hz);

Next I tried using KronProd but for random Kx, Ky, Mx and My. Dont know how to get these for my matrices (K_struct and M_struct).

Finally, as said, vector Freq_Hz is needed in ascending order.

For reference I have attached the files.

[Rehan Jamshed]

For reference how I spotted matrices (5x5) and (75x75):

The whole case is a Structural Analysis Problem with Stiffness and Mass Matrices. The size is 375 Degrees of Freedom (DOFs) x 375 Degrees of Freedom for each one.

There are two beam elements each one having 75 DOF x 75 DOF. Each beam has 3 Nodes. So when connected in series, one node is shared and the total nodes are 5. The sharing is reflected in spy(M_struct) as overlapping of the two 75x75 matrices.

[Matt J]

How is it that you don't have Kx, Ky, Mx, and My? Didn't you derive K and M by executing,

   K=kron(Kx,Ky);
   M=kron(Mx,My);

If not, how are you sure that they are Kronecker products?

[Rehan Jamshed]

The process of building K_struct and M_struct is in such a way that Kx, Ky, Mx, and My are not built standalone or explicitly written.

However, Ky and My are visible as two overlapping blocks and are generally different (if I am correct at getting Ky and My).

It would be great if KronProd gives Freq_Hz similar to mine and in process building Kx, Ky, Mx, and My matrices. I could only guess their sizes through your earlier explanation.

[Rehan Jamshed]

In other words, for my problem of 375x375...

There are two Ky's and two My's each 75x75. Depending on problem they can be different and I can take them out of K_struct.

But I dont know how to extract Kx and Mx out of K_struct and M_struct.