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

andreasnoack · 2018-12-07T09:27:46Z

andyferris · 2018-12-07T12:05:12Z

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 · 2018-12-07T12:22:26Z

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 · 2018-12-07T13:28:39Z

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

Don't allow it
Return Matrix for both I(n) and I(m,n). Very inefficient.
Return Matrix only for I(m,n). Slightly inefficient.
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)
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.

StefanKarpinski · 2018-12-07T15:27:18Z

stdlib/LinearAlgebra/src/uniformscaling.jl

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


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.

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

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.

We don't have FillArrays in the stdlib.

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

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".

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

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.

My prediction is one can't predict when dense, dense diagonal, or purely structural¹ 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)

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.

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.

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

I agree 100%

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

New matrix multiplication API may need the O(1) version of I(n). See: https://github.com/JuliaLang/julia/issues/23919#issuecomment-447216076

andyferris · 2018-12-11T02:05:06Z

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 JuliaLang/LinearAlgebra.jl#454 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 · 2018-12-11T08:07:12Z

@andyferris Your concern here was mentioned in JuliaLang/LinearAlgebra.jl#454. 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.

antoine-levitt · 2018-12-16T08:30:36Z

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 · 2019-01-01T22:23:54Z

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 · 2019-01-03T19:28:47Z

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).

KristofferC · 2019-01-03T19:46:41Z

stdlib/LinearAlgebra/src/uniformscaling.jl

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

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


Should be I instead of J?

JeffBezanson · 2019-01-03T19:52:54Z

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 · 2019-01-03T19:53:35Z

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

StefanKarpinski · 2019-01-03T19:54:07Z

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

mbauman · 2019-01-03T20:46:58Z

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 · 2019-01-09T18:34:24Z

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).

Fixes #23156

stdlib/LinearAlgebra/src/uniformscaling.jl

stdlib/LinearAlgebra/test/uniformscaling.jl

StefanKarpinski · 2019-01-18T15:09:43Z

NEWS.md

