Inconsistent handling of `one` and `zero` for special matrix types.

[jakebolewski]

Sometimes one or zero returns a dense matrix, other times a special matrix type. one and zero are not defined for all major special matrix types in LinAlg.

Ex.

julia> A = randn(5,5)
5x5 Array{Float64,2}:
 -0.584781  -0.200855   0.622116   -0.217739     1.1236
  0.134056  -0.136507   2.14805    -0.43222     -1.26277
  0.459692  -0.665234   2.2665     -3.7361      -1.85841
  0.561744  -1.54699    0.0842729  -0.00529171  -0.258032
 -0.108674  -0.538445  -0.290112   -1.05528      0.496317

julia> SymA = A*A'/2
5x5 Array{Float64,2}:
  1.03961      -0.0196852   0.000103685  -0.12706    0.389328
 -0.0196852     3.21606     4.49128       0.397812  -0.367431
  0.000103685   4.49128    11.6015        0.988823   1.33548
 -0.12706       0.397812    0.988823      1.39123    0.312497
  0.389328     -0.367431    1.33548       0.312497   0.872922

julia> BD = Bidiagonal(A, true)
5x5 Base.LinAlg.Bidiagonal{Float64}:
 -0.584781  -0.200855  0.0       0.0          0.0
  0.0       -0.136507  2.14805   0.0          0.0
  0.0        0.0       2.2665   -3.7361       0.0
  0.0        0.0       0.0      -0.00529171  -0.258032
  0.0        0.0       0.0       0.0          0.496317

julia> one(BD)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(BD)
ERROR: MethodError: `similar` has no method matching similar(::Base.LinAlg.Bidiagonal{Float64}, ::Type{Float64}, ::(Int64,Int64))
Closest candidates are:
  similar(::Range{T}, ::Type{T<:Top}, ::(Int64...,))
  similar(::Range{T}, ::Type{T<:Top}, ::(Integer...,))
  similar(::SubArray{T,N,P<:AbstractArray{T,N},I<:(Union(Array{Int64,1},Colon,Int64,Range{Int64})...,),LD}, ::Any, ::(Int64...,))
  ...

 in zero at abstractarray.jl:298

julia> D = Diagonal(A)
5x5 Base.LinAlg.Diagonal{Float64}:
 -0.584781   0.0       0.0      0.0         0.0
  0.0       -0.136507  0.0      0.0         0.0
  0.0        0.0       2.2665   0.0         0.0
  0.0        0.0       0.0     -0.00529171  0.0
  0.0        0.0       0.0      0.0         0.496317

