Inverting the following positive semidefinite matrix

[graphicsRat]

I'm trying to invert the following positive semidefinite matrix, with zeros along its diagonal. Matrices of this sort that arise from scattered data interpolation with radial basis functions.

[ 0.000000 529.831737 1198.292909 529.831737 1.000000 10.000000 0.000000 ]
[ 529.831737 0.000000 529.831737 1198.292909 1.000000 0.000000 10.000000 ]
[ 1198.292909 529.831737 0.000000 529.831737 1.000000 -10.000000 0.000000 ]
[ 529.831737 1198.292909 529.831737 0.000000 1.000000 0.000000 -10.000000 ]
[ 1.000000 1.000000 1.000000 1.000000 0.000000 0.000000 0.000000 ]
[ 10.000000 0.000000 -10.000000 0.000000 0.000000 0.000000 0.000000 ]
[ 0.000000 10.000000 0.000000 -10.000000 0.000000 0.000000 0.000000 ]

The routines DPOTRI and DGETRI fail because the matrix has zero elements along its diagonal. However, GNU Octave successfully computes the inverts the matrix to be

[ 1.8034e-003 -1.8034e-003 1.8034e-003 -1.8034e-003 2.5000e-001 5.0000e-002 -2.2854e-018 ]
[ -1.8034e-003 1.8034e-003 -1.8034e-003 1.8034e-003 2.5000e-001 -1.4014e-017 5.0000e-002 ]
[ 1.8034e-003 -1.8034e-003 1.8034e-003 -1.8034e-003 2.5000e-001 -5.0000e-002 -1.1356e-017 ]
[ -1.8034e-003 1.8034e-003 -1.8034e-003 1.8034e-003 2.5000e-001 -1.2975e-017 -5.0000e-002 ]
[ 2.5000e-001 2.5000e-001 2.5000e-001 2.5000e-001 -5.6449e+002 1.6031e-015 -2.2001e-015 ]
[ 5.0000e-002 1.5053e-017 -5.0000e-002 1.1470e-017 -1.6031e-015 5.9915e+000 2.3211e-016 ]
[ 1.3878e-017 5.0000e-002 5.8566e-018 -5.0000e-002 1.4668e-015 -1.6985e-017 5.9915e+000 ]

I've verified that the product of both matrices is the identify matrix. Unfortunately, there is no way to integrate Octave into my standalone application. As such I'm still looking for a way to compute the inverse of the above matrix.

[admin]

Hi graphicsRat,
You should try dsytrf then dsytri, that will work.
I tried and got your inverse.
Julie

[Julien Langou]

LAPACK matrix inversion subroutines work in two steps.
First, you call the factorization subroutine. Second, you call the inversion subroutine.

Example:

Code: Select all

CALL DGETRF( M, N, A, LDA, IPIV, INFO )
CALL DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )


You cannot call DGETRI (nor DPOTRI) on this matrix right away, it will indeed produce garbage.

Your matrix is symmetric, but it is not positive definite so you can not use DPOTRF+DPOTRI.

Since it's symmetric, you can use DSYTRF+DSYTRI.
(Julie just implemented a faster version of DSYTRI, available in 3.3. You might want to check this out.)

Julien.

[graphicsRat]

Thanks. Here's what my implementation (DSYTRF +DSYTRI) looks like:

Code: Select all

        // _dim is the dimension of the matrix
        // _uplo is a pointer to the character 'L'
        // _A is a pointer the array of matrix data

        int info;
        int* ipiv = new int [ _dim ];
        double* work = new double [ _dim ];

        DSYTRF( _uplo , &_dim , _A , &_dim , ipiv , work , &_dim , &info );

        if( info != 0 )
        {
                cout << "Step 1 | DSYTRF failed ..." << endl;
                return info;
        }

        int lwork = (int) work[0];
        double* nwork = new double [ lwork ];

        DSYTRI( _uplo , &_dim , _A , &_dim , ipiv , nwork , &info );

        return info;


