{"id":5303,"date":"2013-12-17T17:15:04","date_gmt":"2013-12-17T16:15:04","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=5303"},"modified":"2013-12-17T17:15:04","modified_gmt":"2013-12-17T16:15:04","slug":"computing-eigenvalues-without-balancing-in-python-using-the-nag-library","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=5303","title":{"rendered":"Computing Eigenvalues without balancing in Python using the NAG library"},"content":{"rendered":"

In a recent Stack Overflow query<\/a>, someone asked if you could switch off the balancing step when calculating eigenvalues in Python. In the document A case where balancing is harmful<\/a>, David S. Watkins describes the balancing step as ‘the input matrix A is replaced by a rescaled matrix A* = D-1<\/sup>AD, where D is a diagonal matrix chosen so that, for each i, the ith row and the ith column of A* have roughly the same norm.’ \u00a0<\/em><\/p>\n

Such balancing is usually very useful and so is performed by default by software such as \u00a0MATLAB or Numpy. \u00a0There are times, however, when one would like to switch it off.<\/p>\n

In MATLAB, this is easy and the following is taken from the online MATLAB documentation<\/p>\n

A = [ 3.0     -2.0      -0.9     2*eps;\r\n     -2.0      4.0       1.0    -eps;\r\n     -eps\/4    eps\/2    -1.0     0;\r\n     -0.5     -0.5       0.1     1.0];\r\n\r\n[VN,DN] = eig(A,'nobalance')\r\nVN =\r\n\r\n    0.6153   -0.4176   -0.0000   -0.1528\r\n   -0.7881   -0.3261         0    0.1345\r\n   -0.0000   -0.0000   -0.0000   -0.9781\r\n    0.0189    0.8481   -1.0000    0.0443\r\n\r\nDN =\r\n    5.5616         0         0         0\r\n         0    1.4384         0         0\r\n         0         0    1.0000         0\r\n         0         0         0   -1.0000<\/pre>\n
At the time of writing, it is not possible to directly do this in Numpy (as far as I know at least). Numpy’s eig command currently uses the LAPACK routine DGEEV<\/a> to do the heavy lifting for double precision matrices. \u00a0We can see this by looking at the source code of numpy.linalg.eig<\/a>\u00a0where the relevant subsection is<\/p>\n
lapack_routine = lapack_lite.dgeev\r\n        wr = zeros((n,), t)\r\n        wi = zeros((n,), t)\r\n        vr = zeros((n, n), t)\r\n        lwork = 1\r\n        work = zeros((lwork,), t)\r\n        results = lapack_routine(_N, _V, n, a, n, wr, wi,\r\n                                  dummy, 1, vr, n, work, -1, 0)\r\n        lwork = int(work[0])\r\n        work = zeros((lwork,), t)\r\n        results = lapack_routine(_N, _V, n, a, n, wr, wi,\r\n                                  dummy, 1, vr, n, work, lwork, 0)<\/pre>\n
My plan was to figure out how to tell DGEEV not to perform the balancing step and I’d be done. Sadly, however, it turns out that this is not possible. Taking a look at the reference implementation of DGEEV<\/a>, we can see that the balancing step is always performed\u00a0<\/strong>and is not user controllable–here’s the relevant bit of Fortran<\/p>\n
*     Balance the matrix\r\n*     (Workspace: need N)\r\n*\r\n      IBAL = 1\r\n      CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )<\/pre>\n
So, using DGEEV is a dead-end unless we are willing to modifiy and recompile the lapack source — something that’s rarely a good idea in my experience. There is another LAPACK routine that is of use, however, in the form of DGEEVX<\/a> that allows us to control balancing. \u00a0Unfortunately, this routine is not part of the numpy.linalg.lapack_lite<\/strong>\u00a0interface provided by Numpy and I’ve yet to figure out how to add extra routines to it.<\/p>\n
I’ve also discovered that this functionality is an open feature request in Numpy<\/a>.<\/p>\n
Enter the NAG Library<\/strong><\/p>\n
My University has a site license for the commercial Numerical Algorithms Group (NAG) library<\/a>. \u00a0Among other things, NAG offers an interface to all of LAPACK<\/a> along with an interface to Python<\/a>. \u00a0So, I go through the installation and do<\/p>\n
import numpy as np\r\nfrom ctypes import *\r\nfrom nag4py.util import Nag_RowMajor,Nag_NoBalancing,Nag_NotLeftVecs,Nag_RightVecs,Nag_RCondEigVecs,Integer,NagError,INIT_FAIL\r\nfrom nag4py.f08 import nag_dgeevx\r\n\r\neps = np.spacing(1)\r\nnp.set_printoptions(precision=4,suppress=True) \r\n\r\ndef unbalanced_eig(A):\r\n    \"\"\"\r\n    Compute the eigenvalues and right eigenvectors of a square array using DGEEVX via the NAG library.\r\n    Requires the NAG C library and NAG's Python wrappers http:\/\/www.nag.co.uk\/python.asp\r\n    The balancing step that's performed in DGEEV is not performed here.\r\n    As such, this function is the same as the MATLAB command eig(A,'nobalance')\r\n\r\n    Parameters\r\n    ----------\r\n    A : (M, M) Numpy array\r\n        A square array of real elements.\r\n\r\n        On exit: \r\n        A is overwritten and contains the real Schur form of the balanced version of the input matrix .\r\n\r\n    Returns\r\n    -------\r\n    w : (M,) ndarray\r\n        The eigenvalues\r\n    v : (M, M) ndarray\r\n        The eigenvectors\r\n\r\n    Author: Mike Croucher (www.walkingrandomly.com)\r\n    Testing has been mimimal\r\n    \"\"\"\r\n\r\n    order = Nag_RowMajor\r\n    balanc = Nag_NoBalancing\r\n    jobvl = Nag_NotLeftVecs\r\n    jobvr = Nag_RightVecs\r\n    sense = Nag_RCondEigVecs\r\n\r\n    n = A.shape[0]\r\n    pda = n\r\n    pdvl = 1\r\n\r\n    wr = np.zeros(n)\r\n    wi = np.zeros(n)\r\n\r\n    vl=np.zeros(1);\r\n    pdvr = n\r\n    vr = np.zeros(pdvr*n)\r\n\r\n    ilo=c_long(0)\r\n    ihi=c_long(0)\r\n\r\n    scale = np.zeros(n)\r\n    abnrm = c_double(0)\r\n    rconde = np.zeros(n)\r\n    rcondv = np.zeros(n)\r\n\r\n    fail = NagError()\r\n    INIT_FAIL(fail)\r\n\r\n    nag_dgeevx(order,balanc,jobvl,jobvr,sense, \r\n               n, A.ctypes.data_as(POINTER(c_double)), pda, wr.ctypes.data_as(POINTER(c_double)),\r\n               wi.ctypes.data_as(POINTER(c_double)),vl.ctypes.data_as(POINTER(c_double)),pdvl, \r\n               vr.ctypes.data_as(POINTER(c_double)),pdvr,ilo,ihi, scale.ctypes.data_as(POINTER(c_double)),\r\n               abnrm, rconde.ctypes.data_as(POINTER(c_double)),rcondv.ctypes.data_as(POINTER(c_double)),fail)\r\n\r\n    if all(wi == 0.0):\r\n            w = wr\r\n            v = vr.reshape(n,n)\r\n    else:\r\n            w = wr+1j*wi\r\n            v = array(vr, w.dtype).reshape(n,n)\r\n\r\n    return(w,v)<\/pre>\n
Define a test matrix:<\/p>\n
A = np.array([[3.0,-2.0,-0.9,2*eps],\r\n          [-2.0,4.0,1.0,-eps],\r\n          [-eps\/4,eps\/2,-1.0,0],\r\n          [-0.5,-0.5,0.1,1.0]])<\/pre>\n
Do the calculation<\/p>\n
(w,v) = unbalanced_eig(A)<\/pre>\n
which gives<\/p>\n
(array([ 5.5616,  1.4384,  1.    , -1.    ]),\r\n array([[ 0.6153, -0.4176, -0.    , -0.1528],\r\n       [-0.7881, -0.3261,  0.    ,  0.1345],\r\n       [-0.    , -0.    , -0.    , -0.9781],\r\n       [ 0.0189,  0.8481, -1.    ,  0.0443]]))<\/pre>\n
This is exactly what you get by running the MATLAB command eig(A,’nobalance’).<\/strong><\/p>\n
Note that unbalanced_eig(A) changes the input matri<\/strong>x A to<\/p>\n
array([[ 5.5616, -0.0662,  0.0571,  1.3399],\r\n       [ 0.    ,  1.4384,  0.7017, -0.1561],\r\n       [ 0.    ,  0.    ,  1.    , -0.0132],\r\n       [ 0.    ,  0.    ,  0.    , -1.    ]])<\/pre>\n
According to the NAG documentation<\/a>, this is the real Schur form of the balanced version of the input matrix. \u00a0I can’t see how to ask NAG to not do this. I guess that if it’s not what you want unbalanced_eig() to do, \u00a0you’ll need to pass a copy of the input matrix to NAG.<\/p>\n
The IPython notebook<\/strong><\/p>\n
The code for this article is available as an IPython Notebook<\/a><\/p>\n
The future<\/strong><\/p>\n
This blog post was written using Numpy version 1.7.1. There is an enhancement request for the functionality<\/a> discussed in this article open in Numpy’s git repo and so I expect this article to become redundant pretty soon.<\/p>\n","protected":false},"excerpt":{"rendered":"
In a recent Stack Overflow query, someone asked if you could switch off the balancing step when calculating eigenvalues in Python. In the document A case where balancing is harmful, David S. Watkins describes the balancing step as ‘the input matrix A is replaced by a rescaled matrix A* = D-1AD, where D is a […]<\/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,11,28,31],"tags":[],"class_list":["post-5303","post","type-post","status-publish","format-standard","hentry","category-linear-algebra","category-matlab","category-nag-library","category-python"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1nx","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5303","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=5303"}],"version-history":[{"count":5,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5303\/revisions"}],"predecessor-version":[{"id":5310,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5303\/revisions\/5310"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=5303"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=5303"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=5303"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}