Adding `one` for structured matrices that preserves type by mcognetta · Pull Request #29777 · JuliaLang/julia

mcognetta · 2018-10-23T09:57:51Z

As with several other PRs, one falls back to generic methods for structured matrices and it does not maintain type (zero does not have these problems).

This PR adds specialized one methods that preserve input type for structured matrices (Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal).

It should be noted that when the container type is a range, it gets promoted to a Vector, which mimics the behavior or zero. Otherwise, the container type is preserved.

julia> D = Diagonal(1:3)
3×3 Diagonal{Int64,UnitRange{Int64}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

julia> zero(D)
3×3 Diagonal{Int64,Array{Int64,1}}:
 0  ⋅  ⋅
 ⋅  0  ⋅
 ⋅  ⋅  0

julia> one(D)
3×3 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> D = Diagonal(sparsevec([1, 2]))
2×2 Diagonal{Int64,SparseVector{Int64,Int64}}:
 1  ⋅
 ⋅  2

julia> zero(D)
2×2 Diagonal{Int64,SparseVector{Int64,Int64}}:
 0  ⋅
 ⋅  0

julia> one(D)
2×2 Diagonal{Int64,SparseVector{Int64,Int64}}:
 1  ⋅
 ⋅  1

Some times:

v1.0

julia> D = Diagonal(rand(10000))
julia> B = Bidiagonal(rand(10000), rand(9999), 'U')
julia> T = Tridiagonal(rand(9999), rand(10000), rand(9999))
julia> S = SymTridiagonal(rand(10000), rand(9999))

julia> @time one(D)
 0.213831 seconds (6 allocations: 762.940 MiB, 10.34% gc time)
julia> @time one(B)
 0.212062 seconds (6 allocations: 762.940 MiB, 9.95% gc time)
julia> @time one(T)
 0.212362 seconds (6 allocations: 762.940 MiB, 10.02% gc time)
julia> @time one(S)
 0.209232 seconds (6 allocations: 762.940 MiB, 9.72% gc time)

julia> one(D)
10000×10000 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮              ⋱            ⋮                      
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  1.0

After PR

julia> D = Diagonal(rand(10000))
julia> B = Bidiagonal(rand(10000), rand(9999), 'U')
julia> T = Tridiagonal(rand(9999), rand(10000), rand(9999))
julia> S = SymTridiagonal(rand(10000), rand(9999))

julia> @time one(D)
 0.000024 seconds (7 allocations: 78.375 KiB)
julia> @time one(B)
 0.000043 seconds (9 allocations: 156.594 KiB)
julia> @time one(T)
 0.000069 seconds (11 allocations: 234.813 KiB)
julia> @time one(S)
 0.000043 seconds (9 allocations: 156.594 KiB)

julia> one(D)
10000×10000 Diagonal{Float64,Array{Float64,1}}:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅   …   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0      ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
 ⋮                        ⋮              ⋱            ⋮                      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …   ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0

… ranges)

stevengj · 2018-10-23T19:53:31Z

+one(A::Diagonal{T}) where T = Diagonal(fill!(similar(A.diag), one(T)))
+one(A::Bidiagonal{T}) where T = Bidiagonal(fill!(similar(A.dv), one(T)), zero(A.ev), A.uplo)
+one(A::Tridiagonal{T}) where T = Tridiagonal(zero(A.du), fill!(similar(A.d), one(T)), zero(A.dl))
+one(A::SymTridiagonal{T}) where T = SymTridiagonal(fill!(similar(A.dv), one(T)), zero(A.ev))


similar(foo) in all of these cases is wrong. It should be similar(foo, typeof(one(T)))

The issue is that, if T is a dimensionful type, one returns a dimensionless type.

You could also use fill(one(T), size(foo)), which might be simpler.

zero(A.dl) etcetera is wrong for a related reason, because zero returns a dimensionful value. You should instead use fill(zero(one(T)), size(A.dl)) or similar.

You could also use fill(one(T), size(foo)), which might be simpler.

No, this would always return Vector.

Would be good to test the dimensionful case … see the tests with Furlongs in the test/triangular.jl file.

@fredrikekre is correct about (not) using fill(one(T), size(foo)). This would not preserve the underlying container type (as an example: sparsevector).

As for the dimensionless type, using similar should an abstractvector with the given eltype, which then promotes the result of one to the appropriate dimensionful type.

For example:

