implement pinv with regularization for Diagonal matrices?

[daanhb]

It seems that Diagonal matrices in LinearAlgebra do not fully implement the interface of pinv. There is an implementation of pinv(D::Diagonal) here and of pinv(D::Diagonal, rtol) just below. There is no implementation of a pseudo-inverse with regularization, as in pinv(D::Diagonal; atol = ..., rtol = ...).

On the other hand, the generic routine for dense arrays explicitly checks for diagonal matrices and implements a custom pinv for that case with thresholds here. However, this code does not return a diagonal matrix:

julia> using LinearAlgebra

julia> D = Diagonal((1/10).^(1:5))
5×5 Diagonal{Float64, Vector{Float64}}:
 0.1   ⋅     ⋅      ⋅       ⋅
  ⋅   0.01   ⋅      ⋅       ⋅
  ⋅    ⋅    0.001   ⋅       ⋅
  ⋅    ⋅     ⋅     0.0001   ⋅
  ⋅    ⋅     ⋅      ⋅      1.0e-5

julia> pinv(D; atol = 1e-3)
5×5 Matrix{Float64}:
 10.0    0.0     0.0  0.0  0.0
  0.0  100.0     0.0  0.0  0.0
  0.0    0.0  1000.0  0.0  0.0
  0.0    0.0     0.0  0.0  0.0
  0.0    0.0     0.0  0.0  0.0

Might anyone else be interested in having a more complete implementation of pinv for Diagonal matrices which returns a Diagonal?

Combining the implementation in dense.jl with the one in diagonal.jl, it could look something like this:

function LinearAlgebra.pinv(D::Diagonal{T}; atol::Real = 0.0, rtol::Real = (eps(real(float(oneunit(T))))*length(D.diag))*iszero(atol)) where T
    Di = similar(D.diag, typeof(inv(oneunit(T))))
    if !isempty(D.diag)
        maxabsA = maximum(abs, D.diag)
        abstol = max(rtol * maxabsA, atol)
        for i in 1:length(D.diag)
            if abs(D.diag[i]) > abstol
                invD = inv(D.diag[i])
                if isfinite(invD)
                    Di[i] = invD
                    continue
                end
            end
            # fallback
            Di[i] = zero(T)
        end
    end
    Diagonal(Di)
end

There has been quite a discussion about default choices of the regularization parameters (e.g. JuliaLang/LinearAlgebra.jl#652, JuliaLang/LinearAlgebra.jl#444 and issues referenced there). The implementation above does not change that. It is merely about returning a Diagonal matrix to get this behaviour:

julia> pinv(D; atol=1e-3)
5×5 Diagonal{Float64, Vector{Float64}}:
 10.0     ⋅       ⋅    ⋅    ⋅
   ⋅   100.0      ⋅    ⋅    ⋅
   ⋅      ⋅   1000.0   ⋅    ⋅
   ⋅      ⋅       ⋅   0.0   ⋅
   ⋅      ⋅       ⋅    ⋅   0.0

[daanhb]

On the same topic: for consistency it may be worth adding atol and rtol keyword arguments to pinv(x::Number) as well, returning the same as it would for the corresponding 1-by-1 matrix.