1import numpy as np 2from numpy.linalg import svd 3 4 5def rank_est(A, atol=1e-13, rtol=0): 6 """ 7 Estimate the rank (i.e. the dimension of the nullspace) of a matrix. 8 9 The algorithm used by this function is based on the singular value 10 decomposition of `A`. 11 12 Parameters 13 ---------- 14 A : array_like(float, ndim=1 or 2) 15 A should be at most 2-D. A 1-D array with length n will be 16 treated as a 2-D with shape (1, n) 17 atol : scalar(float), optional(default=1e-13) 18 The absolute tolerance for a zero singular value. Singular 19 values smaller than `atol` are considered to be zero. 20 rtol : scalar(float), optional(default=0) 21 The relative tolerance. Singular values less than rtol*smax are 22 considered to be zero, where smax is the largest singular value. 23 24 Returns 25 ------- 26 r : scalar(int) 27 The estimated rank of the matrix. 28 29 Note: If both `atol` and `rtol` are positive, the combined tolerance 30 is the maximum of the two; that is: 31 32 tol = max(atol, rtol * smax) 33 34 Note: Singular values smaller than `tol` are considered to be zero. 35 36 See also 37 -------- 38 numpy.linalg.matrix_rank 39 matrix_rank is basically the same as this function, but it does 40 not provide the option of the absolute tolerance. 41 42 """ 43 44 A = np.atleast_2d(A) 45 s = svd(A, compute_uv=False) 46 tol = max(atol, rtol * s[0]) 47 rank = int((s >= tol).sum()) 48 49 return rank 50 51 52def nullspace(A, atol=1e-13, rtol=0): 53 """ 54 Compute an approximate basis for the nullspace of A. 55 56 The algorithm used by this function is based on the singular value 57 decomposition of `A`. 58 59 Parameters 60 ---------- 61 A : array_like(float, ndim=1 or 2) 62 A should be at most 2-D. A 1-D array with length k will be 63 treated as a 2-D with shape (1, k) 64 atol : scalar(float), optional(default=1e-13) 65 The absolute tolerance for a zero singular value. Singular 66 values smaller than `atol` are considered to be zero. 67 rtol : scalar(float), optional(default=0) 68 The relative tolerance. Singular values less than rtol*smax are 69 considered to be zero, where smax is the largest singular value. 70 71 Returns 72 ------- 73 ns : array_like(float, ndim=2) 74 If `A` is an array with shape (m, k), then `ns` will be an array 75 with shape (k, n), where n is the estimated dimension of the 76 nullspace of `A`. The columns of `ns` are a basis for the 77 nullspace; each element in numpy.dot(A, ns) will be 78 approximately zero. 79 80 Note: If both `atol` and `rtol` are positive, the combined tolerance 81 is the maximum of the two; that is: 82 83 tol = max(atol, rtol * smax) 84 85 Note: Singular values smaller than `tol` are considered to be zero. 86 87 """ 88 89 A = np.atleast_2d(A) 90 u, s, vh = svd(A) 91 tol = max(atol, rtol * s[0]) 92 nnz = (s >= tol).sum() 93 ns = vh[nnz:].conj().T 94 95 return ns 96