isposdef rewrite with help of cholfact

[fredrikekre]

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

[andreasnoack]

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]

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]

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]

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]

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

[tkelman]

why the promotion with Float32 in such a general method?

[fredrikekre]

Copied from here:

julia/base/linalg/cholesky.jl

Lines 279 to 280 in 2d8c0f9

	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?

[fredrikekre]

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.

[tkelman]

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?

[fredrikekre]

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

[andreasnoack]

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.

[Sacha0]

Perhaps fix the odd spacing? :)

[fredrikekre]

Fixed