Implement normalize and normalize!

NEWS.md

@@ -43,18 +43,28 @@ Library improvements
 
   * `cov` and `cor` don't use keyword arguments anymore and are therefore now type stable ([#13465]).
 
-  * New method for generic QR with column pivoting ([#13480]).
+  * Linear algebra:
 
-  * A new `SparseVector` type allows for one-dimensional sparse arrays. Slicing
-    and reshaping sparse matrices now return vectors when appropriate. The
-    `sparsevec` function returns a one-dimensional sparse vector instead of a
-    one-column sparse matrix.
+    * New `normalize` and `normalize!` convenience functions for normalizing
+      vectors ([#13681]).
+
+    * QR
+
+      * New method for generic QR with column pivoting ([#13480]).
+
+      * New method for polar decompositions of `AbstractVector`s ([#13681]).
+
+    * A new `SparseVector` type allows for one-dimensional sparse arrays.
+      Slicing and reshaping sparse matrices now return vectors when
+      appropriate. The `sparsevec` function returns a one-dimensional sparse
+      vector instead of a one-column sparse matrix. ([#13440])
 
 Deprecated or removed
 ---------------------
 
-  * The method `A_ldiv_B!(SparseMatrixCSC, StrideVecOrMat)` has been deprecated in favor
-    of versions that require the matrix to in factored form ([#13496]).
+  * The method `A_ldiv_B!(SparseMatrixCSC, StrideVecOrMat)` has been deprecated
+    in favor of versions that require the matrix to be in factored form
+    ([#13496]).
 
 Julia v0.4.0 Release Notes
 ==========================

base/exports.jl

@@ -676,6 +676,8 @@ export
     lufact,
     lyap,
     norm,
+    normalize,
+    normalize!,
     nullspace,
     ordschur!,
     ordschur,

base/linalg.jl

@@ -95,6 +95,8 @@ export
     lufact!,
     lyap,
     norm,
+    normalize,
+    normalize!,
     nullspace,
     ordschur!,
     ordschur,

base/linalg/generic.jl

@@ -8,7 +8,7 @@ scale(s::Number, X::AbstractArray) = s*X
 # For better performance when input and output are the same array
 # See https://github.com/JuliaLang/julia/issues/8415#issuecomment-56608729
 function generic_scale!(X::AbstractArray, s::Number)
-    for I in eachindex(X)
+    @simd for I in eachindex(X)
         @inbounds X[I] *= s
     end
     X
@@ -529,3 +529,71 @@ function isapprox{T<:Number,S<:Number}(x::AbstractArray{T}, y::AbstractArray{S};
     d = norm(x - y)
     return isfinite(d) ? d <= atol + rtol*max(norm(x), norm(y)) : x == y
 end
+
+"""
+    normalize!(v, [p=2])
+
+Normalize the vector `v` in-place with respect to the `p`-norm.
+
+# Inputs
+
+- `v::AbstractVector` - vector to be normalized
+- `p::Real` - The `p`-norm to normalize with respect to. Default: 2
+
+# Output
+
+- `v` - A unit vector being the input vector, rescaled to have norm 1.
+        The input vector is modified in-place.
+
+# See also
+
+`normalize`, `qr`
+
+"""
+function normalize!(v::AbstractVector, p::Real=2)
+    nrm = norm(v, p)
+    __normalize!(v, nrm)
+end
+
+@inline function __normalize!{T<:AbstractFloat}(v::AbstractVector{T}, nrm::T)
+    #The largest positive floating point number whose inverse is less than
+    #infinity
+    const δ = inv(prevfloat(typemax(T)))
+
+    if nrm ≥ δ #Safe to multiply with inverse
+        invnrm = inv(nrm)
+        scale!(v, invnrm)
+
+    else #Divide by norm; slower but more correct
+        #Note 2015-10-19: As of Julia 0.4, @simd does not vectorize floating
+        #point division, although vectorized intrinsics like DIVPD exist. I
+        #will leave the @simd annotation in, hoping that someone will improve
+        #the @simd macro in the future. - cjh
+        @inbounds @simd for i in eachindex(v)
+            v[i] /= nrm
+        end
+    end
+
+    v
+end
+
+"""
+    normalize(v, [p=2])
+
+Normalize the vector `v` with respect to the `p`-norm.
+
+# Inputs
+
+- `v::AbstractVector` - vector to be normalized
+- `p::Real` - The `p`-norm to normalize with respect to. Default: 2
+
+# Output
+
+- `v` - A unit vector being a copy of the input vector, scaled to have norm 1
+
+# See also
+
+`normalize!`, `qr`
+"""
+normalize(v::AbstractVector, p::Real=2) = v/norm(v, p)
+

base/linalg/qr.jl

@@ -108,6 +108,49 @@ function _qr(A::Union{Number, AbstractMatrix}, ::Type{Val{true}}; thin::Bool=tru
     full(getq(F), thin=thin), F[:R]::Matrix{eltype(F)}, F[:p]::Vector{BlasInt}
 end
 
+"""
+    qr(v::AbstractVector)
+
+Computes the polar decomposition of a vector.
+
+# Input
+- `v::AbstractVector` - vector to normalize
+
+# Outputs
+- `w` - A unit vector in the direction of `v`
+- `r` - The norm of `v`
+
+# See also
+
+`normalize`, `normalize!`, `qr!`
+"""
+function qr(v::AbstractVector)
+    nrm = norm(v)
+    v/nrm, nrm
+end
+
+"""
+    qr!(v::AbstractVector)
+
+Computes the polar decomposition of a vector.
+
+# Input
+- `v::AbstractVector` - vector to normalize
+
+# Outputs
+- `w` - A unit vector in the direction of `v`
+- `r` - The norm of `v`
+
+# See also
+
+`normalize`, `normalize!`, `qr`
+"""
+function qr!(v::AbstractVector)
+    nrm = norm(v)
+    __normalize!(v, nrm), nrm
+end
+
+
 convert{T}(::Type{QR{T}},A::QR) = QR(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ))
 convert{T}(::Type{Factorization{T}}, A::QR) = convert(QR{T}, A)
 convert{T}(::Type{QRCompactWY{T}},A::QRCompactWY) = QRCompactWY(convert(AbstractMatrix{T}, A.factors), convert(AbstractMatrix{T}, A.T))

test/linalg/generic.jl

@@ -150,3 +150,24 @@ let
     @test LinAlg.axpy!(α, x, deepcopy(y)) == x .* Matrix{Int}[α]
     @test LinAlg.axpy!(α, x, deepcopy(y)) != Matrix{Int}[α] .* x
 end
+
+let
+    v = [3.0, 4.0]
+    @test norm(v) === 5.0
+    w = normalize(v)
+    @test w == [0.6, 0.8]
+    @test norm(w) === 1.0
+    @test norm(normalize!(v) - w, Inf) < eps()
+end
+
+#Test potential overflow in normalize!
+let
+    δ = inv(prevfloat(typemax(Float64)))
+    v = [δ, -δ]
+
+    @test norm(v) === 7.866824069956793e-309
+    w = normalize(v)
+    @test w ≈ [1/√2, -1/√2]
+    @test norm(w) === 1.0
+    @test norm(normalize!(v) - w, Inf) < eps()
+end

test/linalg/qr.jl

@@ -162,3 +162,11 @@ let
     Q=full(qrfact(A)[:Q])
     @test vecnorm(A-Q) < eps()
 end
+
+let
+    debug && println("qr on AbstractVector")
+
+    v = [3.0, 4.0]
+    @test qr(v) == ([0.6, 0.8], 5.0)
+end
+