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.

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

isposdef rewrite with help of cholfact #22245

Merged
merged 1 commit into from
Jun 13, 2017

Conversation

fredrikekre
Copy link
Member

Finally implementing #21559 (comment) after #21595, #21976 and #22188 🎉

I removed the isposdef[!](::Matrix, ::Symbol) methods (similarly as for #22188). They were not documented, and I got the feeling they were more of internal methods to dispatch to LAPACK. This is potentially breaking, since there are no real replacement. The closest would be to deprecate it for isposdef(Hermitian(A, uplo)) but that may throw an ArgumentError in the Hermitian constructor for Matrix{Complex}. Should we:

  1. throw an error saying it is discontinued?
  2. copy the old code to deprecated.jl and throw a depwarn?
  3. keep it?

Fix #20004

@kshyatt kshyatt added the linear algebra Linear algebra label Jun 6, 2017
@andreasnoack
Copy link
Member

Great. I think deprecation for isposdef(Hermitian(A, uplo)) is the best solution

but that may throw an ArgumentError in the Hermitian constructor for Matrix{Complex}

Argh. The Hermitian constructor exception again.

@fredrikekre
Copy link
Member Author

Test failure is caused by this...

julia> typeof(Hermitian(Symmetric(rand(2, 2))))
Hermitian{Float64,Symmetric{Float64,Array{Float64,2}}}

Any objections to make that Hermitian{Float64, Array{Float64,2}} similarly as in #21622? It simplifies methods like this:

isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(Hermitian(A)))

if A would be Symmetric. Otherwise we need a special method for isposdef(::Symmetric).

@fredrikekre
Copy link
Member Author

fredrikekre commented Jun 7, 2017

It is possible to skip the check in the constructor though, so this will work:

@deprecate isposdef(A::AbstractMatrix, UL::Symbol) isposdef(Hermitian{eltype(A), typeof(A)}(A, char_uplo(UL)))

assuming LAPACK.potrf! will return a non-zero info for the case with non-real diagonal elements. EDIT: It doesn't

@fredrikekre
Copy link
Member Author

Actually that method is bugged (essentially the same problem as #22187), so we should probably go with option 1 above?

julia> A = complex.(rand(2, 2), rand(2, 2)); A[2, 1] = conj(A[1, 2]); A
2×2 Array{Complex{Float64},2}:
 0.826553+0.997144im  0.281465+0.162411im
 0.281465-0.162411im  0.146851+0.549507im

julia> ishermitian(A)
false

julia> isposdef(A)
false

julia> isposdef(A, :U)
true

julia> isposdef(Hermitian{eltype(A), typeof(A)}(A, 'U'))
true

Deprecate it as in #22245 (comment) would at least result in the same bugged behavior but feels weird to keep supporting a bugged function ....

@fredrikekre fredrikekre changed the title [WIP] isposdef rewrite with help of cholfact isposdef rewrite with help of cholfact Jun 7, 2017
@fredrikekre
Copy link
Member Author

fredrikekre commented Jun 7, 2017

Rebased after #22264 and added deprecations, should be good to go now.

@Sacha0 Sacha0 added the deprecation This change introduces or involves a deprecation label Jun 7, 2017
isposdef!(S == T ? copy(A) : convert(AbstractMatrix{S}, A))
end
isposdef(A::AbstractMatrix{T}) where {T} =
isposdef!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)))
Copy link
Contributor

Choose a reason for hiding this comment

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

why the promotion with Float32 in such a general method?

Copy link
Member Author

Choose a reason for hiding this comment

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

Copied from here:

cholfact(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}, ::Type{Val{false}}=Val{false}) =
cholfact!(copy_oftype(A, promote_type(typeof(chol(one(eltype(A)))),Float32)))

Essentially to be able to call cholfact inplace.

The alternative would be to define:
isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(A))
that might actually be cleaner, and avoids the copy if A is not Hermitian.
Thoughts?

Copy link
Member Author

Choose a reason for hiding this comment

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

So with isposdef2 as in the comment above we have the following:

julia> for n in (10, 100, 1000)
           A = rand(n,n); B = A'A
           println("non-hermitian $n x $n")
           @btime isposdef($A)
           @btime isposdef2($A)
           println("hermitian $n x $n")
           @btime isposdef($B)
           @btime isposdef2($B)
       end
non-hermitian 10 x 10
  69.209 ns (2 allocations: 912 bytes)
  7.838 ns (0 allocations: 0 bytes)
hermitian 10 x 10
  487.292 ns (8 allocations: 1.05 KiB)
  500.427 ns (9 allocations: 1.08 KiB)
non-hermitian 100 x 100
  4.054 μs (3 allocations: 78.22 KiB)
  7.653 ns (0 allocations: 0 bytes)
hermitian 100 x 100
  37.041 μs (9 allocations: 78.38 KiB)
  36.770 μs (10 allocations: 78.41 KiB)
non-hermitian 1000 x 1000
  808.359 μs (3 allocations: 7.63 MiB)
  7.651 ns (0 allocations: 0 bytes)
hermitian 1000 x 1000
  6.166 ms (9 allocations: 7.63 MiB)
  6.129 ms (10 allocations: 7.63 MiB)

So I will definitely update to the definition in the comment above.

Copy link
Contributor

@tkelman tkelman Jun 8, 2017

Choose a reason for hiding this comment

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

doesn't need to be addressed in this PR, but does that indicate that maybe the other instance of this Float32 promotion should get the same treatment?

Copy link
Member Author

Choose a reason for hiding this comment

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

It's from 186287a but from the comment it seems like it could have been removed in 5e88074?

Copy link
Member

Choose a reason for hiding this comment

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

I'm a bit rusty on the promotion logic here but there is a difference between pivoted and non-pivoted since we don't yet have a generic version for the latter.

@testset "non-symmetric eigen decomposition" begin
d,v = eig(a)
Copy link
Member

Choose a reason for hiding this comment

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

Perhaps fix the odd spacing? :)

Copy link
Member Author

Choose a reason for hiding this comment

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

Fixed

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
deprecation This change introduces or involves a deprecation linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

isposdef() is incorrect
6 participants