[Numpy-discussion] np.linalg.norm terrible behavior

[Xavier Gnata]

Kahan summation algorithm?
pairwise summation?
Looks like BLAS does not use these algorithms.
I'm curious to know which algo are coded in matlab/scilab.

Another issue is that sum is slow.

Xavier

> for the sum function part,
> It's not a good way to fix but... if you want more accuracy
>
> x=(np.array([1]*10000 + [1e4], dtype=np.float32))
> np.sum(x*x)
> 1.0001e+08
>
> You can sort x from small numbers to bigger numbers before you call sum.
>
> -Kibeom Kim
>
> On Sat, Mar 5, 2011 at 6:27 PM, Xavier Gnata <xavier.gnata at gmail.com 
> <mailto:xavier.gnata at gmail.com>> wrote:
>
>     Hi,
>
>     I got this problem in a real life code and it took me some time to
>     figure out that np.linalg.norm has a terrible numerical behavior.
>     The problem is nicely described here
>     http://fseoane.net/blog/2011/computing-the-vector-norm/
>
>     numpy/linalg/linalg.py claims to be a "high-level Python interface to
>     the LAPACK library".
>     That's great but the low level norm function is coded in pure python :
>         x = asarray(x)
>         if ord is None: # check the default case first and handle it
>     immediately
>             return sqrt(add.reduce((x.conj() * x).ravel().real))
>
>     Is there a way to use *by default* the blas function to compute
>     the norm?
>
>     ok....looks like sum as the same terrible numerical behavior
>     x=(np.array([1e4] + [1]*10000, dtype=np.float32))
>     np.sum(x*x)
>     returns 1e+08 instead of 1.0001e+08
>
>     Moreover, np.linalg.norm is slow compare to blas.
>
>     Is there a way/plan to fix it?
>
>     Xavier
>     _______________________________________________
>     NumPy-Discussion mailing list
>     NumPy-Discussion at scipy.org <mailto:NumPy-Discussion at scipy.org>
>     http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion