Pardiso: Schur complement for an SPD matrix

[Zohar]

I have the following system: Ax = b, where A is a 2x2 block matrix:

A = [H   A1  
     A1' Z ]

H is SPD (symmetric positive definite) and Z is made of zeros.

(It's a Hessian with linear constraints.)

I'm solving the system using Schur complement (iparm[35]=2).

It's fine with mtype=-2 (symm, indefinite), but the factorization fails (-3, bad input) when mtype=2 (SPD).

Just to test the code, I tried it on an SPD A, which is solved fine with mtype=2.

It seems that mtype=2 requires that A is SPD, but for Schur it should only require H to be SPD?

[Zohar]

For a refresher on Schur complement, see section 3.5 in

https://igl.ethz.ch/projects/sparse-cholesky-update/sparse-cholesky-update-paper.pdf

My system is similar to eq. (20).

Only A is factorized, and it's SPD (mtype=2).

B is decomposed using lapacke (see the example pardiso_schur.c).

I can do the whole thing by myself: i) factorize A as SPD, ii) solve for B, iii) solve a small dense system with B. The bottleneck, however, is (ii). Forward/backward substitution for a single RHS is fast. For a matrix RHS, it's slow--which pardiso overcomes somehow.

[Shiquan_Su]

Hi, Zohar:

Thanks for reporting these findings. Would you please provide more details of your test?

1, what is your OS?

2, what is your hardware?

3, what is your programming language?

4, would you please provide the source code of your test and the instruction of how you build your source code with Intel MKL?

[Zohar]

1. win11.

2. pc.

3. c++.

4. Use the mkl example "pardiso_schur.c". Feed to it the matrix:

A =
     1     0     1
     0     1     0
     1     0     0

where H = eye(2) and A1' = [1 0].
H is SPD and A is indefinte.

[Fengrui]

Let M=[A B;C D], mtype is for the type of the whole matrix M. oneMKL Pardiso will decide internally how to compute for A^(-1)B when calculating the Schur complement of A which is D-C*A^(-1)*B. The user doesn't need to (actually can not) specify the type of A.

[Zohar]

Yeah, that was the point. M has a special structure where A is SPD while M isn't. Computing inv(A) using Cholesky (for SPD) is more efficient than L'DL (for symmetric only). Thus, it's a design question: shouldn't you let the user specify the type of A, and guessing the algorithm you use, it would be more efficient?

Specifically, if A is inverted, then why care about M definiteness at all?

[Shiquan_Su · Accepted solution]

The user can use mtype=2 (SPD) in this case for Schur complement, but then iparm[35]=1 has to be set.

In this approach, the user will then solve the Schur part on his own.

If the user is using iparm[35]=2 then they cannot specify mtype=2, because PARDISO does the complete factorization and at one point it cannot proceed without pivoting for these kinds of matrices.

Also note PARDISO is not going to do LU factorization for symmetric indefinite matrices (mtype=-2) but an LDL^{T} factorization.

So, it might not be as expensive as it seems.