julia> x = Furlong(3) Furlong{1,Int64}(3) julia> D = Diagonal([x, x, x]) 3×3 Diagonal{Furlong{1,Int64},Array{Furlong{1,Int64},1}}: Furlong{1,Int64}(3) ⋅ ⋅ ⋅ Furlong{1,Int64}(3) ⋅ ⋅ ⋅ Furlong{1,Int64}(3) julia> one(D) 3×3 Diagonal{Furlong{1,Int64},Array{Furlong{1,Int64},1}}: Furlong{1,Int64}(1) ⋅ ⋅ ⋅ Furlong{1,Int64}(1) ⋅ ⋅ ⋅ Furlong{1,Int64}(1) julia> y = Furlong{2}(3) Furlong{2,Int64}(3) julia> E = Diagonal([y, y, y]) 3×3 Diagonal{Furlong{2,Int64},Array{Furlong{2,Int64},1}}: Furlong{2,Int64}(3) ⋅ ⋅ ⋅ Furlong{2,Int64}(3) ⋅ ⋅ ⋅ Furlong{2,Int64}(3) julia> one(E) 3×3 Diagonal{Furlong{2,Int64},Array{Furlong{2,Int64},1}}: Furlong{2,Int64}(1) ⋅ ⋅ ⋅ Furlong{2,Int64}(1) ⋅ ⋅ ⋅ Furlong{2,Int64}(1)

Is this not the behavior that we want?

one should return a multiplicative identity, so it should be unitless, e.g.

julia> x = Furlong(3) Furlong{1,Int64}(3) julia> one(x) 1 julia> x * one(x) == x true

and thus,

julia> D = Diagonal([x, x, x]) 3×3 Diagonal{Furlong{1,Int64},Array{Furlong{1,Int64},1}}: Furlong{1,Int64}(3) ⋅ ⋅ ⋅ Furlong{1,Int64}(3) ⋅ ⋅ ⋅ Furlong{1,Int64}(3) julia> one(D) 3×3 Diagonal{Int64,Array{Int64,1}}: <-- dimensionless 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ 1

should be true.

Ah, I interpreted @stevengj 's comment the opposite way. I will update this and add some tests with types that have a dimension.

This would not preserve the underlying container type (as an example: sparsevector).

Who cares? Why does the underlying container type matter? An array of ones is not sparse anyway.

This way, it would be consistent with zero. Not that I think making it a vector is a bad idea though (especially considering the behavior of zero and one for these matrices backed by ranges).

julia> x=sprand(Float64, 10, 10, .1) 10×10 SparseMatrixCSC{Float64,Int64} with 5 stored entries: [6 , 1] = 0.999066 [4 , 4] = 0.45223 [3 , 5] = 0.187222 [6 , 6] = 0.87619 [8 , 9] = 0.926811 julia> zero(x) 10×10 SparseMatrixCSC{Float64,Int64} with 5 stored entries: [6 , 1] = 0.0 [4 , 4] = 0.0 [3 , 5] = 0.0 [6 , 6] = 0.0 [8 , 9] = 0.0 julia> one(x) 10×10 SparseMatrixCSC{Float64,Int64} with 10 stored entries: [1 , 1] = 1.0 [2 , 2] = 1.0 [3 , 3] = 1.0 [4 , 4] = 1.0 [5 , 5] = 1.0 [6 , 6] = 1.0 [7 , 7] = 1.0 [8 , 8] = 1.0 [9 , 9] = 1.0 [10, 10] = 1.0 julia> UpperTriangular(x) 10×10 UpperTriangular{Float64,SparseMatrixCSC{Float64,Int64}}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ 0.0 0.0 0.187222 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ 0.45223 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ 0.87619 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.926811 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 julia> zero(UpperTriangular(x)) 10×10 UpperTriangular{Float64,SparseMatrixCSC{Float64,Int64}}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0

chriscoey · 2018-10-29T20:07:22Z

This is great! Is it ready?

mcognetta · 2018-10-30T12:50:33Z

Sorry, I've been afk working on my thesis. @stevengj approved it but I'd like to know the consensus on preserving container type in these sorts of operations.

StefanKarpinski · 2018-10-30T14:29:43Z

Seems great to me! Unlike zeros, which has been deprecated, zero and one for matrices clearly fall into the functional, value-based class of linear algebra APIs and it's unusual to call them and then mutate the result. Given that, I can't see why we wouldn't want this.

