Add generic pivoted QR and a bit of restructuring of the qrfacts

NEWS.md

@@ -28,6 +28,8 @@ 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]).
+
 Deprecated or removed
 ---------------------
 

base/linalg/dense.jl

@@ -433,7 +433,7 @@ function factorize{T}(A::Matrix{T})
  end
  return lufact(A)
  end
- qrfact(A,typeof(zero(T)/sqrt(zero(T) + zero(T)))<:BlasFloat?Val{true}:Val{false}) # Generic pivoted QR not implemented yet
+ qrfact(A, Val{true})
 end
 
 (\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)

base/linalg/qr.jl

@@ -24,8 +24,7 @@ immutable QRPivoted{T,S<:AbstractMatrix} <: Factorization{T}
 end
 QRPivoted{T}(factors::AbstractMatrix{T}, τ::Vector{T}, jpvt::Vector{BlasInt}) = QRPivoted{T,typeof(factors)}(factors, τ, jpvt)
 
-function qrfact!{T}(A::AbstractMatrix{T}, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false})
- pivot==Val{true} && warn("pivoting only implemented for Float32, Float64, Complex64 and Complex128")
+function qrfactUnblocked!{T}(A::AbstractMatrix{T})
  m, n = size(A)
  τ = zeros(T, min(m,n))
  for k = 1:min(m - 1 + !(T<:Real), n)
@@ -36,12 +35,66 @@ function qrfact!{T}(A::AbstractMatrix{T}, pivot::Union{Type{Val{false}}, Type{Va
  end
  QR(A, τ)
 end
-qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{false}} = Val{false}) = QRCompactWY(LAPACK.geqrt!(A, min(minimum(size(A)), 36))...)
-qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Val{true}}) = QRPivoted(LAPACK.geqp3!(A)...)
-qrfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}) = qrfact!(copy(A), pivot)
-copy_oftype{T}(A::StridedMatrix{T}, ::Type{T}) = copy(A)
-copy_oftype{T,S}(A::StridedMatrix{T}, ::Type{S}) = convert(AbstractMatrix{S}, A)
-qrfact{T}(A::StridedMatrix{T}, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}) = qrfact!(copy_oftype(A, typeof(one(T)/norm(one(T)))), pivot)
+
+# Find index for columns with largest two norm
+function indmaxcolumn(A::StridedMatrix)
+ mm = norm(slice(A, :, 1))
+ ii = 1
+ for i = 2:size(A, 2)
+ mi = norm(slice(A, :, i))
+ if abs(mi) > mm
+ mm = mi
+ ii = i
+ end
+ end
+ return ii
+end
+
+function qrfactPivotedUnblocked!(A::StridedMatrix)
+ m, n = size(A)
+ piv = collect(UnitRange{BlasInt}(1,n))
+ τ = Array(eltype(A), min(m,n))
+ for j = 1:min(m,n)
+
+ # Find column with maximum norm in trailing submatrix
+ jm = indmaxcolumn(slice(A, j:m, j:n)) + j - 1
+
+ if jm != j
+ # Flip elements in pivoting vector
+ tmpp = piv[jm]
+ piv[jm] = piv[j]
+ piv[j] = tmpp
+
+ # Update matrix with
+ for i = 1:m
+ tmp = A[i,jm]
+ A[i,jm] = A[i,j]
+ A[i,j] = tmp
+ end
+ end
+
+ # Compute reflector of columns j
+ x = slice(A, j:m, j)
+ τj = LinAlg.reflector!(x)
+ τ[j] = τj
+
+ # Update trailing submatrix with reflector
+ LinAlg.reflectorApply!(x, τj, sub(A, j:m, j+1:n))
+ end
+ return LinAlg.QRPivoted{eltype(A), typeof(A)}(A, τ, piv)
+end
+
+# LAPACK version
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}) = QRCompactWY(LAPACK.geqrt!(A, min(minimum(size(A)), 36))...)
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{true}}) = QRPivoted(LAPACK.geqp3!(A)...)
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}) = qrfact!(A, Val{false})
+qrfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = qrfact!(copy(A), args...)
+
+# Generic fallbacks
+qrfact!(A::StridedMatrix, ::Type{Val{false}}) = qrfactUnblocked!(A)
+qrfact!(A::StridedMatrix, ::Type{Val{true}}) = qrfactPivotedUnblocked!(A)
+qrfact!(A::StridedMatrix) = qrfact!(A, Val{false})
+qrfact{T}(A::StridedMatrix{T}, args...) = qrfact!(copy_oftype(A, typeof(zero(T)/norm(one(T)))), args...)
 qrfact(x::Number) = qrfact(fill(x,1,1))
 
 qr(A::Union{Number, AbstractMatrix}, pivot::Union{Type{Val{false}}, Type{Val{true}}}=Val{false}; thin::Bool=true) =
@@ -336,7 +389,7 @@ function A_ldiv_B!{T}(A::QR{T}, B::StridedMatrix{T})
  m, n = size(A)
  minmn = min(m,n)
  mB, nB = size(B)
- Ac_mul_B!(A[:Q], sub(B, 1:m, 1:nB))
+ Ac_mul_B!(A[:Q], sub(B, 1:m, :))
  R = A[:R]
  @inbounds begin
  if n > m # minimum norm solution
@@ -389,6 +442,7 @@ function A_ldiv_B!(A::QRPivoted, B::StridedMatrix)
  B[1:size(A.factors, 2),:] = sub(B, 1:size(A.factors, 2), :)[invperm(A.jpvt),:]
  B
 end
+
 function \{TA,Tb}(A::Union{QR{TA},QRCompactWY{TA},QRPivoted{TA}}, b::StridedVector{Tb})
  S = promote_type(TA,Tb)
  m,n = size(A)

test/linalg/lapack.jl

@@ -523,8 +523,7 @@ for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
  end
 
  ## QR
- FJulia = invoke(qrfact!, Tuple{AbstractMatrix, Type{Val{false}}},
- copy(A), Val{false})
+ FJulia = LinAlg.qrfactUnblocked!(copy(A))
  FLAPACK = Base.LinAlg.LAPACK.geqrf!(copy(A))
  @test_approx_eq FJulia.factors FLAPACK[1]
  @test_approx_eq FJulia.τ FLAPACK[2]

test/linalg/qr.jl

@@ -73,21 +73,20 @@ debug && println("Thin QR decomposition (without pivoting)")
  @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
  end
 
-debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
- if eltya <: BlasFloat
- @inferred qrfact(a, Val{true})
- @inferred qr(a, Val{true})
- end
+debug && println("(Automatic) Fat (pivoted) QR decomposition")
+ @inferred qrfact(a, Val{true})
+ @inferred qr(a, Val{true})
+
  qrpa = factorize(a[1:n1,:])
  q,r = qrpa[:Q], qrpa[:R]
  @test_throws KeyError qrpa[:Z]
- if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
+ p = qrpa[:p]
  @test_approx_eq q'*full(q, thin=false) eye(n1)
  @test_approx_eq q*full(q, thin=false)' eye(n1)
  @test_approx_eq (UpperTriangular(eye(eltya,size(q,2)))*q')*full(q, thin=false) eye(n1)
  @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:n1,p] : a[1:n1,:]
- @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:n1,:]
- @test_approx_eq isa(qrpa, QRPivoted) ? q*r*qrpa[:P].' : q*r a[1:n1,:]
+ @test_approx_eq q*r[:,invperm(p)] a[1:n1,:]
+ @test_approx_eq q*r*qrpa[:P].' a[1:n1,:]
  @test_approx_eq_eps a[1:n1,:]*(qrpa\b[1:n1]) b[1:n1] 5000ε
  @test_approx_eq full(qrpa) a[1:5,:]
  @test_throws DimensionMismatch q*b[1:n1 + 1]
@@ -96,15 +95,15 @@ debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is onl
  @test eye(eltyb,n1)*q ≈ convert(AbstractMatrix{tab},q)
  end
 
-debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
+debug && println("(Automatic) Thin (pivoted) QR decomposition")
  qrpa = factorize(a[:,1:n1])
  q,r = qrpa[:Q], qrpa[:R]
  @test_throws KeyError qrpa[:Z]
- if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia
+ p = qrpa[:p]
  @test_approx_eq q'*full(q, thin=false) eye(n)
  @test_approx_eq q*full(q, thin=false)' eye(n)
- @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:n1]
- @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:n1]
+ @test_approx_eq q*r a[:,p]
+ @test_approx_eq q*r[:,invperm(p)] a[:,1:n1]
  @test_approx_eq full(qrpa) a[:,1:5]
  @test_throws DimensionMismatch q*b[1:n1 + 1]
  @test_throws DimensionMismatch b[1:n1 + 1]*q'
@@ -124,10 +123,8 @@ debug && println("Matmul with QR factorizations")
  @test_approx_eq A_mul_Bc!(full(q,thin=false),q) eye(n)
  @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
  @test_throws BoundsError size(q,-1)
- if eltya == BigFloat
- @test_throws DimensionMismatch Base.LinAlg.A_mul_B!(q,zeros(eltya,n1+1))
- @test_throws DimensionMismatch Base.LinAlg.Ac_mul_B!(q,zeros(eltya,n1+1))
- end
+ @test_throws DimensionMismatch Base.LinAlg.A_mul_B!(q,zeros(eltya,n1+1))
+ @test_throws DimensionMismatch Base.LinAlg.Ac_mul_B!(q,zeros(eltya,n1+1))
 
  qra = qrfact(a[:,1:n1], Val{false})
  q,r = qra[:Q], qra[:R]
@@ -165,8 +162,3 @@ let
  Q=full(qrfact(A)[:Q])
  @test vecnorm(A-Q) < eps()
 end
-
-qrpa = qrfact(areal,Val{true})
-qrpa = convert(QRPivoted{BigFloat},qrpa)
-@test qrpa \ zeros(BigFloat,n) ≈ zeros(BigFloat,n)
-@test qrpa \ zeros(BigFloat,n,n) ≈ zeros(BigFloat,n,n)