pivoted cholesky(I(2)) not implemented about julia

[palday]

Is the fix as easy as defining

julia> using LinearAlgebra

julia> function LinearAlgebra.cholesky!(A::Diagonal{T}, ::RowMaximum; tol=0.0, check::Bool=true) where {T <: Real}
           A .= T.(sqrt.(A))
           rank = size(A, 1)
           uplo = 'L'
           info = 0
           return CholeskyPivoted(A, uplo, 1:rank, rank, tol, info)
       end

julia> function LinearAlgebra.cholesky(A::Diagonal, ::RowMaximum; tol=0.0, check::Bool=true)
           A = sqrt.(A)
           rank = size(A, 1)
           uplo = 'L'
           info = 0
           return CholeskyPivoted(A, uplo, 1:rank, rank, tol, info)
       end

julia> cholesky(2 * I(2))
Cholesky{Float64, Diagonal{Float64, Vector{Float64}}}
U factor:
2×2 Diagonal{Float64, Vector{Float64}}:
 1.41421   ⋅ 
  ⋅       1.41421

julia> cholesky(2 * I(2), RowMaximum())
CholeskyPivoted{Float64, Diagonal{Float64, Vector{Float64}}, UnitRange{Int64}}
L factor with rank 2:
2×2 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.41421   ⋅ 
 0.0      1.41421
permutation:
1:2

julia> cholesky!(2 * I(2), RowMaximum())
ERROR: InexactError: Int64(1.4142135623730951)
Stacktrace:
 [1] Int64
   @ ./float.jl:900 [inlined]
 [2] _broadcast_getindex_evalf
   @ ./broadcast.jl:683 [inlined]
 [3] _broadcast_getindex
   @ ./broadcast.jl:656 [inlined]
 [4] copyto!(dest::Diagonal{Int64, Vector{Int64}}, bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal{Int64, Vector{Int64}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, Type{Int64}, Tuple{Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal{Int64, Vector{Int64}}}, Nothing, typeof(sqrt), Tuple{Diagonal{Int64, Vector{Int64}}}}}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.9.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/structuredbroadcast.jl:160
 [5] materialize!
   @ ./broadcast.jl:884 [inlined]
 [6] materialize!
   @ ./broadcast.jl:881 [inlined]
 [7] #cholesky!#7
   @ ./REPL[5]:2 [inlined]
 [8] cholesky!(A::Diagonal{Int64, Vector{Int64}}, ::RowMaximum)
   @ Main ./REPL[5]:1
 [9] top-level scope
   @ REPL[9]:1

julia> cholesky!(4 * I(2), RowMaximum())
CholeskyPivoted{Int64, Diagonal{Int64, Vector{Int64}}, UnitRange{Int64}}
L factor with rank 2:
2×2 LowerTriangular{Int64, Diagonal{Int64, Vector{Int64}}}:
 2  ⋅
 0  2
permutation:
1:2

julia> cholesky!(2.0 * I(2), RowMaximum())
CholeskyPivoted{Float64, Diagonal{Float64, Vector{Float64}}, UnitRange{Int64}}
L factor with rank 2:
2×2 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.41421   ⋅ 
 0.0      1.41421
permutation:
1:2

julia> cholesky(2 * I(2), RowMaximum())
CholeskyPivoted{Float64, Diagonal{Float64, Vector{Float64}}, UnitRange{Int64}}
L factor with rank 2:
2×2 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.41421   ⋅ 
 0.0      1.41421
permutation:
1:2

julia> cholesky(2.0 * I(2), RowMaximum())
CholeskyPivoted{Float64, Diagonal{Float64, Vector{Float64}}, UnitRange{Int64}}
L factor with rank 2:
2×2 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 1.41421   ⋅ 
 0.0      1.41421
permutation:
1:2

?

I'm unsure how I feel about

function LinearAlgebra.cholesky!(A::Diagonal{T}, ::RowMaximum; tol=0.0, check::Bool=true) where {T <: Real}
    A .= T.(sqrt.(A))
    rank = size(A, 1)
    uplo = 'L'
    info = 0
    return CholeskyPivoted(A, uplo, 1:rank, rank, tol, info)
end

It means that sometimes it works for Int and sometimes doesn't. I guess we could be more restrictive and do

function LinearAlgebra.cholesky!(A::Diagonal{<:AbstractFloat}, ::RowMaximum; tol=0.0, check::Bool=true)
    A .= sqrt.(A)
    rank = size(A, 1)
    uplo = 'L'
    info = 0
    return CholeskyPivoted(A, uplo, 1:rank, rank, tol, info)
end

thoughts @dkarrasch ? The other question is whether we should call collect(1:rank) or just return the range. Let me know and I'll open up a PR.

[palday]

Hmmm, I guess I need to also handle the case where one of the diagonal elements is zero. In that case, I would call collect(1:rank) for consistency when there's an actual pivot.

[dkarrasch]

I'm not sure how exactly the pivoting works, but I thought that pivoted Cholesky is rank revealing in the sense that the (almost) zero diagonal entries come at the end:

julia> using LinearAlgebra

julia> D = Diagonal(collect(0:9));

julia> cholesky(Matrix(D), RowMaximum(), check = false)
CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}
U factor with rank 9:
10×10 UpperTriangular{Float64, Matrix{Float64}}:
 3.0  0.0      0.0      0.0      0.0      0.0  0.0      0.0      0.0  0.0
  ⋅   2.82843  0.0      0.0      0.0      0.0  0.0      0.0      0.0  0.0
  ⋅    ⋅       2.64575  0.0      0.0      0.0  0.0      0.0      0.0  0.0
  ⋅    ⋅        ⋅       2.44949  0.0      0.0  0.0      0.0      0.0  0.0
  ⋅    ⋅        ⋅        ⋅       2.23607  0.0  0.0      0.0      0.0  0.0
  ⋅    ⋅        ⋅        ⋅        ⋅       2.0  0.0      0.0      0.0  0.0
  ⋅    ⋅        ⋅        ⋅        ⋅        ⋅   1.73205  0.0      0.0  0.0
  ⋅    ⋅        ⋅        ⋅        ⋅        ⋅    ⋅       1.41421  0.0  0.0
  ⋅    ⋅        ⋅        ⋅        ⋅        ⋅    ⋅        ⋅       1.0  0.0
  ⋅    ⋅        ⋅        ⋅        ⋅        ⋅    ⋅        ⋅        ⋅   0.0
permutation:
10-element Vector{Int64}:
 10
  9
  8
  7
  6
  5
  4
  3
  2
  1

ADDENDUM: Also, in this simple case with a clear zero we need to set check = false because otherwise the matrix method throws a RankDeficientException(1).