Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC:QR decomposition in julia #5526

Merged
merged 1 commit into from
Jan 31, 2014
Merged

RFC:QR decomposition in julia #5526

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,8 @@ Library improvements

* Balancing options for eigenvector calculations for general matrices ([#5428]).

* Mutating linear algebra functions no longer promote ([#5526]).

* Sparse linear algebra

* Faster sparse `kron` ([#4958]).
Expand Down Expand Up @@ -147,6 +149,8 @@ Library improvements

* LU factorization ([#5381] and [#5430])

* QR factorization ([#5526])

* New function `deleteat!` deletes a specified index or indices and
returns the updated collection

Expand Down Expand Up @@ -177,6 +181,8 @@ Deprecated or removed

* `myindexes` has been renamed to `localindexes` ([#5475])

* `factorize!` is deprecated in favor of `factorize`. ([#5526])

[#4042]: https://github.com/JuliaLang/julia/issues/4042
[#5164]: https://github.com/JuliaLang/julia/issues/5164
[#4026]: https://github.com/JuliaLang/julia/issues/4026
Expand Down Expand Up @@ -220,6 +226,7 @@ Deprecated or removed
[#5025]: https://github.com/JuliaLang/julia/pull/5025
[#4888]: https://github.com/JuliaLang/julia/pull/4888
[#5475]: https://github.com/JuliaLang/julia/pull/5475
[#5526]: https://github.com/JuliaLang/julia/pull/5526

Julia v0.2.0 Release Notes
==========================
Expand Down
1 change: 1 addition & 0 deletions base/deprecated.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,6 +178,7 @@ export PipeString
@deprecate cholpfact(A) cholfact(A, :U, pivot=true)
@deprecate symmetrize!(A) Base.LinAlg.copytri!(A, 'U')
@deprecate symmetrize!(A, uplo) Base.LinAlg.copytri!(A, uplo)
@deprecate factorize!(A) factorize(A)

deprecated_ls() = run(`ls -l`)
deprecated_ls(args::Cmd) = run(`ls -l $args`)
Expand Down
1 change: 0 additions & 1 deletion base/exports.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -607,7 +607,6 @@ export
eigvecs,
expm,
eye,
factorize!,
factorize,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good riddance.

givens,
hessfact!,
Expand Down
1 change: 0 additions & 1 deletion base/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,6 @@ export
sqrtm,
eye,
factorize,
factorize!,
givens,
gradient,
hessfact,
Expand Down
8 changes: 4 additions & 4 deletions base/linalg/bunchkaufman.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
## LD for BunchKaufman, UL for CholeskyDense, LU for LUDense and
## define size methods for Factorization types using it.

type BunchKaufman{T<:BlasFloat} <: Factorization{T}
immutable BunchKaufman{T} <: Factorization{T}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was initially ambivalent about making all the special matrices and factorizations immutable, but I'm increasingly convinced that this is the right thing to do.

LD::Matrix{T}
ipiv::Vector{BlasInt}
uplo::Char
Expand All @@ -18,10 +18,10 @@ function bkfact!{T<:BlasComplex}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric
LD, ipiv = (symmetric ? LAPACK.sytrf! : LAPACK.hetrf!)(string(uplo)[1] , A)
BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
end
bkfact!(A::StridedMatrix, args...) = bkfact!(float(A), args...)
bkfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = bkfact!(copy(A), args...)
bkfact(A::StridedMatrix, args...) = bkfact(float(A), args...)
bkfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(copy(A), uplo, symmetric)
bkfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(convert(Matrix{promote_type(Float32,typeof(sqrt(one(T))))},A),uplo,symmetric)

convert{T}(::Type{BunchKaufman{T}},B::BunchKaufman) = BunchKaufman(convert(Matrix{T},B.LD),B.ipiv,B.uplo,B.symmetric)
size(B::BunchKaufman) = size(B.LD)
size(B::BunchKaufman,d::Integer) = size(B.LD,d)
issym(B::BunchKaufman) = B.symmetric
Expand Down
81 changes: 38 additions & 43 deletions base/linalg/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -287,29 +287,21 @@ function rcswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
end
end

function sqrtm{T<:Real}(A::StridedMatrix{T}, cond::Bool)
issym(A) && return sqrtm(Symmetric(A), cond)

function sqrtm{T<:Real}(A::StridedMatrix{T})
issym(A) && return sqrtm(Symmetric(A))
n = chksquare(A)
SchurF = schurfact!(complex(A))
SchurF = schurfact(complex(A))
R = full(sqrtm(Triangular(SchurF[:T])))
retmat = SchurF[:vectors]*R*SchurF[:vectors]'
retmat2= all(imag(retmat) .== 0) ? real(retmat) : retmat
cond ? (retmat2, norm(R)^2/norm(SchurF[:T])) : retmat2
all(imag(retmat) .== 0) ? real(retmat) : retmat
end
function sqrtm{T<:Complex}(A::StridedMatrix{T}, cond::Bool)
ishermitian(A) && return sqrtm(Hermitian(A), cond)

function sqrtm{T<:Complex}(A::StridedMatrix{T})
ishermitian(A) && return sqrtm(Hermitian(A))
n = chksquare(A)
SchurF = schurfact(A)
R = full(sqrtm(Triangular(SchurF[:T])))
retmat = SchurF[:vectors]*R*SchurF[:vectors]'
cond ? (retmat, norm(R)^2/norm(SchurF[:T])) : retmat
SchurF[:vectors]*R*SchurF[:vectors]'
end

sqrtm{T<:Integer}(A::StridedMatrix{T}, cond::Bool) = sqrtm(float(A), cond)
sqrtm{T<:Integer}(A::StridedMatrix{Complex{T}}, cond::Bool) = sqrtm(complex128(A), cond)
sqrtm(A::StridedMatrix) = sqrtm(A, false)
sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
sqrtm(a::Complex) = sqrt(a)

Expand All @@ -327,7 +319,7 @@ function inv(A::Matrix)
return inv(lufact(A))
end

function factorize!{T}(A::Matrix{T})
function factorize{T}(A::Matrix{T})
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I die a little every time I have to read this polyalgorithm. "A convenient factorization" belies a monstrous Kaonashi.

m, n = size(A)
if m == n
if m == 1 return A[1] end
Expand Down Expand Up @@ -377,9 +369,9 @@ function factorize!{T}(A::Matrix{T})
end
if utri1
if (herm & (T <: Complex)) | sym
return ldltd!(SymTridiagonal(diag(A), diag(A, -1)))
return ldltd(SymTridiagonal(diag(A), diag(A, -1)))
end
return lufact!(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
return lufact(Tridiagonal(diag(A, -1), diag(A), diag(A, 1)))
end
end
if utri
Expand All @@ -388,32 +380,36 @@ function factorize!{T}(A::Matrix{T})
if herm
if T <: BlasFloat
C, info = LAPACK.potrf!('U', copy(A))
elseif typeof(one(T)/one(T)) <: BlasFloat
C, info = LAPACK.potrf!('U', float(A))
else
error("Unable to factorize hermitian $(typeof(A)). Try converting to other element type or use explicit factorization.")
else
S = typeof(one(T)/one(T))
if S <: BlasFloat
C, info = LAPACK.potrf!('U', convert(Matrix{S}, A))
else
C, info = S <: Real ? LAPACK.potrf!('U', complex128(A)) : LAPACK.potrf!('U', complex128(A))
end
end
if info == 0 return Cholesky(C, 'U') end
return factorize!(Hermitian(A))
return factorize(Hermitian(A))
end
if sym
if T <: BlasFloat
C, info = LAPACK.potrf!('U', copy(A))
elseif eltype(one(T)/one(T)) <: BlasFloat
C, info = LAPACK.potrf!('U', float(A))
else
error("Unable to factorize symmetric $(typeof(A)). Try converting to other element type or use explicit factorization.")
S = eltype(one(T)/one(T))
if S <: BlasFloat
C, info = LAPACK.potrf!('U', convert(Matrix{S},A))
else
C, info = S <: Real ? LAPACK.potrf!('U', float64(A)) : LAPACK.potrf!('U', complex(A))
end
end
if info == 0 return Cholesky(C, 'U') end
return factorize!(Symmetric(A))
return factorize(Symmetric(A))
end
return lufact!(A)
return lufact(A)
end
return qrfact!(A,pivot=true)
qrfact(A,pivot=T<:BlasFloat)
end

factorize(A::AbstractMatrix) = factorize!(copy(A))

(\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)
function (\)(A::StridedMatrix, B::StridedVecOrMat)
m, n = size(A)
Expand All @@ -424,32 +420,31 @@ function (\)(A::StridedMatrix, B::StridedVecOrMat)
istriu(A) && return \(Triangular(A, :U),B)
return \(lufact(A),B)
end
return qrfact(A,pivot=true)\B
return qrfact(A,pivot=eltype(A)<:BlasFloat)\B
end

## Moore-Penrose inverse
function pinv{T<:BlasFloat}(A::StridedMatrix{T})
m, n = size(A)
(m == 0 || n == 0) && return Array(T, n, m)
SVD = svdfact(A, true)
Sinv = zeros(T, length(SVD[:S]))
index = SVD[:S] .> eps(real(one(T)))*max(m,n)*maximum(SVD[:S])
Sinv[index] = 1.0 ./ SVD[:S][index]
function pinv{T}(A::StridedMatrix{T})
SVD = svdfact(A, thin=true)
S = eltype(SVD[:S])
m, n = size(A)
(m == 0 || n == 0) && return Array(S, n, m)
Sinv = zeros(S, length(SVD[:S]))
index = SVD[:S] .> eps(real(float(one(T))))*max(m,n)*maximum(SVD[:S])
Sinv[index] = one(S) ./ SVD[:S][index]
return SVD[:Vt]'scale(Sinv, SVD[:U]')
end
pinv{T<:Integer}(A::StridedMatrix{T}) = pinv(float(A))
pinv(a::StridedVector) = pinv(reshape(a, length(a), 1))
pinv(x::Number) = one(x)/x

## Basis for null space
function null{T<:BlasFloat}(A::StridedMatrix{T})
function null{T}(A::StridedMatrix{T})
m, n = size(A)
(m == 0 || n == 0) && return eye(T, n)
SVD = svdfact(A, false)
SVD = svdfact(A, thin=false)
indstart = sum(SVD[:S] .> max(m,n)*maximum(SVD[:S])*eps(eltype(SVD[:S]))) + 1
return SVD[:V][:,indstart:end]
end
null{T<:Integer}(A::StridedMatrix{T}) = null(float(A))
null(a::StridedVector) = null(reshape(a, length(a), 1))

function cond(A::StridedMatrix, p::Real=2)
Expand Down
Loading