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.

Sign up for GitHub

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
Merged

Adding one for structured matrices that preserves type #29777

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
38713b6
added one for structured types. not very elegant
mcognetta Oct 19, 2018
e507488
fix diagonal fill
mcognetta Oct 19, 2018
392c194
adding tests and making the container types stay the same (except for…
mcognetta Oct 23, 2018
3d86c23
removing small enough to put into separate PR
mcognetta Oct 23, 2018
58b5b31
fixing the case of types with dimension and adding tests
mcognetta Oct 24, 2018
2928f39
fixing merge conflict after merge from master
mcognetta Jan 25, 2019
429620a
missed merge conflict in the test file
mcognetta Jan 25, 2019
36921f8
fixing error in test
mcognetta Jan 25, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions stdlib/LinearAlgebra/src/special.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,3 +148,8 @@ function fill!(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}, x)
throw(ArgumentError("array of type $(typeof(A)) and size $(size(A)) can
not be filled with $x, since some of its entries are constrained."))
end

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))
Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Member

@stevengj stevengj Oct 23, 2018 •
edited
Loading

Choose a reason for hiding this comment

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

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

Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Member

Choose a reason for hiding this comment

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

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

No, this would always return Vector.

Copy link
Member

Choose a reason for hiding this comment

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

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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

57 changes: 57 additions & 0 deletions stdlib/LinearAlgebra/test/special.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,4 +252,61 @@ end
@test isa((@inferred vcat(Float64[], spzeros(1))), SparseVector)
end

@testset "zero and one for structured matrices" begin
for elty in (Int64, Float64, ComplexF64)
D = Diagonal(rand(elty, 10))
Bu = Bidiagonal(rand(elty, 10), rand(elty, 9), 'U')
Bl = Bidiagonal(rand(elty, 10), rand(elty, 9), 'L')
T = Tridiagonal(rand(elty, 9),rand(elty, 10), rand(elty, 9))
S = SymTridiagonal(rand(elty, 10), rand(elty, 9))
mats = [D, Bu, Bl, T, S]
for A in mats
@test iszero(zero(A))
@test isone(one(A))
end

@test zero(D) isa Diagonal
@test one(D) isa Diagonal

@test zero(Bu) isa Bidiagonal
@test one(Bu) isa Bidiagonal
@test zero(Bl) isa Bidiagonal
@test one(Bl) isa Bidiagonal
@test zero(Bu).uplo == one(Bu).uplo == Bu.uplo
@test zero(Bl).uplo == one(Bl).uplo == Bl.uplo

@test zero(T) isa Tridiagonal
@test one(T) isa Tridiagonal
@test zero(S) isa SymTridiagonal
@test one(S) isa SymTridiagonal
end

# ranges
D = Diagonal(1:10)
Bu = Bidiagonal(1:10, 1:9, 'U')
Bl = Bidiagonal(1:10, 1:9, 'L')
T = Tridiagonal(1:9, 1:10, 1:9)
S = SymTridiagonal(1:10, 1:9)
mats = [D, Bu, Bl, T, S]
for A in mats
@test iszero(zero(A))
@test isone(one(A))
end

@test zero(D) isa Diagonal
@test one(D) isa Diagonal

@test zero(Bu) isa Bidiagonal
@test one(Bu) isa Bidiagonal
@test zero(Bl) isa Bidiagonal
@test one(Bl) isa Bidiagonal
@test zero(Bu).uplo == one(Bu).uplo == Bu.uplo
@test zero(Bl).uplo == one(Bl).uplo == Bl.uplo

@test zero(T) isa Tridiagonal
@test one(T) isa Tridiagonal
@test zero(S) isa SymTridiagonal
@test one(S) isa SymTridiagonal
end

end # module TestSpecial