Inverting the following positive semidefinite matrix

Open discussion regarding features, bugs, issues, vendors, etc.

Inverting the following positive semidefinite matrix

Postby graphicsRat » Mon Apr 04, 2011 9:36 am

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.
graphicsRat
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby admin » Mon Apr 04, 2011 11:00 am

Hi graphicsRat,
You should try dsytrf then dsytri, that will work.
I tried and got your inverse.
Julie
admin
Site Admin
 
Posts: 486
Joined: Wed Dec 08, 2004 7:07 pm

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Mon Apr 04, 2011 11:02 am

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.
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Tue Apr 05, 2011 2:46 pm

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.
graphicsRat
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Tue Apr 05, 2011 3:45 pm

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 )
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Tue Apr 05, 2011 4:22 pm

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
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Tue Apr 05, 2011 4:54 pm

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
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Tue Apr 05, 2011 5:02 pm

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.
graphicsRat
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Tue Apr 05, 2011 11:22 pm

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
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Tue Apr 05, 2011 11:24 pm

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
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Tue Apr 05, 2011 11:37 pm

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.
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Wed Apr 06, 2011 9:29 am

Julien Langou wrote:
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
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby graphicsRat » Wed Apr 06, 2011 9:56 am

Julien Langou wrote:
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).
graphicsRat
 
Posts: 84
Joined: Wed Mar 25, 2009 3:08 pm

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Wed Apr 06, 2011 10:04 am

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
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Inverting the following positive semidefinite matrix

Postby Julien Langou » Wed Apr 06, 2011 10:06 am

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.
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 2 guests