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 
  33     See also
  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
-------------------------