On Sun, Jan 30, 2011 at 9:15 AM, Sturla Molden <stu...@molden.no> wrote: > Den 30.01.2011 07:28, skrev Algis Kabaila: >> Why not simply numply.linalg.cond? This gives the condition >> number directly (and presumably performs the inspection of >> sv's). Or do you think that sv's give more useful information? > > You can use the singular value decomposition to invert the matrix, solve > linear systems and solve least squares problems. Looking at the topic > you don't just want to compute condition numbers, but invert the > ill-conditioned (nearly singular) matrix. > > Say you want to do: > > invA = np.linalg.inv(a) > > With matrix a nearly singular, you can proceed like this: > > U, s, Vh = np.linalg.svd(a, full_matrices=False) > > Edit small singular values. This will add a small bias but reduce > rounding error: > > singular = s < threshold > invS = 1/s > invS[singular] = 0 # truncate inf to 0 actually helps... > > Et voilá: > > invA = np.dot(np.dot(U,np.diag(invS)),Vh) > > I hope this helps :) > > There is a chapter on SVD in "Numerical Receipes" and almost any linear > algebra textbook. > > Just to verify: > > >>> a = np.diag([1,2,3]) > >>> np.linalg.inv(a) > array([[ 1. , 0. , 0. ], > [ 0. , 0.5 , 0. ], > [ 0. , 0. , 0.33333333]]) > >>> U, s, Vh = np.linalg.svd(a, full_matrices=False) > >>> np.dot(np.dot(U,np.diag(1/s)),Vh) > array([[ 1. , 0. , 0. ], > [ 0. , 0.5 , 0. ], > [ 0. , 0. , 0.33333333]]) > > > Sturla > > If the matrix is not full rank then it is not invertible (http://en.wikipedia.org/wiki/Matrix_inverse) and (matrix) rank (http://en.wikipedia.org/wiki/Matrix_rank) can be computed from the above code. You do have to beware that you can get a generalized inverse (numpy.linalg provides pinv) when your system has an infinite number of solutions. (A generalized inverse is not the problem, the problem is when you expect a unique solution and do not get it.)
Bruce _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
