This is an archival dump of old wiki content --- see scipy.org for current material.

# Rank and nullspace of a matrix

The following module, rank_nullspace.py, provides the functions rank() and nullspace(). (Note that NumPy already provides the function matrix_rank(); the function given here allows an absolute tolerance to be specified along with a relative tolerance.)

rank_nullspace.py

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

Here's a demonstration script.

```   1 import numpy as np
2
3 from rank_nullspace import rank, nullspace
4
5
6 def checkit(a):
7     print "a:"
8     print a
9     r = rank(a)
10     print "rank is", r
11     ns = nullspace(a)
12     print "nullspace:"
13     print ns
14     if ns.size > 0:
15         res = np.abs(np.dot(a, ns)).max()
16         print "max residual is", res
17
18
19 print "-"*25
20
21 a = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
22 checkit(a)
23
24 print "-"*25
25
26 a = np.array([[0.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
27 checkit(a)
28
29 print "-"*25
30
31 a = np.array([[0.0, 1.0, 2.0, 4.0], [1.0, 2.0, 3.0, 4.0]])
32 checkit(a)
33
34 print "-"*25
35
36 a = np.array([[1.0,   1.0j,   2.0+2.0j],
37               [1.0j, -1.0,   -2.0+2.0j],
38               [0.5,   0.5j,   1.0+1.0j]])
39 checkit(a)
40
41 print "-"*25
```

And here is the output of the script.

```-------------------------
a:
[[ 1.  2.  3.]
[ 4.  5.  6.]
[ 7.  8.  9.]]
rank is 2
nullspace:
[[-0.40824829]
[ 0.81649658]
[-0.40824829]]
max residual is 4.4408920985e-16
-------------------------
a:
[[ 0.  2.  3.]
[ 4.  5.  6.]
[ 7.  8.  9.]]
rank is 3
nullspace:
[]
-------------------------
a:
[[ 0.  1.  2.  4.]
[ 1.  2.  3.  4.]]
rank is 2
nullspace:
[[-0.48666474 -0.61177492]
[-0.27946883  0.76717915]
[ 0.76613356 -0.15540423]
[-0.31319957 -0.11409267]]
max residual is 3.88578058619e-16
-------------------------
a:
[[ 1.0+0.j   0.0+1.j   2.0+2.j ]
[ 0.0+1.j  -1.0+0.j  -2.0+2.j ]
[ 0.5+0.j   0.0+0.5j  1.0+1.j ]]
rank is 1
nullspace:
[[ 0.00000000-0.j         -0.94868330-0.j        ]
[ 0.13333333+0.93333333j  0.00000000-0.10540926j]
[ 0.20000000-0.26666667j  0.21081851-0.21081851j]]
max residual is 1.04295984227e-15
-------------------------```

SciPy: Cookbook/RankNullspace (last edited 2015-10-24 17:48:23 by anonymous)