Here's the the lower matrix after DSYTRF
[ 0.000000 ]
[ 529.831737 0.000000 ]
[ 1198.292909 529.831737 0.000000 ]
[ 529.831737 1198.292909 529.831737 0.000000 ]
[ 1.000000 1.000000 1.000000 1.000000 0.000000 ]
[ 10.000000 0.000000 -10.000000 0.000000 0.000000 0.000000 ]
[ 0.000000 10.000000 0.000000 -10.000000 0.000000 0.000000 0.000000 ]

And here's the result after DSYTRI

[ 0.000000 ]
[ -0.000000 0.000000 ]
[ -0.000002 0.000007 0.000000 ]
[ 0.001891 0.000833 529.831737 0.000000 ]
[ -0.000000 -0.001891 1.000000 -0.000000 0.000000 ]
[ 0.000000 0.000000 -0.000368 0.000000 0.000000 -0.000000 ]
[ -0.000000 -0.000000 -0.000002 -10.000000 0.001895 0.000000 -0.000000 ]

It looks a very suspect. Is it correct? ...

And how do I get the inverse from this lower triangular matrix? My matrix theory isn't very strong.

[Julien Langou]

Instead of

Code: Select all

        int lwork = (int) work[0];
        double* nwork = new double [ lwork ];


write:

Code: Select all

        double* nwork = new double [ _dim ];


( Note: I did not try the change, I just read the lines of your code so I am not sure this is it )

[graphicsRat]

Thanks. I've made the change and I got the same result.

How do I go from the output of the DSYTRI to the inverse matrix?

[graphicsRat]

Hooray, DGETRF + DGETRI works!!! ... The inverse is:

[ 0.001803 -0.001803 0.001803 -0.001803 0.250000 0.050000 -0.000000 ]
[ -0.001803 0.001803 -0.001803 0.001803 0.250000 -0.000000 0.050000 ]
[ 0.001803 -0.001803 0.001803 -0.001803 0.250000 -0.050000 -0.000000 ]
[ -0.001803 0.001803 -0.001803 0.001803 0.250000 -0.000000 -0.050000 ]
[ 0.250000 0.250000 0.250000 0.250000 -564.489096 0.000000 0.000000 ]
[ 0.050000 0.000000 -0.050000 0.000000 -0.000000 5.991465 0.000000 ]
[ 0.000000 0.050000 0.000000 -0.050000 0.000000 0.000000 5.991465 ]

[graphicsRat]

Here's related question. Would DSPSV be the correct routine to use in solving linear systems Ax = b, where A is a positive semidefinite matrix -- the sort I'm trying to invert, that has zeros along its diagonal e.g. in the original post.

[Julien Langou]

Hello

1) DPOTRF + DPOTRI will not work on that matrix. DPOTRF will fail and return with INFO = 1, this is because your matrix is symmetric but it is not symmetric positive definite. (Since A(1,1) = 0, e_1^T A e_1 = 0, so the matrix is not definite. )

2) DGETRF + DGETRI should work. I am glad you got that working.

3) DSYTRF + DSYTRI should work as well, it works just fine on my machine, see quickly written C code below.

4) we should use the C interface to LAPACK when working in C. See:
http://www.netlib.org/lapack/#_standard_c_language_apis_for_lapack
This would avoid the use for writing the prototypes (there is .h file to include ready for use)
This would avoid the use for the worskpace query and workspace allocation.
All in all that would be nicer. I try to send a following post later. Feel free to do it yourself and post the code for others to look at, if you have time on your end.

Julien

Code: Select all

#include <stdlib.h>
#include <stdio.h>

extern void dsytrf_( char *uplo, int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info );
extern void dsytri_( char *uplo, int *n, double *a, int *lda, int *ipiv, double *work, int *info );