julia> one(D)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(D)
5x5 Base.LinAlg.Diagonal{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> S = Symmetric(SymA)
5x5 Base.LinAlg.Symmetric{Float64,Array{Float64,2}}:
  1.03961      -0.0196852   0.000103685  -0.12706    0.389328
 -0.0196852     3.21606     4.49128       0.397812  -0.367431
  0.000103685   4.49128    11.6015        0.988823   1.33548
 -0.12706       0.397812    0.988823      1.39123    0.312497
  0.389328     -0.367431    1.33548       0.312497   0.872922

julia> one(S)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(S)
ERROR: MethodError: `similar` has no method matching similar(::Base.LinAlg.Symmetric{Float64,Array{Float64,2}}, ::Type{Float64}, ::(Int64,Int64))
Closest candidates are:
  similar(::Range{T}, ::Type{T<:Top}, ::(Int64...,))
  similar(::Range{T}, ::Type{T<:Top}, ::(Integer...,))
  similar(::SubArray{T,N,P<:AbstractArray{T,N},I<:(Union(Array{Int64,1},Colon,Int64,Range{Int64})...,),LD}, ::Any, ::(Int64...,))
  ...

 in zero at abstractarray.jl:298

julia> H = Hermitian(SymA)
5x5 Base.LinAlg.Hermitian{Float64,Array{Float64,2}}:
  1.03961      -0.0196852   0.000103685  -0.12706    0.389328
 -0.0196852     3.21606     4.49128       0.397812  -0.367431
  0.000103685   4.49128    11.6015        0.988823   1.33548
 -0.12706       0.397812    0.988823      1.39123    0.312497
  0.389328     -0.367431    1.33548       0.312497   0.872922

julia> one(H)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(H)
ERROR: MethodError: `similar` has no method matching similar(::Base.LinAlg.Hermitian{Float64,Array{Float64,2}}, ::Type{Float64}, ::(Int64,Int64))
Closest candidates are:
  similar(::Range{T}, ::Type{T<:Top}, ::(Int64...,))
  similar(::Range{T}, ::Type{T<:Top}, ::(Integer...,))
  similar(::SubArray{T,N,P<:AbstractArray{T,N},I<:(Union(Array{Int64,1},Colon,Int64,Range{Int64})...,),LD}, ::Any, ::(Int64...,))
  ...

 in zero at abstractarray.jl:298

julia> UT = UpperTriangular(A)
5x5 Base.LinAlg.UpperTriangular{Float64,Array{Float64,2}}:
 -0.584781  -0.200855  0.622116  -0.217739     1.1236
  0.0       -0.136507  2.14805   -0.43222     -1.26277
  0.0        0.0       2.2665    -3.7361      -1.85841
  0.0        0.0       0.0       -0.00529171  -0.258032
  0.0        0.0       0.0        0.0          0.496317

julia> one(UT)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(UT)
5x5 Base.LinAlg.UpperTriangular{Float64,Array{Float64,2}}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> ST = SymTridiagonal(SymA)
5x5 Base.LinAlg.SymTridiagonal{Float64}:
  1.03961    -0.0196852   0.0       0.0       0.0
 -0.0196852   3.21606     4.49128   0.0       0.0
  0.0         4.49128    11.6015    0.988823  0.0
  0.0         0.0         0.988823  1.39123   0.312497
  0.0         0.0         0.0       0.312497  0.872922

julia> one(ST)
5x5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(ST)
ERROR: MethodError: `similar` has no method matching similar(::Base.LinAlg.SymTridiagonal{Float64}, ::Type{Float64}, ::(Int64,Int64))

[eschnett]

Is there a documentation that specifies what functions a matrix type needs to implement? Of course, this should really be a concept, but until Julia has these, documentation and maybe a generic test function would catch this:

function test_matrix_type(M::Type)
    z = zero(M)
    @test isa(z, M)
    e = one(M)
    @test isa(e, M)
    @test isa(z + e, M)
    ... more tests ...
end

[jiahao]

Is there a documentation that specifies what functions a matrix type needs to implement?

Not formally. The general interface question is discussed in JuliaLang/julia#6975, and I had dabbled in similar testing ideas in my formal algebra of types notebook. Currently the thinking is that one and zero are the multiplicative and additive identities respectively, so where possible they should be defined so that the types are closed under the algebraic ring (+, *).

[eschnett]

Right, the notebook! Cool stuff.

For matrix types (that are receiving a lot of development at the moment), a more stringent requirement may make sense, and things could probably also be tested more easily.

[mbauman]

These all work now, but there is a little disagreement between one and zero:

zero{T}(x::AbstractArray{T}) = fill!(similar(x), zero(T))
one{T}(x::AbstractMatrix{T}) = eye(T, size(x, 1)) # … after ensuring that the matrix is square

They should probably both use similar.

[mbauman]

Hrm, this issue just got a little bit worse due to the above abstract definitions:

julia> zero(Diagonal([1,2,3]))
┌ Warning: `fill!(D::Diagonal, x)` is deprecated, use `LinearAlgebra.fillstored!(D, x)` instead.
│   caller = zero(::LinearAlgebra.Diagonal{Int64,Array{Int64,1}}) at abstractarray.jl:776
└ @ Base abstractarray.jl:776
3×3 LinearAlgebra.Diagonal{Int64,Array{Int64,1}}:
 0  ⋅  ⋅
 ⋅  0  ⋅
 ⋅  ⋅  0

[Sacha0]

fill!(D::Diagonal, x) should come back in 1.0, alleviating this particular difficulty :). (In full: fill!(D::Diagonal, x) didn't conform to fill!'s general semantics, instead behaving like fillstored! née fillslots!, hence the deprecation in 0.7 to the latter. The plan is to reintroduce fill!(D::Diagonal, x) with appropriate semantics in 1.0.)

[brenhinkeller]

Looks like all is well on this front now in the future; I'm going to close as completed but feel free to undo if that's a mistake for some reason I'm missing!

