From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 22089 invoked by alias); 1 Nov 2002 14:48:02 -0000 Mailing-List: contact gsl-discuss-help@sources.redhat.com; run by ezmlm Precedence: bulk List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sources.redhat.com Received: (qmail 22003 invoked from network); 1 Nov 2002 14:48:01 -0000 Received: from unknown (HELO xsmtp.ethz.ch) (129.132.97.6) by sources.redhat.com with SMTP; 1 Nov 2002 14:48:01 -0000 Received: from xfe2.d.ethz.ch ([192.168.36.11]) by xsmtp.ethz.ch with Microsoft SMTPSVC(5.0.2195.5329); Fri, 1 Nov 2002 15:47:59 +0100 Received: from localhost.localdomain ([129.132.17.66]) by xfe2.d.ethz.ch with Microsoft SMTPSVC(5.0.2195.5329); Fri, 1 Nov 2002 15:47:59 +0100 Subject: cholesky decomposition bug? From: Pedro Gonnet To: gsl-discuss@sources.redhat.com Content-Type: text/plain Content-Transfer-Encoding: 7bit Date: Fri, 01 Nov 2002 23:41:00 -0000 Message-Id: <1036162079.11311.48.camel@selonia> Mime-Version: 1.0 X-OriginalArrivalTime: 01 Nov 2002 14:47:59.0552 (UTC) FILETIME=[AC0EF400:01C281B5] X-SW-Source: 2002-q4/txt/msg00096.txt.bz2 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