Skip to content

Commit

Permalink
Provides specialized Bidiagonal linear solver and eigensystem solver
Browse files Browse the repository at this point in the history 
Addresses #3688

- New tests of Float16 and BigFloat
- Miscellaneous improvements to Bidiagonal tests
- Clean up bidiagonal methods, particularly to reduct overtyping
- Adds BigFloat(x::Float16)
  • Loading branch information
@jiahao
jiahao committed Jan 7, 2014
1 parent 4954ac5 commit d66c9c4
Show file tree
Hide file tree
Showing 4 changed files with 160 additions and 103 deletions.
153 changes: 116 additions & 37 deletions base/linalg/bidiag.jl
View file
Original file line number Diff line number Diff line change
@@ -1,34 +1,21 @@
#### Specialized matrix types ####

## Bidiagonal matrices


# Bidiagonal matrices
type Bidiagonal{T} <: AbstractMatrix{T}
dv::Vector{T} # diagonal
ev::Vector{T} # sub/super diagonal
dv::AbstractVector{T} # diagonal
ev::AbstractVector{T} # sub/super diagonal
isupper::Bool # is upper bidiagonal (true) or lower (false)
function Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)
function Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)
length(ev)==length(dv)-1 ? new(dv, ev, isupper) : throw(DimensionMismatch(""))
end
end

Bidiagonal{T<:BlasFloat}(dv::Vector{T}, ev::Vector{T}, isupper::Bool)=Bidiagonal{T}(copy(dv), copy(ev), isupper)
Bidiagonal{T}(dv::Vector{T}, ev::Vector{T}) = error("Did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument.")
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)=Bidiagonal{T}(copy(dv), copy(ev), isupper)
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}) = error("Did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument.")

#Convert from BLAS uplo flag to boolean internal
function Bidiagonal{T<:BlasFloat}(dv::Vector{T}, ev::Vector{T}, uplo::BlasChar)
if uplo=='U'
isupper = true
elseif uplo=='L'
isupper = false
else
error("Bidiagonal can only be upper 'U' or lower 'L' but you said '$uplo''")
end
Bidiagonal{T}(copy(dv), copy(ev), isupper)
end
Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::BlasChar) = Bidiagonal(copy(dv), copy(ev), uplo=='U' ? true : uplo=='L' ? false : error("Bidiagonal can only be upper 'U' or lower 'L' but you said '$uplo''"))

function Bidiagonal{Td<:Number,Te<:Number}(dv::Vector{Td}, ev::Vector{Te}, isupper::Bool)
T = promote(Td,Te)
function Bidiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te}, isupper::Bool)
T = promote_type(Td,Te)
Bidiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev), isupper)
end

Expand All @@ -37,7 +24,6 @@ Bidiagonal(A::AbstractMatrix, isupper::Bool)=Bidiagonal(diag(A), diag(A, isupper
#Converting from Bidiagonal to dense Matrix
full{T}(M::Bidiagonal{T}) = convert(Matrix{T}, M)
convert{T}(::Type{Matrix{T}}, A::Bidiagonal{T})=diagm(A.dv) + diagm(A.ev, A.isupper?1:-1)
promote_rule{T}(::Type{Matrix{T}}, ::Type{Bidiagonal{T}})=Matrix{T}
promote_rule{T,S}(::Type{Matrix{T}}, ::Type{Bidiagonal{S}})=Matrix{promote_type(T,S)}

#Converting from Bidiagonal to Tridiagonal
Expand All @@ -46,9 +32,19 @@ function convert{T}(::Type{Tridiagonal{T}}, A::Bidiagonal{T})
z = zeros(T, size(A)[1]-1)
A.isupper ? Tridiagonal(A.ev, A.dv, z) : Tridiagonal(z, A.dv, A.ev)
end
promote_rule{T}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{T}})=Tridiagonal{T}
promote_rule{T,S}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{S}})=Tridiagonal{promote_type(T,S)}

#################
# BLAS routines #
#################

#Singular values
svdvals{T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
svd {T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))

