DSYTRF question

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

DSYTRF question

Postby kdallen » Tue Apr 22, 2008 4:45 pm

Hi, I'm trying to get L and D from DSYTRF but can't decipher the explanation in the comments within the source code. Is there somewhere else that explains how L and D can be extracted?
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Postby kdallen » Wed Apr 23, 2008 2:02 am

Golub has the following example:

Code: Select all
( 10 20  30 )   ( 1 0 0 ) (10 0 0 ) ( 1 2 3 )
( 20 45  80 ) = ( 2 1 0 ) ( 0 5 0 ) ( 0 1 4 )
( 30 80 171 )   ( 3 4 1 ) ( 0 0 1 ) ( 0 0 1 )


And DSYTRF yields the following factorization (0.25 and 0.062 are rounded for display):

Code: Select all
( 10    0     0 )
(  2   81     0 )
(  3 0.25 0.062 )

( 1 3 3 ) = IPIV


Can somone explain how to get the above 3 matrices from the DSYTRF factorization? Thanks.
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Postby sven » Wed Apr 23, 2008 9:47 am

The short answer is that you cannot produce your three matrices from DSYTRF. This is because DSYTRF has performed a permutation at the second step, whereas your example is based upon no permutations.

In exact arithmetic DSYTRF would produce the factorization

P A P^T = L D L^T,

where

Code: Select all
   A = ( 10  20  30  )
       ( 20  45  80  )
       ( 30  80 171  )

   P = ( 1  0  0 ),  L = ( 1   0     0 ),  D = ( 10   0  0    ).
       ( 0  0  1 )       ( 3   1     0 )       (  0  81  0    )
       ( 0  1  0 )       ( 2  20/81  1 )       (  0   0  5/81 )


The vector IPIV = ( 1 3 3 ) means that no interchange (permutation) took place at steps 1 and 3, but at step 2, a permutation in the (2,3) plane took place.

Unless you wish to do something special, you generally should not need the explicit factors as routines such as DSYTRS take the output from DSYTRF.

You can see more information in Section 4.4 of the third edition of Golub and Van Loan. Hope this helps.

Sven Hammarling.
sven
 
Posts: 144
Joined: Wed Dec 22, 2004 4:28 am

Postby kdallen » Wed Apr 23, 2008 5:11 pm

Sven,

Thanks, I understand it now based on your explanation. Now the hard part. How to construct D when IPIV has negative values?

The reason I'm trying to construct L, D, (and P) is for handling constraints in orthogonal least squares problems. In ordinary LS problems w/ constraints A = J^T J, so things work out nicely. In orthogonal least squares, A is symmetric, sometimes positive definite, but not always, so LDL^T factorization comes in real handy.
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Postby sven » Thu Apr 24, 2008 12:15 pm

Two consecutive negative values in IPIV simply mean that D has a two by two block instead of a one by one block. The actual element in L that is used to store the off-diagonal element of D is zero.

If you have access to MATLAB, the a good way to see what is going on is to use a command such as

[L,D,P] = ldl(A)

which returns L, D and P such that

P^T A P = L D L^T.

You can compare the output from MATLAB with that from LAPACK.

Sven.
sven
 
Posts: 144
Joined: Wed Dec 22, 2004 4:28 am

Postby kdallen » Fri Apr 25, 2008 4:13 pm

Sven,

Thanks, almost there. No access to MATLAB, unfortunately, but an example in Golub w/ pivoting helped. I managed to factor a 12x12 matrix (w/ all 2x2 block diagonals for D), but 2 factorizations were required, one to determine P and the 2nd to factor the permutation.

Is it possible to determine P w/o factoring the original matrix?
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Postby Julien Langou » Fri Apr 25, 2008 4:28 pm

Is it possible to determine P w/o factoring the original matrix?


Nope, the easiest (known) way to find P is to factor the original matrix.
If you find a less expensive way to get P, let us know !!!!! :wink:
Julien Langou
 
Posts: 733
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Postby kdallen » Fri Apr 25, 2008 4:52 pm

Julien,

Thanks. I was afraid of that. I did manage to find a link that explains the gory details:

http://www-curri.u-strasbg.fr/documenta ... #m86Plchri

Example 3-2. Symmetric indefinite matrix factorization

But, the method looks REAL painful. I guess the choices are:

1. Perform two factorizations. Straightforward, but 2x computations needed for the factorization.
2. Use the L = P1*L1*P2*L2, etc. method, though L is really PL. Looks very ugly
3. Forget about using LDL^T and just calculate A^-1
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Postby sven » Sun Apr 27, 2008 4:00 pm

You can use the result that
Code: Select all
L(r) P(p,q) = P(p,q) M(r),  r < p,q

where M(r) is the unit lower triangular matrix L(r) with elements l(p,r) and l(q,r) interchanged. Moving the permutations matrices in this way is what allows us to express the factorization as
Code: Select all
P^T A P = L D L^T

The same idea is used to express Gaussian elimination with partial pivoting as
Code: Select all
P A = L U

Of course here we have the minor complication of two by two blocks, but I hope this helps you to see how to readily extract the factorization from the information returned by LAPACK.

Good luck,

Sven.
sven
 
Posts: 144
Joined: Wed Dec 22, 2004 4:28 am

Postby kdallen » Wed Apr 30, 2008 1:07 am

Sven,

Thanks for all your help. I was able to come up with a concise algorithm for forming P, L, D based on element swapping of DSYTRF output. :D
kdallen
 
Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

Re:

Postby holyjoly » Mon Jul 16, 2012 6:39 pm

kdallen wrote:I was able to come up with a concise algorithm for forming P, L, D based on element swapping of DSYTRF output. :D


Hi,

I realise that it's been a while since you did this, but would you mind outlining the algorithm you came up with?

Many thanks,
Joe
holyjoly
 
Posts: 1
Joined: Mon Jul 16, 2012 5:24 pm


Return to User Discussion

Who is online

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