Search code examples
scipyorthogonal

What is the "automatic cutoff" in scipy.linalg.orth()?


http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.orth.html

And how can I tune it? I am getting fewer basis than expected. Google gives me no obviously useful result on the first few pages.


Solution

  • orth uses the singular value decomposition. The definition of orth currently lives in https://github.com/scipy/scipy/blob/master/scipy/linalg/decomp_svd.py, and as I write this, the complete code for orth (dropping the docstring) is:

    def orth(A):
        u, s, vh = svd(A, full_matrices=False)
        M, N = A.shape
        eps = numpy.finfo(float).eps
        tol = max(M, N) * numpy.amax(s) * eps
        num = numpy.sum(s > tol, dtype=int)
        Q = u[:, :num]
        return Q
    

    The result is that singular values less than tol are considered to be 0, and those directions are not considered part of the range of A. tol is a relative tolerance: it is set to max(M, N)*eps times the maximum singular value of A, where eps is the floating point machine epsilon.

    orth provides no argument to control how tol is computed, but as you can see, the function is only a few lines. If you want to use a different method for deciding which singular values to ignore, you can use orth as a starting point for writing your own function.