####################
# Generic routines #
####################

function show(io::IO, M::Bidiagonal)
println(io, summary(M), ":")
Expand All @@ -62,11 +58,10 @@ size(M::Bidiagonal) = (length(M.dv), length(M.dv))
size(M::Bidiagonal, d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(M.dv) : 1)

#Elementary operations
copy(M::Bidiagonal) = Bidiagonal(copy(M.dv), copy(M.ev), copy(M.isupper))
round(M::Bidiagonal) = Bidiagonal(round(M.dv), round(M.ev), M.isupper)
iround(M::Bidiagonal) = Bidiagonal(iround(M.dv), iround(M.ev), M.isupper)
for func in (conj, copy, round, iround)
func(M::Bidiagonal) = Bidiagonal(func(M.dv), func(M.ev), M.isupper)
end

conj(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), M.isupper)
transpose(M::Bidiagonal) = Bidiagonal(M.dv, M.ev, !M.isupper)
ctranspose(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), !M.isupper)

Expand All @@ -92,17 +87,101 @@ end
/(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper)
==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)

SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal)
SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular)
*(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B)

# solver uses tridiagonal gtsv!
function \{T<:BlasFloat}(M::Bidiagonal{T}, rhs::StridedVecOrMat{T})
stride(rhs, 1)==1 || solve(M, rhs) #generic fallback
z = zeros(size(M, 1) - 1)
apply(LAPACK.gtsv!, M.isupper ? (z, copy(M.dv), copy(M.ev), copy(rhs)) : (copy(M.ev), copy(M.dv), z, copy(rhs)))
#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
@eval begin
($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B)
#($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B)
end
end


#Linear solvers
\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), b::AbstractVector{T}) = naivesub!(A, b, similar(b))
\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...)
A_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(A,b)
At_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(transpose(A),b)
Ac_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(ctranspose(A),b)
for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
function ($func)(A::Bidiagonal, B::AbstractMatrix)
for i=1:size(B,2)
($func)(A, B[:,i])
end
B
end
end end
for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
function ($func)(A::Bidiagonal, B::AbstractMatrix)
for i=1:size(B,1)
($func)(A, B[i,:])
end
B
end
end end

function inv(A::Union(Bidiagonal, Triangular))
B = eye(eltype(A), size(A, 1))
for i=1:size(B,2)
naivesub!(A, B[:,i])
end
B
end

#Generic solver using naive substitution
function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector=b)
N = size(A, 2)
N==length(b)==length(x) || throw(DimensionMismatch(""))

if !A.isupper #do forward substitution
for j = 1:N
x[j] = b[j]
j>1 && (x[j] -= A.ev[j-1] * x[j-1])
x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
end
else #do backward substitution
for j = N:-1:1
x[j] = b[j]
j<N && (x[j] -= A.ev[j] * x[j+1])
x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
end
end
x
end

#Wrap bdsdc to compute singular values and vectors
svdvals{T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
svd {T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))
# Eigensystems
eigvals(M::Bidiagonal) = M.dv
function eigvecs{T}(M::Bidiagonal{T})
n = length(M.dv)
Q=Array(T, n, n)
blks = [0; find(x->x==0, M.ev); n]
if M.isupper
for idx_block=1:length(blks)-1, i=blks[idx_block]+1:blks[idx_block+1] #index of eigenvector
v=zeros(T, n)
v[blks[idx_block]+1] = one(T)
for j=blks[idx_block]+1:i-1 #Starting from j=i, eigenvector elements will be 0
v[j+1] = (M.dv[i]-M.dv[j])/M.ev[j] * v[j]
end
Q[:,i] = v/norm(v)
end
else
for idx_block=1:length(blks)-1, i=blks[idx_block+1]:-1:blks[idx_block]+1 #index of eigenvector
v=zeros(T, n)
v[blks[idx_block+1]] = one(T)
for j=(blks[idx_block+1]-1):-1:max(1,(i-1)) #Starting from j=i, eigenvector elements will be 0
v[j] = (M.dv[i]-M.dv[j+1])/M.ev[j] * v[j+1]
end
Q[:,i] = v/norm(v)
end
end
Q #Actually Triangular
end
eigfact(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))