int main (){

double *A;
double *work;
int *ipiv;
int info, lwork;
int n=7;

A = (double *) malloc(49 * sizeof (double));
ipiv = (int *) malloc(7 * sizeof (int));

A[ 0]=   0.000000; A[ 7]= 0xDEADBEEF; A[14]= 0xDEADBEEF; A[21]= 0xDEADBEEF; A[28]=0xDEADBEEF; A[35]= 0xDEADBEEF; A[42]= 0xDEADBEEF;
A[ 1]= 529.831737; A[ 8]=   0.000000; A[15]= 0xDEADBEEF; A[22]= 0xDEADBEEF; A[29]=0xDEADBEEF; A[36]= 0xDEADBEEF; A[43]= 0xDEADBEEF;
A[ 2]=1198.292909; A[ 9]= 529.831737; A[16]=   0.000000; A[23]= 0xDEADBEEF; A[30]=0xDEADBEEF; A[37]= 0xDEADBEEF; A[44]= 0xDEADBEEF;
A[ 3]= 529.831737; A[10]=1198.292909; A[17]= 529.831737; A[24]=   0.000000; A[31]=0xDEADBEEF; A[38]= 0xDEADBEEF; A[45]= 0xDEADBEEF;
A[ 4]=   1.000000; A[11]=   1.000000; A[18]=   1.000000; A[25]=   1.000000; A[32]=  0.000000; A[39]= 0xDEADBEEF; A[46]= 0xDEADBEEF;
A[ 5]=  10.000000; A[12]=   0.000000; A[19]= -10.000000; A[26]=   0.000000; A[33]=  0.000000; A[40]=   0.000000; A[47]= 0xDEADBEEF;
A[ 6]=   0.000000; A[13]=  10.000000; A[20]=   0.000000; A[27]= -10.000000; A[34]=  0.000000; A[41]=   0.000000; A[48]=   0.000000;

lwork = -1;
work = (double *) malloc(1 * sizeof (double));
dsytrf_( "L", &n, A, &n, ipiv, work, &lwork, &info );
lwork = (int) work[0];
free(work);
work = (double *) malloc(lwork * sizeof (double));
dsytrf_( "L", &n, A, &n, ipiv, work, &lwork, &info );
free(work);
work = (double *) malloc(n * sizeof (double));
dsytri_( "L", &n, A, &n, ipiv, work, &info );
free(work);

int i,j;

for (i=0;i<7;i++) { for (j=0;j<=i;j++) printf("%14.7f ",A[i+j*7]); printf("\n"); }

free(ipiv);
free(A);

return 0;

}


[Julien Langou]

By the way why is the title of your post : "Inverting the following positive semidefinite matrix" ???
Do you realize that your matrix is indefinite? It has negative and positive eigenvalues.
(As consequence any "PO" routine will not work.)
Julien.

[Julien Langou]

Here's related question. Would DSPSV be the correct routine to use in solving linear systems Ax = b, where A is a positive semidefinite matrix -- the sort I'm trying to invert, that has zeros along its diagonal e.g. in the original post.

1) Your matrix is not symmetric positive semidefinite. (See previous post.)

2) We have a factorization subroutine for symmetric positive semidefinite. I do not think this concerns you but I mention this here for completeness. The name of the subroutine is DPSTRF. There is no solve associated to it since it's not clear what solving a symmetric semidefinite positive systems means ... (There are uses for the factorization subroutine however, thus it's presence in the LAPACK library.)

3) For your matrix, which is symmetric indefinite, you can either use DSYTRF+DSYTRS (and the combo driver is named DSYSV, so DSYSV=DSYTRF+DSYTRS). This is when you want to store your matrix in full format. I.e.: storing the symmetric matrix requires n^2 data storage. (See the code in the previous post where I needed 49 entries for a 7x7 symmetric matrix.) If you want to save space, you can use the packed storage. The storage is then the minimum n*(n+1)/2 (so about n^2/2, for 7x7 would be 28 instead of 49). (See data storage in the header of the subroutine.) The name of the subroutines are now DSPTRF+DSPTRS (and combo driver is named DSPSV, so DSPSV=DSPTRF+DSPTRS). You can also compute inverses in both of these formats. Note: SP routines save half the space but are NOTICEABLY slower than their SY counterparts. So if space is not an issue, stick with full format and the SY subroutines.

Hope this helps some,
Good luck,
J.

[graphicsRat]

Code: Select all

...
A[ 0]=   0.000000; A[ 7]= 0xDEADBEEF; A[14]= 0xDEADBEEF; A[21]= 0xDEADBEEF; A[28]=0xDEADBEEF; A[35]= 0xDEADBEEF; A[42]= 0xDEADBEEF;
A[ 1]= 529.831737; A[ 8]=   0.000000; A[15]= 0xDEADBEEF; A[22]= 0xDEADBEEF; A[29]=0xDEADBEEF; A[36]= 0xDEADBEEF; A[43]= 0xDEADBEEF;
A[ 2]=1198.292909; A[ 9]= 529.831737; A[16]=   0.000000; A[23]= 0xDEADBEEF; A[30]=0xDEADBEEF; A[37]= 0xDEADBEEF; A[44]= 0xDEADBEEF;
A[ 3]= 529.831737; A[10]=1198.292909; A[17]= 529.831737; A[24]=   0.000000; A[31]=0xDEADBEEF; A[38]= 0xDEADBEEF; A[45]= 0xDEADBEEF;
A[ 4]=   1.000000; A[11]=   1.000000; A[18]=   1.000000; A[25]=   1.000000; A[32]=  0.000000; A[39]= 0xDEADBEEF; A[46]= 0xDEADBEEF;
A[ 5]=  10.000000; A[12]=   0.000000; A[19]= -10.000000; A[26]=   0.000000; A[33]=  0.000000; A[40]=   0.000000; A[47]= 0xDEADBEEF;
A[ 6]=   0.000000; A[13]=  10.000000; A[20]=   0.000000; A[27]= -10.000000; A[34]=  0.000000; A[41]=   0.000000; A[48]=   0.000000;
...


So this is how you specify a lower/upper matrix i.e. you create the full matrix but only initialize the lower or upper part!

I've been creating a matrix of from an array whose length equal to the number of initialized elements!!! (This explains the triangular monstrosities in my earlier post.)

Thanks.

[graphicsRat]

Here's related question. Would DSPSV be the correct routine to use in solving linear systems Ax = b, where A is a positive semidefinite matrix -- the sort I'm trying to invert, that has zeros along its diagonal e.g. in the original post.

For your matrix, which is symmetric indefinite, you can either use DSYTRF+DSYTRS (and the combo driver is named DSYSV, so DSYSV=DSYTRF+DSYTRS). This is when you want to store your matrix in full format. I.e.: storing the symmetric matrix requires n^2 data storage. (See the code in the previous post where I needed 49 entries for a 7x7 symmetric matrix.) If you want to save space, you can use the packed storage. The storage is then the minimum n*(n+1)/2 (so about n^2/2, for 7x7 would be 28 instead of 49). (See data storage in the header of the subroutine.) The name of the subroutines are now DSPTRF+DSPTRS (and combo driver is named DSPSV, so DSPSV=DSPTRF+DSPTRS). You can also compute inverses in both of these formats. Note: SP routines save half the space but are NOTICEABLY slower than their SY counterparts. So if space is not an issue, stick with full format and the SY subroutines.

Thanks!

Just to confirm, DGESV would be the driver to use in order to solve an un-symmetric matrix (having zeros along its diagonal).

[Julien Langou]

Just to confirm, DGESV would be the driver to use in order to solve an un-symmetric matrix (having zeros along its diagonal).

correct

[Julien Langou]

I've been creating a matrix of from an array whose length equal to the number of initialized elements!!! (This explains the triangular monstrosities in my earlier post.)

You could have used DSPTRF+DSPFTRI in that case. (Packed format.) As I wrote earlier, those subroutines are really slow for large matrices.