LU Decomposition error

LU Decomposition error

[Peter Roche]

Dear All,

I hope this isn't a repeat of an already dealt with problem, but I
couldn't see anything in the archives. I have encountered some problems
with the gsl_linalg_LU_decomp rountine.  The routine ran without exiting
with an error message, but I found that if I multiply the decomposed
Lower and Upper triangular matrices I do not get the original matrix I
started with before decomposing.  I have included a section of my code
and also a section of the original matrix (sym_denom) and the result of
matrix multiplying the decomposed lower and upper triangular matrices
(LU check).

Some of the rows are the same as the original but are not in the correct
row. For example row 4 in "LU check" has the same elements as row 2
"sym_denom"!

Has anyone else encountered similar problems or can help me with this
one.

Cheers, Peter

------------------------------------------------------------------------

  sym_denom_LU = gsl_matrix_calloc(num_ind,num_ind);
  sym_denom_Lower = gsl_matrix_calloc(num_ind,num_ind);
  sym_denom_Upper = gsl_matrix_calloc(num_ind,num_ind);

  perm = gsl_permutation_calloc(num_ind);
  gsl_matrix_memcpy(sym_denom_LU,sym_denom);
  gsl_linalg_LU_decomp(sym_denom_LU,perm,&sign);

  gsl_matrix_memcpy(sym_denom_Lower,sym_denom_LU);
  gsl_matrix_memcpy(sym_denom_Upper,sym_denom_LU);

  /* check on LU decomp */
  for(ina=0;ina<num_ind;ina++){
    for(inb=ina;inb<num_ind;inb++){
      if(inb > ina ) gsl_matrix_set(sym_denom_Lower,ina,inb,0.0);
    }
    for(inb=ina+1;inb<num_ind;inb++){
      gsl_matrix_set(sym_denom_Upper,inb,ina,0.0);
    }
  }

  /* Multiply lower/upper matrices to get original sym_denom */

gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,1.0,sym_denom_Lower,sym_denom_Upper);

------------------------------------------------------------------------


sym_denom
 1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20 -3.066652e-19   4.391243e-19
 2    -1.652853e-19  -3.923216e-20  -3.608558e-20  -1.378546e-20 -7.772215e-20   9.616245e-20
 3    -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20   7.866459e-20
 4    -6.114344e-20  -1.378546e-20  -2.635082e-20  -9.419878e-21 -1.627012e-20   3.311913e-20
 5    -3.066652e-19  -7.772215e-20  -3.128613e-20  -1.627012e-20 -2.676652e-19   2.252427e-19
 6     4.391243e-19   9.616245e-20   7.866459e-20   3.311913e-20 2.252427e-19  -3.018037e-19
 7    -1.748050e-20  -3.759400e-21  -2.671663e-20  -8.956211e-21 2.023082e-20   4.257951e-21
 8     2.556328e-19   6.230603e-20   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19
 9     4.908470e-20   1.223947e-20   3.696976e-20   1.238438e-20 3.418982e-21  -2.760516e-20


LU check - should be the same as sym_denom
 1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20 -3.066652e-19   4.391243e-19
 2     3.491029e-19   8.709227e-20   5.217473e-20   2.250474e-20 2.709399e-19  -2.325323e-19
 3     2.556328e-19   6.230603e-20   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19
 4    -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20   7.866459e-20
 5    -1.652853e-19  -3.923216e-20  -3.608558e-20  -1.378546e-20 -7.772215e-20   9.616245e-20
 6     2.847849e-19   6.186536e-20   4.510521e-20   2.035696e-20 1.565680e-19  -2.154436e-19
 7    -1.748050e-20  -3.759400e-21  -2.671663e-20  -8.956211e-21 2.023082e-20   4.257951e-21
 8    -1.493627e-19  -4.009839e-20  -1.313729e-19  -4.382593e-20 -5.583881e-21   8.188773e-20
 9     1.981133e-20   5.659666e-21   2.816777e-20   9.600430e-21 -8.449448e-21  -8.121479e-21

*********************************************************

Peter Roche
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Silver Street, Cambridge, CB3 9EW

email: p.j.p.roche@damtp.cam.ac.uk

*********************************************************

Re: LU Decomposition error

[Slaven Peles]

