From ed86c0467de2f8903d9a18c28fa478412cd9a52e Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 30 Oct 2016 16:27:45 -0400 Subject: [PATCH] Add the condition_number() function to the matrices module. --- dunshire/matrices.py | 56 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/dunshire/matrices.py b/dunshire/matrices.py index 66c2176..2d4bb17 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 @@ -388,3 +388,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 -- 2.44.2