+  * `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]).


Are these changes to NEWS unrelated?

StefanKarpinski · 2019-01-22T21:53:31Z

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

andreasnoack · 2019-01-22T21:58:23Z

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

NEWS.md

StefanKarpinski · 2019-01-22T22:03:47Z

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

andreasnoack added linear algebra Linear algebra minor change Marginal behavior change acceptable for a minor release feature Indicates new feature / enhancement requests labels Dec 7, 2018

andreasnoack force-pushed the an/In branch from b73bc32 to b103c35 Compare December 7, 2018 09:31

StefanKarpinski reviewed Dec 7, 2018

View reviewed changes

andreasnoack added the triage This should be discussed on a triage call label Dec 11, 2018

tkf mentioned this pull request Dec 14, 2018

Matrix Multiplication API JuliaLang/LinearAlgebra.jl#473

Closed

KristofferC reviewed Jan 3, 2019

View reviewed changes

stevengj mentioned this pull request Jan 3, 2019

print arrays of Bool as 0/1 #30575

Merged

andreasnoack and others added 2 commits January 17, 2019 13:21

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

c4e4308

Fixes #23156

rename arg and update version compat

11d006a

JeffBezanson force-pushed the an/In branch from 43353c1 to 11d006a Compare January 17, 2019 18:25

stevengj reviewed Jan 17, 2019

View reviewed changes

stdlib/LinearAlgebra/src/uniformscaling.jl Outdated Show resolved Hide resolved

stevengj reviewed Jan 17, 2019

View reviewed changes

stdlib/LinearAlgebra/test/uniformscaling.jl Show resolved Hide resolved

more updates

55d1930

StefanKarpinski added this to the 1.2 milestone Jan 17, 2019

StefanKarpinski removed the triage This should be discussed on a triage call label Jan 17, 2019

JeffBezanson removed the minor change Marginal behavior change acceptable for a minor release label Jan 17, 2019

stevengj approved these changes Jan 18, 2019

View reviewed changes

StefanKarpinski reviewed Jan 18, 2019

View reviewed changes

stdlib/LinearAlgebra/test/uniformscaling.jl Outdated Show resolved Hide resolved

StefanKarpinski reviewed Jan 18, 2019

View reviewed changes

StefanKarpinski self-assigned this Jan 22, 2019

StefanKarpinski added 2 commits January 22, 2019 17:01

fix NEWS merge artifact

7cee3b5

Update stdlib/LinearAlgebra/test/uniformscaling.jl

b723718

StefanKarpinski reviewed Jan 22, 2019

View reviewed changes

NEWS.md Show resolved Hide resolved

Update NEWS.md

3d2cd38

StefanKarpinski merged commit 7c955c4 into master Jan 23, 2019

StefanKarpinski deleted the an/In branch January 23, 2019 05:26

AzamatB mentioned this pull request Feb 5, 2019

Standard MvNormal, Inconsistency with Normal, and I JuliaStats/Distributions.jl#826

Closed

Lecrapouille mentioned this pull request Nov 26, 2024

Matrix identity not calling one() ? JuliaLang/LinearAlgebra.jl#659

Closed

dkarrasch mentioned this pull request Sep 3, 2019

add Kronecker products for LinearMaps JuliaLinearAlgebra/LinearMaps.jl#61

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

andreasnoack commented Dec 7, 2018

andyferris commented Dec 7, 2018

fredrikekre commented Dec 7, 2018

andreasnoack commented Dec 7, 2018 •

edited

Loading

StefanKarpinski Dec 7, 2018 •

edited

Loading

mbauman Dec 7, 2018

jlapeyre Dec 7, 2018

mbauman Dec 7, 2018

jlapeyre Dec 7, 2018 •

edited

Loading

jlapeyre Dec 8, 2018

StefanKarpinski Dec 10, 2018

andreasnoack Dec 11, 2018

jlapeyre Dec 11, 2018 •

edited

Loading

tkf Dec 14, 2018

andyferris commented Dec 11, 2018

andreasnoack commented Dec 11, 2018

antoine-levitt commented Dec 16, 2018

stevengj commented Jan 1, 2019 •

edited

Loading

StefanKarpinski commented Jan 3, 2019

KristofferC Jan 3, 2019

JeffBezanson commented Jan 3, 2019

JeffBezanson commented Jan 3, 2019

StefanKarpinski commented Jan 3, 2019

mbauman commented Jan 3, 2019 •

edited

Loading

stevengj commented Jan 9, 2019 •

edited

Loading

StefanKarpinski Jan 18, 2019

StefanKarpinski commented Jan 22, 2019

andreasnoack commented Jan 22, 2019

StefanKarpinski commented Jan 22, 2019

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

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

Conversation

andreasnoack commented Dec 7, 2018

andyferris commented Dec 7, 2018

fredrikekre commented Dec 7, 2018

andreasnoack commented Dec 7, 2018 • edited Loading

StefanKarpinski Dec 7, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jlapeyre Dec 7, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jlapeyre Dec 11, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andyferris commented Dec 11, 2018

andreasnoack commented Dec 11, 2018

antoine-levitt commented Dec 16, 2018

stevengj commented Jan 1, 2019 • edited Loading

StefanKarpinski commented Jan 3, 2019

Choose a reason for hiding this comment

JeffBezanson commented Jan 3, 2019

JeffBezanson commented Jan 3, 2019

StefanKarpinski commented Jan 3, 2019

mbauman commented Jan 3, 2019 • edited Loading

stevengj commented Jan 9, 2019 • edited Loading

Choose a reason for hiding this comment

StefanKarpinski commented Jan 22, 2019

andreasnoack commented Jan 22, 2019

StefanKarpinski commented Jan 22, 2019

andreasnoack commented Dec 7, 2018 •

edited

Loading

StefanKarpinski Dec 7, 2018 •

edited

Loading

jlapeyre Dec 7, 2018 •

edited

Loading

jlapeyre Dec 11, 2018 •

edited

Loading

stevengj commented Jan 1, 2019 •

edited

Loading

mbauman commented Jan 3, 2019 •

edited

Loading

stevengj commented Jan 9, 2019 •

edited

Loading