topic Re: Pardiso: Schur complement for an SPD matrix in Intel® oneAPI Math Kernel Library

Pardiso: Schur complement for an SPD matrix

Zohar — Mon, 27 Jan 2025 11:32:30 GMT

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?

Re: Pardiso: Schur complement for an SPD matrix

Zohar — Tue, 28 Jan 2025 21:31:03 GMT

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.

Re:Pardiso: Schur complement for an SPD matrix

Shiquan_Su — Thu, 30 Jan 2025 16:11:37 GMT

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?

Re: Re:Pardiso: Schur complement for an SPD matrix

Zohar — Thu, 30 Jan 2025 21:21:56 GMT

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.

Re: Pardiso: Schur complement for an SPD matrix

Fengrui — Sat, 15 Feb 2025 01:35:45 GMT

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.

Re: Pardiso: Schur complement for an SPD matrix

Zohar — Sun, 16 Feb 2025 00:05:10 GMT

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?

Re:Pardiso: Schur complement for an SPD matrix

Shiquan_Su — Mon, 24 Mar 2025 21:38:25 GMT

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.