X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=dunshire%2Fmatrices.py;h=bcf83778d62752436f7894b5a6fbad4cc3e1012e;hb=040374ca134b2f3d962b91a9dac97a7600032685;hp=6ed4f851487235dba117ef1521151bd432db9076;hpb=bdb596b84a06d0c97e39d42586a51fc36ba44186;p=dunshire.git diff --git a/dunshire/matrices.py b/dunshire/matrices.py index 6ed4f85..bcf8377 100644 --- a/dunshire/matrices.py +++ b/dunshire/matrices.py @@ -5,7 +5,7 @@ class:`cvxopt.base.matrix` class). from math import sqrt from cvxopt import matrix -from cvxopt.lapack import gees, syevr +from cvxopt.lapack import gees, gesdd, syevr from . import options @@ -140,7 +140,10 @@ def eigenvalues(symmat): domain_dim = symmat.size[0] eigs = matrix(0, (domain_dim, 1), tc='d') - syevr(symmat, eigs) + + # Create a copy of ``symmat`` here because ``syevr`` clobbers it. + dummy = matrix(symmat, symmat.size) + syevr(dummy, eigs) return list(eigs) @@ -209,7 +212,7 @@ def eigenvalues_re(anymat): return [eig.real for eig in eigs] -def identity(domain_dim): +def identity(domain_dim, typecode='i'): """ Create an identity matrix of the given dimensions. @@ -219,6 +222,10 @@ def identity(domain_dim): domain_dim : int The dimension of the vector space on which the identity will act. + typecode : {'i', 'd', 'z'}, optional + The type code for the returned matrix, defaults to 'i' for integers. + Can also be 'd' for real double, or 'z' for complex double. + Returns ------- @@ -247,7 +254,7 @@ def identity(domain_dim): entries = [int(i == j) for i in range(domain_dim) for j in range(domain_dim)] - return matrix(entries, (domain_dim, domain_dim)) + return matrix(entries, (domain_dim, domain_dim), tc=typecode) def inner_product(vec1, vec2): @@ -384,3 +391,57 @@ def vec(mat): """ return matrix(mat, (len(mat), 1)) + + +def condition_number(mat): + """ + Return the condition number of the given matrix. + + The condition number of a matrix quantifies how hard it is to do + numerical computation with that matrix. It is usually defined as + the ratio of the norm of the matrix to the norm of its inverse, and + therefore depends on the norm used. One way to compute the condition + number with respect to the 2-norm is as the ratio of the matrix's + largest and smallest singular values. Since we have easy access to + those singular values, that is the algorithm we use. + + The larger the condition number is, the worse the matrix is. + + Parameters + ---------- + mat : matrix + The matrix whose condition number you want. + + Returns + ------- + + float + The nonnegative condition number of ``mat``. + + Examples + -------- + + >>> condition_number(identity(1, typecode='d')) + 1.0 + >>> condition_number(identity(2, typecode='d')) + 1.0 + >>> condition_number(identity(3, typecode='d')) + 1.0 + + >>> A = matrix([[2,1],[1,2]], tc='d') + >>> abs(condition_number(A) - 3.0) < options.ABS_TOL + True + + >>> A = matrix([[2,1j],[-1j,2]], tc='z') + >>> abs(condition_number(A) - 3.0) < options.ABS_TOL + True + + """ + num_eigs = min(mat.size) + eigs = matrix(0, (num_eigs, 1), tc='d') + gesdd(mat, eigs) + + if len(eigs) > 0: + return eigs[0]/eigs[-1] + else: + return 0