GiggleLiu · 2018-10-30T16:05:11Z

Shouldn't zero and one always maintain type (#29839)?
Type preserving can make programs more predictable.

StefanKarpinski · 2018-10-30T21:11:29Z

No, that's been claimed over there but it's far from concluded.

StefanKarpinski · 2018-10-30T21:12:45Z

Also, maintaining the input type is what these definitions do... or were you just supporting the change?

GiggleLiu · 2018-10-31T02:34:05Z

I have given an example in
#29839
This example shows how type non-preserving can affect program design. This is exactly the bug I encountered in Flux.

Please have a look. @StefanKarpinski

StefanKarpinski · 2018-10-31T02:49:15Z

Yeah, I saw it, it's not really definitive.

mcognetta · 2018-11-11T06:28:38Z

Sorry to ping, but I believe this is ready for a final review. @stevengj @fredrikekre @mschauer

fredrikekre · 2018-11-11T08:21:55Z

LGTM, but are we allowed to merged such (technically) breaking changes? Will we then revert if it breaks some package?

StefanKarpinski · 2018-11-12T02:00:48Z

Yes, mark it as “minor change”. We’ll run PkgEval before releasing 1.1.

mcognetta · 2019-01-25T23:48:41Z

Bump on this one. The CI failure seems erroneous.

stevengj · 2019-02-06T14:22:27Z

Yes, macOS travis failure seems to be a network glitch. (No point in restarting, since nowadays Travis macOS is failing for another reason.)

StefanKarpinski · 2019-02-07T15:10:31Z

@mcognetta: would you mind making a PR to add a NEWS item about this?

mcognetta · 2019-02-07T15:33:45Z

Yes, I'll add it shortly.

mcognetta added 4 commits October 19, 2018 18:20

added one for structured types. not very elegant

38713b6

fix diagonal fill

e507488

adding tests and making the container types stay the same (except for…

392c194

… ranges)

removing small enough to put into separate PR

3d86c23

fredrikekre added linear algebra Linear algebra minor change Marginal behavior change acceptable for a minor release labels Oct 23, 2018

stevengj reviewed Oct 23, 2018

View reviewed changes

fixing the case of types with dimension and adding tests

58b5b31

GiggleLiu mentioned this pull request Oct 29, 2018

zero should be type preserving #29839

Closed

stevengj approved these changes Oct 30, 2018

View reviewed changes

mschauer approved these changes Nov 11, 2018

View reviewed changes

mcognetta added 3 commits January 25, 2019 11:50

fixing merge conflict after merge from master

2928f39

missed merge conflict in the test file

429620a

fixing error in test

36921f8

stevengj merged commit fb9f1fd into JuliaLang:master Feb 6, 2019

mcognetta deleted the add_one_structured_matrix branch February 6, 2019 18:02

StefanKarpinski added the needs news A NEWS entry is required for this change label Feb 7, 2019

mcognetta added a commit to mcognetta/julia that referenced this pull request Feb 7, 2019

added news item for JuliaLang#29777

70b71d0

ararslan pushed a commit that referenced this pull request Feb 7, 2019

NEWS.md item for #29777 (#30997)

abb09f8

Uh oh!

Uh oh!

Conversation

mcognetta commented Oct 23, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

stevengj Oct 23, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

chriscoey commented Oct 29, 2018

Uh oh!

mcognetta commented Oct 30, 2018

Uh oh!

StefanKarpinski commented Oct 30, 2018

Uh oh!

GiggleLiu commented Oct 30, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

StefanKarpinski commented Oct 30, 2018

Uh oh!

StefanKarpinski commented Oct 30, 2018

Uh oh!

GiggleLiu commented Oct 31, 2018

Uh oh!

StefanKarpinski commented Oct 31, 2018

Uh oh!

mcognetta commented Nov 11, 2018

Uh oh!

fredrikekre commented Nov 11, 2018

Uh oh!

StefanKarpinski commented Nov 12, 2018

Uh oh!

mcognetta commented Jan 25, 2019

Uh oh!

stevengj commented Feb 6, 2019

Uh oh!

StefanKarpinski commented Feb 7, 2019

Uh oh!

mcognetta commented Feb 7, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

7 participants

stevengj Oct 23, 2018 •

edited

Loading

GiggleLiu commented Oct 30, 2018 •

edited

Loading