[SciPy-dev] small bug(s) in scipy.linalg

[Alexander Schmolck]

I think scipy.linalg.basic has a few issues with 1D vectors vs row and column
vectors (respective shapes: (n,), (1,n), (n,1)). It would seem clear to me,
for example that `norm` ought to produce exactly the same with a 1D vector as
with a row vector or column vector. Currently it doesn't, and the very common
case of the 2 norm for a vector is, I think, handled inefficiently.

I've created a "fixed" the `norm` function according to these criteria, but
before I submit a properly tested patch, I wanted to make sure that this
behavior is indeed regarded as buggy.

Here is my modified and completely untested version, just to give you the
picture:

UNTESTED CODE

def norm(x, ord=2):
    """matrix and vector norm.

    Inputs:

      x -- a rank-1 (vector) or rank-2 (matrix) array
      ord -- the order of norm.

     Comments:

       For vectors ord can be any real number including Inf or -Inf.
         ord = Inf, computes the maximum of the magnitudes
         ord = -Inf, computes minimum of the magnitudes
         ord is finite, computes sum(abs(x)**ord)**(1.0/ord)

       For matrices ord can only be + or - 1, 2, Inf.
         ord = 2 computes the largest singular value
         ord = -2 computes the smallest singular value
         ord = 1 computes the largest row sum
         ord = -1 computes the smallest row sum
         ord = Inf computes the largest column sum
         ord = -Inf computes the smallest column sum
    """
    x = asarray(x)
    nd = len(x.shape)
    Inf = scipy_base.Inf
    if nd > 2:
        raise ValueError, "Improper number of dimensions to norm."
    if nd == 2:
        # a 'real' matrix, i.e. not a column or row vector (Nx1 or 1xN)
        if min(self.me.shape) != 1:
            if ord == 2:
                return scipy_base.amax(decomp.svd(x)[1])
            elif ord == -2:
                return scipy_base.amin(decomp.svd(x)[1])
            elif ord == 1:
                return scipy_base.amax(scipy_base.sum(abs(x)))
            elif ord == Inf:
                return scipy_base.amax(scipy_base.sum(abs(x),axis=1))
            elif ord == -1:
                return scipy_base.amin(scipy_base.sum(abs(x)))
            elif ord == -Inf:
                return scipy_base.amin(scipy_base.sum(abs(x),axis=1))
            else:
                raise ValueError, "Invalid norm order for matrices."
        # Nx1 or 1xN to 1D vector
        else:
            x = ravel(x)
    # a vector
    if ord == Inf:
        return scipy_base.amax(abs(x))
    elif ord == -Inf:
        return scipy_base.amin(abs(x))
    elif ord == 2:
        return scipy_base.sqrt(dot(a, a))
    else:
        return scipy_base.sum(abs(x)**ord)**(1.0/ord)