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

Add some (matrix) functions for UniformScaling #28872

Merged
merged 25 commits into from
Feb 19, 2020

Conversation

ZeppLu
Copy link
Contributor

@ZeppLu ZeppLu commented Aug 24, 2018

Notations used here:

  • I: identity matrix, known as UniformScaling in julia
  • x: scalar
  • v: column vector
  • M: matrix

Methods list:

  • x .\ I
  • I ^ x, I .^ x
  • exp, log, (a)sin/cos/tan/csc/sec/cot(h)
  • sqrt
  • v * I, v / I
  • real, imag
  • isposdef
  • tr when iszero(J.λ)
  • matrix decompositions
  • broadcast functions of which 0 is a fixed point (so that matrix structure is preserved) (through BroadcastStyle?):
    • real, imag, abs, abs2, angle, conj, sqrt, cbrt, ...
    • float, complex, big, rationalize, ... and other casting functions
  • or more?

@mbauman
Copy link
Member

mbauman commented Aug 24, 2018

This is great — welcome! The failures look to be due to missing imports (or Base. qualifications) of the matrix functions from base.

@mbauman mbauman added the linear algebra Linear algebra label Aug 24, 2018
In fact none of other matrices in LinearAlgebra
specializes cbrt(), I saw sqrt() and then added
cbrt() as well without ever thinking.
changes: M    = ...
     to: M(x) = ...
in order to make such conversion reusable in
following tests
I forgot it...
@ZeppLu
Copy link
Contributor Author

ZeppLu commented Aug 25, 2018

I'm having trouble understanding

\(A::Union{Bidiagonal{T},AbstractTriangular{T}}, J::UniformScaling) where {T<:Number} =
rmul!(inv(A), J.λ)
\(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
\(A::AbstractMatrix, J::UniformScaling) = rmul!(inv(A), J.λ)

IIUC, line 175-176 is redundant, since Union{Bidiagonal{T},AbstractTriangular{T}} has already been covered by AbstractMatrix. Or is there any special cases that should be paid attention to?

@andreasnoack
Copy link
Member

My guess is that they are there to avoid ambiguities. You can try to comment them out and see what happens.

`J` is used by following tests, but I mistakenly changed its value
I made v a scalar, forgive me
`Union{Bidiagonal,AbstractTriangular} \ UniformScaling`
seems redundant because there's another method
`AbstractMatrix \ UniformScaling`
@ZeppLu
Copy link
Contributor Author

ZeppLu commented Aug 27, 2018

Just did a git log -L to find out the first commit that adds this line, which is done by @andreasnoack ! There surely are times a programmer forgets why he wrote a piece of code, and I come across the same feeling now and then, hahaha

@andreasnoack
Copy link
Member

Hi ZeppLu. It's great that you are working on this. However, it would be great if you could push your changes a little less frequently since each time you push changes, it triggers new CI runs.

@ZeppLu
Copy link
Contributor Author

ZeppLu commented Aug 30, 2018

This is my first time working with a CI, apologize for any troubles I've made and thanks for your advice!

@dalum
Copy link
Contributor

dalum commented Aug 31, 2018

Looks good! Maybe you could also add ^(J::UniformScaling, M::Matrix) = exp!(log(J.λ)*M). On that note, maybe ^(x::Number, M::Matrix) = exp!(log(x)*M) should probably be added first? Just a thought. 🙂

@ZeppLu
Copy link
Contributor Author

ZeppLu commented Sep 2, 2018

Adding ^(x::Number, M::Matrix) and exp! to stdlib is beyond the purpose of this PR, which only focuses on UniformScaling. I'd suggest you open another PR.

BTW, I'm going to be quite busy this month (Sep), and will be back here as soon as I'm free again.

*(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ*J2.λ)
*(B::BitArray{2}, J::UniformScaling) = *(Array(B), J::UniformScaling)
*(J::UniformScaling, B::BitArray{2}) = *(J::UniformScaling, Array(B))
*(A::AbstractMatrix, J::UniformScaling) = A*J.λ
*(A::AbstractVecOrMat, J::UniformScaling) = 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.

I don't think this is well-defined for vectors. If we think of J being the λ-fold of the identity matrix, the multiplying a vertical vector from the left is meaningless.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sorry my linear algebra class was not taught in English, I can't get what "the λ-fold of the identity matrix" means, could you explain it in detail?

Copy link
Member

Choose a reason for hiding this comment

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

It's just saying that it's some identity matrix multiplied by λ.

This operation, multiplying a vec*mat is frequently an error. I typically think of an identity matrix as being larger than 1x1, so I initially agreed. But this is well-defined if J is considered to be 1x1... and we do say it automatically infers its size based upon the operation at hand.

I do think some extra consideration is still needed here, though, as multiplying rand(5) * fill(2, 1, 1) yields a 5×1 matrix and not a 5-element vector (as it would return as it's currently defined).

Copy link
Member

Choose a reason for hiding this comment

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

Yes, the issue is the distinction between a simple vector v and the linear map x \mapsto v*x. In the first case, v should be a vector as we know it. In the second, v should be an nx1 matrix. Take your example rand(5) * fill(2, 1, 1). As it stands, it is non-sensical. Checking with @which what is actually done leads to

(*)(a::AbstractVector, B::AbstractMatrix) = reshape(a,length(a),1)*B

You see, to perform the product you first need to reshape the vector into a matrix, and only then multiplication works as expected. This and the lines above the quote are somewhat strange "convenience" functions, I would say. Anyway, if we wanted to be consistent with what we already have, we should have

*(A::AbstractMatrix, J::UniformScaling) = A*J.λ
*(a::AbstractVector, J::UniformScaling) = reshape(a*J.λ, length(a), 1)

Copy link
Member

@mbauman mbauman Mar 28, 2019

Choose a reason for hiding this comment

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

You'll probably enjoy the discussion in #4774 if you've not discovered that one yet. In short, we've decided to allow vectors to participate in the linear algebra of matrices — and in many cases that means they behave as though they have an implicit trailing ×1 dimension. This particular case is covered in #20423. Sure it's implemented as a reshape, but the important part is the behavior — and this is the behavior we've decided to have.

Copy link
Member

@dkarrasch dkarrasch Mar 29, 2019

Choose a reason for hiding this comment

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

Thanks @mbauman, for the link. I remember reading it some long time ago. I agree it's a design decision, and I imagined that behavior is desired. Would

*(a::AbstractVector, J::UniformScaling) = reshape(a, length(a), 1) * J

be then more appropriate? We would leave the multiplication handling to more generic parts of the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks you two! I also consider *(a::AbstractVector, J::UniformScaling) = reshape(a, length(a), 1) * J more reasonable

*(J::UniformScaling, A::AbstractVecOrMat) = J.λ*A
*(x::Number, J::UniformScaling) = UniformScaling(x*J.λ)
*(J::UniformScaling, x::Number) = UniformScaling(J.λ*x)

/(J1::UniformScaling, J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ)
/(J::UniformScaling, A::AbstractMatrix) = lmul!(J.λ, inv(A))
/(A::AbstractMatrix, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ
/(A::AbstractVecOrMat, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : 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.

I think this should remain AbstractMatrix for the same reason as above.

@@ -79,11 +81,15 @@ istril(::UniformScaling) = true
issymmetric(::UniformScaling) = true
ishermitian(J::UniformScaling) = isreal(J.λ)

isposdef(J::UniformScaling) = isreal(J.λ) && J.λ > 0
Copy link
Member

Choose a reason for hiding this comment

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

How about

isposdef(J::UniformScaling) = isposdef(J.λ)

?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

you're right!

- *(v::AbstractVector, J::UniformScaling)
- /(v::AbstractVector, J::UniformScaling)
- isposdef(J::UniformScaling)
@chriscoey
Copy link

hey @ZeppLu, thanks for doing this, any chance you could finish it off?

@chriscoey
Copy link

also, could this fix the broadcasting with I issue at #23197?

@dkarrasch
Copy link
Member

Anyone care to take a quick look? @jagot @JeffreySarnoff ?

@jagot
Copy link
Contributor

jagot commented Feb 19, 2020

I think it looks fine!

@dkarrasch dkarrasch merged commit 4188b2c into JuliaLang:master Feb 19, 2020
@dkarrasch dkarrasch changed the title [WIP] Complete missing methods for UniformScaling Add some (matrix) functions for UniformScaling Feb 19, 2020
@vtjnash
Copy link
Member

vtjnash commented Feb 19, 2020

This introduced some build issues:

WARNING: Method definition isposdef(LinearAlgebra.UniformScaling{T} where T<:Number) in module LinearAlgebra at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:109 overwritten at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:111.
WARNING: Method definition ^(LinearAlgebra.UniformScaling{T} where T<:Number, Number) in module LinearAlgebra at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:118 overwritten at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:269.
WARNING: Method definition broadcasted(typeof(Base.:(^)), LinearAlgebra.UniformScaling{T} where T<:Number, Number) in module LinearAlgebra at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:267 overwritten at /data/vtjnash/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/uniformscaling.jl:272.

dkarrasch added a commit that referenced this pull request Feb 19, 2020
@dkarrasch
Copy link
Member

Sorry! Fix in #34819.

vtjnash pushed a commit that referenced this pull request Feb 20, 2020
birm pushed a commit to birm/julia that referenced this pull request Feb 22, 2020
Co-authored-by: Daniel Karrasch <Daniel.Karrasch@gmx.de>
birm pushed a commit to birm/julia that referenced this pull request Feb 22, 2020
KristofferC pushed a commit that referenced this pull request Apr 11, 2020
Co-authored-by: Daniel Karrasch <Daniel.Karrasch@gmx.de>
KristofferC pushed a commit that referenced this pull request Apr 11, 2020
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.

9 participants