Retrieve inverse of the Cholesky factor from MKL pardiso

[Gwz]

Hello,

I'm trying to use MKL pardiso (version MKL 2021.2) to obtain the inverse of Cholesky factor L^-1 by solving Lx = I (dentity) using phase 331 (here A = L * L' and A is well-conditioned).

But I found what I retrieve seems to be not matching up with the true L^-1 (despite ordering P). It's a bit confusing and I'm wondering if there is some parameter resulting in this phenomenon.

( The procedure was originally intended for a Lanczos algorithm where the largest eigenvalue of a matrix like L^-1 * X * L^-T is needed )

I've attached an example of 4 * 4 matrix with code reproducing this issue. Any help would be appreciated. Thank you in advance !

Best Regards.

int testpardiso(void) {
    
    int          n = 4;
    
    /*
     A =

        0.331829000000000                   0                   0  -0.000157649000000
                        0   0.033182900000000                   0  -0.022505800000000
                        0                   0   0.331829000000000   0.022470700000000
       -0.000157649000000  -0.022505800000000   0.022470700000000   0.580879000000000
    */
    int      Ap[5]  = {0, 4, 7, 9, 10};
    int      Ai[10] = {0, 1, 2, 3, 1, 2, 3, 2, 3, 3};
    double   Ax[10] = {0.33182890474490068, 0, 0, -0.00015764871901342671, 0.33182890470071952,
                      0, -0.02250575799193312, 0.33182890465891607, 0.022470695607439466,
                      0.58087859335405345,};
    double  eye[16] = {1.0, 0.0, 0.0, 0.0,
                       0.0, 1.0, 0.0, 0.0,
                       0.0, 0.0, 1.0, 0.0,
                       0.0, 0.0, 0.0, 1.0};
    
    double  aux[16] = {1.0, 0.0, 0.0, 0.0,
                       0.0, 1.0, 0.0, 0.0,
                       0.0, 0.0, 1.0, 0.0,
                       0.0, 0.0, 0.0, 1.0};
    double Linv[16] = {0.0};
    int    perm[4] = {0, 1, 2, 3};
    
    int maxfct  = 1; // Maximum number of factors
    int mnum    = 1; // The matrix used for the solution phase
    int mtype   = 2; // Real and symmetric positive definite
    int msglvl  = 0; // Print information
    int idummy  = 0; // Dummy variable for taking up space
    int phs12   = 12;
    int phs331  = 331;
    int phsfree = -1;
    int error   = 0;
    
    void *pdsWorker[64];
    
    memset(pdsWorker, 0, sizeof(void *) * 64);
    
    int params[64] = {
        
        1, /* Non-default value */ 3, /* P Nested dissection */ 0, /* Reserved          */
        0, /* No CG             */ 2, /* No user permutation */ 0, /* No overwriting    */
        0, /* Refinement report */ 0, /* Auto ItRef step     */ 0, /* Reserved          */
        8, /* Perturb           */ 0, /* Disable scaling     */ 0, /* No transpose      */
        0, /* Disable matching  */ 0, /* Report on pivots    */ 0, /* Output            */
        0, /* Output            */ 0, /* Output              */-1, /* No report         */
        0, /* No report         */ 0, /* Output              */ 0, /* Pivoting          */
        0, /* nPosEigVals       */ 0, /* nNegEigVals         */ 0, /* Classic factorize */
        0,                         0,                           1, /* Matrix checker    */
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         1, /* 0-based solve       */ 0,
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         0,                           0,
        0,                         0, /* No diagonal         */ 0,
        0,                         0,                           0,
        0,                         0,                           0,
        0
    };
    
    pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phs12, &n,
            Ax, Ap, Ai, perm, &idummy, params, &msglvl,
            eye, aux, &error);
    
    printf("Permutation: \n");
    for (int i = 0; i < n; ++i) {
        printf("%d, ", perm[i]);
    }
    printf("\n");
    
    if (error) {
        printf("Error code: %d \n", error);
        goto cleanup;
    }
    
    pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phs331, &n,
            Ax, Ap, Ai, &idummy, &n, params, &msglvl,
            eye, Linv, &error);
    
    if (error) {
        printf("Error code: %d \n", error);
        goto cleanup;
    }
    
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            printf("%10.6e, ", Linv[i + n * j]);
        }
        printf("\n");
    }
    
cleanup:
    pardiso(pdsWorker, &maxfct, &mnum, &mtype, &phsfree, &n,
            NULL, NULL, NULL, &idummy, &idummy, params, &msglvl,
            NULL, NULL, &error);
    
    return error;
}

Testing environment:

MacBook Pro Apple Silicon 8GB RAM

Pardiso from MKL 2021.2 built using Rosetta

[VidyalathaB_Intel]

Hi,

Thanks for reaching out to us.

As per the documentation, Ax[] Array should contain the non-zero elements but in the provided code we could see that the array contains 0 entries as well which might be the issue with the mismatch of the results.

So please refer to the below link for more details regarding the arguments of the pardiso call.

https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso.html

You can also look at the sample codes which are present in the installed directory of mkl for better understanding.

Regards,

Vidya.

[Gwz]

Hello,

Thank you for your prompt reply. I believe that the answer is helpful but I still have a few concerns.

The zero elements in Ax[] array are numerical zeros due to cancellation instead of structural ones (there are consecutive solves where Ax[] is all non-zero). I have also tried perturbing the zeros by 1e-10 and the result seems unchanged.

I'm looking more closely into the documentation and will be grateful if there is any more comment available.

Best regards,

Gwz

[VidyalathaB_Intel]

Hi,

>>what I retrieve seems to be not matching up with the true L^-1

Could you please share with us the expected results and the results that you are getting so that we could verify them from our end?

Regards,

Vidya.

[Gwz]

Hello,

Thank you again for your help and for clarity I'm using the output format from MATLAB. 

The full data matrix A is given by

A =

    0.3318         0         0   -0.0002
         0    0.0332         0   -0.0225
         0         0    0.3318    0.0225
   -0.0002   -0.0225    0.0225    0.5809

and the output L^-1 I extracted from Pardiso ( by Phase 331 ) is 

Linv =

    1.7360         0         0    0.0005
    0.0000    1.7383         0    0.0673
   -0.0000   -0.0046    1.7383   -0.0674
         0         0         0    1.3121

From my understanding I should have A *  transpose(Linv) * Linv = I but

>> A * Linv' * Linv

ans =

    1.0000    0.0000    0.0000    0.0000
   -0.0000    0.0976    0.0024   -0.0351
    0.0000   -0.0000    1.0000   -0.0000
   -0.0000   -0.0000    0.0000    1.0000

seems still a bit far from identity.

Best Regards,

Gwz

[VidyalathaB_Intel]

Hi,

Thanks for your patience.

We are able to reproduce your issue. we are working on it and we will get back to you soon.

Regards,

Vidya.

[Ruqiu_C_Intel · Accepted solution]

Hello,

It seems to be an issue on user side. The A matrix in the reproducer (lines 61-64) does not match the A matrix actually used in the code (lines 6-11). In particular, the (2,2) element is wrong. It is 0.033182900000000, while in the code it is 0.33182890470071952.

If use the same A as in the code and then multiply that by Linv, we do get essentially the identity matrix as expected - the (2,2) element of that matrix is now 1 instead of 0.0976.

Regards,

Ruqiu

[Gwz]

Thanks for your reply!

I have noticed the issue and there's no problem now.

Best regards,

Gwz