Avoid overflow in p norms

[andreasnoack]

Use division instead of multiplication with reciprocal when scaling elements in vecnorm methods. This fixes JuliaLang/LinearAlgebra.jl#222 and supersedes #11789. The logic used here is a little different from the other generic_vecnorms. I try to use the element type of the array when the input is an AbstractVector{T<:Number} and for general iterables I have to require that they are non-empty.

cc: @stevengj

[stevengj]

Didn't the old implementation work for Array{T} for any T that defines norm(x::T)? See #11043. (I should have added a test case for that, grrr.) The new implementation loses this because it requires zero(T) (not to mention abs instead of norm); I don't think that should be necessary.

[andreasnoack]

The new version doesn't require zero(T). It only uses zero(T) when T<:Number. For all other element types it uses zero(x[1]).

[andreasnoack]

...but I'm fine with using norm instead of abs. I'll change that now.

[stevengj]

Can you add a test case for the vecnorm of an array of arrays, if we don't have one yet?

[stevengj]

It's kind of annoying to compute norm(x[1]) twice in this case. Maybe add an i0 argument to generic_vecnorm2 and loop from i = i0:n. That way you can pass i0 == 2 to avoid recomputing the first element, which can be passed via ssq.

[andreasnoack]

I don't think norm will actually be computed here. E.g. for a Vector{Complex{Float64}} I get that

julia> @code_llvm zero(norm(x2[1]))

define double @julia_zero_21182(double) {
top:
  ret double 0.000000e+00
}

[yuyichao]

@code_llvm compute the argument before getting its type.

[yuyichao]

Even with inbounds,

julia> a = Complex128[1, 2, 3]
3-element Array{Complex{Float64},1}:                                              
 1.0+0.0im                                                                        
 2.0+0.0im                                                                        
 3.0+0.0im                                                                        

julia> f(x) = (@inbounds v = zero(norm(x[1])); v)
f (generic function with 1 method)

julia> @code_llvm f(a)

define double @julia_f_21134(%jl_value_t*) {
top:
  %1 = alloca %Complex.12, align 8
  %2 = bitcast %jl_value_t* %0 to %Complex.12**
  %3 = load %Complex.12** %2, align 8
  %4 = load %Complex.12* %3, align 8
  store %Complex.12 %4, %Complex.12* %1, align 8
  %5 = call double @julia_norm_21121(%Complex.12* %1)
  ret double 0.000000e+00
}

I don't think type inference or llvm is smart enough to figure out norm(::Complex) is effect free. There's even more junk emitted without @inbounds.

[stevengj]

I'm mostly worried about it for something like an array of matrices, where norm is expensive.

[StefanKarpinski]

@stevengj, should we just merge this and fix the issues in a separate change?

[andreasnoack]

Hold on a bit. I'll add the tests @stevengj asked for and there might also be few other adjustments based on the Slack discussion this evening.

[staticfloat]

@andreasnoack bumping this since it's slated for 0.4.

[andreasnoack]

I'm on it. I've been reading up on why BLAS implementation and some critique of it. Scaling all elements even when it's not necessary is expensive. I might just fix the overflow problem problem for now and then we can reconsider the algorithm in 0.5.

[andreasnoack]

This one works now. I've patched the existing versions to handle underflow without being dead slow because of the division, i.e. the value are only scaled when necessary. This is also the reason to avoid the BLAS implementation because it scales even when it is not necessary and therefore it is very slow.

cc: @stevengj

[stevengj]

LGTM.

[stevengj]

Oh, darn it, I forgot to check whether you added a test for the array-of-arrays case (you didn't) like I had suggested.

It looks like you broke norm for this case, because you call one, which isn't defined for a general Banach space. e.g. norm(Vector{Int}[[1,2], [3,4]]) no longer works. This is a regression, and needs to be fixed.

[andreasnoack]

I'll look into it.

Den torsdag den 3. september 2015 skrev Steven G. Johnson <
notifications@github.com>:

Oh, darn it, I forgot to check whether you added a test for the
array-of-arrays case (you didn't).

It looks like you broke norm for this case, because you call one, which
isn't defined for a general Banach space. e.g. norm(Vector{Int}[[1,2],
[3,4]]) is broken.

—
Reply to this email directly or view it on GitHub
#12571 (comment).

[stevengj]

I have PR almost ready.