{"id":5092,"date":"2013-09-16T15:16:44","date_gmt":"2013-09-16T14:16:44","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=5092"},"modified":"2013-09-16T15:16:44","modified_gmt":"2013-09-16T14:16:44","slug":"numpy-matlab-and-singular-matrices","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=5092","title":{"rendered":"Numpy, MATLAB and singular matrices"},"content":{"rendered":"

Last week I gave a live demo of the IPython notebook<\/a> to a group of numerical analysts and one of the computations we attempted to do was to solve the following linear system using Numpy’s solve command.<\/p>\n

\"<\/p>\n

Now, the matrix shown above is singular<\/a> and so we expect that we might have problems. Before looking at how Numpy deals with this computation, lets take a look at what happens if you ask MATLAB<\/a> to do it<\/p>\n

>> A=[1 2 3;4 5 6;7 8 9];\r\n>> b=[15;15;15];\r\n>> x=A\\b\r\nWarning: Matrix is close to singular\r\nor badly scaled. Results may be\r\ninaccurate. RCOND =  1.541976e-18. \r\nx =\r\n  -39.0000\r\n   63.0000\r\n  -24.0000<\/pre>\n
MATLAB gives us a warning that the input matrix is close<\/strong> to being singular (note that it didn’t actually recognize that it is<\/strong> singular) along with an estimate of the reciprocal of the condition number<\/a>. It tells us that the results may be inaccurate and we’d do well to check. So, lets check:<\/p>\n
>> A*x\r\n\r\nans =\r\n   15.0000\r\n   15.0000\r\n   15.0000\r\n\r\n>> norm(A*x-b)\r\n\r\nans =\r\n2.8422e-14<\/pre>\n
We seem to have dodged the bullet since, despite the singular nature of our matrix, MATLAB has able to find a valid solution. MATLAB was right to have warned us though…in other cases we might not have been so lucky.<\/p>\n
Let’s see how Numpy deals with this using the IPython notebook:<\/p>\n
In [1]:\r\nimport numpy\r\nfrom numpy import array\r\nfrom numpy.linalg import solve\r\nA=array([[1,2,3],[4,5,6],[7,8,9]])\r\nb=array([15,15,15])\r\nsolve(A,b)\r\n\r\nOut[1]:\r\n\r\narray([-39.,  63., -24.])<\/pre>\n
It gave the same result as MATLAB [See note 1], presumably because it’s using the exact same LAPACK routine<\/a>, but there was no warning of the singular nature of the matrix. \u00a0During my demo, it was generally felt by everyone in the room that a warning should have been given, particularly when working in an interactive setting.<\/p>\n
If you look at the documentation for Numpy’s solve<\/a> command you’ll see that it is supposed to throw an exception when the matrix is singular but it clearly didn’t do so here. The exception is sometimes thrown though:<\/p>\n
In [4]:\r\n\r\nC=array([[1,1,1],[1,1,1],[1,1,1]])\r\nx=solve(C,b)\r\n\r\n---------------------------------------------------------------------------\r\nLinAlgError                               Traceback (most recent call last)\r\n in ()\r\n      1 C=array([[1,1,1],[1,1,1],[1,1,1]])\r\n----> 2 x=solve(C,b)\r\n\r\nC:\\Python32\\lib\\site-packages\\numpy\\linalg\\linalg.py in solve(a, b)\r\n    326     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)\r\n    327     if results['info'] > 0:\r\n--> 328         raise LinAlgError('Singular matrix')\r\n    329     if one_eq:\r\n    330         return wrap(b.ravel().astype(result_t))\r\n\r\nLinAlgError: Singular matrix<\/pre>\n
It seems that Numpy is somehow checking for exact singularity but this will rarely be detected due to rounding errors. Those I’ve spoken to consider that MATLAB’s approach of estimating the condition number and warning when that is high would be better behavior since it alerts the user to the fact that the matrix is badly conditioned.<\/p>\n
Thanks to Nick Higham<\/a> and David Silvester<\/a> for useful discussions regarding this post.<\/p>\n
Notes<\/strong>
\n[1] – The results really are identical which you can see by rerunning the calculation after evaluating format long<\/strong> in MATLAB and numpy.set_printoptions(precision=15)<\/strong> in Python<\/p>\n","protected":false},"excerpt":{"rendered":"
Last week I gave a live demo of the IPython notebook to a group of numerical analysts and one of the computations we attempted to do was to solve the following linear system using Numpy’s solve command. Now, the matrix shown above is singular and so we expect that we might have problems. Before looking […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[73,4,11,31],"tags":[],"class_list":["post-5092","post","type-post","status-publish","format-standard","hentry","category-linear-algebra","category-math-software","category-matlab","category-python"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1k8","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5092","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=5092"}],"version-history":[{"count":23,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5092\/revisions"}],"predecessor-version":[{"id":5116,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5092\/revisions\/5116"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=5092"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=5092"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=5092"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}