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

Preserve types when adding/subtracting Herm/Sym/UniformScaling #29500

Merged
merged 5 commits into from
Dec 13, 2018

Conversation

dalum
Copy link
Contributor

@dalum dalum commented Oct 3, 2018

This PR implements some basic type–preserving methods to +, -, real and imag for Hermitian and Symmetric matrices, which are true by definition, and I believe were probably just left out due to oversight. In addition to making Hermitian/Symmetric matrices form a group under addition (geeky yay) it is also slightly faster (more yay!). The original reason for this PR, however, was due to some head–scratching over why I couldn't add a complex UniformScaling to a Hermitian matrix, which I believe should be pretty generically implemented here (using the AbstractArray constructor). I actually found that there was a test to check that it didn't work (added in this PR: #19228), but couldn't find a good argument why adding a complex UniformScaling to a Hermitian should throw, except that, well ... it did. So I also removed those two tests.

This should probably be two PRs, functionality wise, but hopefully the changes aren't too plentiful to review.

@fredrikekre fredrikekre added breaking This change will break code linear algebra Linear algebra labels Oct 3, 2018
Copy link
Member

@andreasnoack andreasnoack left a comment

Choose a reason for hiding this comment

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

Thanks for adding these. I've added a few comments. Please also add some tests for real and imag.

# matrix breaks the hermiticity, if the UniformScaling is non-real.
# However, to preserve type stability, we do not special-case a
# UniformScaling{<:Complex} that happens to be real.
function (+)(A::Hermitian, J::UniformScaling{<:Complex})
Copy link
Member

Choose a reason for hiding this comment

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

Couldn't this just be

Suggested change
function (+)(A::Hermitian, J::UniformScaling{<:Complex})
(+)(A::Hermitian{T,S}, J::UniformScaling{<:Complex}) where {T,S} = S(A) + J

Copy link
Member

Choose a reason for hiding this comment

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

That would result in two copies though. Maybe ... = (R = S(A); R .+= J; R).

Copy link
Member

Choose a reason for hiding this comment

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

Or I guess

R = S(A); R[diagind(R)] .+= I.λ; R

since we can't broadcast UniformScaling.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@andreasnoack @fredrikekre I don't think this would be an ideal solution. We need to promote to a complex type before adding to the diagonal, so we would have to special case real Hermitian matrices. But the copytri! may be unnecessary, it was just the fastest conversion I could find.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually, the tests fail for S(A) when S is a SubArray, so copytri! is probably necessary here.

Copy link
Member

Choose a reason for hiding this comment

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

Hm. This was what full used to do and we don't have a complete replacement. It was assumed that using type parameters or explicit constructor would be sufficient but your example is a case where a generic function would be better.

stdlib/LinearAlgebra/src/uniformscaling.jl Outdated Show resolved Hide resolved
stdlib/LinearAlgebra/src/symmetric.jl Outdated Show resolved Hide resolved
@dalum
Copy link
Contributor Author

dalum commented Nov 3, 2018

@andreasnoack Would you mind taking another look at this? I'm not sure if there is a better solution to the conversion issue. 🙂

@StefanKarpinski
Copy link
Member

Yes, please!

Copy link
Member

@andreasnoack andreasnoack left a comment

Choose a reason for hiding this comment

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

Just a small comment but except for that it looks ready. I would have been nice to avoid the calls to _return_type and instead rely on existing promotion but I couldn't come up with a better solution (i.e. one that didn't allocate more).

stdlib/LinearAlgebra/src/uniformscaling.jl Outdated Show resolved Hide resolved
@kshyatt
Copy link
Contributor

kshyatt commented Dec 12, 2018

Can this be merged?

@StefanKarpinski
Copy link
Member

Rerunning CI. @andreasnoack, please merge if CI passes and you feel this is in good shape.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
breaking This change will break code linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants