This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

cholesky decomposition bug?


hi!

i'm using the cholesky decomposition for a newton-raphson iteration and
seem to have found a bug: a matrix that gsl_cholesky_decomp chokes on
yet is symmetric and positive definate, since maple can decompose it
fine (with linalg[cholesky]).

the matrix in question is:

H := [[ 1.105615e+04 , 2.984318e+02 , 2.136758e+02 , 1.687062e+02 ,
2.186212e+02 , 1.984495e+02 , 1.949485e+02 , 2.085723e+02 , 1.620816e+02
, 1.416917e+02 , 1.682196e+02 , 2.297207e+02 , 1.707851e+02 ,
1.971006e+02 , 1.953821e+02 , 2.589295e+02 , 1.942593e+02 , 1.321812e+02
, 1.810242e+02 , 2.226683e+02 ],
[ 2.984318e+02 , 9.637381e+03 , 2.134579e+02 , 2.119008e+02 ,
1.634247e+02 , 1.770826e+02 , 1.727033e+02 , 1.912182e+02 , 1.510900e+02
, 1.682738e+02 , 1.924433e+02 , 3.008754e+02 , 1.851525e+02 ,
1.588338e+02 , 2.090053e+02 , 2.382960e+02 , 2.375855e+02 , 9.384388e+01
, 1.980848e+02 , 1.904123e+02 ],
[ 2.136758e+02 , 2.134579e+02 , 1.767244e+04 , 1.871965e+02 ,
2.105550e+02 , 2.048553e+02 , 2.041612e+02 , 2.065549e+02 , 1.940995e+02
, 1.805209e+02 , 1.899688e+02 , 1.996041e+02 , 1.923128e+02 ,
2.060599e+02 , 1.990169e+02 , 2.141627e+02 , 1.946252e+02 , 1.744385e+02
, 1.947292e+02 , 2.109173e+02 ],
[ 1.687062e+02 , 2.119008e+02 , 1.871965e+02 , 1.699480e+04 ,
1.396089e+02 , 1.718975e+02 , 1.715848e+02 , 1.737547e+02 , 1.877004e+02
, 2.379468e+02 , 2.256384e+02 , 2.420733e+02 , 2.135947e+02 ,
1.561675e+02 , 2.049513e+02 , 1.603504e+02 , 2.356509e+02 , 1.318139e+02
, 2.131295e+02 , 1.592888e+02 ],
[ 2.186212e+02 , 1.634247e+02 , 2.105550e+02 , 1.396089e+02 ,
5.998519e+03 , 2.336708e+02 , 2.337985e+02 , 2.313066e+02 , 2.049457e+02
, 1.431296e+02 , 1.606709e+02 , 1.421151e+02 , 1.740419e+02 ,
2.583737e+02 , 1.861049e+02 , 2.467230e+02 , 1.522068e+02 , 2.375777e+02
, 1.755065e+02 , 2.557145e+02 ],
[ 1.984495e+02 , 1.770826e+02 , 2.048553e+02 , 1.718975e+02 ,
2.336708e+02 , 1.481967e+04 , 2.169281e+02 , 2.145591e+02 , 2.071013e+02
, 1.756634e+02 , 1.837965e+02 , 1.693388e+02 , 1.907383e+02 ,
2.262182e+02 , 1.947694e+02 , 2.144640e+02 , 1.773474e+02 , 2.085101e+02
, 1.906065e+02 , 2.230084e+02 ],
[ 1.949485e+02 , 1.727033e+02 , 2.041612e+02 , 1.715848e+02 ,
2.337985e+02 , 2.169281e+02 , 1.457797e+04 , 2.144316e+02 , 2.087695e+02
, 1.767615e+02 , 1.840878e+02 , 1.666899e+02 , 1.912437e+02 ,
2.270293e+02 , 1.944093e+02 , 2.124257e+02 , 1.762793e+02 , 2.118772e+02
, 1.906624e+02 , 2.226032e+02 ],
[ 2.085723e+02 , 1.912182e+02 , 2.065549e+02 , 1.737547e+02 ,
2.313066e+02 , 2.145591e+02 , 2.144316e+02 , 1.557088e+04 , 2.018461e+02
, 1.732062e+02 , 1.834598e+02 , 1.783635e+02 , 1.895151e+02 ,
2.224422e+02 , 1.959291e+02 , 2.194959e+02 , 1.812970e+02 , 1.974257e+02
, 1.907073e+02 , 2.229237e+02 ],
[ 1.620816e+02 , 1.510900e+02 , 1.940995e+02 , 1.877004e+02 ,
2.049457e+02 , 2.071013e+02 , 2.087695e+02 , 2.018461e+02 , 1.500642e+04
, 2.020212e+02 , 1.974530e+02 , 1.643343e+02 , 2.016612e+02 ,
2.121628e+02 , 1.945283e+02 , 1.817208e+02 , 1.825294e+02 , 2.141151e+02
, 1.972919e+02 , 2.005543e+02 ],
[ 1.416917e+02 , 1.682738e+02 , 1.805209e+02 , 2.379468e+02 ,
1.431296e+02 , 1.756634e+02 , 1.767615e+02 , 1.732062e+02 , 2.020212e+02
, 1.536298e+04 , 2.244253e+02 , 2.085502e+02 , 2.154934e+02 ,
1.640777e+02 , 1.993281e+02 , 1.462429e+02 , 2.194752e+02 , 1.568557e+02
, 2.106538e+02 , 1.577937e+02 ],
[ 1.682196e+02 , 1.924433e+02 , 1.899688e+02 , 2.256384e+02 ,
1.606709e+02 , 1.837965e+02 , 1.840878e+02 , 1.834598e+02 , 1.974530e+02
, 2.244253e+02 , 1.805351e+04 , 2.145629e+02 , 2.091938e+02 ,
1.745101e+02 , 2.014848e+02 , 1.684338e+02 , 2.167138e+02 , 1.570624e+02
, 2.075122e+02 , 1.732744e+02 ],
[ 2.297207e+02 , 3.008754e+02 , 1.996041e+02 , 2.420733e+02 ,
1.421151e+02 , 1.693388e+02 , 1.666899e+02 , 1.783635e+02 , 1.643343e+02
, 2.085502e+02 , 2.145629e+02 , 1.493801e+04 , 2.018336e+02 ,
1.494437e+02 , 2.095509e+02 , 1.940979e+02 , 2.479355e+02 , 1.012540e+02
, 2.090713e+02 , 1.684941e+02 ],
[ 1.707851e+02 , 1.851525e+02 , 1.923128e+02 , 2.135947e+02 ,
1.740419e+02 , 1.907383e+02 , 1.912437e+02 , 1.895151e+02 , 2.016612e+02
, 2.154934e+02 , 2.091938e+02 , 2.018336e+02 , 1.810595e+04 ,
1.852501e+02 , 1.998485e+02 , 1.751521e+02 , 2.070869e+02 , 1.703398e+02
, 2.042293e+02 , 1.823287e+02 ],
[ 1.971006e+02 , 1.588338e+02 , 2.060599e+02 , 1.561675e+02 ,
2.583737e+02 , 2.262182e+02 , 2.270293e+02 , 2.224422e+02 , 2.121628e+02
, 1.640777e+02 , 1.745101e+02 , 1.494437e+02 , 1.852501e+02 ,
1.067627e+04 , 1.904392e+02 , 2.228286e+02 , 1.629951e+02 , 2.335320e+02
, 1.845050e+02 , 2.368922e+02 ],
[ 1.953821e+02 , 2.090053e+02 , 1.990169e+02 , 2.049513e+02 ,
1.861049e+02 , 1.947694e+02 , 1.944093e+02 , 1.959291e+02 , 1.945283e+02
, 1.993281e+02 , 2.014848e+02 , 2.095509e+02 , 1.998485e+02 ,
1.904392e+02 , 1.879300e+04 , 1.939897e+02 , 2.059799e+02 , 1.642005e+02
, 2.008211e+02 , 1.932790e+02 ],
[ 2.589295e+02 , 2.382960e+02 , 2.141627e+02 , 1.603504e+02 ,
2.467230e+02 , 2.144640e+02 , 2.124257e+02 , 2.194959e+02 , 1.817208e+02
, 1.462429e+02 , 1.684338e+02 , 1.940979e+02 , 1.751521e+02 ,
2.228286e+02 , 1.939897e+02 , 1.238953e+04 , 1.795780e+02 , 1.714368e+02
, 1.816107e+02 , 2.372049e+02 ],
[ 1.942593e+02 , 2.375855e+02 , 1.946252e+02 , 2.356509e+02 ,
1.522068e+02 , 1.773474e+02 , 1.762793e+02 , 1.812970e+02 , 1.825294e+02
, 2.194752e+02 , 2.167138e+02 , 2.479355e+02 , 2.070869e+02 ,
1.629951e+02 , 2.059799e+02 , 1.795780e+02 , 1.823611e+04 , 1.299884e+02
, 2.094123e+02 , 1.711028e+02 ],
[ 1.321812e+02 , 9.384388e+01 , 1.744385e+02 , 1.318139e+02 ,
2.375777e+02 , 2.085101e+02 , 2.118772e+02 , 1.974257e+02 , 2.141151e+02
, 1.568557e+02 , 1.570624e+02 , 1.012540e+02 , 1.703398e+02 ,
2.335320e+02 , 1.642005e+02 , 1.714368e+02 , 1.299884e+02 , 1.892783e+02
, 1.638115e+02 , 2.070998e+02 ],
[ 1.810242e+02 , 1.980848e+02 , 1.947292e+02 , 2.131295e+02 ,
1.755065e+02 , 1.906065e+02 , 1.906624e+02 , 1.907073e+02 , 1.972919e+02
, 2.106538e+02 , 2.075122e+02 , 2.090713e+02 , 2.042293e+02 ,
1.845050e+02 , 2.008211e+02 , 1.816107e+02 , 2.094123e+02 , 1.638115e+02
, 1.862820e+04 , 1.846298e+02 ],
[ 2.226683e+02 , 1.904123e+02 , 2.109173e+02 , 1.592888e+02 ,
2.557145e+02 , 2.230084e+02 , 2.226032e+02 , 2.229237e+02 , 2.005543e+02
, 1.577937e+02 , 1.732744e+02 , 1.684941e+02 , 1.823287e+02 ,
2.368922e+02 , 1.932790e+02 , 2.372049e+02 , 1.711028e+02 , 2.070998e+02
, 1.846298e+02 , 1.254316e+04 ]];


the decomposition maple makes of it is:

H_decomp := linalg[cholesky](H);
H_decomp := [
[105.1482287 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0
, 0 , 0 , 0 , 0],
[2.838200925 , 98.12912724 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0],
[2.032138845 , 2.116499834 , 132.9053454 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0
, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.604460694 , 2.113001756 , 1.350313222 , 130.3301101 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[2.079171496 , 1.605268465 , 1.526893322 , 1.003752945 , 77.38397921 , 0
, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.887330889 , 1.750000030 , 1.484636349 , 1.251950716 , 2.887082926 ,
121.6591086 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.854035036 , 1.706335118 , 1.480618291 , 1.250710876 , 2.890628896 ,
1.630238334 , 120.6517760 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,
0 , 0],
[1.983602602 , 1.891266563 , 1.493702741 , 1.262631264 , 2.850696839 ,
1.606761019 , 1.598619621 , 124.6847197 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0
, 0 , 0 , 0 , 0],
[1.541458206 , 1.495122152 , 1.413055408 , 1.382335423 , 2.530182159 ,
1.565376089 , 1.572074261 , 1.442547287 , 122.4111775 , 0 , 0 , 0 , 0 ,
0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.347542434 , 1.675844964 , 1.310974969 , 1.768382079 , 1.729827167 ,
1.323641212 , 1.326898236 , 1.235063345 , 1.493532365 , 123.8678684 , 0
, 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.599832941 , 1.914850951 , 1.374398712 , 1.666303963 , 1.944841904 ,
1.378314916 , 1.374753521 , 1.303700391 , 1.443971457 , 1.643141506 ,
134.2704793 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[2.184732000 , 3.002927875 , 1.420625516 , 1.767085918 , 1.664547393 ,
1.239802449 , 1.213152876 , 1.245708351 , 1.161415429 , 1.503078647 ,
1.400425240 , 122.0908593 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ],
[1.624231831 , 1.839847237 , 1.392856686 , 1.574619177 , 2.119354991 ,
1.432651624 , 1.430558357 , 1.348316500 , 1.474235764 , 1.568470812 ,
1.370422137 , 1.419360714 , 134.4491452 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
[1.874502333 , 1.564403862 , 1.496851687 , 1.134297577 , 3.211787701 ,
1.701702625 , 1.700690243 , 1.583901843 , 1.531743839 , 1.135508126 ,
1.098497559 , 0.9829025604 , 1.150644419 , 103.1554013 , 0 , 0 , 0 , 0 ,
0 , 0],
[1.858158739 , 2.076156976 , 1.435959382 , 1.501142285 , 2.264155308 ,
1.455551739 , 1.446313209 , 1.388872279 , 1.406499897 , 1.430793981 ,
1.306847290 , 1.473425421 , 1.261602271 , 1.525203360 , 136.9592402 , 0
, 0 , 0 , 0 , 0],
[2.462518895 , 2.357168387 , 1.536202871 , 1.145892758 , 3.028059609 ,
1.588321716 , 1.564733985 , 1.545708165 , 1.272828449 , 0.9825575887 ,
1.044628849 , 1.328362376 , 1.064501759 , 1.809354494 , 1.141910256 ,
111.1150532 , 0 , 0 , 0 , 0],
[1.847480480 , 2.367716759 , 1.398435785 , 1.732488173 , 1.818083010 ,
1.316982967 , 1.302708527 , 1.279199191 , 1.317001453 , 1.598130994 ,
1.423083398 , 1.785419956 , 1.314944896 , 1.267837628 , 1.249322484 ,
1.273238816 , 134.8990011 , 0 , 0 , 0],
[1.257093930 , 0.9199714433 , 1.278630221 , 0.9677463603 , 2.979472977 ,
1.584885363 , 1.605255264 , 1.415202404 , 1.577247465 , 1.103975140 ,
0.9964763051 , 0.6281958678 , 1.072501292 , 1.966654671 , 0.9707070461 ,
1.224718391 , 0.7232103829 , 12.45930095 , 0 , 0],
[1.721609600 , 1.968819365 , 1.407495265 , 1.567608335 , 2.132791696 ,
1.427777070 , 1.422056942 , 1.354132976 , 1.435476074 , 1.526905906 ,
1.355647244 , 1.474599758 , 1.297844876 , 1.474526076 , 1.214092737 ,
1.292258828 , 1.274377980 , 10.39392870 , 135.9463624 , 0],
[2.117660970 , 1.879176528 , 1.524668865 , 1.149861781 , 3.163610293 ,
1.667662718 , 1.656928525 , 1.580730966 , 1.431538848 , 1.080101266 ,
1.084712121 , 1.129558133 , 1.123157322 , 1.949123684 , 1.141320935 ,
1.755782669 , 0.9763374964 , 13.51905074 , 0.02297779832 ,
110.9713066]];

the c-code that does the decomposition is, abreviated:

	A = gsl_matrix_view_array(H,n,n);
	b = gsl_vector_view_array(g,n);
	x = gsl_vector_view_array(x_delta,n);
	gsl_set_error_handler_off();
	if ((res = gsl_linalg_cholesky_decomp(&A.matrix)) != 0) {
		...
		}
	gsl_linalg_cholesky_solve(&A.matrix,&b.vector,&x.vector);

this stuff works sometimes, so i don't think it's a coding problem on my
behalf.

cheers and many thanks!
pedro gonnet


	


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]