Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

PARDISO. SCHUR Complement matrix

marcsolal
Beginner
1,855 Views

I would like to use PARDISO to compute the SCHUR complement. I am trying to understand the documentation and I have questions:

- I understand that the Schur complement matrix is obtained in the solution vector. I had a look on PARDISO 5.0 (not the Intel software) documentation and the SCHUR complement is returned as a sparse matrix. Is MKL PARDISO returning a full matrix or a sparse matrix ? Also which element is where in the matrix.

- I will probably need to access the full factorization ie the Schur complement but also the other matrices (L11, L21, U11,U12 from the documentation). I understand it is possible to compute these matrices, but how can we access those and in which form.

- do you have a code example of PARDISO for computation of the SCHUR complement and access to all the matrices.

Thanks a lot.

Marc

0 Kudos
26 Replies
Alexander_K_Intel2
1,594 Views

Hi,

Please see my comments below

Thanks,

Alex

- I understand that the Schur complement matrix is obtained in the solution vector. I had a look on PARDISO 5.0 (not the Intel software) documentation and the SCHUR complement is returned as a sparse matrix. Is MKL PARDISO returning a full matrix or a sparse matrix ? Also which element is where in the matrix.

[akalinki] Current version of MKL Pardiso support schur complement matrix in dense format only

- I will probably need to access the full factorization ie the Schur complement but also the other matrices (L11, L21, U11,U12 from the documentation). I understand it is possible to compute these matrices, but how can we access those and in which form.

[akalinki] You are correct, in current version of MKL you can handle this matrix via solving step but cannot access it's internal representation

- do you have a code example of PARDISO for computation of the SCHUR complement and access to all the matrices.

[akalinki] file pardiso_schur_c in general mkl example folder show hot to calcluate schur complement of small matrix 

Thanks a lot.

Marc

0 Kudos
marcsolal
Beginner
1,594 Views

I am sorry.  I am confused. I understand the Schur complement is a full matrix. This is fine, but I do not understand your answer about triangular matrices. You tell it can be used in the solving step. But what is the solving step doing when the Schur complement is used ? Solving for the full matrix A.X=B or something else. If the solving is done for the full matrix, it means that you need to factorize the Schur complement matrix in the solving step. Am I correct? Also, there is a remark in the documentation about using LAPACK to compute the Schur complement vector for large matrices. What do you call Schur complement vector? Is this the triangular matrices? Or the solution of the system or something else.

Sorry for asking may be naive questions.

Thanks

 

0 Kudos
Alexander_K_Intel2
1,594 Views

marcsolal wrote:

I am sorry.  I am confused. I understand the Schur complement is a full matrix. This is fine, but I do not understand your answer about triangular matrices. You tell it can be used in the solving step. But what is the solving step doing when the Schur complement is used ? 

In therms from this paper phase 331 solve system with first, lower-triangular matrix, phase 333 with third, upper-triangular matrix

marcsolal wrote:

 If the solving is done for the full matrix, it means that you need to factorize the Schur complement matrix in the solving step. Am I correct? 

Yes and no. If you set iparm[35] to 1 then only schur complement will be computed. If you set it to 2 then pardiso on factorization step will compute all factorized matrix and also will factorized Schur complement matrix. 

marcsolal wrote:

Also, there is a remark in the documentation about using LAPACK to compute the Schur complement vector for large matrices. What do you call Schur complement vector? Is this the triangular matrices? Or the solution of the system or something else.

It is part of vector full size vector that correspond to Schut complement matrix. Actually that's the element of array x for which correspond positions in perm array set to 1.

0 Kudos
marcsolal
Beginner
1,594 Views

Thanks for the answers. I think I need to be more specific. I have a (sparse) linear system with a 4 block matrix and I am trying to eliminate some dof from the system.

Let say my matrix is   A=[A11 A12;A21 A22]

My right hand side vector is B=[B1;B2]

And my system is:   A11 X1+A12X2=B1

                               A21 X1+A22X2=B2

I would like to express the system only in function of X2

I get: ( A22- A21A11-1A12)X2=B2-A21A11-1B1

The matrix on the left is the Schur complement matrix and PARDISO should give this. My problem is to be able to update the right hand side. A21A11-1B1 can be find from the triangular matrices. But if I do not have access to the triangular matrices, I need to find another way.

I hope I made myself clear.

0 Kudos
Alexander_K_Intel2
1,594 Views

Hi. Looks like there is a trick how it could be calculated without internal matrix data. I attach the photo with my calculations, and I really sorry for quality of this photo - i have only my telephone and paper from my daughter album for drawing :) 

0 Kudos
marcsolal
Beginner
1,594 Views

Now I get it. It was exactly the meaning of my original question. I wanted to know if I can use the triangular matrix independently of the Schur complement matrix. I did not understand the meaning of forward and backward substitution. So the Schur complement matrix is not factorized in steps 331 and 333, which is exactly what I need. I am assuming step 332 in case of Schur complement matrix both factorizes and solves the Schur complement matrix which is a full matrix. It is why the documentation recommends to use Lapack routines to do that.

Thanks I think I can do what I need now.

 

0 Kudos
Konstantin_K_2
Beginner
1,594 Views

Hi there,

I am trying to calculate large Schur complement for a symmetric sparse matrix using MKL PARDISO pardiso_schur_c.c program.

My parameters are:

MKL_INT n = 604034;

MKL_INT nonz = 2063839;

If I use the last 500 elements in perm vector = 1 (dimensions of S), then the S matrix is symmetric. If I increase the number of elements, say 1000, I get unsymmetric S.

What is the content of S? Does it contain the Schur complement, or the upper triangle is different from the lower triangle? Can anybody help me to understand this? Thanks in advance.

Konstantin

0 Kudos
MariaZh
Employee
1,594 Views

Konstantin K. wrote:

Hi there,

I am trying to calculate large Schur complement for a symmetric sparse matrix using MKL PARDISO pardiso_schur_c.c program.

My parameters are:

MKL_INT n = 604034;

MKL_INT nonz = 2063839;

If I use the last 500 elements in perm vector = 1 (dimensions of S), then the S matrix is symmetric. If I increase the number of elements, say 1000, I get unsymmetric S.

What is the content of S? Does it contain the Schur complement, or the upper triangle is different from the lower triangle? Can anybody help me to understand this? Thanks in advance.

Konstantin

Hello Konstantin,

Please take a look at the related topic http://software.intel.com/en-us/articles/intel-mkl-support-to-new-functionality-schur-complement, I suppose it'll be helpful for understanding how Schur complement approach is implemented in MKL Pardiso.
Regarding your question, on the output Schur complement is returned as a general dense matrix.

Best regards,
Maria

 

0 Kudos
Konstantin_K_2
Beginner
1,594 Views

Hi Maria,

 

Thank you very much for your response. I understand that the Schur matrix is a dense matrix with dimensions (n_schur,n_schur). I’ve been to the site you’ve suggested but my problem is quite different.

I have a large symmetric sparse matrix (upper triangle only). The task is to compute Schur complement on A22 , where A22 has different size. I have a check for symmetry of the Schur matrix at the end of the program. I’ve noticed that when the dimension of A22 is 500 x 500, the resulting Schur matrix is symmetric. If I increase the dimension, say 800 x 800, the Schur matrix I get is non-symmetrical. My understanding (and experience) is that if A is symmetric and can be partitioned as

|A11   A12|, A11 and A22 are invertible, then the Schur matrix, S = A22 - A21(A11) -1A12 is always

|A21   A22 |  symmetric.

My puzzle is why am I getting unsymmetrical result when the size is greater than 500? Is this a sign that something went wrong? Any thoughts? Thanks.

 

Kind Regards,

Konstantin

0 Kudos
MariaZh
Employee
1,594 Views

Hi Konstantin,

Can you please provide a reproducer, so I'll be able to investigate your case carefully?

Best regards,
Maria

0 Kudos
Konstantin_K_2
Beginner
1,594 Views

Hi Maria,

 

Thank you very much for your response and willingness to help. Much appreciated. Let me give you some background of my story so you can understand the problem better.

