From d1638062ec3a28dc64b91f308858211db243db7b Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Thu, 13 Oct 2016 08:32:19 -0400 Subject: [PATCH] Add eigenvalues_re() to compute real parts of eigenvalues of any matrix. --- src/dunshire/matrices.py | 72 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 71 insertions(+), 1 deletion(-) diff --git a/src/dunshire/matrices.py b/src/dunshire/matrices.py index bf8ae03..384b11f 100644 --- a/src/dunshire/matrices.py +++ b/src/dunshire/matrices.py @@ -3,9 +3,10 @@ Utility functions for working with CVXOPT matrices (instances of the class:`cvxopt.base.matrix` class). """ +from copy import copy from math import sqrt from cvxopt import matrix -from cvxopt.lapack import syevr +from cvxopt.lapack import gees, syevr import options @@ -95,6 +96,10 @@ def eigenvalues(symmat): """ Return the eigenvalues of the given symmetric real matrix. + On the surface, this appears redundant to the :func:`eigenvalues_re` + function. However, if we know in advance that our input is + symmetric, a better algorithm can be used. + Parameters ---------- @@ -140,6 +145,71 @@ def eigenvalues(symmat): return list(eigs) +def eigenvalues_re(anymat): + """ + Return the real parts of the eigenvalues of the given square matrix. + + Parameters + ---------- + + anymat : matrix + The square matrix whose eigenvalues you want. + + Returns + ------- + + list of float + A list of the real parts (in no particular order) of the + eigenvalues of ``anymat``. + + Raises + ------ + + TypeError + If the input matrix is not square. + + Examples + -------- + + This is symmetric and has two real eigenvalues: + + >>> A = matrix([[2,1],[1,2]], tc='d') + >>> sorted(eigenvalues_re(A)) + [1.0, 3.0] + + But this rotation matrix has eigenvalues `i` and `-i`, both of whose + real parts are zero: + + >>> A = matrix([[0,-1],[1,0]]) + >>> eigenvalues_re(A) + [0.0, 0.0] + + If the input matrix is not square, it doesn't have eigenvalues:: + + >>> A = matrix([[1,2],[3,4],[5,6]]) + >>> eigenvalues_re(A) + Traceback (most recent call last): + ... + TypeError: input matrix must be square + + """ + if not anymat.size[0] == anymat.size[1]: + raise TypeError('input matrix must be square') + + domain_dim = anymat.size[0] + eigs = matrix(0, (domain_dim, 1), tc='z') + + # Create a copy of ``anymat`` here for two reasons: + # + # 1. ``gees`` clobbers its input. + # 2. We need to ensure that the type code of ``dummy`` is 'd' or 'z'. + # + dummy = matrix(anymat, anymat.size, tc='d') + + gees(dummy, eigs) + return [eig.real for eig in eigs] + + def identity(domain_dim): """ Create an identity matrix of the given dimensions. -- 2.44.2