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

inv of Symmetric/Hermitian with more eltypes than BlasFloat #22882

Merged
merged 1 commit into from
Aug 4, 2017

Conversation

fredrikekre
Copy link
Member

#22686 and #22767 implemented the infrastructure to invert Symmetric / Hermitian matrices with any element type that supports LU factorization (and backsolve). The relaxations here are necessary as well. For instance, before this inv(::Symmetric{Int}) returned a Matrix rather than Symmetric and inv(::Symmetric{Int}) MethodErrord.

@fredrikekre fredrikekre added the linear algebra Linear algebra label Jul 20, 2017
@tkelman
Copy link
Contributor

tkelman commented Jul 20, 2017

this makes #22686 (comment) more of an issue

@fredrikekre
Copy link
Member Author

I find it unlikely that a matrix which is a result of lufact can not also represent the inverse though?

julia/base/linalg/lu.jl

Lines 317 to 318 in 8cd3c3f

inv!(A::LU{T,<:StridedMatrix}) where {T} =
@assertnonsingular A_ldiv_B!(A.factors, copy(A), eye(T, size(A, 1))) A.info

@test inv(Symmetric(asym))::Symmetric ≈ inv(asym)
@test inv(Hermitian(aherm))::Hermitian ≈ inv(aherm)
for uplo in (:U, :L)
@test inv(Symmetric(asym, uplo))::Symmetric ≈ inv(full(Symmetric(asym, uplo)))
Copy link
Member

Choose a reason for hiding this comment

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

Perhaps Array instead of full? :)

@fredrikekre
Copy link
Member Author

The previous fallback assumes that the inverse will have the same element type as the factorization, so this should not be a problem.

@fredrikekre fredrikekre force-pushed the fe/symmetric-inv branch 2 times, most recently from 56c5122 to b324d0d Compare July 23, 2017 09:46
@fredrikekre
Copy link
Member Author

Do you still think this is problematic @tkelman, or anyone else?

@tkelman
Copy link
Contributor

tkelman commented Jul 25, 2017

Main concern for me is whether this requires additional methods to be defined for it to work on inv of symmetric/hermitian wrappers of user defined array and/or element types

@fredrikekre
Copy link
Member Author

Okay, yea, I did some testing, and for a user defined number type, this does not change anything; if you define all that's needed for lufact you have enough for the backsolve. For a custom array, it is enough with getindex / setindex!. The difference with this PR is that the custom array is used as right hand side in the backsolve. Pre this PR a Matrix is used as right hand side and then converted back to the user defined array type when creating the Symmetric return object.

@@ -364,8 +364,8 @@ function _inv(A::HermOrSym)
end
B
end
inv(A::Hermitian{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Hermitian{T,S}(_inv(A), A.uplo)
inv(A::Symmetric{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Symmetric{T,S}(_inv(A), A.uplo)
inv(A::Hermitian{<:Any,<:StridedMatrix}) = Hermitian(_inv(A), Symbol(A.uplo))
Copy link
Member Author

@fredrikekre fredrikekre Jul 26, 2017

Choose a reason for hiding this comment

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

Re my comment above regarding user define array types, should this be changed to:

inv(A::Hermitian{<:Any,S}) where {S<:StridedMatrix} = (Ai = _inv(A); Hermitian{eltype(Ai), S}(Ai, A.uplo)

such that we enforce a conversion back to the same array type? I noticed that if you don't implement getindex to return your user defined array you get back an Array from the backsolve.
Edit: No, this does not work, it is not always possible to convert back to the same type, e.g. with Ints. So I think the current implementation is the best.

Copy link
Member

@Sacha0 Sacha0 left a comment

Choose a reason for hiding this comment

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

👍! :)

@fredrikekre
Copy link
Member Author

Restarting Travis since there was some time this went through. Hoping to get more than 0 out of 3 this time :trollface:

@tkelman
Copy link
Contributor

tkelman commented Aug 4, 2017

Nope. None your fault though.

@tkelman tkelman merged commit 09f059a into JuliaLang:master Aug 4, 2017
@fredrikekre fredrikekre deleted the fe/symmetric-inv branch August 4, 2017 06:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants