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

Overload UniformScaling to make I(n) create a Diagonal matrix. #30298

Merged
merged 6 commits into from Jan 23, 2019

Conversation

@andreasnoack
Copy link
Member

commented Dec 7, 2018

Fixes #23156

@andyferris

This comment has been minimized.

Copy link
Contributor

commented Dec 7, 2018

This is... interesting. I’m not sure I love punning on call like this (though I admit toying with the idea from time to time, I generally find it to be the wrong solution...).

Fixes

That issue isn’t a bug... ;)

@fredrikekre

This comment has been minimized.

Copy link
Member

commented Dec 7, 2018

Should this be I(m, n)? We decided to remove the single integer Matrix(I, n) constructors, #24438 (comment), #24470. But then we can't return a Diagonal... How about I(m, n) but require m == n, and revisit the m != n case if we get some general banded matrix structure in future versions of LinearAlgebra?

@andreasnoack

This comment has been minimized.

Copy link
Member Author

commented Dec 7, 2018

That issue isn’t a bug... ;)

It's just to automatically close the issue but you probably knew that. See "fixed" issue for the discussion.

Should this be I(m, n)

I think I(n) is fine for the square case as it matches the notation for identity matrices with a specific size so requiring two size arguments seems redundant. We are in special territory here anyway. It is an open question what to do with the rectangular case, though. I think the options are

  1. Don't allow it
  2. Return Matrix for both I(n) and I(m,n). Very inefficient.
  3. Return Matrix only for I(m,n). Slightly inefficient.
  4. Define it in SparseArrays and make it sparse. Maybe together with 3) such that the behavior will depend on loading SparseArrays. (I suspect some people will hate this)
  5. Modify Diagonal such that it can be rectangular. It's pretty convenient to know that a type is square but could probably work. This would also be a bit handly for the SVD where S is often rectangular.

But in any case, I'll suggest that we postpone the decision regarding I(m,n) and focus on the square case for now since it's the more common case.

⋅ ⋅ 0.7
```
"""
(J::UniformScaling)(n::Integer) = Diagonal(fill(J.λ, n))

This comment has been minimized.

Copy link
@StefanKarpinski

StefanKarpinski Dec 7, 2018

Member

It's tempting to use something O(1) sized that acts like a vector filled with a single value here. That would make 5I(n) just as efficient as (5I)(n) which, with this implementation are different efficiencies. I was trying to figure out which type in the sea of range types would be the right one but couldn't.

This comment has been minimized.

Copy link
@mbauman

mbauman Dec 7, 2018

Member

StepRangeLen(J.λ, 0, n), but probably a step too cute.

This comment has been minimized.

Copy link
@jlapeyre

jlapeyre Dec 7, 2018

Contributor

StepRangeLen carries redundant data.

There is

using FillArrays
(J::UniformScaling)(n::Integer) = Diagonal(Fill(J.λ, n))

This carries redundant data, too. typeof(1I(3)) is Diagonal{Int64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}.

And the big pile of right brackets is kinda scary.

This comment has been minimized.

Copy link
@mbauman

mbauman Dec 7, 2018

Member

We don't have FillArrays in the stdlib.

This comment has been minimized.

Copy link
@jlapeyre

jlapeyre Dec 7, 2018

Contributor

Yes, I was thinking more long term. Of course an implementation with Eye can't go in v1.1

This comment has been minimized.

Copy link
@jlapeyre

jlapeyre Dec 8, 2018

Contributor

Why not? We can add new features in minor releases. Not sure if we should, but we can.

1.1 branch date approaching: Dec. 6

It seems to me that using FillArrays would involve committing at the last minute to substantially more new design than the current PR. That's why I conflated "should" with "can".

This comment has been minimized.

Copy link
@StefanKarpinski

StefanKarpinski Dec 10, 2018

Member

Ah yes, timing-wise 1.1 is too soon. I thought you were saying in terms of compatibility.

This comment has been minimized.

Copy link
@andreasnoack

andreasnoack Dec 11, 2018

Author Member

I expect that I(n) will mostly be used in places where performance isn't critical so I'm not sure it's worth being too cute here to avoid the extra allocations. Sure, if we had the functionality of FillArrays available then we should use it here but maybe this application isn't sufficient motivation to migrate the functionality to base.

This comment has been minimized.

Copy link
@jlapeyre

jlapeyre Dec 11, 2018

Contributor
  1. My prediction is one can't predict when dense, dense diagonal, or purely structural1 will be convenient or efficient. Offer a few discoverable options and let the market decide. FillArrays already has typeof(copy(Eye(n))) == Diagonal{Float64,Array{Float64,1}} and typeof(Array(Eye(n))) == Array{Float64,2}. (I don't know if there are alternatives to FillArrays)

  2. If Eye(n) is reliable (i.e. works as a dropin replacement for all the code that uses eye), then recommending use FillArrays and Eye(n) doesn't seem overly burdensome.

  3. In order to investigate use patterns and efficiency and the question in (2) above (but mainly to test a related PR): I tried representing identity matrices with each of 1) Eye (Diagonal{Float64,Ones{Float64,1,Tuple{Base.OneTo{Int64}}}}), 2) Diagonal{Float64,Array{Float64,1}}, and 3) Array{Float64,2} in an existing package. The result for the test suite (with it's choice of sizes, types, etc.) is that Array{Float64,2} is fastest, both for jit and execution. This includes allocating every time eye is used. IIRC the time changes by about a factor of two. Maybe it's not too important for most code. I guess some of the gaps in performance will be closed. But, the ability to do linear indexing on small matrices that stay in cache might be hard to beat.

  4. this application isn't sufficient motivation to migrate the functionality to base

I agree 100%

1 "dense", "allocated", "packed"? None are precise.

This comment has been minimized.

Copy link
@tkf

tkf Dec 14, 2018

Contributor

New matrix multiplication API may need the O(1) version of I(n). See: #23919 (comment)

@andyferris

This comment has been minimized.

Copy link
Contributor

commented Dec 11, 2018

I still feel a slight concern on having an object that is both a container and a "callable". There are some functions where a positional argument might either be a function to be called or a container to use (e.g. @ararslan is on a roll with #30323, #30328, etc - and this seems to be a recurring pattern). My concern is at some point for one or more of these methods it might become ambiguous whether an input argument of I should called, indexed or iterated.

Apart from that, I is better than eye, but I'm not sure by how much. I've read the recent posts on #23156 but still feel a bit unconvinced on how much clarity this brings (and in what situations).

It's just to automatically close the issue but you probably knew that. See "fixed" issue for the discussion.

Sorry @andreasnoack - I merely was amused you didn't use "resolves" or "closes" - https://help.github.com/articles/closing-issues-using-keywords/ 😄

@andreasnoack andreasnoack added the triage label Dec 11, 2018

@andreasnoack

This comment has been minimized.

Copy link
Member Author

commented Dec 11, 2018

@andyferris Your concern here was mentioned in #23156. I don't think it will be much of an issue in practice but I can't offer much more than my gut feeling. Allowing I(n) provides some convenience for a situation that seems to be pretty common at the cost of being potentially ambiguous in a situation I think will be pretty uncommon. I've added the triage label such that a decision can be made. I think we know the benefits and drawbacks at this point.

@tkf tkf referenced this pull request Dec 14, 2018
@antoine-levitt

This comment has been minimized.

Copy link
Contributor

commented Dec 16, 2018

I just want to mention that having I(n) return bools is counter intuitive wrt zeros(n, n) for instance. It feels like exposing an implementation detail (bool are the lowest common denominator that works for multiplication, but not usually what you want in linear algebra). Of course it's the only sensible choice given how UniformScaling is implemented, but still. From what I can see, the point of having a way to create simply identity matrices is as a starting point to modify them afterwards, or for teaching purposes. Having it be non mutable boolean seems to me likely to be more confusing than anything.

@stevengj

This comment has been minimized.

Copy link
Member

commented Jan 1, 2019

Related to this, should we change the display of Bool arrays so that the entries display as 1 and 0?

The type information is still printed at the top, the display is more compact, it is more intuitive for algebraic uses, and reflects the fact that Bool <: Number in Julia.

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Jan 3, 2019

should we change the display of Bool arrays so that the entries display as 1 and 0?

I think we should at least try this on master and see how it is for a while. It's easy to revert if we don't like it after a while. People who think they'll hate it may be surprised (or people who think they'll like it may also be surprised and hate it).

@@ -48,6 +48,31 @@ julia> [1 2im 3; 1im 2 3] * I
"""
const I = UniformScaling(true)

"""
(J::UniformScaling)(n::Integer)

This comment has been minimized.

Copy link
@KristofferC

KristofferC Jan 3, 2019

Contributor

Should be I instead of J?

@JeffBezanson

This comment has been minimized.

Copy link
Member

commented Jan 3, 2019

Seems ok to me --- if we do this now we get the maximum amount of time to experiment before releasing it in 1.2.

@JeffBezanson

This comment has been minimized.

Copy link
Member

commented Jan 3, 2019

I also really like the idea of printing Bools as 1 and 0 when possible; much easier to read.

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Jan 3, 2019

Triage seems to be on board with this—let's go for it.

@mbauman

This comment has been minimized.

Copy link
Member

commented Jan 3, 2019

Yes, I'm in favor here, too, especially with the time for experimentation.

I'll be curious if folks will generally be satisfied with a Diagonal return type (and not a plain-old Matrix a la eye). E.g., if you can't use I because it lacks a size, how frequently do you also need to modify an off-diagonal element or pass to C or an ill-considered ::Matrix signature or somesuch?

@stevengj

This comment has been minimized.

Copy link
Member

commented Jan 9, 2019

Now that #30575 is merged, I think it would be good to merge this I(n) PR. Returning a Diagonal rather than Matrix seems like a good idea to me because of the sparse case (e.g. I(n) will appear frequently in kronecker products of sparse matrices).

I'd also be in favor of a future PR extending Diagonal to allow non-square matrices and using that for I(m,n).

@StefanKarpinski StefanKarpinski added this to the 1.2 milestone Jan 17, 2019

@StefanKarpinski StefanKarpinski removed the triage label Jan 17, 2019

NEWS.md Outdated
* `mul!`, `rmul!` and `lmul!` methods for `UniformScaling` ([#29506]).
* `Symmetric` and `Hermitian` matrices now preserve the wrapper when scaled with a number ([#29469]).
* Exponentiation operator `^` now supports raising a `Irrational` to an `AbstractMatrix` power ([#29782]).
* `UniformScaling` instances are now callable such that e.g. `I(3)` will produce a `Diagonal` matrix ([#30298]).

This comment has been minimized.

Copy link
@StefanKarpinski

StefanKarpinski Jan 18, 2019

Member

Are these changes to NEWS unrelated?

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Jan 22, 2019

Bump, @andreasnoack: what's up with the NEWS changes? Merge artifact?

@andreasnoack

This comment has been minimized.

Copy link
Member Author

commented Jan 22, 2019

Yes that must have happened during a rebasing. I'll fix it tomorrow.

@StefanKarpinski StefanKarpinski self-assigned this Jan 22, 2019

NEWS.md Show resolved Hide resolved
@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Jan 22, 2019

I fixed it up on GitHub. If this passes CI, please squash-merge.

@StefanKarpinski StefanKarpinski merged commit 7c955c4 into master Jan 23, 2019

2 of 3 checks passed

continuous-integration/appveyor/pr AppVeyor build failed
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
julia freebsd ci Build done
Details

@StefanKarpinski StefanKarpinski deleted the an/In branch Jan 23, 2019

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.