#Singular values
function svdfact(M::Bidiagonal, thin::Bool=true)
U, S, V = svd(M)
SVD(U, S, V')
end
30 changes: 0 additions & 30 deletions base/linalg/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,25 +102,6 @@ for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
end
end

A_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(A,b)
At_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(transpose(A),b)
Ac_ldiv_B!(A::Triangular, b::AbstractVector)=naivesub!(ctranspose(A),b)
for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
function ($func)(A::Triangular, B::AbstractMatrix)
for i=1:size(B,2)
($func)(A, B[:,i])
end
B
end
end end
for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
function ($func)(A::Triangular, B::AbstractMatrix)
for i=1:size(B,1)
($func)(A, B[i,:])
end
B
end
end end
#Generic solver using naive substitution
function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b)
N = size(A, 2)
Expand Down Expand Up @@ -152,17 +133,6 @@ function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b)
x
end

\{T<:Number}(A::Triangular{T}, b::AbstractVector{T}) = naivesub!(A, b, similar(b))
\{T<:Number}(A::Triangular{T}, B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...)

function inv(A::Triangular)
B = eye(eltype(A), size(A, 1))
for i=1:size(B,2)
naivesub!(A, B[:,i])
end
B
end

#Generic eigensystems
eigvals(A::Triangular) = diag(A.UL)
det(A::Triangular) = prod(eigvals(A))
Expand Down
2 changes: 1 addition & 1 deletion base/mpfr.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ BigFloat(x::Integer) = BigFloat(BigInt(x))
BigFloat(x::Union(Bool,Int8,Int16,Int32)) = BigFloat(convert(Clong,x))
BigFloat(x::Union(Uint8,Uint16,Uint32)) = BigFloat(convert(Culong,x))

BigFloat(x::Float32) = BigFloat(float64(x))
BigFloat(x::Union(Float16,Float32)) = BigFloat(float64(x))
BigFloat(x::Rational) = BigFloat(num(x)) / BigFloat(den(x))

convert(::Type{Rational}, x::BigFloat) = convert(Rational{BigInt}, x)
Expand Down
78 changes: 43 additions & 35 deletions test/linalg.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import Base.LinAlg
import Base.LinAlg: BlasFloat
import Base.LinAlg: BlasComplex, BlasFloat, BlasReal

n = 10
srand(1234321)
Expand Down Expand Up @@ -550,15 +550,14 @@ for elty in (Float32, Float64, Complex64, Complex128)
end

#Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences
function test_approx_eq_vecs(a, b)
n = size(a)[1]
@test n==size(b)[1]
elty = typeof(a[1])
@test elty==typeof(b[1])
function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing)
n = size(a, 1)
@test n==size(b,1) && size(a,2)==size(b,2)
if error==nothing error=n^2*(eps(S)+eps(T)) end
for i=1:n
ev1, ev2 = a[:,i], b[:,i]
deviation = min(abs(norm(ev1-ev2)),abs(norm(ev1+ev2)))
@test_approx_eq_eps deviation 0.0 n^2*eps(abs(convert(elty, 1.0)))
@test_approx_eq_eps deviation 0.0 error
end
end

Expand Down Expand Up @@ -605,7 +604,6 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt
end

#SymTridiagonal (symmetric tridiagonal) matrices
n=5
Ainit = randn(n)
Binit = randn(n-1)
for elty in (Float32, Float64)
Expand All @@ -632,46 +630,56 @@ for elty in (Float32, Float64)
test_approx_eq_vecs(v, evecs)
end


