Skip to content

Add det(Symmetric/Hermitian) based on bkfact. #13166

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

Merged
merged 1 commit into from
Sep 18, 2015
Merged

Add det(Symmetric/Hermitian) based on bkfact. #13166

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
25 changes: 25 additions & 0 deletions base/linalg/bunchkaufman.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,31 @@ function A_ldiv_B!{T<:BlasComplex}(B::BunchKaufman{T}, R::StridedVecOrMat{T})
(issym(B) ? LAPACK.sytrs! : LAPACK.hetrs!)(B.uplo, B.LD, B.ipiv, R)
end

function det(F::BunchKaufman)
M = F.LD
p = F.ipiv
n = size(F.LD, 1)

i = 1
d = one(eltype(F))
while i <= n
if p[i] > 0
# 1x1 pivot case
d *= M[i,i]
i += 1
else
# 2x2 pivot case. Make sure not to square before the subtraction by scaling with the off-diagonal element. This is safe because the off diagonal is always large for 2x2 pivots.
if F.uplo == 'U'
d *= M[i, i + 1]*(M[i,i]/M[i, i + 1]*M[i + 1, i + 1] - (issym(F) ? M[i, i + 1] : conj(M[i, i + 1])))
else
d *= M[i + 1,i]*(M[i, i]/M[i + 1, i]*M[i + 1, i + 1] - (issym(F) ? M[i + 1, i] : conj(M[i + 1, i])))
end
i += 2
end
end
return d
end

## reconstruct the original matrix
## TODO: understand the procedure described at
## http://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf
10 changes: 9 additions & 1 deletion base/linalg/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,16 @@ A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{
*(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
*(A::StridedMatrix, B::HermOrSym) = A*full(B)

factorize(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A))
bkfact(A::HermOrSym) = bkfact(A.data, symbol(A.uplo), issym(A))
factorize(A::HermOrSym) = bkfact(A)
Copy link
Member

Choose a reason for hiding this comment

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

Should this not try cholfact first @andreasnoack ? Since this is the generic factorization method. Currently:

julia> A = rand(Complex128, 4, 4); A = A'A; H = Hermitian(A);

julia> typeof(factorize(A))
Base.LinAlg.Cholesky{Complex{Float64},Array{Complex{Float64},2}}

julia> typeof(factorize(H))
Base.LinAlg.BunchKaufman{Complex{Float64},Array{Complex{Float64},2}}

and then change this and this to bkfact(Hermitian(A)) and bkfact(Symmetric(A)) ?

Copy link
Contributor

Choose a reason for hiding this comment

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

cholfact can fail, bunch kaufman won't, so why would starting with cholfact be better?

Copy link
Member

Choose a reason for hiding this comment

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

This type of check is already happening in factorize? Since factorize is a general method which essentially should return the "best" factorization I think the same thing should apply if you have a Hermitian input matrix?

Copy link
Member

Choose a reason for hiding this comment

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

(And since #21976 cholfact does not throw)

Copy link
Contributor

Choose a reason for hiding this comment

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

it doesn't throw, but it gives you an answer that you can't use for many indefinite inputs

Copy link
Member

Choose a reason for hiding this comment

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

That is already what is happening (unless the matrix is annotated with Hermitian/Symmetric). Are you suggesting we should not even try the Cholesky factorization in factorize and do this:

diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index ab2dac2..c39940f 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -768,12 +768,7 @@ function factorize(A::StridedMatrix{T}) where T
             return UpperTriangular(A)
         end
         if herm
-            cf = cholfact(A)
-            if cf.info == 0
-                return cf
-            else
-                return factorize(Hermitian(A))
-            end
+            return factorize(Hermitian(A))
         end
         if sym
             return factorize(Symmetric(A))

Copy link
Contributor

Choose a reason for hiding this comment

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

it's a bit of a judgement call on how frequent you expect indefinite input to be and worth benchmarking the cost of trying cholesky before bunch kaufman vs going to bunch kaufman first. The latter has a faster worst case but can be a bit slower on inputs that cholesky works for, which may be more common in some fields

Copy link
Member

Choose a reason for hiding this comment

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

Okay, what I am really after is that my interpretation of the documentation is that this should return a Cholesky factorization since H is positive definite:

julia> A = rand(3, 3); A = A'A; H = Hermitian(A);

julia> isposdef(H)
true

julia> typeof(factorize(H))
Base.LinAlg.BunchKaufman{Float64,Array{Float64,2}}

Copy link
Member

Choose a reason for hiding this comment

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

My original suggestion here would also make it consistent with sparse matrices https://github.com/JuliaLang/julia/blob/master/base/sparse/linalg.jl#L927-L942

Copy link
Contributor

Choose a reason for hiding this comment

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

yeah, part of the problem is cholmod can't do sparse bkfact, there are other libraries but it would take some thought to figure out how to make them sanely coexist


# Is just RealHermSymComplexHerm, but type alias seems to be broken
det{T<:Real,S}(A::Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}) = real(det(bkfact(A)))
det{T<:Real}(A::Symmetric{T}) = det(bkfact(A))
det(A::Symmetric) = det(bkfact(A))

\{T,S<:StridedMatrix}(A::HermOrSym{T,S}, B::StridedVecOrMat) = \(bkfact(A.data, symbol(A.uplo), issym(A)), B)

inv{T<:BlasFloat,S<:StridedMatrix}(A::Hermitian{T,S}) = Hermitian{T,S}(inv(bkfact(A.data, symbol(A.uplo))), A.uplo)
inv{T<:BlasFloat,S<:StridedMatrix}(A::Symmetric{T,S}) = Symmetric{T,S}(inv(bkfact(A.data, symbol(A.uplo), true)), A.uplo)

Expand Down
10 changes: 10 additions & 0 deletions test/linalg/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,16 @@ let n=10
# cond
@test cond(Hermitian(asym)) ≈ cond(asym)

# det
@test det(asym) ≈ det(Hermitian(asym, :U))
@test det(asym) ≈ det(Hermitian(asym, :L))
if eltya <: Real
@test det(asym) ≈ det(Symmetric(asym, :U))
@test det(asym) ≈ det(Symmetric(asym, :L))
end
@test det(a + a.') ≈ det(Symmetric(a + a.', :U))
@test det(a + a.') ≈ det(Symmetric(a + a.', :L))

# rank
let A = a[:,1:5]*a[:,1:5]'
@test rank(A) == rank(Hermitian(A))
Expand Down