My example comes from genetics. We have a matrix, say A that contains coefficients of relationships between animals in a given pedigree (probabilities of common genes coming from parents). This matrix is not easy to compute because it is very large (10M x10M). However, we have a number of animals that are genotyped (their genome is scanned) and we are interested only in those animals.

Therefore, we can partition A as

|A11  A12|

|A21  A22| , where A22 contains the animals of interest. We have an algorithm that can calculate A22 relatively inexpensive. However, we need the inverse of A22. Therefore the work to get the A22_inverse is doubled. Luckily, we have a fast algorithm to compute the inverse of A for millions of animals. However, we need only the portion of A_inverse that contains the genotyped animas, e.g. A22_inverse. That’s where the Schur complement from Pardiso comes into the picture. I have the entire A_inverse and I want to calculate A22 part only. Then the test A22 * A22_inverse = I, where I is identity matrix will tell me whether I am doing the right thing.  

My problem is that I get Schur = A22_inverse as an unsymmetrical matrix, which contradicts the theory.

There is another puzzle I found yesterday. I copied the lower triangle of the Schur matrix to the upper part and I got the required result (A22 * A22_inverse = I). This suggests to me that the program keeps the Schur result in the lower part of the matrix, which contradicts to the description of the matrix in the manual.

Anyway, attached please find the following files:

Test_schur.cpp – testing program, so you can reproduce the results.

Hol.ain – the matrix A_inverse in row compressed sparse format. The first part of the file contains IA and the second JA and A.

You can disable the use of Blitz array library. It is there because I have a version of the test program that uses Blitz arrays instead of C arrays.

Sorry for the long story, but I am very excited about using this program. It will save us lots of computing time. However, I need a proof that I am doing the right thing. Once again, thanks.

Kind regards,

Konstantin

0 Kudos
MariaZh
Employee
1,594 Views

Konstantin,

Thank you very much for the answer!
I was able to reproduce your problem and I can see that there is indeed problem with a lower triangular part of the matrix, so I'm working on it.
As a workaround I suggest you to use only upper triangular part of the Schur complement matrix in your computations.
Thanks for the issue.

Best regards,
Maria

0 Kudos
Konstantin_K_2
Beginner
1,594 Views

Hi Maria,

I am glad that you’ve confirmed my findings. Just a comment: I believe that the problem is with the upper triangle, not the lower. In my other test program I copied the lower to the upper triangle and the test A22*A22_inv = I worked.

Kind regards,

Konstantin

0 Kudos
MariaZh
Employee
1,594 Views

Hello Konstantin,

You're right, of course. My apologies for the mix up.

Best regards,
Maria.

 

0 Kudos
Jens_E_
New Contributor I
1,594 Views

Hi Maria and Konstantin,

I'm also having issues with the Schur complement matrix returned from Pardiso not being symmetric when the original (sparse) matrix is symmetric. It seems that one of of the triangular parts of the returned matrix is correct. I'm only passing the one triangle of the original matrix to Pardiso.

Is it guaranteed that the triangular (seemingly correct) part of the Schur complement is indeed correct?

I'm using mkl 11.3 update 3. Is this issue fixed in Update 4?

Best and thanks,

Jens

 

0 Kudos
Alexander_K_Intel2
1,594 Views

Hi Jens,

You are right, currently only one triangular returned by Schur complement is correct and the issue still exist in last version of MKL - MKL2017u1

Thanks,

Alex

0 Kudos
marcsolal
Beginner
1,594 Views

I am willing to use the Schur functionality of PARADISO on symmetric matrices. Can you please confirm which triangle is the correct one. Am I understanding correctly that the lower triangle is the one to use?

Thanks,

Marc 

0 Kudos
MariaZh
Employee
1,594 Views

Hi Marc,

You are right, lower triangle of the output matrix is correct and should be used.

Best regards,
Maria

 

0 Kudos
marcsolal
Beginner
1,594 Views

Thanks. I am assuming Schur is not working with DSS handle. Right?

Thanks,

Marc

0 Kudos
Gennady_F_Intel
Moderator
1,448 Views

the problem is fixed in MKL 2017 u2. This update is released at Feb 22. Please check how this fix works on your side. 

0 Kudos
Reply