#Bidiagonal matrices
dv = randn(n)
ev = randn(n-1)
for elty in (Float32, Float64, Complex64, Complex128)
if (elty == Complex64)
dv += im*randn(n)
ev += im*randn(n-1)
for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
dv = convert(Vector{elty}, randn(n))
ev = convert(Vector{elty}, randn(n-1))
b = convert(Vector{elty}, randn(n))
if (elty <: Complex)
dv += im*convert(Vector{elty}, randn(n))
ev += im*convert(Vector{elty}, randn(n-1))
b += im*convert(Vector{elty}, randn(n))
end
for isupper in (true, false) #Test upper and lower bidiagonal matrices
T = Bidiagonal{elty}(dv, ev, isupper)
T = Bidiagonal(dv, ev, isupper)

@test size(T, 1) == n
@test size(T, 1) == size(T, 2) == n
@test size(T) == (n, n)
@test full(T) == diagm(dv) + diagm(ev, isupper?1:-1)
@test Bidiagonal(full(T), isupper) == T
z = zeros(elty, n)

# idempotent tests
@test conj(conj(T)) == T
@test transpose(transpose(T)) == T
@test ctranspose(ctranspose(T)) == T
for func in (conj, transpose, ctranspose)
@test func(func(T)) == T
end

if (elty <: Real)
Tfull = full(T)
#Linear solver
Tfull = full(T)
condT = cond(complex128(Tfull))
x = T \ b
tx = Tfull \ b
@test norm(x-tx,Inf) <= 4*condT*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf))

#Test eigenvalues/vectors
d1, v1 = eig(T)
d2, v2 = eig((elty<:Complex?complex128:float64)(Tfull))
@test_approx_eq isupper?d1:reverse(d1) d2
if elty <: Real
test_approx_eq_vecs(v1, isupper?v2:v2[:,n:-1:1])
end

if (elty <: BlasReal)
#Test singular values/vectors
@test_approx_eq svdvals(Tfull) svdvals(T)
u1, d1, v1 = svd(Tfull)
u2, d2, v2 = svd(T)
@test_approx_eq d1 d2
test_approx_eq_vecs(u1, u2)
test_approx_eq_vecs(v1, v2)

#Test eigenvalues/vectors
#d1, v1 = eig(Tfull)
#d2, v2 = eigvals(T), eigvecs(T)
#@test_approx_eq d1 d2
#test_approx_eq_vecs(v1, v2)
#@test_approx_eq_eps 0 norm(v1 * diagm(d1) * inv(v1) - Tfull) eps(elty)*n*(n+1)
#@test_approx_eq_eps 0 norm(v2 * diagm(d2) * inv(v2) - Tfull) eps(elty)*n*(n+1)
if elty <: Real
test_approx_eq_vecs(u1, u2)
test_approx_eq_vecs(v1, v2)
end
@test_approx_eq_eps 0 normfro(u2*diagm(d2)*v2'-Tfull) n*max(n^2*eps(relty), normfro(u1*diagm(d1)*v1'-Tfull))
end
end
end
Expand All @@ -692,8 +700,8 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt
@test_approx_eq_eps D*v DM*v n*eps(relty)*(elty<:Complex ? 2:1)
@test_approx_eq_eps D*U DM*U n^2*eps(relty)*(elty<:Complex ? 2:1)
if relty != BigFloat
@test_approx_eq_eps D\v DM\v n*eps(relty)*(elty<:Complex ? 2:1)
@test_approx_eq_eps D\U DM\U n^2*eps(relty)*(elty<:Complex ? 2:1)
@test_approx_eq_eps D\v DM\v 2n^2*eps(relty)*(elty<:Complex ? 2:1)
@test_approx_eq_eps D\U DM\U 2n^3*eps(relty)*(elty<:Complex ? 2:1)
end
for func in (det, trace)
@test_approx_eq_eps func(D) func(DM) n^2*eps(relty)
Expand All @@ -703,7 +711,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt
@test_approx_eq_eps func(D) func(DM) n^2*eps(relty)
end
end
if elty <: Complex && relty <:BlasFloat
if elty <: BlasComplex
for func in (logdet, sqrtm)
@test_approx_eq_eps func(D) func(DM) n^2*eps(relty)
end
Expand Down

0 comments on commit d66c9c4

Please sign in to comment.