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

Adding `one` for structured matrices that preserves type #29777

Merged
merged 8 commits into from Feb 6, 2019

Conversation

@mcognetta
Copy link
Contributor

commented Oct 23, 2018

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

This comment has been minimized.

Copy link
@stevengj

stevengj Oct 23, 2018

Member

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.

This comment has been minimized.

Copy link
@stevengj

stevengj Oct 23, 2018

Member

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

This comment has been minimized.

Copy link
@stevengj

stevengj Oct 23, 2018

Member

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.

This comment has been minimized.

Copy link
@fredrikekre

fredrikekre Oct 23, 2018

Member

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

No, this would always return Vector.

This comment has been minimized.

Copy link
@stevengj

stevengj Oct 23, 2018

Member

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

This comment has been minimized.

Copy link
@mcognetta

mcognetta Oct 24, 2018

Author Contributor

@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?

This comment has been minimized.

Copy link
@fredrikekre

fredrikekre Oct 24, 2018

Member

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.

This comment has been minimized.

Copy link
@mcognetta

mcognetta Oct 24, 2018

Author Contributor

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

This comment has been minimized.

Copy link
@stevengj

stevengj Oct 24, 2018

Member

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 comment has been minimized.

Copy link
@mcognetta

mcognetta Oct 24, 2018

Author Contributor

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

This comment has been minimized.

Copy link

commented Oct 29, 2018

This is great! Is it ready?

@mcognetta

This comment has been minimized.

Copy link
Contributor Author

commented Oct 30, 2018

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

This comment has been minimized.

Copy link
Member

commented Oct 30, 2018

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

This comment has been minimized.

Copy link
Contributor

commented Oct 30, 2018

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

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Oct 30, 2018

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

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Oct 30, 2018

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

@GiggleLiu

This comment has been minimized.

Copy link
Contributor

commented Oct 31, 2018

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

This comment has been minimized.

Copy link
Member

commented Oct 31, 2018

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

@mcognetta

This comment has been minimized.

Copy link
Contributor Author

commented Nov 11, 2018

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

@fredrikekre

This comment has been minimized.

Copy link
Member

commented Nov 11, 2018

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

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Nov 12, 2018

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

@mcognetta

This comment has been minimized.

Copy link
Contributor Author

commented Jan 25, 2019

Bump on this one. The CI failure seems erroneous.

@stevengj

This comment has been minimized.

Copy link
Member

commented Feb 6, 2019

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

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

2 of 3 checks passed

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

@mcognetta mcognetta deleted the mcognetta:add_one_structured_matrix branch Feb 6, 2019

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Feb 7, 2019

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

@mcognetta

This comment has been minimized.

Copy link
Contributor Author

commented Feb 7, 2019

Yes, I'll add it shortly.

mcognetta added a commit to mcognetta/julia that referenced this pull request Feb 7, 2019
ararslan added a commit that referenced this pull request Feb 7, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
7 participants
You can’t perform that action at this time.