On Wednesday 13 November 2002 02:39 pm, Peter Roche wrote:
> couldn't see anything in the archives. I have encountered some problems
> with the gsl_linalg_LU_decomp rountine.  The routine ran without exiting
> with an error message, but I found that if I multiply the decomposed
> Lower and Upper triangular matrices I do not get the original matrix I
> started with before decomposing.  I have included a section of my code
> and also a section of the original matrix (sym_denom) and the result of
> matrix multiplying the decomposed lower and upper triangular matrices
> (LU check).
>
> Some of the rows are the same as the original but are not in the correct
> row. For example row 4 in "LU check" has the same elements as row 2
> "sym_denom"!
>
> Has anyone else encountered similar problems or can help me with this
> one.
>
> Cheers, Peter
>

Have you checked if your matrix is singular? It was reported before that LU 
decomposition does not exit with error when the input matrix is singular 
matrix. I've just had a look at function gsl_linalg_LU_decomp in lu.c in GSL 
1.2. It does check for zeros at the main diagonal after the decomposition, 
and proceeds with calculation in case there are none. But it seems to me that 
it does nothing if it finds a zero, i.e. when the input matrix is singular. I 
copy that piece of gsl_linalg_LU_decomp code below.  

======== from gsl_linalg_LU_decomp in lu.c ==================

	  ajj = gsl_matrix_get (A, j, j);

	  if (ajj != 0.0)
	    {
	      for (i = j + 1; i < N; i++)
		{
		  REAL aij = gsl_matrix_get (A, i, j) / ajj;
		  gsl_matrix_set (A, i, j, aij);

		  for (k = j + 1; k < N; k++)
		    {
		      REAL aik = gsl_matrix_get (A, i, k);
		      REAL ajk = gsl_matrix_get (A, j, k);
		      gsl_matrix_set (A, i, k, aik - aij * ajk);
		    }
		}
	    }
	}

      return GSL_SUCCESS;
    }
}
==========================================

I guess there should be an else statement at the end calling GSL_ERROR. 

Cheers,
Slaven


>   sym_denom_LU = gsl_matrix_calloc(num_ind,num_ind);
>   sym_denom_Lower = gsl_matrix_calloc(num_ind,num_ind);
>   sym_denom_Upper = gsl_matrix_calloc(num_ind,num_ind);
>
>   perm = gsl_permutation_calloc(num_ind);
>   gsl_matrix_memcpy(sym_denom_LU,sym_denom);
>   gsl_linalg_LU_decomp(sym_denom_LU,perm,&sign);
>
>   gsl_matrix_memcpy(sym_denom_Lower,sym_denom_LU);
>   gsl_matrix_memcpy(sym_denom_Upper,sym_denom_LU);
>
>   /* check on LU decomp */
>   for(ina=0;ina<num_ind;ina++){
>     for(inb=ina;inb<num_ind;inb++){
>       if(inb > ina ) gsl_matrix_set(sym_denom_Lower,ina,inb,0.0);
>     }
>     for(inb=ina+1;inb<num_ind;inb++){
>       gsl_matrix_set(sym_denom_Upper,inb,ina,0.0);
>     }
>   }
>
>   /* Multiply lower/upper matrices to get original sym_denom */
>
> gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,1.0,sym_denom_Lo
>wer,sym_denom_Upper);
>
> ------------------------------------------------------------------------
>
>
> sym_denom
>  1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20
> -3.066652e-19   4.391243e-19 2    -1.652853e-19  -3.923216e-20 
> -3.608558e-20  -1.378546e-20 -7.772215e-20   9.616245e-20 3   
> -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20  
> 7.866459e-20 4    -6.114344e-20  -1.378546e-20  -2.635082e-20 
> -9.419878e-21 -1.627012e-20   3.311913e-20 5    -3.066652e-19 
> -7.772215e-20  -3.128613e-20  -1.627012e-20 -2.676652e-19   2.252427e-19 6 
>    4.391243e-19   9.616245e-20   7.866459e-20   3.311913e-20 2.252427e-19 
> -3.018037e-19 7    -1.748050e-20  -3.759400e-21  -2.671663e-20 
> -8.956211e-21 2.023082e-20   4.257951e-21 8     2.556328e-19   6.230603e-20
>   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19 9    
> 4.908470e-20   1.223947e-20   3.696976e-20   1.238438e-20 3.418982e-21 
> -2.760516e-20
>
>
> LU check - should be the same as sym_denom
>  1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20
> -3.066652e-19   4.391243e-19 2     3.491029e-19   8.709227e-20  
> 5.217473e-20   2.250474e-20 2.709399e-19  -2.325323e-19 3     2.556328e-19 
>  6.230603e-20   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19 4  
>  -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20  
> 7.866459e-20 5    -1.652853e-19  -3.923216e-20  -3.608558e-20 
> -1.378546e-20 -7.772215e-20   9.616245e-20 6     2.847849e-19  
> 6.186536e-20   4.510521e-20   2.035696e-20 1.565680e-19  -2.154436e-19 7   
> -1.748050e-20  -3.759400e-21  -2.671663e-20  -8.956211e-21 2.023082e-20  
> 4.257951e-21 8    -1.493627e-19  -4.009839e-20  -1.313729e-19 
> -4.382593e-20 -5.583881e-21   8.188773e-20 9     1.981133e-20  
> 5.659666e-21   2.816777e-20   9.600430e-21 -8.449448e-21  -8.121479e-21
>
> *********************************************************
>
> Peter Roche
> Department of Applied Mathematics and Theoretical Physics,
> University of Cambridge, Silver Street, Cambridge, CB3 9EW
>
> email: p.j.p.roche@damtp.cam.ac.uk
>
> *********************************************************

