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