julia> versioninfo()
Julia Version 1.8.2
Commit 36034abf260 (2022-09-29 15:21 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.3.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores

julia> A = randn(5,5)
5×5 Matrix{Float64}:
 -2.46511     0.351087  -0.40178   -0.395841  -1.8384
  0.0284581   0.627275  -0.542904  -0.374328  -1.14128
  0.916458    1.35424   -0.285203  -0.501903   0.838456
  1.27476    -0.104871   0.622806  -0.240977  -0.349744
 -0.415512   -0.602154  -0.178554  -0.423585  -0.355635

julia> using LinearAlgebra

julia> BD = Bidiagonal(A, :U)
5×5 Bidiagonal{Float64, Vector{Float64}}:
 -2.46511  0.351087    ⋅          ⋅          ⋅
   ⋅       0.627275  -0.542904    ⋅          ⋅
   ⋅        ⋅        -0.285203  -0.501903    ⋅
   ⋅        ⋅          ⋅        -0.240977  -0.349744
   ⋅        ⋅          ⋅          ⋅        -0.355635

julia> one(BD)
5×5 Bidiagonal{Float64, Vector{Float64}}:
 1.0  0.0   ⋅    ⋅    ⋅
  ⋅   1.0  0.0   ⋅    ⋅
  ⋅    ⋅   1.0  0.0   ⋅
  ⋅    ⋅    ⋅   1.0  0.0
  ⋅    ⋅    ⋅    ⋅   1.0

julia> zero(BD)
5×5 Bidiagonal{Float64, Vector{Float64}}:
 0.0  0.0   ⋅    ⋅    ⋅
  ⋅   0.0  0.0   ⋅    ⋅
  ⋅    ⋅   0.0  0.0   ⋅
  ⋅    ⋅    ⋅   0.0  0.0
  ⋅    ⋅    ⋅    ⋅   0.0

julia> SymA = A*A'/2
5×5 Matrix{Float64}:
  4.94891   1.30725    -1.50593   -1.34556     0.853041
  1.30725   1.06583     0.130683   0.0608651   0.135917
 -1.50593   0.130683    1.85507    0.33816    -0.615462
 -1.34556   0.0608651   0.33816    1.10215    -0.17564
  0.853041  0.135917   -0.615462  -0.17564     0.436511

julia> S = Symmetric(SymA)
5×5 Symmetric{Float64, Matrix{Float64}}:
  4.94891   1.30725    -1.50593   -1.34556     0.853041
  1.30725   1.06583     0.130683   0.0608651   0.135917
 -1.50593   0.130683    1.85507    0.33816    -0.615462
 -1.34556   0.0608651   0.33816    1.10215    -0.17564
  0.853041  0.135917   -0.615462  -0.17564     0.436511

julia> one(S)
5×5 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(S)
5×5 Symmetric{Float64, Matrix{Float64}}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> H = Hermitian(SymA)
5×5 Hermitian{Float64, Matrix{Float64}}:
  4.94891   1.30725    -1.50593   -1.34556     0.853041
  1.30725   1.06583     0.130683   0.0608651   0.135917
 -1.50593   0.130683    1.85507    0.33816    -0.615462
 -1.34556   0.0608651   0.33816    1.10215    -0.17564
  0.853041  0.135917   -0.615462  -0.17564     0.436511

julia> one(H)
5×5 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

julia> zero(H)
5×5 Hermitian{Float64, Matrix{Float64}}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> ST = SymTridiagonal(SymA)
5×5 SymTridiagonal{Float64, Vector{Float64}}:
 4.94891  1.30725    ⋅          ⋅         ⋅
 1.30725  1.06583   0.130683    ⋅         ⋅
  ⋅       0.130683  1.85507    0.33816    ⋅
  ⋅        ⋅        0.33816    1.10215  -0.17564
  ⋅        ⋅         ⋅        -0.17564   0.436511

julia> one(ST)
5×5 SymTridiagonal{Float64, Vector{Float64}}:
 1.0  0.0   ⋅    ⋅    ⋅
 0.0  1.0  0.0   ⋅    ⋅
  ⋅   0.0  1.0  0.0   ⋅
  ⋅    ⋅   0.0  1.0  0.0
  ⋅    ⋅    ⋅   0.0  1.0

julia> zero(ST)
5×5 SymTridiagonal{Float64, Vector{Float64}}:
 0.0  0.0   ⋅    ⋅    ⋅
 0.0  0.0  0.0   ⋅    ⋅
  ⋅   0.0  0.0  0.0   ⋅
  ⋅    ⋅   0.0  0.0  0.0
  ⋅    ⋅    ⋅   0.0  0.0

julia> zero(Diagonal([1,2,3]))
3×3 Diagonal{Int64, Vector{Int64}}:
 0  ⋅  ⋅
 ⋅  0  ⋅
 ⋅  ⋅  0