Re: LU Decomposition error

[Brian Gough]

Peter Roche writes:
 > I hope this isn't a repeat of an already dealt with problem, but I
 > couldn't see anything in the archives. I have encountered some problems
 > with the gsl_linalg_LU_decomp rountine.  The routine ran without exiting
 > with an error message, but I found that if I multiply the decomposed
 > Lower and Upper triangular matrices I do not get the original matrix I
 > started with before decomposing.  I have included a section of my code
 > and also a section of the original matrix....

Please send a complete example program which reproduces the
problem. Thanks.

Re: LU Decomposition error

[Peter Roche]

On more careful reflection the LU decomposition seems OK.  I had simply
not taken into consideration that the marix product of the Upper and Lower
triangular matrices gives a permutation of the original matrix, that is
P.A = L.U, where A is the original matrix, L the lower triangular matrix,
U the Upper triangular and P a vector containing a record of the row
permutations carried out in the decomposition.

Cheers, Peter

*********************************************************

Peter Roche
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Silver Street, Cambridge, CB3 9EW

email: p.j.p.roche@damtp.cam.ac.uk

*********************************************************

On Wed, 13 Nov 2002, Slaven Peles wrote:

> On Wednesday 13 November 2002 02:39 pm, Peter Roche wrote:
> > couldn't see anything in the archives. I have encountered some problems
> > with the gsl_linalg_LU_decomp rountine.  The routine ran without exiting
> > with an error message, but I found that if I multiply the decomposed
> > Lower and Upper triangular matrices I do not get the original matrix I
> > started with before decomposing.  I have included a section of my code
> > and also a section of the original matrix (sym_denom) and the result of
> > matrix multiplying the decomposed lower and upper triangular matrices
> > (LU check).
> >
> > Some of the rows are the same as the original but are not in the correct
> > row. For example row 4 in "LU check" has the same elements as row 2
> > "sym_denom"!
> >
> > Has anyone else encountered similar problems or can help me with this
> > one.
> >
> > Cheers, Peter
> >
>
> Have you checked if your matrix is singular? It was reported before that LU
> decomposition does not exit with error when the input matrix is singular
> matrix. I've just had a look at function gsl_linalg_LU_decomp in lu.c in GSL
> 1.2. It does check for zeros at the main diagonal after the decomposition,
> and proceeds with calculation in case there are none. But it seems to me that
> it does nothing if it finds a zero, i.e. when the input matrix is singular. I
> copy that piece of gsl_linalg_LU_decomp code below.
>
> ======== from gsl_linalg_LU_decomp in lu.c ==================
>
> 	  ajj = gsl_matrix_get (A, j, j);
>
> 	  if (ajj != 0.0)
> 	    {
> 	      for (i = j + 1; i < N; i++)
> 		{
> 		  REAL aij = gsl_matrix_get (A, i, j) / ajj;
> 		  gsl_matrix_set (A, i, j, aij);
>
> 		  for (k = j + 1; k < N; k++)
> 		    {
> 		      REAL aik = gsl_matrix_get (A, i, k);
> 		      REAL ajk = gsl_matrix_get (A, j, k);
> 		      gsl_matrix_set (A, i, k, aik - aij * ajk);
> 		    }
> 		}
> 	    }
> 	}
>
>       return GSL_SUCCESS;
>     }
> }
> ==========================================
>
> I guess there should be an else statement at the end calling GSL_ERROR.
>
> Cheers,
> Slaven
>
>
> >   sym_denom_LU = gsl_matrix_calloc(num_ind,num_ind);
> >   sym_denom_Lower = gsl_matrix_calloc(num_ind,num_ind);
> >   sym_denom_Upper = gsl_matrix_calloc(num_ind,num_ind);
> >
> >   perm = gsl_permutation_calloc(num_ind);
> >   gsl_matrix_memcpy(sym_denom_LU,sym_denom);
> >   gsl_linalg_LU_decomp(sym_denom_LU,perm,&sign);
> >
> >   gsl_matrix_memcpy(sym_denom_Lower,sym_denom_LU);
> >   gsl_matrix_memcpy(sym_denom_Upper,sym_denom_LU);
> >
> >   /* check on LU decomp */
> >   for(ina=0;ina<num_ind;ina++){
> >     for(inb=ina;inb<num_ind;inb++){
> >       if(inb > ina ) gsl_matrix_set(sym_denom_Lower,ina,inb,0.0);
> >     }
> >     for(inb=ina+1;inb<num_ind;inb++){
> >       gsl_matrix_set(sym_denom_Upper,inb,ina,0.0);
> >     }
> >   }
> >
> >   /* Multiply lower/upper matrices to get original sym_denom */
> >
> > gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,1.0,sym_denom_Lo
> >wer,sym_denom_Upper);
> >
> > ------------------------------------------------------------------------
> >
> >
> > sym_denom
> >  1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20
> > -3.066652e-19   4.391243e-19 2    -1.652853e-19  -3.923216e-20
> > -3.608558e-20  -1.378546e-20 -7.772215e-20   9.616245e-20 3
> > -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20
> > 7.866459e-20 4    -6.114344e-20  -1.378546e-20  -2.635082e-20
> > -9.419878e-21 -1.627012e-20   3.311913e-20 5    -3.066652e-19
> > -7.772215e-20  -3.128613e-20  -1.627012e-20 -2.676652e-19   2.252427e-19 6
> >    4.391243e-19   9.616245e-20   7.866459e-20   3.311913e-20 2.252427e-19
> > -3.018037e-19 7    -1.748050e-20  -3.759400e-21  -2.671663e-20
> > -8.956211e-21 2.023082e-20   4.257951e-21 8     2.556328e-19   6.230603e-20
> >   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19 9
> > 4.908470e-20   1.223947e-20   3.696976e-20   1.238438e-20 3.418982e-21
> > -2.760516e-20
> >
> >
> > LU check - should be the same as sym_denom
> >  1    -7.920023e-19  -1.652853e-19  -1.560040e-19  -6.114344e-20
> > -3.066652e-19   4.391243e-19 2     3.491029e-19   8.709227e-20
> > 5.217473e-20   2.250474e-20 2.709399e-19  -2.325323e-19 3     2.556328e-19
> >  6.230603e-20   1.710786e-19   5.691356e-20 3.972437e-20  -1.485747e-19 4
> >  -1.560040e-19  -3.608558e-20  -8.714679e-20  -2.635082e-20 -3.128613e-20
> > 7.866459e-20 5    -1.652853e-19  -3.923216e-20  -3.608558e-20
> > -1.378546e-20 -7.772215e-20   9.616245e-20 6     2.847849e-19
> > 6.186536e-20   4.510521e-20   2.035696e-20 1.565680e-19  -2.154436e-19 7
> > -1.748050e-20  -3.759400e-21  -2.671663e-20  -8.956211e-21 2.023082e-20
> > 4.257951e-21 8    -1.493627e-19  -4.009839e-20  -1.313729e-19
> > -4.382593e-20 -5.583881e-21   8.188773e-20 9     1.981133e-20
> > 5.659666e-21   2.816777e-20   9.600430e-21 -8.449448e-21  -8.121479e-21
> >
> > *********************************************************
> >
> > Peter Roche
> > Department of Applied Mathematics and Theoretical Physics,
> > University of Cambridge, Silver Street, Cambridge, CB3 9EW
> >
> > email: p.j.p.roche@damtp.cam.ac.uk
> >
> > *********************************************************
>