[SciPy-dev] small bug(s) in scipy.linalg
Alexander Schmolck
a.schmolck at gmx.net
Sun Feb 9 18:13:28 CST 2003
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)
More information about the Scipy-dev
mailing list